Optical Vortex Modelling

Numerically modelling the evolution of a vortex lattice in an optical beam and comparison with measured results. Possible extension to include a lattice with defects in a linear medium and a lattice in a nonlinear medium.

• http://www.homepages.ucl.ac.uk/~ucapphj/tweezers/people.htm

Current Status

23 August - Convolving Gaussian Example

http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/Progress#Gaussian_Convolving_Function

Need to integrate the following algorithm with the Param2.py model.

23 August - Gird Size Affecting the model.

Testing the model with different grid sizes it is apparent that the -g 20 and -g 10 give different results. This implies that the FFT is being impacted by the grid size and in particular non periodic signals are causing a frequency signal spread and amplitude reduction.

Windowing

There are 2 examples below one in matlab and another in python that illustrate how to apply windowing function in the frequency domain.

Matlab

http://www.mathworks.com/access/helpdesk/help/toolbox/signal/index.html?/access/helpdesk/help/toolbox/signal/f4-9659.html

Python

http://ftp.osc.edu/supercomputing/training/python/python_pt2_0702.pdf

x = r_[0:1:512j]; h01 = hamming(512)
theta = 2*pi*x; F=fftfreq(1024)
y01 = 2*sin(10*theta) + 3*cos(20*theta) + sin(30*theta) + 2*cos(45*theta);
ys= y01*h01; 
Ys= fft(ys, 1024); Y01 = fft(y01, 1024)
Ysn= fftshift(abs(Ys)/max(abs(Ys)))
Y01n = fftshift(abs(Y01)/max(abs(Y01)))
subplot(2,1,1);plot(F,Ysn, linewidth=2);title('FFT of Smoothed Signal');axis([400, 620, -0.01, 1.01])
subplot(2,1,2);plot(F,Y01n,linewidth=2);title('FFT of Unmodified

Simple Example

• apply the technique to a simple example below
• integrate with the numerical model

1. retest model with different grid sizes -g 10 -g 20

http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/VortexLattice#Scalar_Diffraction_Theory

Leakage Explained

When the FFT of a non-periodic signal is computed then the resulting frequency spectrum suffers from leakage. Leakage results in the signal energy smearing out over a wide frequency range in the FFT when it should be in a narrow frequency range. Figure 3 illustrates the effect of leakage. The left-top graph shows a 10 Hz sine wave with amplitude 1.0 that is periodic in the time frame. The resulting FFT (bottom-left) shows a narrow peak at 10 Hz in the frequency axis with a height of 1.0 as expected. Note the dB scale is used to highlight the shape of the FFT at low levels. The right-top graph shows a sine wave that is not periodic in the time frame resulting in leakage in the FFT (bottom-right). The amplitude is less than the expected 1.0 value and the signal energy is more dispersed. The dispersed shape of the FFT makes it more difficult to identify the frequency content of the measured signal.

http://www.lds-group.com/docs/site_documents/AN014%20Understanding%20FFT%20Windows.pdf

21 August - Different Sized Charged Cores

Based on the following parameterized model, with different charged core sizes we can see a different number of similar size core vortex are generated. Whilst it is clear to see the vortex count incrementing between 4 & 5 vortex centers, the count get slightly trickier to determine between the 12 and 15 vortex cores.

With higher numbers of approx 20 the model output cannot be easily determined, similarly with a very low core charge of 1 there is no break up into smaller vortex centers.

charged core 4pi

• bakera@treebeard:~/inglethorpe/summer-project/results$ python Param2.py --steps 35 -g 20 -r 6 -i color -p 4

charged core 5 pi

• bakera@treebeard:~/inglethorpe/summer-project/results$ python Param2.py --steps 35 -g 20 -r 6 -i color -p 5

charge core 12 pi

• bakera@treebeard:~/inglethorpe/summer-project/results$ python Param2.py --steps 35 -g 20 -r 6 -i color -p 12

charged core 15 pi

• bakera@treebeard:~/inglethorpe/summer-project/results$ python Param2.py --steps 35 -g 20 -r 6 -i color -p 15

# Date : 29th July
# Author : Alexander Baker

import sys
from pylab import figure, show, clf, savefig, cm
from numpy.fft import fft, fft2, ifft, ifft2, fftshift
from optparse import OptionParser
from numpy import *

def generateResolution(n):
   return 2**n, (2**n)-1

def store_value(option, opt_str, value, parser):
    setattr(parser.values, option.dest, value)
   
def main():

   parser = OptionParser()

   parser.add_option("-r", "--resolution", action="callback", callback=store_value, type="int", nargs=1, dest="resolution", help="resolution of the grid parameter")   
   parser.add_option("-g", "--grid", action="callback", callback=store_value, type="int", nargs=1, dest="grid", help="grid size parameter")   
   parser.add_option("-s", "--steps", action="callback", callback=store_value, type="int", nargs=1, dest="steps", help="number of steps")
   parser.add_option("-b", "--beam", action="callback", callback=store_value, type="int", nargs=1, dest="beam", help="beam size")
   parser.add_option("-c", "--core", action="callback", callback=store_value, type="float", nargs=1, dest="core", help="beam size")
   parser.add_option("-p", "--phase", action="callback", callback=store_value, type="int", nargs=1, dest="phase", help="phase size parameter")   
   parser.add_option("-i", "--image", action="callback", callback=store_value, type="string", nargs=1, dest="image", help="phase size parameter")   

   (options, args) = parser.parse_args()
   
   grid_resolution = 6       
   grid_size = 8
   steps = 35
   beam = 8
   core = 0.3
   phase = 15
   image = 'bw'
   
   if options.resolution:
      grid_resolution = options.resolution   
      
   if options.grid:
      grid_size = options.grid       

   if options.steps:
      steps = options.steps

   if options.beam:
      beam = options.beam

   if options.core:
      core = options.core
      
   if options.phase:
      phase = options.phase

   if options.image:
      image = options.image

   print '\n######################################################'
   print '# Numerical Lattice Model'
   print '# Alexander Baker - August 2007'
   print '######################################################\n'

   #grid_size = generate(6)[0] * 1.
   #grid_size_1 = generate(6)[1] * 1.
      
   grid_resolution, grid_resolution_1 = generateResolution(grid_resolution)
   #grid_resolution_1 = generateResolution(grid_resolution)[1] * 1.

   #x1 = arange(-8.,-8+63.* 16/64, 8./64)
   #x2 = arange(-8.,-8+63.* 16/64, 8./64)

   #x1 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)
   #x2 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)

   # The grid size effects the number of steps we apply to get to the upper limit

   print "grid_size : [%d]\ngrid_resolution [%d]\ngrid_resolution-1 [%d]" %(grid_size, grid_resolution, grid_resolution_1)

   x1 = arange(-1 * 1.0* grid_size,(-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, (grid_size * 1.0)/grid_resolution)

   x2 = arange(-1 * 1.0 * grid_size,(-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, (grid_size * 1.0)/grid_resolution)

   #xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
   xmin, xmax, ymin, ymax = -1 * 1.0 * grid_size, (-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, -1 * 1.0 * grid_size, (-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution
   extent = xmin, xmax, ymin, ymax

   print "\ngrid extent X-[%f][%f] Y-[%f][%f]" % (xmin, xmax, ymin, ymax)
   print "shape x1", shape(x1)
   print "shape x2", shape(x2)
   print "step [%f]" % ((grid_size * 1.0)/grid_resolution)
   print "start [%f], end[%f]" % (-1 * 1.0* grid_size, grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution)

   [x1,x2] = meshgrid(x1,x2)

   rbeam = beam
   wcore = core
   vpos = 1
   _del = 1


   print "\nrbeam [%d]\nwcore [%f]\nvpos [%f]\n_del [%f]" % (rbeam, wcore, vpos, _del)

   #
   # Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
   #

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
     exp(1j*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos))) * \
          tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
        exp(1j*(arctan2(x2-vpos,x1-vpos) + \
          arctan2(x2+vpos,x1+vpos))) * \
             tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
           exp(1j*(+arctan2(x2,x1) + \
           arctan2(x2+vypos,x1+vxpos))) * \
           tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1+vxpos)**2 + (x2+vypos)**2), 0.5)/wcore)'''


   y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)*exp(1j* phase *(+arctan2(x2,x1))) * tanh(pow((x1**2 + x2**2), 0.5)/wcore)          

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
       exp(1j*15*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos) + \
       arctan2(x2+vpos,x1-vpos) + \
       arctan2(x2-vpos,x1+vpos) )) * \
            tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
            tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
   tanh(pow(((x1-vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''


   fig1 = figure(1)
   fig1.clf()
   ax1a = fig1.add_subplot(121)

   if image == 'bw':
   ax1a.imshow(angle((y)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   else:
   ax1a.imshow(angle((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax1a.set_title(r'Angle')
   ax1b = fig1.add_subplot(122)

   if image == 'bw':
       ax1b.imshow(abs((y)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   else:
       ax1b.imshow(abs((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)

   ax1b.set_title(r'Amplitude')
   savefig('big_start_' + str(wcore)+'_' + str(vpos) +'.png')  

   #print "\nbeam power (start) - [%f]" % (sum(sum(abs(y**2))))

   #u1 = arange(-1.,-1+63.* 2/64, 1./64)
   #u2 = arange(-1.,-1+63.* 2/64, 1./64)

   #u1 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)
   #u2 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)

   # the grid resolution is different since we need to same number of points as the outer
   # grid but over a smaller range... the grid_resolution in this case is dependent on the
   # grid_size and resolution of the whole grid

   u1 = arange(-1.0,-1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), 1.0/grid_resolution)
   u2 = arange(-1.0,-1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), 1.0/grid_resolution)

   u1min, u1max, u2min, u2max = -1.0, -1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), -1.0, -1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution)

   print "\npropagation grid X-[%f][%f] Y-[%f][%f]" % (u1min, u1max, u2min, u2max)

   print "shape u1", shape(u1)
   print "shape u2", shape(u2)

   print "\nbeam power (start) - [%f]" % (sum(sum(abs(y**2))))

   [u1,u2] = meshgrid(u1,u2)

   t = exp(2*pi*1j*(u1**2 + u2**2)*_del)
   w = fftshift(t)

   for i in arange(100,100+steps, 1):
     z = fft2(y)
     zp = z * w   
     yp = ifft2(zp)
     p = (exp(+0.01*pi*1j*(x1**2 + x2**2)*_del + 0.05*pi*1j*y*conj(y))*_del); 
     yp = yp * p
     y = yp
     zp = fft2(yp)
     zp = zp * w
     yp = ifft2(zp)

     fig3 = figure(3)
     fig3.clf()
     ax3 = fig3.add_subplot(111)

     if image == 'bw':   
        ax3.imshow(abs((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
     else:
        ax3.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
     #savefig('split_step_amp_' + str(i) + '.png')
     ax3 = fig3.add_subplot(111)
     if image == 'bw':   
        ax3.imshow(angle((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
     else :
        ax3.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)    
     #savefig('split_step_phase_' + str(i) + '.png')  

     print sum(sum(abs(yp**2))), i-100
   print "beam power (end) - [%f]" % (sum(sum(abs(yp**2))))

   fig2 = figure(2)
   fig2.clf()
   ax2a = fig2.add_subplot(121)
   if image == 'bw':   
      ax2a.imshow(angle((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   else:
      ax2a.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax2a.set_title(r'Angle')
   ax2b = fig2.add_subplot(122)
   if image == 'bw':   
      ax2b.imshow(abs((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   else:
      ax2b.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax2b.set_title(r'Amplitude')
   savefig('big_end_' + str(wcore)+'_' + str(vpos) +'.png')  

   #show()

   print '\ndone. ok'


if __name__ == "__main__":   
   main()      

17th August - Highly Chared Core Images

The following link seems to show the results that we were trying to repicate with the numerical model, where a single highly charged core can break up into a number of separate cores.

http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/VortexLatticeResults

17th August - Higher Resolution and Grid size images

These images were generated from the parameterised code below, these images are generated in a larger grid and with a might higher number of data points, 10^6 per step.

• end:
  

• end good:
  

# Date : 29th July
# Author : Alexander Baker

import sys
from pylab import figure, show, clf, savefig, cm
from numpy.fft import fft, fft2, ifft, ifft2, fftshift
from optparse import OptionParser
from numpy import *

def generateResolution(n):
   return 2**n, (2**n)-1

def store_value(option, opt_str, value, parser):
    setattr(parser.values, option.dest, value)
   
def main():

   parser = OptionParser()

   parser.add_option("-r", "--resolution", action="callback", callback=store_value, type="int", nargs=1, dest="resolution", help="resolution of the grid parameter")   
   parser.add_option("-g", "--grid", action="callback", callback=store_value, type="int", nargs=1, dest="grid", help="grid size parameter")   
   parser.add_option("-s", "--steps", action="callback", callback=store_value, type="int", nargs=1, dest="steps", help="number of steps")
   parser.add_option("-b", "--beam", action="callback", callback=store_value, type="int", nargs=1, dest="beam", help="beam size")
   parser.add_option("-c", "--core", action="callback", callback=store_value, type="float", nargs=1, dest="core", help="beam size")


   (options, args) = parser.parse_args()
   
   grid_resolution = 6       
   grid_size = 8
   steps = 35
   beam = 8
   core = 0.3
   
   if options.resolution:
      grid_resolution = options.resolution   
      
   if options.grid:
      grid_size = options.grid       

   if options.steps:
      steps = options.steps

   if options.beam:
      beam = options.beam

   if options.core:
      core = options.core
      
   print '\n######################################################'
   print '# Numerical Lattice Model'
   print '# Alexander Baker - August 2007'
   print '######################################################\n'

   #grid_size = generate(6)[0] * 1.
   #grid_size_1 = generate(6)[1] * 1.
      
   grid_resolution, grid_resolution_1 = generateResolution(grid_resolution)
   #grid_resolution_1 = generateResolution(grid_resolution)[1] * 1.

   #x1 = arange(-8.,-8+63.* 16/64, 8./64)
   #x2 = arange(-8.,-8+63.* 16/64, 8./64)

   #x1 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)
   #x2 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)

   # The grid size effects the number of steps we apply to get to the upper limit

   print "grid_size : [%d]\ngrid_resolution [%d]\ngrid_resolution-1 [%d]" %(grid_size, grid_resolution, grid_resolution_1)

   x1 = arange(-1 * 1.0* grid_size,(-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, (grid_size * 1.0)/grid_resolution)

   x2 = arange(-1 * 1.0 * grid_size,(-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, (grid_size * 1.0)/grid_resolution)

   #xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
   xmin, xmax, ymin, ymax = -1 * 1.0 * grid_size, (-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution, -1 * 1.0 * grid_size, (-1 * grid_size) + grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution
   extent = xmin, xmax, ymin, ymax

   print "\ngrid extent X-[%f][%f] Y-[%f][%f]" % (xmin, xmax, ymin, ymax)
   print "shape x1", shape(x1)
   print "shape x2", shape(x2)
   print "step [%f]" % ((grid_size * 1.0)/grid_resolution)
   print "start [%f], end[%f]" % (-1 * 1.0* grid_size, grid_resolution_1 * 1.0 * (2 * grid_size)/grid_resolution)

   [x1,x2] = meshgrid(x1,x2)

   rbeam = beam
   wcore = core
   vpos = 1
   _del = 1


   print "\nrbeam [%d]\nwcore [%f]\nvpos [%f]\n_del [%f]" % (rbeam, wcore, vpos, _del)

   #
   # Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
   #

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
     exp(1j*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos))) * \
          tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
        exp(1j*(arctan2(x2-vpos,x1-vpos) + \
          arctan2(x2+vpos,x1+vpos))) * \
             tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
           exp(1j*(+arctan2(x2,x1) + \
           arctan2(x2+vypos,x1+vxpos))) * \
           tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1+vxpos)**2 + (x2+vypos)**2), 0.5)/wcore)'''


   y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)*exp(1j*15*(+arctan2(x2,x1))) * tanh(pow((x1**2 + x2**2), 0.5)/wcore)          

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
       exp(1j*15*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos) + \
       arctan2(x2+vpos,x1-vpos) + \
       arctan2(x2-vpos,x1+vpos) )) * \
            tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
            tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
   tanh(pow(((x1-vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''


   fig1 = figure(1)
   fig1.clf()
   ax1a = fig1.add_subplot(121)
   ax1a.imshow(angle((y)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   ax1a.set_title(r'Angle')
   ax1b = fig1.add_subplot(122)
   ax1b.imshow(abs((y)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   ax1b.set_title(r'Amplitude')
   savefig('big_start_' + str(wcore)+'_' + str(vpos) +'.png')  

   #print "\nbeam power (start) - [%f]" % (sum(sum(abs(y**2))))

   #u1 = arange(-1.,-1+63.* 2/64, 1./64)
   #u2 = arange(-1.,-1+63.* 2/64, 1./64)

   #u1 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)
   #u2 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)

   # the grid resolution is different since we need to same number of points as the outer
   # grid but over a smaller range... the grid_resolution in this case is dependent on the
   # grid_size and resolution of the whole grid

   u1 = arange(-1.0,-1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), 1.0/grid_resolution)
   u2 = arange(-1.0,-1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), 1.0/grid_resolution)

   u1min, u1max, u2min, u2max = -1.0, -1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution), -1.0, -1+1.0*grid_resolution_1* ((2 * 1.0)/grid_resolution)

   print "\npropagation grid X-[%f][%f] Y-[%f][%f]" % (u1min, u1max, u2min, u2max)

   print "shape u1", shape(u1)
   print "shape u2", shape(u2)

   print "\nbeam power (start) - [%f]" % (sum(sum(abs(y**2))))

   [u1,u2] = meshgrid(u1,u2)

   t = exp(2*pi*1j*(u1**2 + u2**2)*_del)
   w = fftshift(t)

   for i in arange(1,steps, 1):
     z = fft2(y)
     zp = z * w   
     yp = ifft2(zp)
     p = (exp(+0.01*pi*1j*(x1**2 + x2**2)*_del + 0.05*pi*1j*y*conj(y))*_del); 
     yp = yp * p
     y = yp
     zp = fft2(yp)
     zp = zp * w
     yp = ifft2(zp)

     fig3 = figure(3)
     fig3.clf()
     ax3 = fig3.add_subplot(111)
     ax3.imshow(abs((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
     #savefig('split_step_amp_' + str(i) + '.png')
     ax3 = fig3.add_subplot(111)
     ax3.imshow(angle((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)  
     #savefig('split_step_phase_' + str(i) + '.png')  

     #print sum(sum(abs(yp**2)))
   print "beam power (end) - [%f]" % (sum(sum(abs(yp**2))))

   fig2 = figure(2)
   fig2.clf()
   ax2a = fig2.add_subplot(121)
   ax2a.imshow(angle((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   ax2a.set_title(r'Angle')
   ax2b = fig2.add_subplot(122)
   ax2b.imshow(abs((yp)), cmap=cm.gist_gray, alpha=.9, interpolation='bilinear', extent=extent)
   ax2b.set_title(r'Amplitude')
   savefig('big_end_' + str(wcore)+'_' + str(vpos) +'.png')  

   #show()

   print '\ndone. ok'


if __name__ == "__main__":   
   main()      

10th August - Grey Scale Images

Tweeking the parameters for the model slightly, with a highly charged core function, the pattern in the lattice appears to form a regular square pattern.

• start highly charged center 8pi with 25 steps:
  

• end highly charged center 8pi with 25 step:
  

9th August - Parameterised Model

Have introduced the concept of parameterisation to the model

• --grid (default 127 elements)
• --resolution (default 64)

These parameters is it hoped can be used to tune the model so that we can get a highly charged core breaking up into a triangular lattice pattern. It is hoped that this will be confirmation of the model.

# Date : 29th July
# Author : Alexander Baker

import sys
from pylab import figure, show, clf, savefig, cm
from numpy.fft import fft, fft2, ifft, ifft2, fftshift
from optparse import OptionParser
from numpy import *

def generateResolution(n):
   return 2**n, (2**n)-1

def store_value(option, opt_str, value, parser):
    setattr(parser.values, option.dest, value)
   
def main():

   parser = OptionParser()

   parser.add_option("-r", "--resolution", action="callback", callback=store_value, type="int", nargs=1, dest="resolution", help="resolution of the grid parameter")   
   parser.add_option("-g", "--grid", action="callback", callback=store_value, type="int", nargs=1, dest="grid", help="grid size parameter")   

   (options, args) = parser.parse_args()
   
   grid_resolution = 6       
   grid_size = 
   
   if options.resolution:
      grid_resolution = options.resolution   
      
   if options.grid:
      
      
      
   #grid_size = generate(6)[0] * 1.
   #grid_size_1 = generate(6)[1] * 1.
      
   grid_resolution = generateResolution(grid_resolution)[0] * 1.
   grid_resolution_1 = generateResolution(grid_resolution)[1] * 1.

   #x1 = arange(-8.,-8+63.* 16/64, 8./64)
   #x2 = arange(-8.,-8+63.* 16/64, 8./64)

   x1 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)
   x2 = arange(-8.,-8+grid_resolution_1* 16/grid_resolution, 8./grid_resolution)

   xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
   extent = xmin, xmax, ymin, ymax

   [x1,x2] = meshgrid(x1,x2)

   rbeam = 8;
   wcore = 0.3
   vpos = 1
   _del = 0.2

   #
   # Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
   #

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
     exp(1j*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos))) * \
          tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
        exp(1j*(arctan2(x2-vpos,x1-vpos) + \
          arctan2(x2+vpos,x1+vpos))) * \
             tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
          tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore)'''

   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
           exp(1j*(+arctan2(x2,x1) + \
           arctan2(x2+vypos,x1+vxpos))) * \
           tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
          tanh(pow(((x1+vxpos)**2 + (x2+vypos)**2), 0.5)/wcore)'''


   '''y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)*exp(1j*4*(+arctan2(x2,x1))) * tanh(pow((x1**2 + x2**2), 0.5)/wcore)'''          

   y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
       exp(1j*15*(+arctan2(x2,x1) + \
       arctan2(x2-vpos,x1-vpos) + \
       arctan2(x2+vpos,x1+vpos) + \
       arctan2(x2+vpos,x1-vpos) + \
       arctan2(x2-vpos,x1+vpos) )) * \
            tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
            tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
            tanh(pow(((x1+vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
   tanh(pow(((x1-vpos)**2 + (x2+vpos)**2), 0.5)/wcore)


   fig1 = figure(1)
   fig1.clf()
   ax1a = fig1.add_subplot(121)
   ax1a.imshow(angle((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax1a.set_title(r'Angle')
   ax1b = fig1.add_subplot(122)
   ax1b.imshow(abs((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax1b.set_title(r'Amplitude')
   savefig('big_start_' + str(wcore)+'_' + str(vpos) +'.png')  


   print sum(sum(abs(y**2)))

   #u1 = arange(-1.,-1+63.* 2/64, 1./64)
   #u2 = arange(-1.,-1+63.* 2/64, 1./64)

   u1 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)
   u2 = arange(-1.,-1+grid_resolution_1* 2/grid_resolution, 1./grid_resolution)

   [u1,u2] = meshgrid(u1,u2)

   t = exp(2*math.pi*1j*(u1**2 + u2**2)*_del)
   w = fftshift(t)

   for i in arange(1, 80, 1):
     z = fft2(y)
     zp = z * w   
     yp = ifft2(zp)
     p = (exp(+0.001*math.pi*1j*(x1**2 + x2**2) + 0.005*math.pi*1j*y*conj(y))); 
     yp = yp * p
     y = yp
     zp = fft2(yp)
     zp = zp * w
     yp = ifft2(zp)

     fig3 = figure(3)
     fig3.clf()
     ax3 = fig3.add_subplot(111)
     ax3.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
     #savefig('split_step_amp_' + str(i) + '.png')
     ax3 = fig3.add_subplot(111)
     ax3.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)  
     #savefig('split_step_phase_' + str(i) + '.png')  

     #print sum(sum(abs(yp**2)))
   print sum(sum(abs(yp**2)))

   fig2 = figure(2)
   fig2.clf()
   ax2a = fig2.add_subplot(121)
   ax2a.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax2a.set_title(r'Angle')
   ax2b = fig2.add_subplot(122)
   ax2b.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
   ax2b.set_title(r'Amplitude')
   savefig('big_end_' + str(wcore)+'_' + str(vpos) +'.png')  

   #show()


if __name__ == "__main__":   
   main()      

7th August - Plots

Have produced some further plots.

• start:
  

• middle:
  

• end:
  

# Date : 7th August
# Author : Alexander Baker

import sys
from pylab import figure, show, clf, savefig, cm
from numpy.fft import fft, fft2, ifft, ifft2, fftshift
from numpy import *

x1 = arange(-8.,-8+63.* 16/64, 8./64)
x2 = arange(-8.,-8+63.* 16/64, 8./64)

xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
extent = xmin, xmax, ymin, ymax

[x1,x2] = meshgrid(x1,x2)

rbeam = 8
wcore = 0.05
vpos = 0.8
_del = 0.445 

#
# Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
#

y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
  exp(1j*(+arctan2(x2,x1) + \
    arctan2(x2-vpos,x1-vpos) + \
    arctan2(x2+vpos,x1+vpos) + \
    arctan2(x2+vpos,x1-vpos) + \
    arctan2(x2-vpos,x1+vpos) + \
    arctan2(x2-2*vpos,x1) + \
    arctan2(x2+2*vpos,x1) + \
    arctan2(x2,x1-2*vpos) + \
    arctan2(x2,x1+2*vpos))) * \
       tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
       tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1+vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1-vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1)**2 + (x2-2*vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1)**2 + (x2+2*vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1-2*vpos)**2 + (x2)**2), 0.5)/wcore) * \
       tanh(pow(((x1+2*vpos)**2 + (x2)**2), 0.5)/wcore)

fig1 = figure(1)
fig1.clf()
ax1a = fig1.add_subplot(121)
ax1a.imshow(angle((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax1a.set_title(r'Angle')
ax1b = fig1.add_subplot(122)
ax1b.imshow(abs((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax1b.set_title(r'Amplitude')


print sum(sum(abs(y**2)))

u1 = arange(-1.,-1+63.* 2/64, 1./64)
u2 = arange(-1.,-1+63.* 2/64, 1./64)

[u1,u2] = meshgrid(u1,u2)

t = exp(2 * math.pi*1j*(u1**2 + u2**2) * _del)
w = fftshift(t)

for i in arange(200, 380, 1):
  z = fft2(y)
  zp = z * w   
  yp = ifft2(zp)
  p = (exp(+0.0001*math.pi*1j*(x1**2 + x2**2)*_del + 0.0001*math.pi*1j*y*conj(y))*_del); 
  yp = yp * p
  y = yp
  zp = fft2(yp)
  zp = zp * w
  yp = ifft2(zp)
  
  fig3 = figure(3)
  fig3.clf()
  ax3 = fig3.add_subplot(111)
  ax3.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
  savefig('split_step_amp_' + str(i) + '.png')
  #ax3 = fig3.add_subplot(111)
  #ax3.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)  
  #savefig('split_step_phase_' + str(i) + '.png')  

  print sum(sum(abs(yp**2)))
print sum(sum(abs(yp**2)))

fig2 = figure(2)
fig2.clf()
ax2a = fig2.add_subplot(121)
ax2a.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax2a.set_title(r'Angle')
ax2b = fig2.add_subplot(122)
ax2b.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax2b.set_title(r'Amplitude')

show()

2nd August - Profiling

There were some concerns earlier about the need to increase the grid size and the speed of the algorithm. I did some profiling with hotshot and kcachegrind, in the way described here. It turns out that most of the time is spent creating and saving the images. When the savefig calls are turned off the calculation gets much faster. Kcachegrind gives you nice overviews of where time is spent in %, as well as time per line of source code and much more.

* without savefig:

• with savefig:
  

One comment: hotshot seems to be less than perfect when a lot of time is spent inside C extension code. There is a new module cProfile (Python 2.5 only) that should be better. The script to get it to work can be found here (haven't tried it).

http://amath.colorado.edu/faculty/fperez/python/profiling/

Ralf Gommers

2nd August - Where next?

I think that putting a number of vortices in and getting them to self-organise would be a good test of whether the code is right. Another is that a single highly charged vortex (i.e. a phase change of 2*pi*m, where m > 1) should break up on propagation into m singly charged vortices (which
presumably arrange into a triangular lattice). There's also more complex phenomena like "vortex pinning", but that's probably a little further in the future.

• self-organising vortex patttern
• single highly charged vortex breaking up into smaller singly charged vortices (triangular?)
• vortex pinning

clear
nv = 4; %number of vortices along one side the hexagon
Nvort = 3*nv^2 - 3*nv + 1 % total vortices in the hexagonal array.
%-------calculate all the vortex co-ordinates-------------
l = 0;
for k = 1:nv;
    for j = 1:(2*nv-k);
        l=l+1;
        V(l,1) = -nv+j + 0.5*(k-1);
        V(l,2) = 0 + (sqrt(3)/2)*(k-1);
    end
end
for m = 1:1.5*nv^2 - 2.5*nv+1;
    W(m,1) = V(2*nv-1+m,1);
    W(m,2) = -V(2*nv-1+m,2);
end
C = cat(1,V,W);
%plot(C(:,1),C(:,2))

%-----------------define spatial grid---------------------
x1 = [-8.:8./64.:-8+63.*16/64.];
x2 = [-8.:8./64.:-8+63.*16/64.];
s1 = size(x1);
s2 = size(x2);
M1 = repmat(x1,s1(2),1);
M2 = repmat(-x2',1,s2(2));
%--------------beam and vortex parameters----------------
rbeam = 8;              % radius of gaussian background
rlens = 4;              % radius of optional lens function to focus beam
wcore = 0.3;            % core size parameter
C1 = vpos*C(:,1);       % vortex x-co-ordinates
C2 = vpos*C(:,2);       % vortex y-coordinates
L = zeros(Nvort,1)-1;   % define array of vortex charges
L(nv) = 0*L(nv);        % manipulate charge of single vortex (e.g. wrong sign, wrong magnitude or zero)
C1(nv) = C1(nv) + 0.0;  % displace a vortex
C2(nv) = C2(nv) + 0.0;  % displace a vortex

del = 0.2;              % step size
n1 = 0.001;             % linear inhomogeneity
n2 = 0.001;             % nonlinearity

%----------------total phase function------------------
yphase = zeros(s1(2)) + 1;
for m = 1:Nvort;
    y = exp(i*atan2(M2-C2(m),M1-C1(m)).*L(m));
    yphase = yphase.*y;
end
%figure(1), pcolor(angle(yphase)), shading interp, colormap gray, axis square
%----------------total core function-------------------
ycore = zeros(s1(2)) + 1;
for n = 1:Nvort;
    if L(n) == 0;       % checks to see if any vortex has zero charge, ie is missing and if so makes core function 1
        y = zeros(s1(2)) + 1;
    else
        y = tanh(sqrt((M1-C1(n)).^2 + (M2-C2(n)).^2 )./wcore);
    end
    ycore = ycore.*y;
end
%figure(2), pcolor(ycore), shading interp, colormap gray, axis square
y = 1.0*exp(-2*(M1.^2 + M2.^2)./rbeam^2).*yphase.*ycore.*exp(0*2*pi*i*(M1.^2 + M2.^2)/rlens);
energy = sum(sum(abs(y.^2)))
%figure(1), surf(y.*conj(y)), axis square, shading interp, colormap gray

%-----------define spatial frequency grid--------------
u1 = [-1.:1./64.:-1+63.*2/64.];
u2 = [-1.:1./64.:-1+63.*2/64.];
t1 = size(u1);
t2 = size(u2);
P1 = repmat(u1,t1(2),1);
P2 = repmat(-u2',1,t2(2));
T = exp(2*pi*i*(P1.^2 + P2.^2)*del);
w = fftshift(T);
%--------------split step propagation------------------
for j = 1:100;
    z = (1/del)*fft2(y);
    zp = z.*w;
    yp = ifft2(zp);
        p = (exp(n1*2*pi*i*(M1.^2 + M2.^2)*del + n2*2*pi*i*y.*conj(y))*del); 
    yp = yp.*p;
    y = yp;
    %-----this bit makes the movie using getframe----
    pcolor(y.*conj(y))
    axis square
    colormap gray
    shading interp
    M(j) = getframe;
end
movie(M,1,10)
%figure(5), surf(y.*conj(y)), axis square, shading interp, colormap gray
%figure(4), pcolor(angle(y)), axis square, shading interp   %if included his plots out the phase function
energy_p = sum(sum(abs(yp.^2)))

1 August - Example Outputs from Model

• start 1:
  

• end 1:
  

• start 2:
  

• end 2:
  

• start 3:
  

• end 3:
  

31 July - Vectorized Functions

Thanks to Ralf Gommers for help on the plots and for pointing me to the vectorized numpy.tanh and numpy.arctan. Looks like some good results coming through at last.

• 9 start:
  

• 9 - mid:
  

• 8 final:
  

• 8 start phase:
  

• 8 final phase:
  

# Date : 29th July
# Author : Alexander Baker

import sys
from pylab import figure, show, clf, savefig, cm
from numpy.fft import fft, fft2, ifft, ifft2, fftshift
from numpy import *

x1 = arange(-8.,-8+63.* 16/64, 8./64)
x2 = arange(-8.,-8+63.* 16/64, 8./64)

xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
extent = xmin, xmax, ymin, ymax

[x1,x2] = meshgrid(x1,x2)

rbeam = 8;
wcore = 0.05
vpos = 0.8

#
# Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
#

y = 1.0*exp(-2*(x1**2 + x2**2)/rbeam**2)* \
  exp(1j*(+arctan2(x2,x1) + \
    arctan2(x2-vpos,x1-vpos) + \
    arctan2(x2+vpos,x1+vpos) + \
    arctan2(x2+vpos,x1-vpos) + \
    arctan2(x2-vpos,x1+vpos) + \
    arctan2(x2-2*vpos,x1) + \
    arctan2(x2+2*vpos,x1) + \
    arctan2(x2,x1-2*vpos) + \
    arctan2(x2,x1+2*vpos))) * \
       tanh(pow((x1**2 + x2**2), 0.5)/wcore) * \
       tanh(pow(((x1-vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1+vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1+vpos)**2 + (x2-vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1-vpos)**2 + (x2+vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1)**2 + (x2-2*vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1)**2 + (x2+2*vpos)**2), 0.5)/wcore) * \
       tanh(pow(((x1-2*vpos)**2 + (x2)**2), 0.5)/wcore) * \
       tanh(pow(((x1+2*vpos)**2 + (x2)**2), 0.5)/wcore)

fig1 = figure(1)
fig1.clf()
ax1a = fig1.add_subplot(121)
ax1a.imshow(angle((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax1a.set_title(r'Angle')
ax1b = fig1.add_subplot(122)
ax1b.imshow(abs((y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax1b.set_title(r'Amplitude')


print sum(sum(abs(y**2)))

u1 = arange(-1.,-1+63.* 2/64, 1./64)
u2 = arange(-1.,-1+63.* 2/64, 1./64)

[u1,u2] = meshgrid(u1,u2)

t = exp(1j*(u1**2 + u2**2))
w = fftshift(t)

for i in arange(1, 80, 1):
  z = fft2(y)
  zp = z * w   
  yp = ifft2(zp)
  p = (exp(+0.01*1j*(x1**2 + x2**2) + 0.05*1j*y*conj(y))); 
  yp = yp * p
  y = yp
  zp = fft2(yp)
  zp = zp * w
  yp = ifft2(zp)
  
  fig3 = figure(3)
  fig3.clf()
  ax3 = fig3.add_subplot(111)
  ax3.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
  savefig('split_step_amp_' + str(i) + '.png')
  ax3 = fig3.add_subplot(111)
  ax3.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)  
  savefig('split_step_phase_' + str(i) + '.png')  

  print sum(sum(abs(yp**2)))
print sum(sum(abs(yp**2)))

fig2 = figure(2)
fig2.clf()
ax2a = fig2.add_subplot(121)
ax2a.imshow(angle((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax2a.set_title(r'Angle')
ax2b = fig2.add_subplot(122)
ax2b.imshow(abs((yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
ax2b.set_title(r'Amplitude')

show()

31 July - Image Directories

Please see these for the full set of images for amplitude and phase propagation of range of vortex configurations.

• http://www.alexfb.com/summer-project/images-eight-vortex/
• http://www.alexfb.com/summer-project/images-four-vortex/
• http://www.alexfb.com/summer-project/images-single-vortex/
• http://www.alexfb.com/summer-project/images-none-vortex/

29th July - Animated Images

• download directory

29th July - Revised Algorithm

Based on a change to the algorithm (see below) we have generated a number of movies with the following script to output some animated results. please see download directory for range of results for :

• none, 1, 4 and 8 vorticies
• movies generated at 12 and 6 frames per second.

See link below for details on script to generate movie files.

• http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/VortexMovie

# Date : 29th July
# Author : Alexander Baker

import math, sys

from pylab import *
from numpy import *

from numpy.fft import *
from numpy.fft import fft
from numpy.fft import ifft
from numpy.fft import fftshift

import matplotlib.backends
import matplotlib.mathtext

x1 = arange(-8.,-8+63.* 16/64, 8./64)
x2 = arange(-8.,-8+63.* 16/64, 8./64)
y = reshape(zeros(15876, dtype=Float), (126,126))

rbeam = 8;
wcore = 0.05
vpos = 0.8

#
# Step 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
#

for yvalue in x2:
   for xvalue in x1:
       oPoint = 1.0*exp(-2*(xvalue**2 + yvalue**2)/rbeam**2)* \
          exp(1j*(+math.atan2(yvalue,xvalue) + \
                 math.atan2(yvalue-vpos,xvalue-vpos) + \
                 math.atan2(yvalue+vpos,xvalue+vpos) + \
                 math.atan2(yvalue+vpos,xvalue-vpos) + \
                 math.atan2(yvalue-vpos,xvalue+vpos) + \
                 math.atan2(yvalue-2*vpos,xvalue) + \
                 math.atan2(yvalue+2*vpos,xvalue) + \
                 math.atan2(yvalue,xvalue-2*vpos) + \
                 math.atan2(yvalue,xvalue+2*vpos))) * \
               math.tanh(pow((xvalue**2 + yvalue**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue-vpos)**2 + (yvalue-vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue+vpos)**2 + (yvalue+vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue+vpos)**2 + (yvalue-vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue-vpos)**2 + (yvalue+vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue)**2 + (yvalue-2*vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue)**2 + (yvalue+2*vpos)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue-2*vpos)**2 + (yvalue)**2), 0.5)/wcore) * \
               math.tanh(pow(((xvalue+2*vpos)**2 + (yvalue)**2), 0.5)/wcore)
       y[yvalue][xvalue] = oPoint

xmin, xmax, ymin, ymax = min(x1), max(x1), min(x2), max(x2)
extent = xmin, xmax, ymin, ymax

[x1,x2] = meshgrid(x1,x2)

im2 = imshow(abs(fftshift(y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
show()
im2 = imshow(angle(fftshift(y)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
show()


print sum(sum(abs(y**2)))

u1 = arange(-1.,-1+63.* 2/64, 1./64)
u2 = arange(-1.,-1+63.* 2/64, 1./64)

[u1,u2] = meshgrid(u1,u2)

t = exp(1j*(u1**2 + u2**2))
w = fftshift(t)

for i in arange(1, 80, 1):
  z = fft2(y)
  zp = z * w   
  yp = ifft2(zp)
  p = (exp(+0.01*1j*(x1**2 + x2**2) + 0.05*1j*y*conj(y))); 
  yp = yp * p
  y = yp
  zp = fft2(yp)
  zp = zp * w
  yp = ifft2(zp)
  
  im2 = imshow(abs(fftshift(yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
  savefig('/home/bakera/inglethorpe/summer-project/images-eight-vortex/split_step_amp_' + str(i) + '.png')
  im2 = imshow(angle(fftshift(yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)  
  savefig('/home/bakera/inglethorpe/summer-project/images-eight-vortex/split_step_phase_' + str(i) + '.png')  

  print sum(sum(abs(yp**2)))
print sum(sum(abs(yp**2)))

im2 = imshow(angle(fftshift(yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
show()
im2 = imshow(abs(fftshift(yp)), cmap=cm.jet, alpha=.9, interpolation='bilinear', extent=extent)
show()

27th July - Revised Algorithm From Phil

• http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/VortexLatticeMatlab

Set of Images for 4 rotating vortex

• 4 vortex start

• 4 vortex focused

• 4 vortex free space

Start and End profile for 2 vortex

• start profile

• end profile

25th July - Grid Size

Questions around the shape and noise in the grid results, we are thinking that increasing the precision of the numerical grid might give us better results. The current results are based on a grid size of 3969 points, however we are thinking that increasing the number of points above this level. We are not sure what level we need to increase too, however the grid sizes as can be seen below can grow substantially. This might prohibit that amount of numerical modelling that we can perform.

>>> for a in range(0, 15, 1):
...     print ("min[%f], 2^n-1[%f], step[%f], max[%f], grid size[%d]") % (-4, pow(2, a)-1, (8.0)/pow(2,a), (pow(2,a)-1)
* ((8.0)/pow(2,a)), (pow(2,a)-1) * (pow(2,a)-1) )
...
min[-4.000000], 2^n-1[0.000000], step[8.000000], max[0.000000], grid size[0]
min[-4.000000], 2^n-1[1.000000], step[4.000000], max[4.000000], grid size[1]
min[-4.000000], 2^n-1[3.000000], step[2.000000], max[6.000000], grid size[9]
min[-4.000000], 2^n-1[7.000000], step[1.000000], max[7.000000], grid size[49]
min[-4.000000], 2^n-1[15.000000], step[0.500000], max[7.500000], grid size[225]
min[-4.000000], 2^n-1[31.000000], step[0.250000], max[7.750000], grid size[961]

*** CURRENT ***
min[-4.000000], 2^n-1[63.000000], step[0.125000], max[7.875000], grid size[3969]

***POSSIBLE EXTENSIONS***
min[-4.000000], 2^n-1[127.000000], step[0.062500], max[7.937500], grid size[16129]
min[-4.000000], 2^n-1[255.000000], step[0.031250], max[7.968750], grid size[65025]
min[-4.000000], 2^n-1[511.000000], step[0.015625], max[7.984375], grid size[261121]
min[-4.000000], 2^n-1[1023.000000], step[0.007813], max[7.992188], grid size[1046529]
min[-4.000000], 2^n-1[2047.000000], step[0.003906], max[7.996094], grid size[4190209]
min[-4.000000], 2^n-1[4095.000000], step[0.001953], max[7.998047], grid size[16769025]
min[-4.000000], 2^n-1[8191.000000], step[0.000977], max[7.999023], grid size[67092481]
min[-4.000000], 2^n-1[16383.000000], step[0.000488], max[7.999512], grid size[268402689]
>>>

20th July - Early Results

Early results from the split step beam propagation method, using the PoonBanjerjee^? algorithm the following plots were obtained Early Results

17th July

Implemented the PoonBanerjee algorithm for split step propagation see link below.

http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/PoonBanerjee

15th July - First Python Cut

First rough cut of a split beam propagation method in 2-d for 1 central vortex. This is a very basic model and am not sure about the symmetry and need to extend to 3-d to be able to simulate a phase variation into the model.

http://www.alexfb.com/cgi-bin/twiki/view/PtPhysics/VortexAlgorithm

9th July - Test Scripts

Have developed approach and some tools to plot the gaussian beam. The put together script to translate test function to frequency domain.

Email from Phil

I've been reading up a bit on the famous split-step beam propagation method in a few books I found in the library. One of them lays it out in a particularly simple flow diagram form:

• 1. Take your initial beam profile (gaussian, a vortex etc etc) at z = 0
• 2. Fourier transform it
• 3. Propagate it a (small!) distance dz, which is done my multiplying by a
• simple exponential factor
• 4. Take the inverse Fourier transform
• 5. Apply the nonlinearity (if any), which again is just multiplying by an exponential
• 6. Go back to step 2 and repeat as many time as you want.

For propagation in a linear medium we can just cut out step 5, and in this case all the propagation can be done in one big step from the initial beam profile.

For the nonlinearity we just have to work out the form, whether it's a change in the refractive index that depends on the intensity, or a spatial variation of the refractive index as can be found in some optical fibres. Ultimately I think it would be good to include both, as one should compensate the other to stop the beam either expanding uncontrollably or collapsing to a filament.

Next Steps

• include the tanh() core function and the exponential phase functions in the z=0 model
• work on how to multiply the angular spectrum from fft with appropriate phase function
• continue on to ifft() function post phase adustments back into time domain.

Numerical Modelling Approach

The approach to modelling a collection of vortices is as follows

• model one vortex at z=0, model the gaussian envelope and the core function
• plot this in 2 dimensions
• develop a scalar diffraction computation that will allow the propagation of beam model
• extend model to multiple vortices
• extend scalar diffraction to support propagation of multi vortex beam

Feel free to post something in here?

This approach is based on the work from the following papers.

• Rozas.PRL.97.pdf: Rozas
• Roux

Gaussian 3d plots

These are some plots using matplotlib and python script below that show different sized gaussian beams.

gaussian plots

# simple script to calculate the underlying 
# gaussian beam function
# Alexander Baker

import math, sys

# using standard python extension libs
# for maths and plotting, scipt generates
# a range of images for different size
# gaussians.

import pylab as p
import matplotlib.axes3d as p3

from numpy import *

font = {'fontname'   : 'Courier',
        'color'      : 'r',
        'fontweight' : 'bold',
        'fontsize'   : 11}

x = arange(-12, 12, 0.1)   
y = arange(-12, 12, 0.1)
[x,y] = p.meshgrid(x,y)

variance = arange(1, 50, 0.1)

fig=p.figure()
ax = p3.Axes3D(fig)

for value in variance:
   z = pow(math.e, -((x**2)/(value**2)+ (y**2)/(value**2)))
   #ax.plot_wireframe(x,y,z)
   ax.plot_surface(x,y,z)
   ax.set_xlabel('X')
   ax.set_ylabel('Y')
   ax.set_zlabel('Z')
   #p.show()
   #save('data', z)
   p.savefig('gaussian' + str(value) + '.png')

Scalar Diffraction Theory

• Basic plot from time to frequency domain

# simple script to calculate the fft for function
# Alexander Baker

import math, sys

# using standard python extension libs
# for maths and plotting

from pylab import *
from numpy import *

from numpy.fft import fft

x = arange(0, 100, 0.01)
y = sin(x/3) + cos(x/10)
subplot(2,1,1)
ylabel('sin(x/3)+cos(x/10)')
grid(True)
plot(x,y)

subplot(2,1,2)
ylabel('fft(sin(x/3)+cos(x/10))')
grid(True)
plot(x[0:50],abs(fft(y)[0:50]))

show()

savefig('fft.png')

• start-amplitude.bmp: amplitude start

Early Results

From Poon Split Beam Propagation Method.

• start amplitude:
  

• phase 0:
  

• phase 21:
  

• phase 125:
  

• amplitude 0:
  

Useful Links

• http://www.scipy.org/Cookbook/Matplotlib/mplot3D
• http://stsdas.stsci.edu/numarray/Doc/node60.html
• http://www.scipy.org/Cookbook/BuildingArrays

• start:
  

• end: