Program for generating a 2D array of magnetic grid points.

The program is called prog22.f90 and is supplied with the CPO3D package together with a compiled version that can be activated by the command ‘prog22’. The output data will be put into the file ‘tempout.dat’.

The program uses, as an example, an axial magnetic field Bz given by the Glaser formula:

Bz = B0/(1 + (z/w)**2)

with B0 = 0.01 Tesla, w = 10mm. The focal length of the lens is then given by

f = w/sin(pi/sqrt(1+K))

where

K = e*(w*B0)**2/(8*m*V)

and where eV is the electron energy, taken here as 1000eV. Therefore f = 34.153mm.

The '2D grid' option is used here for the magnetic field. The grid of field components is held in the file tempout.dat, reproduced below in the file xmpl3d46.dat.

The triggering letter used in xmpl3d46.dat is 'y', for an incomplete array of 2D grid points. In fact the grid points in tempout.dat are complete and ordered, but the results are unchanged if unnecessary points are removed and/or the order of the points is changed (in any way).

Alternatively, the output file can be generated by using a modification of the program given as appendix 2 of cpo3d.dat. The axial field is calculated using the above formula and then the program is used to generate the off-axis components.

Other methods could also be used -for example a user-produced program for the field.

Program used to generate the data file (in Fortran 95):

! Program for generating an input file for use in xmpl3d46.dat,

! giving a grid of magnetic field values.

! An axial magnetic field Bz given by the Glaser formula:

! Bz = B0/(1 + (z/w)**2)

! with B0 = 0.01 Tesla, w = 10mm.

! The 2D grid is relevant for an cylindrically symmetrically system, in the

! rho,z plane.

DIMENSION b(5000) , db(5000) , ddb(5000) , dddb(5000) ,

& ddddb(5000)

! -for a maximum of 5000 axial points

LOGICAL cyl

! =.true. for 2D cylindrically symmetric system, when required output is

! a 2D array of points

cyl = .TRUE.

! cyl = .FALSE.

dz = 0.5 !=spacing of grid points

! Set coordinates maxima and minima:

IF (cyl) THEN

rmax = 1.0

xmax = rmax !(which substitutes for rmax below)

nr = (xmax/dz + 1.E-6) + 1 !=number of radial grid points

ELSE

xmin = -1.0

xmax = 1.0

ymin = -1.0

ymax = 1.0

ENDIF

zmin = -64.

zmax = 64.

zmin = zmin - 12*dz

zmax = zmax + 12*dz

! -to restore ends that will be removed below.

! Find integer min and max coordinates:

IF (xmin.LT.0.) xmin = xmin - dz

nxmin = (xmin/dz + 0.0001)

IF (xmax.LT.0.) xmax = xmax - dz

nxmax = (xmax/dz + 0.0001)

IF (ymin.LT.0.) ymin = ymin - dz

nymin = (ymin/dz + 0.0001)

IF (ymax.LT.0.) ymax = ymax - dz

nymax = (ymax/dz + 0.0001)

nz = (zmax - zmin + 1.E-10)/dz + 1

diff7 = 1./(60.*dz) ! -needed for 7-point differentiation

OPEN (UNIT=2,FILE='tempout.dat',STATUS='UNKNOWN')

! =results file, with name and format for use by cpo3d

n1 = 1

! Axial values:

DO j = 1 , nz

z = dz*real(j-1) + zmin

b(j) = 0.01/(1. + (z/10.)**2)

! =Glaser formula

ENDDO

n1 = 1

n2 = nz

! =present number of points in z direction

! 1st differential:

n1 = n1 + 3

n2 = n2 - 3

! -Cannot include 3 points at each end, which are used for

! differentiation

DO i = n1 , n2

db(i) = ( b(i+3)-b(i-3)

& - 9.*(b(i+2)-b(i-2))

& + 45.*(b(i+1)-b(i-1)) )*diff7

! =1st differential

ENDDO

! 2nd differential:

n1 = n1 + 3

n2 = n2 - 3

! -Cannot include 3 points at each end

DO i = n1 , n2

ddb(i) = ( db(i+3)-db(i-3)

& - 9.*(db(i+2)-db(i-2))

& + 45.*(db(i+1)-db(i-1)) )*diff7

! =2nd differential

ENDDO

! 3rd differential:

n1 = n1 + 3

n2 = n2 - 3

DO i = n1 , n2

dddb(i) = ( ddb(i+3)-ddb(i-3)

& - 9.*(ddb(i+2)-ddb(i-2))

& + 45.*(ddb(i+1)-ddb(i-1)) )*diff7

ENDDO

! 4th differential:

n1 = n1 + 3

n2 = n2 - 3

DO i = n1 , n2

ddddb(i) = ( dddb(i+3)-dddb(i-3)

& - 9.*(dddb(i+2)-dddb(i-2))

& + 45.*(dddb(i+1)-dddb(i-1)) )*diff7

ENDDO

!

nz = n2 - n1 + 1 !=new number of points in z direction

zmaxx = zmin + real(n2 - 1)*dz

zminn = zmin + real(n1 - 1)*dz

! =z limits with differentiation points excluded

IF (cyl) THEN

! Cylindrically symmetric system, 2D

rmin = 0.

WRITE (2,'(2F8.3)') rmin , rmax , zminn , zmaxx

nn = (nxmax + 1)*nz

WRITE (2,'(I8, F8.3)') nn , dz !=total number of grid points, spacing

ELSE

! 3D system

xmin = real(nxmin)*dz

xmax = real(nxmax)*dz

ymin = real(nymin)*dz

ymax = real(nymax)*dz

WRITE (2,'(2F8.3)') xmin , xmax , ymin , ymax , zminn , zmaxx

nn = (nxmax - nxmin + 1)*(nymax - nymin + 1)*nz

WRITE (2,'(I8, F8.3)') nn , dz !=total number of grid points, spacing

ENDIF

! =coordinates of origin point, grid spacing

! (and note that the range of z is smaller than the original zmin to

! zmax)

IF (cyl) THEN

! Cylindrically symmetric system, 2D

DO ir = 1 , nr

r = real(ir - 1)*dz

r2 = r**2

DO iz = n1 , n2

z = dz*real(iz-1) + zmin

bz = b(iz) - 0.25*r2*ddb(iz) + (r2**2/64.)*ddddb(iz)

br = -0.5*r*db(iz) + (r*r2/16.)*dddb(iz)

WRITE (2,'(2F8.3,1P2E12.4)') r , z , br , bz

! =required field components

ENDDO

ENDDO

ELSE

! 3D system

DO ix = nxmin , nxmax

x = real(ix)*dz

DO iy = nymin , nymax

y = real(iy)*dz

r2 = x**2 + y**2

r = sqrt(r2)

IF (abs(r) .GT. 1.E-10) THEN

cx = x/r

cy = y/r

! =direction cosines in x and y directions

ELSE

cx = 0.

cy = 0.

ENDIF

DO iz = n1 , n2

z = dz*real(iz-1) + zmin

bz = b(iz) - 0.25*r2*ddb(iz) + (r2**2/64.)*ddddb(iz)

br = -0.5*r*db(iz) + (r*r2/16.)*dddb(iz)

bx = br*cx

by = br*cy

WRITE (2,'(3F8.3,1P3E12.4)') x , y , z , bx , by , bz

! =required field components

ENDDO

ENDDO

ENDDO

ENDIF

STOP

END