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