Programs relevant to xmpl3d27, Simplest unbalanced deflector, for integrating aberration coefficients or for using previously calculated coefficients to compare with direct ray tracing.
These can be copied from the example file, which is supplied with the CPO3D package.
(in Fortran)
C Calculate or test formulas for aberration coefficients for an unbalanced
C deflector.
C Either use axial values output by CPO3D with dflab.dat,
C in the output data tmp3a.dat, to integrate coefficients.
C Or use final ray output parameters to test coefficients.
C Call dflabbat, then it gives data in temp.dat.
DIMENSION z(-500:500) , f1(-500:500) , p(-500:500) ,
& pd(-500:500) , f2(-500:500) , f3(-500:500) , tn(7) ,
& picht(-500:500) , R(-500:500) , ra(-500:500) ,
& rad(-500:500) , radd(-500:500) , sqi(-500:500) ,
& pdp(-500:500) , pdd(-500:500) , f1d(-500:500) ,
& f1dd(-500:500)
LOGICAL integrate_coeffs , weak_field , adjust_D
OPEN (UNIT=1,FILE='tempdat.dat',STATUS='OLD')
C =trajectory output file
OPEN (UNIT=2,FILE='temp.dat',STATUS='UNKNOWN')
C =results file
integrate_coeffs = .TRUE.
c integrate_coeffs = .FALSE.
weak_field = .TRUE.
c weak_field = .FALSE.
adjust_D = .FALSE.
adjust_D = .TRUE.
IF (integrate_coeffs) THEN
dv = 0.1
C = dv of deflector
n = 200
C = number of z values, positive, excluding z = 0 (reflected for negative)
z2 = 5.0
C = z final, for axial field information
dx = 0.025
dy = dx
C = x and y for finding f2,f3
ELSE
dv = 0.1
C = dv of deflector
energy = 20.
C = energy
dv = dv/energy
nrays = 27
ENDIF
z_spec = 8.0
C = z of special plane
C coeffs to be used if they are not recalculated:
C obtained with integration of principal ray (ie not weak field):
D = 2.63440E+00
A = -2.46621E+01
C1= -2.14419E+00
C2= 2.94321E+02
C weak field:
D = 2.64332E+00
A = -2.49471E+01
C1= -1.98914E+00
C2= 3.01404E+02
C
IF (.not.integrate_coeffs) GOTO 100
C -ie do not integrate aberration coefficients
z1 = -z2
C = z initial for axial field information
n2 = n-2
n4 = n-4
dz = z2/real(n)
C read info for along axis:
C potential:
DO i = 0 , n
z(i) = dz*real(i)
z(-i) = -z(i)
read(1,*)d1 , d2 , d3 , p(i)
C d1,2,3 = dummies
p(i) = p(i) + 1.
C -ie constant potential of 1 added
p(-i) = p(i)
ENDDO
C fields:
DO i = 0 , n
read(1,*)d1 , d2 , d3 , ex , d4 , ez
C d1,2,3,4 = dummies
f1(i) = -ex
f1(-i) = f1(i)
pd(i) = -ez
pd(-i) = -pd(i)
ENDDO
IF (.NOT.weak_field) THEN
C Differentiate:
DO i = -n2 , n2
f1d(i) = (8.*(f1(i+1)-f1(i-1))-(f1(i+2)-f1(i-2)))/(12.*dz)
pdd(i) = (8.*(pd(i+1)-pd(i-1))-(pd(i+2)-pd(i-2)))/(12.*dz)
ENDDO
DO i = -n4 , n4
f1dd(i)
& = (8.*(f1d(i+1)-f1d(i-1))-(f1d(i+2)-f1d(i-2)))/(12.*dz)
ENDDO
ENDIF
C read info for along line at (x,y) = (0,dy):
y = dy
x = 0.
DO i = 0 , n
read(1,*) d1 , d2 , d3 , ex , ey , d4
f2(i) = ey/(2.*dy)
f3(i) = (ex+f1(i))/(3.*dy**2)
IF ( .NOT.weak_field .AND. (i.LE.n4) ) THEN
f2(i) = f2(i) - pdd(i)/4.
f3(i) = f3(i) - f1dd(i)/24.
C (the integration of the double differentials gives zero for the
C weak field situation)
ENDIF
f2(-i) = f2(i)
f3(-i) = f3(i)
ENDDO
c
IF (weak_field) THEN
DO i = -n , n
ra(i) = z_spec - z(i)
C = ray that points to z = z_spec
rad(i) = -1.
C with slope -1
radd(i) = 0.
ENDDO
nn = n
C =limiting n
ELSE
C Integrating Picht’s equation:
tn(1) = 5./720.
tn(2) = -33./720.
tn(3) = 93./720.
tn(4) = -142./720.
tn(5) = 63./720.
tn(6) = -657./720.
tn(7) = 49./720.
DO i = -n , n
tsq = (pd(i)/p(i))**2
C = (phidash/phi)**2=T**2
picht(i) = tsq*3.
ENDDO
C Integrate principal ray through lens,
C ray starts at z = -z_spec, distance 2*z_spec from axis, slope -1, so that
C it points to z = z_spec
C It is ray that is converged by a perfect lens.
C Use my own algorith on Picht’s equation, an adaptation of Numerov’s (see
C Heddle’s book), taken from cpo2d1.for
C Use non-relativistic form of Picht’s equation (see Lencova and Lenc):
C d2R/dtsq = (3/16)*tsq*R
C Transformation used is R = r*phi_star**0.25
sqrtps1 = sqrt(p(-n))
trans = sqrt(sqrtps1)
R(-n) = (z_spec+z2)*trans
C =R at z = z1 (= -z2)
R(-n+1) = R(-n) - dz*trans
C -because slope is -1
dzsq = dz**2
a = 0.0625*dzsq/12.
DO i = (-n+2) , (-n+5)
R(i) = (2.*R(i-1)-R(i-2) - a*(picht(i-2)*R(i-2)
& + 10.*picht(i-1)*R(i-1)))/(1.+a*picht(i))
C =Numerov’s algorithm for 3rd points
ENDDO
a = 0.0625*dzsq
DO i = (-n+6) , n
R(i) = 2.*R(i-1)-R(i-2)
DO it = 1 , 6
R(i) = R(i) + a*tn(7-it)*picht(i-it)*R(i-it)
ENDDO
R(i) = R(i)/(1.+a*tn(7)*picht(i))
C = my algorithm for later points
ENDDO
C
DO i = -n , n
trans = p(i)**(-0.25)
ra(i) = R(i)*trans
ENDDO
DO i = -n2 , n2
rad(i) = (8.*(ra(i+1)-ra(i-1))-(ra(i+2)-ra(i-2)))/(12.*dz)
ENDDO
DO i = -n4 , n4
radd(i) = (8.*(rad(i+1)-rad(i-1))-(rad(i+2)-rad(i-2)))
& /(12.*dz)
ENDDO
nn = n4
c = limiting n
C Find final value and slope of r:
trans = p(n)**(-0.25)
endr = R(n)*trans
endr1 = R(n-1)*trans
cc write(2,'(''zi,r,rd= '',1p3e12.4)')z(n) , endr , endr1
ENDIF
C
C integrate aberration coeffs:
D = 0.
A = 0.
C1 = 0.
C2 = 0.
DO i = -nn , nn
C nn = n for weak_field, n4 otherwise
sqi(i) = 1./sqrt(p(i))
pdp(i) = pd(i)/p(i)
ENDDO
v1 = 3./64.
v2 = 3./8.
v3 = 5./16.
v4 = 1./8.
DO i = -nn , nn
u = f1(i)*sqi(i)*ra(i)
D = D + 0.5*u
A = A + f2(i)*sqi(i)*ra(i)**2
C2 = C2 + 3.*f3(i)*sqi(i)*ra(i)**3
C1 = C1 + u*(v1*(pdp(i)*ra(i))**2 + v2*rad(i)**2
& + v3*ra(i)*radd(i) - v4*pdp(i)*ra(i)*rad(i))
ENDDO
v1 = dz*(sqi(nn)/rad(nn))/dv
C (have to use values at nn because they are not known at n, but should
C be the same)
v2 = v1/rad(nn)
v3 = v2/rad(nn)
D = -v1*D
A = -v2*A
C2 = -v3*C2
C1 = v3*C1
write(2,'('' D ='',1pe13.5/
& '' A ='',e13.5/
& '' C1='',e13.5/
& '' C2='',e13.5)')D,A,C1,C2
stop
100 continue
C =jump point for not integrating coefficients
C
C Analyse ray data:
write(2,'(''Using pre-calculated D etc:''/1p4e11.3)')
& D,A,C1,C2
D = D*dv
A = A*dv
C1 = C1*dv
C2 = C2*dv
DO i = 1 , nrays
read(1,*)x1 , y1 , z1 , vx1 , vy1 , vz1 , en , x2 , y2 , z2 ,
& vx2 , vy2 , vz2
C = initial and final ray parameters
xd2 = vx2/vz2
yd2 = vy2/vz2
C =final slopes
C Use coeffs to calc expected x2 and y2:
x2_calc = D + A*xd2 + C1*(3.*xd2**2 + yd2**2)
& + C2*0.5*(xd2**2 - yd2**2)
y2_calc = -A*yd2 + C1*2.*xd2*yd2 - C2*xd2*yd2
write(2,'(''x2,y2,x1,y1= '',1p4e11.3
& /''x2calc,y2calc= '',2e11.3)')
& x2 , y2 , x1 , y1 , x2_calc , y2_calc
IF ( adjust_D .AND. (i.EQ.1) ) THEN
C First ray is the principal ray
C Adjust D, using the principal ray
D = D + x2 - x2_calc
dd = D/dv
write(2,'(''Adjusted D= ''/1pe11.3)') dd
ENDIF
ENDDO
STOP
END