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) ,

& 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.

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

C with slope -1

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

ENDDO

DO i = -n4 , n4

& /(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

ENDDO

C (have to use values at nn because they are not known at n, but should

C be the same)

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