;After compile this file, you can use the functions and procedures
;provided below to calculate electrical field, force and torque
;designed by a holographical phasemask. 

;Simple application: force and torque for planar wave

;1.generate scattering coefficients for normal incident on a sphere
;with radius 0.5um, refraction index 1.5+i*0.1. The incident light has
;default wave length 0.532um in vacumm, particle sits in the default
;medium water. Incident wave is circular polarized.

;IDL>coeff_normal,[0.5,dcomplex(1.5,0.1)],[1,dcomplex(0,1)],coeff=coeff_n

;2.calculate the force and torque since normal incident, there is only
;force and torque in z direction

;IDL> print,coefftoforce(coeff=coeff_n)
; 1.6376768e-16   1.6854324e-16       2.6832659
;IDL> print,coefftotorque(coeff=coeff_n)
; 2.5466212e-15   2.7503788e-15       2.5550933

;3.now we rotate the incident light by Euler angle (0.1,0.5,-0.1), and
;get the scattering coefficients for same particle and medium

;IDL>oldcoeff,coeff_n,0.1,0.5,-0.1,coeff=coeff_i

;4.we can check the torque and force for the new coefficients

;IDL>print,coefftoforce(coeff=coeff_i)
; 1.2799996      0.12842834       2.3547872

;IDL> print,coefftotorque(coeff=coeff_i)
; 1.2188573      0.12229365       2.2423052

;5.it is easy to check the result satisfy geometric relations:
;IDL>print,forcez(coeff_n)*[sin(0.5)*cos(0.1),sin(0.5)*sin(0.1),cos(0.5)]
;IDL>print,coefftoforce(coeff=coeff_i)
; 1.2799995      0.12842833       2.3547873
; 1.2799996      0.12842834       2.3547872
;
;and

;IDL>print,torquez(coeff_n)*[sin(0.5)*cos(0.1),sin(0.5)*sin(0.1),cos(0.5)]
;IDL>print,coefftotorque(coeff=coeff_i)
; 1.2188573      0.12229365       2.2423053
; 1.2188573      0.12229365       2.2423052


;A more advanced example: (in a clean directory, and notice that all
;the default parameters like wave length can be changed by keywords,
;see the comment in the codes for detail):

;1. calculate planar wave spectrum given by a blank phase mask, ie, a
;simple optical tweezer formed in defaut medium water, with defaul
;wave length 0.532um in vacumm.

;IDL>linfo = EMDphasemasktolightinfo(dblarr(800,800))

;2. confirm the spectrum by looking at the intensity profile in water

;IDL>intensity = LinfotoIntensity(linfo,[-2,2,0.1],[-2,2,0.1],[0,1,1])
;IDL>tvscale,bytscl(intensity)

;3. calculate the scattering coefficients when a 0.5um radius silica
;sphere moves along x coordinate within 10 pix range. Particle by
;default sits in water and 1 pix = 0.135um.
;return value is the time used in the calculation

;IDL> time =
; Outfield_multipoints([findgen(1,10),fltarr(1,10),fltarr(1,10)],
;                      [0.5,1.46],linfo)


;4. calculate the optical force from the coefficients which has
;already been saved in to directory './data/'
; force = CoefftoForce()

; MODIFICATION HISTORY:
;
;-
;10/06/08 Bo Sun 1st working version, only fz can be calculated
;10/13/08 Bo Sun added reference and WigerD function
;10/14/08 Bo Sun added FixWd and forcexy function. The force part is
;         completed.
;10/14/08 Bo Sun copy EMD calculation part from old version, this
;         completes the whole calculation from phase mask to force
;10/15/08 Bo Sun added vhot procedure, this will enable automatic
;         compilation
;10/17/08 Bo Sun change loop variable in Outfield_multipoints and
;         linfotocoeff from int to long to prevent potential error.
;10/18/08 Bo Sun add directory test in linfotocoeff to reduce
;         redundant output messages.

;10/19/08 Bo Sun modify force codes in order to store the coefficients 
;         from one batch call of linfotocoeff into a big file, this
;         will increas speed loading the files
;10/29/08 Bo Sun added torque codes, fix an error in Mieexpansionout
;         when particle is absorbing. 
;10/30/08 Bo Sun added examples in the comment part.


; References:
;[BH]     Absorption and Scattering of Light by Small Particles,
;         C.F.Bohren and D.R.Huffman
;[MTL]    Scattering, Absorption and Emission by an Arbitrary 
;         Finite Particle. M.I.Mishchenko, L.D.Travis and A.A.Lacis
;[NLKB]   ToolBox for Calculation of Optical Forces and
;         Torques. T.A.Nieminen, V.L.Y.Loke, G.Knoner and A.M.Branczyk
;[CIGR]   Rapid and Stable Determination of Rotation Matrices
;         Between Spherical Harmonics by Direct Recursion. C.H.Choi,
;         J.Ivanic,M.S.Gordon and K.Ruedenberg.
;[JHCPLM] The Measurable Distinction Between the Spin and 
;         Orbital Angular Momenta of Electromagnetic Radiation,
;         J.H.Crichton and P.L.Marston
;[BSYRDG] Theory of holographic optical trapping, 
;         B. Sun, Y. Roichman and D. G. Grier
;[PRTMPT] P. R. T. Munro and P. Török, Opt. Express 15, 9293 (2007). 
;[EW]     E. Wolf, Proc. Royal Soc. London A 253, 349 (1959). 

;------------------codes----------------------


;****************small utility functions************:

Function coorgenxyz,xrange,yrange,zrange,dim,boundary=boundary

;Generate a 3-d mesh lattice, varing x first.

npts = long(dim[0])*long(dim[1])*long(dim[2])

if(keyword_set(boundary)) then begin
x = double(xrange[0])+(dindgen(dim[0]))*double(xrange[1]-xrange[0])/double(dim[0]-1)
y = double(yrange[0])+(dindgen(dim[1]))*double(yrange[1]-yrange[0])/double(dim[1]-1)
z = double(zrange[0])+(dindgen(dim[2]))*double(zrange[1]-zrange[0])/double(dim[2]-1)
endif

if(not(keyword_set(boundary))) then begin
x = double(xrange[0])+(dindgen(dim[0]))*double(xrange[1]-xrange[0])/double(dim[0])
y = double(yrange[0])+(dindgen(dim[1]))*double(yrange[1]-yrange[0])/double(dim[1])
z = double(zrange[0])+(dindgen(dim[2]))*double(zrange[1]-zrange[0])/double(dim[2])
endif


temp = dblarr(dim[0],dim[1],dim[2])
for i = 0 ,dim[0]-1 do temp[i,*,*]=x[i]
x = reform(temp,1,npts)
temp[*] = 0

for i = 0 ,dim[1]-1 do temp[*,i,*]=y[i]
y = reform(temp,1,npts)
temp[*] = 0

for i = 0 ,dim[2]-1 do temp[*,*,i]=z[i]
z = reform(temp,1,npts)

return,[x,y,z]

end

Function coorgenrtp,rrange,trange,prange,dim,boundary=boundary
;For given [r,t,p] range, generate a mesh of Cartisian coordinate.
convdim = dim
convdim[0]=dim[2]
convdim[1]=dim[1]
convdim[2]=dim[0]

sphcoor = coorgenxyz(prange,trange,rrange,convdim,boundary=boundary)
costheta = cos(sphcoor[1,*])
sintheta = sin(sphcoor[1,*])
cosphi = cos(sphcoor[0,*])
sinphi = sin(sphcoor[0,*])

x = sphcoor[2,*]*sintheta*cosphi
y = sphcoor[2,*]*sintheta*sinphi
z = sphcoor[2,*]*costheta

return,[x,y,z]

end


Function Oldvector,vector,theta,phi
;This rotates a vector by Euler angle(-phi,theta,phi)
;adopted from chapter 2 of [MTL] 
tr = double(-theta)
pr = double(phi)

x0 = vector[0,*]
y0 = vector[1,*]
z0 = vector[2,*]

T00 = (cos(tr)-1.d)*(cos(pr)^2)+1.d
T01 = sin(pr)*cos(pr)*(cos(tr)-1.d)
T02 = -sin(tr)*cos(pr)
T10 = T01
T11 = (sin(pr)^2)*(cos(tr)-1.d)+1.d
T12 = -sin(tr)*sin(pr)
T20 = -T02
T21 = -T12
T22 = cos(tr)

x = T00*x0+T01*y0+T02*z0
y = T10*x0+T11*y0+T12*z0
z = T20*x0+T21*y0+T22*z0
return,[x,y,z]
end


Function tensorproduct, Vec1,Vec2
;Calculate the cross product of Vec1 X Vec2
vx = Vec1[1,*]*Vec2[2,*]-Vec1[2,*]*Vec2[1,*]
vy = Vec1[2,*]*Vec2[0,*]-Vec1[0,*]*Vec2[2,*]
vz = Vec1[0,*]*Vec2[1,*]-Vec1[1,*]*Vec2[0,*]
return, [vx,vy,vz]
end

Function scalorproduct, Vec1,Vec2
;Calculate the scalor product of Vec1 . Vec2
return, Vec1[0,*]*vec2[0,*]+Vec1[1,*]*vec2[1,*]+Vec1[2,*]*vec2[2,*]
end

Function unitvector,Vec
;Normalize a vector
norm = sqrt(Vec[0,*]^2+Vec[1,*]^2+Vec[2,*]^2)
univ = [Vec[0,*]/norm,Vec[1,*]/norm,Vec[2,*]/norm]
return,univ
end

Function scalorvectorproduct, Scalor,Vec
;Multiply a scalor to a vector
return, [Vec[0,*]*Scalor[0,*],Vec[1,*]*Scalor[0,*],Vec[2,*]*Scalor[0,*]]
end

Function fixdir,dir
s = strlen(dir)
if(s eq 0) then return, './'
last = strpos(dir,'/',/reverse_search)
if (last ne (strlen(dir)-1)) then dir = dir+'/'
return, dir
end



;**********Wigner Rotational codes*************
Function WignerD,al,be,ga,nc
;This function will calculate the Wigner rotational matrix for given
;Euler angles al, be, ga. Algorithm adopted from reference [CIGR]
;return value is a (2*nc+1) x (2*nc+1) x (nc+1) dimensional array, and
;D[m+nc,m'+nc,n] contains the value of Wiger-D matrix element
;D^n_{m,m'}
 
D = dcomplexarr(2*nc+1,2*nc+1,nc+1)
;generate the first order
;D^1_{0,0}
D[nc,nc,1] = cos(be)
;D^1_{+,0}
D[nc+1,nc,1] = sin(be)/sqrt(2)*dcomplex(-cos(al),sin(al))
;D^1_{-,0}
D[nc-1,nc,1] = sin(be)/sqrt(2)*dcomplex(cos(al),sin(al))
;D^1_{0-}
D[nc,nc-1,1] = -sin(be)/sqrt(2)*dcomplex(cos(ga),sin(ga))
;D^1_{0+}
D[nc,nc+1,1] = sin(be)/sqrt(2)*dcomplex(cos(ga),-sin(ga))
;D^1_{--}
D[nc-1,nc-1,1] = (1+cos(be))/2*dcomplex(cos(al+ga),sin(al+ga))
;D^1_{-+}
D[nc-1,nc+1,1] = (1-cos(be))/2*dcomplex(cos(al-ga),sin(al-ga))
;D^1_{+-}
D[nc+1,nc-1,1] = (1-cos(be))/2*dcomplex(cos(al-ga),sin(ga-al))
;D^1_{++}
D[nc+1,nc+1,1] = (1+cos(be))/2*dcomplex(cos(al+ga),-sin(al+ga))
;generate the higher order
for n=2d, nc-1 do begin
;for -n< m'<n
indm = findgen(2*n+1)-n
for mp = -n+1d, n-1 do begin
a = sqrt((n+indm)*(n-indm)/(n+mp)/(n-mp))
b1 = sqrt(((n+indm)*(n+indm-1)/(n+mp)/(n-mp)/2.)>0)
b2 = sqrt(((n-indm)*(n-indm-1)/(n+mp)/(n-mp)/2.)>0)
D[nc+indm,mp+nc,n] = a[*]*D[nc,nc,1]*D[indm+nc,mp+nc,n-1]+$
                b1[*]*D[nc+1,nc,1]*D[indm-1+nc,mp+nc,n-1]+$
                b2[*]*D[nc-1,nc,1]*D[indm+1+nc,mp+nc,n-1]
endfor
;for m' = -n

c = sqrt((n+indm)*(n-indm)/n/(2*n-1))
d1 = sqrt(((n+indm)*(n-1+indm)/(2*n)/(2*n-1))>0)
d2 = sqrt(((n-indm)*(n-1-indm)/(2*n)/(2*n-1))>0)
D[nc+indm,-n+nc,n] = c[*]*D[nc,nc-1,1]*D[indm+nc,nc+1-n,n-1]+$
                d1[*]*D[nc+1,nc-1,1]*D[indm-1+nc,nc+1-n,n-1]+$
                d2[*]*D[nc-1,nc-1,1]*D[indm+1+nc,nc+1-n,n-1]

;for m' = n

;c = sqrt((n+indm)*(n-indm)/n/(2*n-1))
;d1 = sqrt(((n+indm)*(n-1+indm)/(2*n)/(2*n-1))>0)
;d2 = sqrt(((n-indm)*(n-1-indm)/(2*n)/(2*n-1))>0)
D[nc+indm,n+nc,n] = c[*]*D[nc,nc+1,1]*D[indm+nc,nc-1+n,n-1]+$
               d1[*]*D[nc+1,nc+1,1]*D[indm-1+nc,nc-1+n,n-1]+$
               d2[*]*D[nc-1,nc+1,1]*D[indm+1+nc,nc-1+n,n-1]

endfor
return,D
end

pro FixWd,libdir,nc,wdx=wdx,wdy=wdy
;This procedure will look into libdir for two files "WDx.gdf" and
;"WDy.gdf". They are the Wigner-D matrix used for calculating force
;and torque in x and y direction. If they exist, and their order are
;higher than nc, then return the two matrices. In all the other cases,
;two new matrices with order nc are generated, returned and stored in
;libdir.
libdir = fixdir(libdir)
redoflag = 1
fwdx = file_info(libdir+'WDx.gdf')
;check the format of wdx

if(fwdx.exists) then begin
wdx = read_gdf(libdir+'WDx.gdf')
swdx = size(wdx)
if((swdx[0] eq 3)and(swdx[3] eq nc+1)and(swdx[1] eq swdx[2])$
   and(swdx[1] eq 2*swdx[3]-1)) then redoflag = 0
endif
if(redoflag) then begin
wdx = WignerD(0.,!dpi/2,0.,nc)
write_gdf,wdx,libdir+'WDx.gdf'
endif

redoflag = 1
fwdy = file_info(libdir+'WDy.gdf')
;check the format of wdy

if(fwdx.exists) then begin
wdy = read_gdf(libdir+'WDy.gdf')
swdy = size(wdy)
;check the format of wdx
if((swdy[0] eq 3)and(swdy[3] eq nc+1)and(swdy[1] eq swdy[2])$
   and(swdy[1] eq 2*swdy[3]-1))  then redoflag = 0
endif
if(redoflag) then begin
wdy = WignerD(-!dpi/2,!dpi/2,0,nc)
write_gdf,wdy,libdir+'WDy.gdf'
endif

end

;*******Effective magnetic dipole calculation*********

Function EMDplanardecomposition,phasemask,scale=scale,kwave=kwave,$
                                beta=beta,kin=kin,Ein=Ein,$
                                prefix=prefix,polar=polar,$
                                fobjective=fobjective
;This function returns planar wave spectrum modulated by the SLM and
;passing through a telescope, see reference [PRTMPT] for detail.
;Ein and kin gives the spectrum goes into the final objective lense,
;they are also saved into files inf keyword prefix is given
;return value has no meaning.

;wave number for 0.532um in vacuum
if (not keyword_set(kwave)) then kwave = 2*!dpi/0.532d 
;magnificatin of telescope, ie the ratio of the two focal length of
;the two lenses that make the telescope
if (not keyword_set(beta)) then beta = 1.d
;the focal lense of the final objective lense, the number here should
;be adjusted according to your working system, defaul value gives
;reasonable intensity profile,one way to determine the value is to
;check the distance between dispaced point traps 
if (not keyword_set(fobjective)) then fobjective = 5400.d
masksize = (size(phasemask))[1:2]
;phasemask resolution , this reflects your SLM point wise distance
if (not keyword_set(scale)) then scale = fobjective/kwave/0.15974999d 
;polar is a 2-elements complex array giving the polarization of
;illuminating beam, ie [1, ci] would be a circular poloarization
if (not keyword_set(polar)) then polar = [1d,0]

ci = dcomplex(0,1)
fftmask = fresnel(phasemask,field=1)


k = kwave
kx = ((dindgen(masksize[0],masksize[1]) mod masksize[0]) - masksize[0]/2.d)
ky = (floor(dindgen(masksize[0],masksize[1])/masksize[0]) - masksize[1]/2.d)

kx = -kx/scale & ky = -ky/scale & rhok = sqrt(kx^2+ky^2)

theta1 = asin(rhok/k) & phi1 = atan(ky,kx)
theta2 = asin(sin(theta1)/beta) & phi2 = phi1+!dpi
Eorigin = [polar[0],polar[1],0]
kx = 0 & ky = 0 & kz = 0 
factor1 = sqrt(cos(theta2)/cos(theta1))*fftmask
E0x = Eorigin[0]/2.*(cos(theta2)+cos(theta1)+(cos(theta2)-cos(theta1))*cos(2*phi2))+$
      Eorigin[1]/2.*(cos(theta2)-cos(theta1))*sin(2*phi2)
E0y = Eorigin[1]/2.*(cos(theta2)+cos(theta1)+(cos(theta1)-cos(theta2))*cos(2*phi2))+$
      Eorigin[0]/2.*(cos(theta2)-cos(theta1))*sin(2*phi2)
E0z = -sin(theta2)*(Eorigin[0]*cos(phi2)+Eorigin[1]*sin(phi2))
E0x = factor1*E0x & E0y = factor1*E0y & E0z = factor1*E0z

E0x = reform(E0x,1,masksize[0]*masksize[1])
E0y = reform(E0y,1,masksize[0]*masksize[1])
E0z = reform(E0z,1,masksize[0]*masksize[1])
theta2 = reform(theta2,1,masksize[0]*masksize[1])
phi2 = reform(phi2,1,masksize[0]*masksize[1])
kin = [sin(theta2)*cos(phi2),sin(theta2)*sin(phi2),cos(theta2)]
Ein = [E0x,E0y,E0z]

if(keyword_set(prefix)) then begin
write_gdf,kin,prefix+'_kin.gdf'
write_gdf,Ein,prefix+'_Ein.gdf'
endif

return,0
end


Function raytracing_all,kin,kout,Ein
;This function will do a ray tracing of light propagating in kin
;direction to kout direction, with electric field Ein
;return value contains the electric field of waves propagating along
;each kout direction, see reference [EW] for detail

kin = double(kin)
kout = double(kout)
mediannorm = tensorproduct(kin,kout) 
threshold = double((mediannorm[0,*]^2+mediannorm[1,*]^2+$
                    mediannorm[2,*]^2) lt 1e-5)
mediannorm = mediannorm+[threshold,threshold,threshold]
mediannorm = unitvector(mediannorm)
rayoutnorm = tensorproduct(kout,mediannorm) 
rayoutnorm = unitvector(rayoutnorm)
rayinnorm = tensorproduct(kin,mediannorm) 
rayinnorm = unitvector(rayinnorm)
medianprojection = scalorproduct(mediannorm,Ein)
rayinprojection = scalorproduct(rayinnorm,Ein)
Eoutt = scalorvectorproduct(rayinprojection,rayoutnorm) +$ 
       scalorvectorproduct(medianprojection,mediannorm)
Eout = scalorvectorproduct(threshold,Ein)+$
       scalorvectorproduct((1-threshold),Eoutt)

return, Eout
end

Function EMDphasemasktolightinfo,phasemask,scale=scale,$
                                 lambda=lambda,beta=beta,$
                                 prefix=prefix,loud=loud,$
                                 bblock=bblock,err=err,$
                                 coneangobj=coneangobj,$
                                 linfodim=linfodim,$
                                 fobjective=fobjective,$
                                 minfo=minfo,polar=polar,$
                                 threashold=threashold
;This function calculate the planar wave spectrum coming out of the
;final objective lense. See referece [BSYRDG] for detail


if (not keyword_set(lambda)) then lambda =0.532d ; green
kwave = 2*!dpi/lambda
;resolution of planar wave spectrum coming out of the final objective
;lense
if (not keyword_set(linfodim)) then linfodim = [100,200]
;refractive index of medium, default for water
if (not keyword_set(minfo)) then minfo = 1.33d
; convergent angle of objective, this has big effect with the z
; trapping
if (not keyword_set(coneangobj)) then coneangobj =!dpi/3.d
; used to filter out noise high frequency modes from phase mask,
; setting non-zero value can substantially incread the speed. The
; safe value should be set by first do fm = fresnel(phasemask), then
; choose a reasonable noise threashold. Default value works fine to
; point, line, ring traps and vortex 
if (not keyword_set(threashold)) then threashold = 1e-5 

lightinfo = dcomplexarr(2,linfodim[0],linfodim[1])
err = dblarr(linfodim[0],linfodim[1])
ci = dcomplex(0,1)

;we are assuming SLM put in vacuum
nothing = EMDplanardecomposition(phasemask,scale=scale,beta=beta,$
          kwave=kwave,kin=kin,Ein=Ein,prefix=prefix,polar=polar,$
          fobjective = fobjective)
                                 
eein = abs(scalorproduct(ein,conj(ein)))

wwgood = where(eein[*] gt threashold)
ein=ein[*,wwgood]
kin=kin[*,wwgood]
thetakin = acos(kin[2,*])
; if a beam block is set to cut central lights
if (keyword_set(bblock)) then begin
wgood = where(thetakin gt bblock)
kin = kin[*,wgood] & Ein = Ein[*,wgood]
thetakin = acos(kin[2,*])
endif

dtheta = (!dpi/2.)/linfodim[0] & dphi = (2.*!dpi/linfodim[1])
ncpt = (size(kin))[2]*((size(kin))[0] eq 2) + ((size(kin))[0] eq 1)
Eout = dcomplexarr(3,ncpt)
phikin = atan(kin[1,*],kin[0,*])

origink = [fobjective*tan(thetakin)*cos(phikin),$
           fobjective*tan(thetakin)*sin(phikin),0*phikin]
for th = 0,floor(coneangobj/dtheta)-1 do begin
theta = (th+0.5d)*dtheta                             
for ph = 0,linfodim[1]-1 do begin
if (keyword_set(loud)) then $
print,'EMD lightinfo calculation ',th,ph,$
       ' of ', floor(coneangobj/dtheta)-1,linfodim[1]
phi = (ph+0.5)*dphi

kout = [sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta)]

koutt = [reform(replicate(kout[0],ncpt),1,ncpt),$
         reform(replicate(kout[1],ncpt),1,ncpt),$
         reform(replicate(kout[2],ncpt),1,ncpt)]
phaseshift = exp(-ci*kwave*minfo*$
  (origink[0,*]*kout[0]+origink[1,*]*kout[1]));+origink[2,*]*kout[2]))
temp1 = scalorproduct(kin,koutt)/minfo
Eoutt = raytracing_all(kin,koutt,Ein)
temp2 = sqrt(temp1*double(temp1 gt 0))*phaseshift

Eout = scalorvectorproduct(temp2,Eoutt)

if (ncpt gt 1) then Edecompose = total(Eout,2)
if (ncpt le 1) then Edecompose = Eout

normfactor = sin(th*dtheta)*dtheta*dphi


E1 = oldvector([1,0,0],theta,phi)
E2 = oldvector([0,1,0],theta,phi)
lightinfo[0,th,ph] = normfactor*scalorproduct(Edecompose,E1)
lightinfo[1,th,ph] = normfactor*scalorproduct(Edecompose,E2)

endfor

endfor

return,lightinfo
end

Function LinfotoIntensity,linfo,xrange,yrange,zrange,$
         field=field,lambda=lambda,loud=loud,minfo=minfo$
         
;This function calculate the field and intensity of a planar wave
;spectrum given by linfo at lattice points defined by xrange, yrange,
;zrange. xrange is a 3-elements array with [xmin, xmax, dx], dx is the
;lattice spacing, and the units are um. The same format for yrange and
;zrange.
;Return value is a 3-d array with the intensity at each point
;If field is set, the calculated electrical field will be passed out

;wave length in vacumm, green light 0.532um as default
if(not(keyword_set(lambda))) then lambda = 0.532d
;refraction index of medium, water as default
if(not(keyword_set(minfo))) then minfo = 1.33d
k = 2*!dpi*minfo/lambda ; units :1/um

dimension=round([(xrange[1]-xrange[0])/xrange[2],$
           (yrange[1]-yrange[0])/yrange[2],$
           (zrange[1]-zrange[0])/zrange[2]])
pos = coorgenxyz(xrange,yrange,zrange,dimension)
dim = (size(linfo))[2:3]
dim = double(dim)
ci = dcomplex(0,1)
npts = n_elements(pos[0,*])
field = dcomplexarr(3,npts)
field[*] = 0

dtheta = !dpi/(2.*dim[0])
dphi = 2*!dpi/dim[1]
for i = 0.d,dim[0]-1 do begin
theta = (i+0.5d)*dtheta
if (keyword_set(loud)) then $
print, 'calculating intensity from lightinfo,', $
        'dealing direction :',theta

for j = 0.d,dim[1]-1 do begin
phi = (j+0.5d)*dphi
if ((abs(linfo[0,i,j])+abs(linfo[1,i,j])) gt 0) then begin
tempf = Oldvector([1,0,0],theta,phi)*linfo[0,i,j]+$
        Oldvector([0,1,0],theta,phi)*linfo[1,i,j]
kx = k*sin(theta)*cos(phi)
ky = k*sin(theta)*sin(phi)
kz = k*cos(theta)
phase = kx*pos[0,*]+ky*pos[1,*]+kz*pos[2,*]
field[0,*] += tempf[0]*exp(ci*phase[0,*])
field[1,*] += tempf[1]*exp(ci*phase[0,*])
field[2,*] += tempf[2]*exp(ci*phase[0,*])
endif
endfor
endfor

intensity = abs(field[0,*])^2+$
            abs(field[1,*])^2+abs(field[2,*])^2
intensity = reform(intensity,dimension[0],$
            dimension[1],dimension[2])
return, intensity
end








;********Mie scattering code*******************


Function Mieexpansionout,ka,np,nm,xiout=xiout,dxiout=dxiout
;This function returns the scattering coefficients from Mie expansion
;method, see ref.[1], chapter 4.
ci = dcomplex(0,1)
m = dcomplex(np/nm)
mx = dcomplex(ka*m)
x = double(ka)
xc = ka + 4.05 * ka^(1./3.) + 8. ; (4.16) convergence condition
nc = floor(xc)+1

a = dcomplexarr(1,nc+1)
b = dcomplexarr(1,nc+1)
D = dcomplexarr(1,nc+1)

;******************* We first calculate An(mx) ***************************************
for ii = 1, nc do begin ;equation (8)
u = ii+0.5d
term = 0 & n = 2 ; start iteration from n = 2
an = 2d*u/mx ; the interation result for n = 1 
fnn = an ;numerator in equation(8)
fdn = 1.d ;denominator in equation(8)
pnn = an ;[an...a1]
pdn = 1e80 ;[an...a2]
jumpnn = 0 ;numerator jump flag, refer to equation (10),(11) and disscusion there in
jumpdn = 0 ;denorminator jump flag, refer to equation (10),(11) and disscusion there in
while (term eq 0) do begin 
an = (-1)^(n+1)*2*(u+n-1)/mx
pnn = an+1.d/pnn ; partial numerator [an...a1]
pdn = an+1.d/pdn ; partial denominator [an...a2]

if (jumpnn eq 1) then begin
pnn = Zn 
fnn = fnn*xsin*pnn ; jump happend last step
endif

if (jumpdn eq 1) then begin
pdn = Zd 
fdn = fdn*xsid*pdn ; jump happend last step
endif
 
if ((abs(pnn) lt 1e-8) and (jumpnn eq 0)) then begin ; partial numeritor is very small, speciall treatment needed
amp1 = (-1.d)^n * 2.d * (u + n) / mx
amp2 = -(-1.d)^n * 2.d * (u + n + 1) / mx
xsin = pnn * amp1 + 1 ;xsi for numerator 
Zn = (amp2*xsin+pnn)/xsin ;Z for numerator, equation(11)
jumpnn = 3 & pnn = 1 ; turn on the flag
endif

if ((abs(pdn) lt 1e-8) and (jumpdn eq 0)) then begin ; partial denorminator is very small, speciall treatment needed
amp1 = (-1.d)^n * 2.d * (u + n) / mx
amp2 = -(-1.d)^n * 2.d * (u + n + 1) / mx
xsid = pdn * amp1 + 1 
Zd = (amp2*xsid+pdn)/xsid
jumpdn = 3 & pdn = 1
endif

if (jumpnn eq 0) then begin ; no jump needed
fnn = fnn*pnn ;numeritor up to order n
endif

if (jumpdn eq 0) then begin ; no jump needed
fdn = fdn*pdn ;numeritor up to order n
endif


term = ((abs(pnn/pdn-1) le 1e-10) and ((jumpdn+jumpnn) eq 0))
n = n+1
jumpnn -= (jumpnn ne 0) & jumpdn -= (jumpdn ne 0)
endwhile

D[ii] = fnn/fdn - ii/mx
endfor

;******************We then calculate Deirmenjian functions psi(x),xsi(x)*************
ASpsi = dcomplexarr(1,nc+1)
ASeta = dcomplexarr(1,nc+1)
dxiout = dcomplexarr(1,nc+1)
xiout = dcomplexarr(1,nc+1)
ASpsi[0] = sin(x) & ASeta[0] = - cos(x) 
xiout[0] = dcomplex(sin(x),-cos(x))
ASxsi = xiout
for ii = 1, nc do begin ;equation (8)
u = ii+0.5d
term = 0 & n = 2 ; start iteration from n = 2
an = 2d*u/x ; the interation result for n = 1 
fnn = an ;numerator in equation(8)
fdn = 1.d ;denominator in equation(8)
pnn = an ;[an...a1]
pdn = 1e80 ;[an...a2]
jumpnn = 0 ;numerator jump flag, refer to equation (10),(11) and disscusion there in
jumpdn = 0 ;denorminator jump flag, refer to equation (10),(11) and disscusion there in
while (term eq 0) do begin 
an = (-1)^(n+1)*2*(u+n-1)/x
pnn = an + 1.d/pnn ; partial numerator [an...a1]
pdn = an + 1.d/pdn ; partial denominator [an...a2]

if (jumpnn eq 1) then begin
pnn = Zn 
fnn = fnn*xsin*pnn ; jump happend last step
endif

if (jumpdn eq 1) then begin
pdn = Zd 
fdn = fdn*xsid*pdn ; jump happend last step
endif

if ((abs(pnn) lt 1e-8) and (jumpnn eq 0)) then begin ; partial numeritor is very small, speciall treatment needed
amp1 = (-1.d)^n * 2.d * (u + n) / x
amp2 = -(-1.d)^n * 2.d * (u + n + 1) / x
xsin = pnn * amp1 + 1 ;xsi for numerator 
Zn = (amp2*xsin+pnn)/xsin ;Z for numerator, equation(11)
jumpnn = 3 & pnn = 1 ; turn on the flag
endif

if ((abs(pdn) lt 1e-8) and (jumpdn eq 0)) then begin ; partial denorminator is very small, speciall treatment needed
amp1 = (-1.d)^n * 2.d * (u + n) / x
amp2 = -(-1.d)^n * 2.d * (u + n + 1) / x
xsid = pdn * amp1 + 1 
Zd = (amp2*xsid+pdn)/xsid
jumpdn = 3 & pdn = 1
endif

if (jumpnn eq 0) then begin ; no jump needed
fnn = fnn*pnn ;numeritor up to order n
endif

if (jumpdn eq 0) then begin ; no jump needed
fdn = fdn*pdn ;numeritor up to order n
endif

term = ((abs(pnn/pdn-1) le 1e-10) and ((jumpdn+jumpnn) eq 0))

n = n+1
jumpnn -= (jumpnn ne 0) & jumpdn -= (jumpdn ne 0)
endwhile

ASpsi[ii] = ASpsi[ii-1]/fnn*fdn
ASeta[ii] = ASeta[ii-1]*fdn/fnn-1/ASpsi[ii-1] ;Abramowitz and Stegun (10.1.31)
ASxsi[ii] = ASpsi[ii]+ci*ASeta[ii]
xiout[ii] = ASxsi[ii]
dxiout[ii] = xiout[ii-1]-ii*xiout[ii]/x
endfor

for n = 1.d, nc do begin
    fac  = D[n]/m + n/x
    a[n] = (fac * ASpsi[n] - ASpsi[n-1])/(fac * ASxsi[n] - ASxsi[n-1])
    fac  = m*D[n] + n/x
    b[n] = (fac * ASpsi[n] - ASpsi[n-1])/(fac * ASxsi[n] - ASxsi[n-1])
endfor

return,[a,b]
end





pro Oldcoeff, newcoeff,al,be,ga,coeff=coeff
;Suppose there is a new coordinate obtained by rotating the old one
;with Euler angles (a, b ,c), and in this new frame we know the
;expansion coefficients in the basis of vectorial spherical
;harmonics. Then this function returns the expansion coefficients in
;terms of spherical harmonics of the original frame. In practice, we
;may have expressed E&M fields as amn*Gmn(new), then we need to know
;what is the expansion interms of a'mn**Gmn(old). Here a, a' are sets
;of expansion coefficients and Gmn are sets of vectorial spherical
;harmonics as of ref.[MTL].
;newcoeff is a 4 x (2*nc+1) x (nc+1) array, an example with max(n) = 1
;is shown below
;              m=-1      m=0      m = 1
;                     a[0,nc,0]                   n=0
;         a[0,nc-1,1] a[0,nc,1]  a[0,nc+1,1]      n=1
;here the first index range from 0 to 3 indicating the coefficients of
; M_mn^1, M_mn^2,N_mn^1, N_mn^2 separately.
;We will take advantage that in the new frame, only a[*,m=1,-1,*] are not
;zero.rotd is the Wigner-D matrix for given euler angles a, b, c
;The basi formula is G_mn(old) = G_in(new) D_im^n, here i =+- 1.
nc = n_elements(newcoeff[0,0,*])-1

coeff = dcomplexarr(4,2*nc+1,nc+1)
;rotation matrix,last index means +- 1 for m, and middle index means 0
;for old, 1 for new, used in the recurrence relations.
D = dcomplexarr(2*nc+1,2,2) 
;generate the first order
D00 = cos(be)
Dp0 = sin(be)/sqrt(2)*dcomplex(-cos(al),sin(al))
Dm0 = sin(be)/sqrt(2)*dcomplex(cos(al),sin(al))
;D^1_{0-}
D[nc,0,0] = -sin(be)/sqrt(2)*dcomplex(cos(ga),sin(ga))
;D^1_{0+}
D[nc,0,1] = sin(be)/sqrt(2)*dcomplex(cos(ga),-sin(ga))
;D^1_{--}
D[nc-1,0,0] = (1+cos(be))/2*dcomplex(cos(al+ga),sin(al+ga))
;D^1_{-+}
D[nc-1,0,1] = (1-cos(be))/2*dcomplex(cos(al-ga),sin(al-ga))
;D^1_{+-}
D[nc+1,0,0] = (1-cos(be))/2*dcomplex(cos(al-ga),sin(ga-al))
;D^1_{++}
D[nc+1,0,1] = (1+cos(be))/2*dcomplex(cos(al+ga),-sin(al+ga))

for i =0,3 do $
coeff[i,*,1] = D[*,0,0]*newcoeff[i,nc-1,1]$
             + D[*,0,1]*newcoeff[i,nc+1,1]
;a = dblarr(2*nc+1,1,2)
;b = a

for n = 2.d, nc-1 do begin

indc = dindgen(2*n+1)-n

D[*,1,*]=0

D[nc+indc,1,0] = sqrt((n+indc)*(n-indc)/(n+1)/(n-1))$
              *D00*D[indc+nc,0,0]$
            + sqrt((n+indc)*(n+indc-1)/2/(n+1)/(n-1))$
              *Dp0*D[indc+nc-1,0,0]$
            + sqrt((n-indc)*(n-indc-1)/2/(n+1)/(n-1))$
              *Dm0*D[indc+nc+1,0,0]

D[nc+indc,1,1] = sqrt((n+indc)*(n-indc)/(n+1)/(n-1))$
              *D00*D[indc+nc,0,1]$
            + sqrt((n+indc)*(n+indc-1)/2/(n+1)/(n-1))$
              *Dp0*D[indc+nc-1,0,1]$
            + sqrt((n-indc)*(n-indc-1)/2/(n+1)/(n-1))$
              *Dm0*D[indc+nc+1,0,1]



coeff[0,*,n] = D[*,1,0]*newcoeff[0,nc-1,n]$
             + D[*,1,1]*newcoeff[0,nc+1,n]
coeff[1,*,n] = D[*,1,0]*newcoeff[1,nc-1,n]$
             + D[*,1,1]*newcoeff[1,nc+1,n]
coeff[2,*,n] = D[*,1,0]*newcoeff[2,nc-1,n]$
             + D[*,1,1]*newcoeff[2,nc+1,n]
coeff[3,*,n] = D[*,1,0]*newcoeff[3,nc-1,n]$
             + D[*,1,1]*newcoeff[3,nc+1,n]
D[*,0,*] = D[*,1,*]
endfor        


end


pro Coeff_Normal,pinfo, linfo,minfo=minfo,$
                 lambda=lambda,coeff=coeff

;This function calculates the scattering coefficients expanded with
;basis as in ref.[2], ie G_mn^y (P=M or N, y=1 or 2)
;We know for normal incident only m=1 and -1 expansion exists,
;and for these VSH, after rotaion by 90 degree, \phi->\phi-!dpi/2,
;so a expansion with polarization in y also easily obtained.

;The return is a 4 x (2*nc+1)*(nc+1) array as for M_mn^1, M_mn^2,
;N_mn^1, N_mn^2 separately
;pinfo   : [radius,refraction index]
;info    : complex amplitude of polarization in x and y, [ax,ay]
;minfo   : refractive index of medium
;lambda  : wave length in vacumm

if(not(keyword_set(minfo))) then minfo = 1.33d
if(not(keyword_set(lambda))) then lambda = 0.532d

ka = dcomplex(2,0)*!dpi*pinfo[0]*minfo/lambda
abcoeff = Mieexpansionout(ka,pinfo[1],minfo)
nc = n_elements(abcoeff[0,*])-1
coeff = dcomplexarr(4,2*nc+1,nc+1)
ci = dcomplex(0,1)
for n = 1d,nc do begin
; E_n/ \gamma_1n, E_n is from ref.[1],page 93, 
;\gamma_1n from ref.[2]
Eg_np = ci^(n mod 4) * sqrt(4*!dpi*(2*n+1))


;x polarization first
coeff[0,nc+1,n] = ci/4.-ci*abcoeff[1,n]/2.;a_1n
coeff[0,nc-1,n] = ci/4.-ci*abcoeff[1,n]/2.;a_{-1,n}
coeff[1,nc+1,n] = ci/4.
coeff[1,nc-1,n] = ci/4.
coeff[2,nc+1,n] = ci/4.-ci*abcoeff[0,n]/2.
coeff[2,nc-1,n] = -ci/4.+ci*abcoeff[0,n]/2.
coeff[3,nc+1,n] = ci/4. 
coeff[3,nc-1,n] = -ci/4.
coeff[*,*,n] *= linfo[0] ;for x polarization

;y polarization contribution
coeff[0,nc+1,n] += linfo[1]*(1/4. -abcoeff[1,n]/2.)
coeff[0,nc-1,n] += linfo[1]*(-1/4.+abcoeff[1,n]/2.)
coeff[1,nc+1,n] += linfo[1]*(1/4.) 
coeff[1,nc-1,n] += linfo[1]*(-1/4.)
coeff[2,nc+1,n] += linfo[1]*(1/4.-abcoeff[0,n]/2.)
coeff[2,nc-1,n] += linfo[1]*(1/4.-abcoeff[0,n]/2.)
coeff[3,nc+1,n] += linfo[1]*(1/4.)
coeff[3,nc-1,n] += linfo[1]*(1/4.)

coeff[*,nc+1,n] *= Eg_np
coeff[*,nc-1,n] *= Eg_np

endfor

end

;****************Force and Torque code********************************

Pro  lightinfotocoeff,ppos,pinfo,linfo,outdir=outdir,$
                         namestart=namestart,mpp=mpp,$
                 loud=loud,minfo=minfo,lambda=lambda,$
                                  smallfile=smallfile
;This function will calculate the total expansion coefficients for a
;given spectrum of incoiming light, for each position ppos[*,i], a
;file is created to store the coefficients. The files are in the
;directory 'outdir', with  a name like posxxxxx.gdf, the numbers in
;filenames start from namestart
;since coeff[*,*,0] is empty, we use coeff[0:2,*,0] to store the
;position of particle.

if(not(keyword_set(minfo))) then minfo = 1.33d
if(not(keyword_set(mpp))) then mpp = 0.135d ;um/pix
if(not(keyword_set(lambda))) then lambda = 0.532d
if(not(keyword_set(outdir))) then outdir = './'
if(not(keyword_set(namestart))) then namestart = 0
k = dcomplex(2,0)*mpp*!dpi*minfo/lambda  ;unit:1/pix
outdir = fixdir(outdir)

;directory to store scattering coefficients of each position that
;particle sits 
dirtest = file_info(outdir)
if(not(dirtest.directory)) then begin
command = 'mkdir '+outdir &spawn,command
endif

dim = double((size(linfo))[2:3])
ci = dcomplex(0,1)
npts = n_elements(ppos[0,*])

;get the expansion for normal incident
Coeff_normal,pinfo,[1,0],minfo=minfo,$
             lambda=lambda,coeff=ncoeffx
Coeff_normal,pinfo,[0,1],minfo=minfo,$
             lambda=lambda,coeff=ncoeffy

nc = n_elements(ncoeffx[0,0,*])-1

expancoef = dcomplexarr(npts,4,2*nc+1,nc+1)
tempcoef = dcomplexarr(4,2*nc+1,nc+1)


for i = 0.d,dim[0]-1 do begin
theta = (i+0.5d)*!dpi/(2.d*dim[0])
if(keyword_set(loud)) then begin
print,i*100/(dim[0]-1),' % finished for ',npts,' points'
endif
for j = 0.d,dim[1]-1 do begin
phi = (j+0.5d)*2.d*!dpi/dim[1]
if ((abs(linfo[0,i,j])+abs(linfo[1,i,j])) gt 0) then begin
kx = sin(theta)*cos(phi)
ky = sin(theta)*sin(phi)
kz = cos(theta)

Oldcoeff,ncoeffx,phi,theta,-phi,coeff=oldcoeffx
Oldcoeff,ncoeffy,phi,theta,-phi,coeff=oldcoeffy

for pt = 0l,npts-1 do begin

phasept = exp(ci*k*(kx*ppos[0,pt]+$
        ky*ppos[1,pt]+kz*ppos[2,pt]))
weight1 = linfo[0,i,j]*phasept
weight2 = linfo[1,i,j]*phasept
expancoef[pt,*,*,*] += oldcoeffx*weight1+oldcoeffy*weight2

endfor
endif
endfor
endfor

if(keyword_set(smallfile)) then begin
for pt = 0l,npts-1 do begin
fname = outdir+'pos'+string(pt+namestart,format='(I06)')$
        +'.gdf'
tempcoef[*,*,*] = expancoef[pt,*,*,*]
tempcoef[0,0:2,0] = ppos[0:2,pt]*mpp
tempcoef[1,0:2,0] = [pinfo,minfo]
tempcoef[2,0:2,0] = [dim,lambda]
write_gdf,tempcoef,fname
endfor
endif else begin
for pt = 0l,npts-1 do begin
expancoef[pt,0,0:2,0] = ppos[0:2,pt]*mpp
expancoef[pt,1,0:2,0] = [pinfo,minfo]
expancoef[pt,2,0:2,0] = [dim,lambda]
endfor
fname = outdir+'pos'+string(namestart,format='(I06)')+'.gdf'
write_gdf,expancoef,fname
endelse
end



Function Outfield_multipoints,ppos,pinfo,linfo,minfo=minfo,$
                              outdir=outdir,loud=loud,step=step,$
                              nameorder=nameorder,lambda=lambda,$
                              mpp=mpp,smallfile=smallfile
;This is a interface to call lightinfotocoeff, given the points you
;would like to calculate (ppos), particle parameters (pinfo), angular
;spectrum of incident light (linfo). A series of posxxxxx.gdf will be
;generated under outdir.
;step is the number of points in one batch
;return value is the total time used for the calculation.
if(not(keyword_set(nameorder))) then nameorder = 0l
if(not(keyword_set(step))) then step = 100
if(not(keyword_set(outdir))) then outdir = './data/'
starttime = systime(1)
npts = n_elements(ppos[0,*])
step = long(step)
nameorder = long(nameorder)
for i = 0l,floor(double(npts-1)/double(step)) do begin
w1 = i*step & w2 = i*step+step-1
if (w2 ge npts) then w2=npts-1

if (keyword_set(loud)) then begin
print, 'points num ',w1+1,' to ',w2+1,' in total of ',npts,' points'
endif

;lightinfotocoeff,ppos[*,w1:w2],pinfo,linfo,outdir=outdir,$
;                 namestart=w1+nameorder,mpp=mpp,$
;                 loud=loud,minfo=minfo,lambda=lambda
namestart = i+nameorder ;defaul using big file to save time
if (keyword_set(smallfile)) then namestart=w1+nameorder
lightinfotocoeff,ppos[*,w1:w2],pinfo,linfo,outdir=outdir,$
                 namestart=namestart,mpp=mpp,$
                 loud=loud,minfo=minfo,lambda=lambda,$
                 smallfile=smallfile

endfor
nameorder = nameorder+i+1l
endtime = systime(1)
return, [double(endtime-starttime),npts]
end


Function forcez,coeff,lambda=lambda

if (not(keyword_set(lambda))) then lambda = 0.532d ;in um
k = 2*!dpi/lambda
nc = n_elements(coeff[0,0,*])-1
result = 0d
ci = dcomplex(0,1)

for n = 1d,nc-1 do begin
m = dindgen(1,2*n+1)-n
a1 = m/(n*(n+1))*imaginary(conj(coeff[0,nc+m,n])*(-ci)*coeff[2,nc+m,n])
b1 = 1./(n+1)*$
    sqrt(n*(n+2)*(n+1-m)*(n+m+1)/(2*n+1)/(2*n+3))*$
    imaginary($
    coeff[0,nc+m,n]*conj(coeff[0,nc+m,n+1])$
    +coeff[2,nc+m,n]*conj(coeff[2,nc+m,n+1]))
a2 = m/(n*(n+1))*imaginary(conj(coeff[1,nc+m,n])*(ci)*coeff[3,nc+m,n])
b2 = 1./(n+1)*$
    sqrt(n*(n+2)*(n+1-m)*(n+m+1)/(2*n+1)/(2*n+3))*$
    imaginary($
    -coeff[1,nc+m,n]*conj(coeff[1,nc+m,n+1])$
    -coeff[3,nc+m,n]*conj(coeff[3,nc+m,n+1]))

result += total(a1+b1+a2+b2)

endfor

return, 2*result/(k^2)
end

Function forcexy,coeff,wd,lambda=lambda
;This function will rotate the coefficients using the matrix wd first,
;then calculate z force using the new coefficients
newcoeff = coeff
nc = n_elements(coeff[0,0,*])-1
for i = 1,nc do begin
newcoeff[0,*,i] = wd[*,*,i]##coeff[0,*,i]
newcoeff[1,*,i] = wd[*,*,i]##coeff[1,*,i]
newcoeff[2,*,i] = wd[*,*,i]##coeff[2,*,i]
newcoeff[3,*,i] = wd[*,*,i]##coeff[3,*,i]
endfor
return,forcez(newcoeff,lambda=lambda)
end




Function CoefftoForce, datadir = datadir,lambda=lambda,$
                       libdir = libdir,coeff=coeff,pos=pos
;datadir gives the directory where previously calculated coefficients
;can be found
;libdir gives the directory where previously calculated Wigner-D
;matrice can be found, it has to contain two gdf files
;WDx.gdf is the Wigner-D matrix with Euler angles (\alpha = 0,\beta =
;\pi/2, \gamma = 0) which will be used to calculate force along x
;direction. 
;WDy.gdf is the Wigner-D matrix with Euler angles (\alpha =
;-\pi/2,\beta = \pi/2, \gamma = 0) which will be used to calculate force along y 
;direction.
;if either of the two are missing, new files will be generated
;if keyword coeff is set, the force corresponding to coeff will be
;returned.
;pos returns the 3-d coordinate of the particle corresponding of the
;returned force field
if(not(keyword_set(libdir))) then libdir ='./'

if(keyword_set(coeff)) then begin
force = dblarr(3)
pos = dblarr(3)
nc = n_elements(coeff[0,0,*])-1
FixWd,libdir,nc,wdx=wdx,wdy=wdy
force[2] = forcez(coeff,lambda=lambda)
force[0] = forcexy(coeff,wdx,lambda=lambda)
force[1] = forcexy(coeff,wdy,lambda=lambda)
pos[*] = coeff[0,0:2,0]
return, force
endif

if(not(keyword_set(datadir))) then datadir = './data/'
datadir = fixdir(datadir)
f = file_search(datadir+'pos*.gdf',count=nfiles)
force = dblarr(3)
nforce = force
pos = dblarr(3)
npos = pos
if(nfiles ge 1) then begin
for nf = 0, nfiles-1 do begin
allcoeff = read_gdf(f[nf])
scoeff = size(allcoeff)
;this is a multipoint coefficient-set file
if(scoeff[0] eq 4) then begin
npts = n_elements(allcoeff[*,0,0,0])
for pt = 0l,npts-1 do begin
coeff = allcoeff[pt,*,*,*]
nc = n_elements(coeff[0,0,0,*])-1
coeff = reform(coeff,4,2*nc+1,nc+1)
FixWd,libdir,nc,wdx=wdx,wdy=wdy
force[2] = forcez(coeff,lambda=lambda)
force[0] = forcexy(coeff,wdx,lambda=lambda)
force[1] = forcexy(coeff,wdy,lambda=lambda)
nforce = [[nforce],[force]]
pos[*] = allcoeff[pt,0,0:2,0]
npos =[[npos],[pos]]
endfor
endif 
;this file contains a single set of coefficients
if(scoeff[0] eq 3) then begin
nc = n_elements(allcoeff[0,0,*])-1
FixWd,libdir,nc,wdx=wdx,wdy=wdy
force[2] = forcez(allcoeff,lambda=lambda)
force[0] = forcexy(allcoeff,wdx,lambda=lambda)
force[1] = forcexy(allcoeff,wdy,lambda=lambda)
nforce = [[nforce],[force]]
pos[*] = allcoeff[0,0:2,0]
npos = [[npos],[pos]]
endif
endfor
endif
force = nforce[*,1:*]
pos = npos[*,1:*]
return, force

end

;This function returns torque in z direction according to coeff

Function torquez,coeff,lambda=lambda

if (not(keyword_set(lambda))) then lambda = 0.532d ;in um
k = 2*!dpi/lambda
nc = n_elements(coeff[0,0,*])-1
result = 0d
ci = dcomplex(0,1)

for n = 1d,nc-1 do begin
m = dindgen(1,2*n+1)-n
a1 = m*abs(coeff[0,m+nc,n])^2
a2 = m*abs(coeff[2,m+nc,n])^2
b1 = m*abs(coeff[1,m+nc,n])^2
b2 = m*abs(coeff[3,m+nc,n])^2
result += total(a1+a2-b1-b2)
endfor
return, -result/k^2
end


Function torquexy,coeff,wd,lambda=lambda
;This function will rotate the coefficients using the matrix wd first,
;then calculate z torque using the new coefficients
newcoeff = coeff
nc = n_elements(coeff[0,0,*])-1
for i = 1,nc do begin
newcoeff[0,*,i] = wd[*,*,i]##coeff[0,*,i]
newcoeff[1,*,i] = wd[*,*,i]##coeff[1,*,i]
newcoeff[2,*,i] = wd[*,*,i]##coeff[2,*,i]
newcoeff[3,*,i] = wd[*,*,i]##coeff[3,*,i]
endfor
return,torquez(newcoeff,lambda=lambda)
end


Function CoefftoTorque, datadir = datadir,lambda=lambda,$
                       libdir = libdir,coeff=coeff,pos=pos
;datadir gives the directory where previously calculated coefficients
;can be found
;libdir gives the directory where previously calculated Wigner-D
;matrice can be found, it has to contain two gdf files
;WDx.gdf is the Wigner-D matrix with Euler angles (\alpha = 0,\beta =
;\pi/2, \gamma = 0) which will be used to calculate torque along x
;direction. 
;WDy.gdf is the Wigner-D matrix with Euler angles (\alpha =
;-\pi/2,\beta = \pi/2, \gamma = 0) which will be used to calculate torque along y 
;direction.
;if either of the two are missing, new files will be generated
;if keyword coeff is set, the torque corresponding to coeff will be
;returned.
;pos returns the 3-d coordinate of the particle corresponding of the
;returned torque field
if(not(keyword_set(libdir))) then libdir ='./'

if(keyword_set(coeff)) then begin
torque = dblarr(3)
pos = dblarr(3)
nc = n_elements(coeff[0,0,*])-1
FixWd,libdir,nc,wdx=wdx,wdy=wdy
torque[2] = torquez(coeff,lambda=lambda)
torque[0] = torquexy(coeff,wdx,lambda=lambda)
torque[1] = torquexy(coeff,wdy,lambda=lambda)
pos[*] = coeff[0,0:2,0]
return, torque
endif

if(not(keyword_set(datadir))) then datadir = './data/'
datadir = fixdir(datadir)
f = file_search(datadir+'pos*.gdf',count=nfiles)
torque = dblarr(3)
ntorque = torque
pos = dblarr(3)
npos = pos
if(nfiles ge 1) then begin
for nf = 0, nfiles-1 do begin
allcoeff = read_gdf(f[nf])
scoeff = size(allcoeff)
;this is a multipoint coefficient-set file
if(scoeff[0] eq 4) then begin
npts = n_elements(allcoeff[*,0,0,0])
for pt = 0l,npts-1 do begin
coeff = allcoeff[pt,*,*,*]
nc = n_elements(coeff[0,0,0,*])-1
coeff = reform(coeff,4,2*nc+1,nc+1)
FixWd,libdir,nc,wdx=wdx,wdy=wdy
torque[2] = torquez(coeff,lambda=lambda)
torque[0] = torquexy(coeff,wdx,lambda=lambda)
torque[1] = torquexy(coeff,wdy,lambda=lambda)
torque = [[ntorque],[torque]]
pos[*] = allcoeff[pt,0,0:2,0]
npos =[[npos],[pos]]
endfor
endif 
;this file contains a single set of coefficients
if(scoeff[0] eq 3) then begin
nc = n_elements(allcoeff[0,0,*])-1
FixWd,libdir,nc,wdx=wdx,wdy=wdy
torque[2] = torquez(allcoeff,lambda=lambda)
torque[0] = torquexy(allcoeff,wdx,lambda=lambda)
torque[1] = torquexy(allcoeff,wdy,lambda=lambda)
ntorque = [[ntorque],[torque]]
pos[*] = allcoeff[0,0:2,0]
npos = [[npos],[pos]]
endif
endfor
endif
torce = ntorque[*,1:*]
pos = npos[*,1:*]
return, torque

end



;If you want to normalize your force and torque, devide your answer by
;the power.
Function power,coeff

s1 = total(coeff[1,*,1:*]*conj(coeff[1,*,1:*]))
s2 = total(coeff[3,*,1:*]*conj(coeff[3,*,1:*]))
return,real_part(s1+s2)
end


pro vhot
print,'Vectorial Calculation for Holographic Optical Tweezer 2.0'
end