Programs for analysing the output data from xmpl3d26.dat, Aberrations due to a localised potential defect on an aperture
These can be copied from the example file, which is supplied with the CPO3D package.
There are two programs (which are rather old and so are in an old version of Fortran).
Program patchfor.for:
c Integrate paraxial aberration coefficients for patch field,
c using axial values output by CPO3D with patch.dat,
c or take output data temp0a.dat and strip it of all except the last lines.
c Call patchbat, then it gives data in temp.dat, which is put into patchres.dat
c and also into temp0a.dat, for calling patchabb, which gives final analysis
c of means, rms. etc, in temp.dat (for putting in patchres.dat).
c (note: for the single defect results the C coefficients were *0.5 immediately after the
c integrations, but for the double defect the *0.5 is applied, correctly, when
c the deflections are being calculated).
dimension p(200),pd(200),pdd(200),z(200),ez(200),ex(200),
& g0(200),g1(200),g2(200),r(200),rd(200),rdd(200),u(200),
& u1(200),xf(50),yf(50),xdf(50),ydf(50),g0d(200),
& g0dd(200),g1d(200),g1dd(200),pddd(200),pdddd(200),
& g3(200),chiy1(200)
c =potential, 1st differential, 2nd, sqrt pot, etc
logical integrate_coeffs,dbl_defect
OPEN (UNIT=1,FILE='temp0a.dat',STATUS='OLD')
C =trajectory output file
OPEN (UNIT=2,FILE='temp.dat',STATUS='UNKNOWN')
C =results file
nrays=56
nrays=106
c (nrays not relevant when integrating coefficients)
integrate_coeffs=.true.
cc integrate_coeffs=.false.
dbl_defect=.true.
dbl_defect=.false.
dv=0.01
dv=1.0
c =dv of defect
zi=-2.50
zi=-1.00
dxy=0.01
dxy=0.005
dxy=0.001
c =x or y for 2nd, 3rd and 4th sets of fields
c coeffs to be used if they are not recalculated:
rad_ap=0.1
c =radius of aperture
if(.not.dbl_defect)then
C0=1.3557E-01
C1=1.7131E-01
C2=1.8634E-01
C3=3.5065E+01
c dxy=0.005:
C0,C1,C2,C3= 1.3557E-01 1.4886E-01 1.8529E-01 7.5930E+00
else
c For double defect
C0=0.
C1=3.5328E-01
C2=0.
C3=8.5721E-01
c dxy=0.005:
C1=3.5267E-01 C3=5.2161E-01
endif
if(.not.integrate_coeffs)goto 100
c -ie do not integrate aberration coefficients
n=101
c =number of z values
zfocus=3.
zf=-zi
n1=n-1
n2=n-2
n3=n-3
n4=n-4
n5=n-5
n6=n-6
n7=n-7
n8=n-8
n9=n-9
n10=n-10
dz=(zf-zi)/real(n1)
c read info for along axis:
do i=1,n
z(i)=zi+dz*real(i-1)
read(1,*)d1,d2,d3,p(i)
c d1,2,3=dummies
enddo
do i=1,n
read(1,*)d1,d2,d3,ex(i),d4,ez(i)
pd(i)=-ez(i)
g0(i)=-ex(i)
enddo
do i=3,n2
pdd(i)=(8.*(pd(i+1)-pd(i-1))-(pd(i+2)-pd(i-2)))/(12.*dz)
g0d(i)=(8.*(g0(i+1)-g0(i-1))-(g0(i+2)-g0(i-2)))/(12.*dz)
enddo
do i=5,n4
g0dd(i)=(8.*(g0d(i+1)-g0d(i-1))-(g0d(i+2)-g0d(i-2)))/(12.*dz)
pddd(i)=(8.*(pdd(i+1)-pdd(i-1))-(pdd(i+2)-pdd(i-2)))/(12.*dz)
enddo
do i=7,n6
pdddd(i)=(8.*(pddd(i+1)-pddd(i-1))-(pddd(i+2)-pddd(i-2)))
& /(12.*dz)
enddo
c read info for along line at (x,y)=(0,dxy):
x=0.
y=dxy
do i=1,n
read(1,*)d1,d2,d3,ex2,ey2
if((i.ge.5).and.(i.le.n4))then
epsx=ex2+g0(i)-0.125*g0dd(i)*y**2
g2(i)=epsx/y**2
endif
enddo
c read info for along line at (x,y)=(dxy,0):
x=dxy
y=0.
do i=1,n
read(1,*)d1,d2,d3,ex2,ey2
if((i.ge.7).and.(i.le.n6))then
epsx=ex2+g0(i)-0.5*pdd(i)*x-0.25*g0dd(i)*3.*x**2
& +0.0625*pdddd(i)*x**3
chiy1(i)=epsx+g2(i)*2.*x**2
endif
enddo
sp=sqrt(1.5)
c read info for along line at (x,y)=(dxy,sp*dxy):
x=dxy
y=sp*dxy
write(2,'(''i,z,g0,1,2,3: '')')
do i=1,n
read(1,*)d1,d2,d3,ex2,ey2
if((i.ge.7).and.(i.le.n6))then
epsy=ey2-0.5*pdd(i)*y-0.25*g0dd(i)*x*y
& +0.0625*pdddd(i)*y*(x**2+y**2)
chiy2=epsy-g2(i)*2.*x*y
g1(i)=-(3.*sp*chiy1(i)+2.*chiy2)/(sp*dxy)
g3(i)=2.*(sp*chiy1(i)+chiy2)/(sp*dxy**3)
write(2,'(i3,f6.2,1p4e11.2)')i,z(i),g0(i),g1(i),
& g2(i),g3(i)
endif
enddo
do i=9,n8
g1d(i)=(8.*(g1(i+1)-g1(i-1))-(g1(i+2)-g1(i-2)))/(12.*dz)
enddo
do i=11,n10
g1dd(i)=(8.*(g1d(i+1)-g1d(i-1))-(g1d(i+2)-g1d(i-2)))/(12.*dz)
enddo
c
c tests:
sh=sqrt(0.5)
write(2,'(''x,y,test,calc,error:''/
& ''1 line for pots, 2 for fields'')')
i=51
c -for z=0
do nrad=1,5
c read info for a radial line at z=0:
c potentials:
do ii=1,11
read(1,*)x,y,d3,phi_test
phi_calc=p(i)-0.25*pdd(i)*(x**2+y**2)+g0(i)*x
& -0.125*g0dd(i)*x*(x**2+y**2)+0.5*g1(i)*(x**2-y**2)
& +g2(i)*x*(x**2-3.*y**2)/3.
& +(1./64.)*pdddd(i)*(x**2+y**2)**2
& +(1./12.)*g1dd(i)*(-3.*x**2*y**2+y**4)
& +g3(i)*(x**4-6.*x**2*y**2+y**4)/4.
error=2.*abs(phi_calc-phi_test)/(phi_calc+phi_test)
write(2,'(2f6.3,3f10.6)')x,y,phi_test,phi_calc,error
enddo
c fields:
do ii=1,11
read(1,*)x,y,d3,ex_test,ey_test
ex_calc=0.5*pdd(i)*x-g0(i)+0.125*g0dd(i)*(3.*x**2+y**2)
& -g1(i)*x-g2(i)*(x**2-y**2)
& -0.0625*pdddd(i)*x*(x**2+y**2)
& +0.5*g1dd(i)*x*y**2
& -g3(i)*x*(x**2-3.*y**2)
uu=ex_calc+ex_test
if(abs(uu).gt.1.e-10)then
error_x=2.*abs(ex_calc-ex_test)/(ex_calc+ex_test)
else
error_x=0.
endif
ey_calc=0.5*pdd(i)*y+0.25*g0dd(i)*x*y
& +g1(i)*y+g2(i)*2.*x*y
& -0.0625*pdddd(i)*y*(x**2+y**2)
& -(1./6.)*g1dd(i)*y*(2.*y**2-3.*x**2)
& -g3(i)*y*(y**2-3.*x**2)
uu=ey_calc+ey_test
if(abs(uu).gt.1.e-10)then
error_y=2.*abs(ey_calc-ey_test)/(ey_calc+ey_test)
else
error_y=0.
endif
write(2,'(2f6.3,3f10.6/12x,3f10.6)')
& x,y,ex_test,ex_calc,error_x,
& ey_test,ey_calc,error_y
enddo
enddo
c
c integrate aberration coeffs:
C0=0.
C1=0.
C2=0.
C3=0.
do i=5,n4
C0=C0+g0(i)
C1=C1+g1(i)
C2=C2+g2(i)
C3=C3+g3(i)
enddo
C0=C0*dz/dv
C1=C1*dz*rad_ap/dv
C2=C2*dz*rad_ap**2/dv
C3=C3*dz*rad_ap**3/dv
c -all scaled for radius of aperture
100 continue
write(2,'(''C0,C1,C2,C3= '',1p4e12.4)')C0,C1,C2,C3
c
c analyse ray data:
C0=C0*0.5*dv
c (was C0=C0*dv for the single defect)
C1=C1*0.5*dv
C2=C2*0.5*dv
C3=C3*0.5*dv
write(2,'(''x,ya, dxd,dyd, using C0 etc:'')')
do i=1,nrays
read(1,*)xi,yi,zi,vxi,vyi,vzi
c =initial ray parameters
xdi=vxi/vzi
ydi=vyi/vzi
c =initial slopes
xa=xi-zi*xdi
ya=yi-zi*ydi
c =coordinates at aperture
xa=xa/rad_ap
ya=ya/rad_ap
c -scaled to radius of the aperture
read(1,*)xf(i),yf(i),zf,vxf,vyf,vzf
c =cpo final values
xdf(i)=vxf/vzf
ydf(i)=vyf/vzf
c =final slopes
dxd=xdf(i)-xdi
dyd=ydf(i)-ydi
c =change in slope
c use coeffs to calc expected dxd,dyd:
dxd_calc=C0+C1*xa+C2*(xa**2-ya**2)+C3*xa*(xa**2-3.*ya**2)
dyd_calc=-C1*ya-C2*2.*xa*ya+C3*ya*(-3.*xa*ya+ya**2)
write(2,'(2f9.5,2x,2f10.6,2x,2f10.6)')
& xa,ya,dxd,dyd,dxd_calc,dyd_calc
enddo
stop
end
Program patchabf.for:
c analyse for dimensionless aberration figure for patch field.
c take output of patchfor.for, in temp.dat, put in temp0a.dat, run patchabb,
c then output, which gives final analysis, is in temp.dat (for putting in
c patchres.dat).
dimension dxd(500),dyd(500)
c =changes in slope
dimension area_strip(100),d_strip(100)
integer hist(500)
logical double
double=.true.
c -for double defect (symmetric in x)
double=.false.
c -for single defect (not symmetric in x)
nrays=200
nrays=31
rescale=100.
c =energy/dv, to give the required dimensionless result
ndisc=20
c =number of disc sizes for convolution, at end
disc_max=1.0
c =max disc diameter
ppi=4.*atan(1.)
OPEN (UNIT=1,FILE='temp0a.dat',STATUS='OLD')
C = file with data
OPEN (UNIT=2,FILE='temp.dat',STATUS='UNKNOWN')
C =results file
ns=100
c =number for histogram
xm=0.
xmin=1.e10
xmax=-1.e10
ymin=1.e10
do i=1,nrays
read(1,*)d1,d2,d3,d4,dxd(i),dyd(i)
c = ray parameters
dxd(i)=rescale*dxd(i)
dyd(i)=rescale*dyd(i)
if(double)dxd(i)=abs(dxd(i))
c (a few are negative)
c delete:
cc write (2,'(''i,dxd,dyd= '',i4,1p2e10.2)')i,dxd(i),dyd(i)
xm=xm+dxd(i)
C (ym=0., all y values are negative)
xmin=min(xmin,dxd(i))
xmax=max(xmax,dxd(i))
ymin=min(ymin,dyd(i))
c (all y values are negative)
enddo
xm=xm/real(nrays)
if(double)then
xm=0.
xmin=0.
endif
ym=0.
ymax=0.
xrms=0.
yrms=0.
xmabs=0.
ymabs=0.
c =mean absolute deviations
do i=1,nrays
xrms=xrms+(dxd(i)-xm)**2
yrms=yrms+(dyd(i)-ym)**2
xmabs=xmabs+abs(dxd(i)-xm)
ymabs=ymabs+abs(dyd(i)-xm)
enddo
xrms=sqrt(xrms/real(nrays))
yrms=sqrt(yrms/real(nrays))
xmabs=xmabs/real(nrays)
ymabs=ymabs/real(nrays)
write(2,'(''mean x,y= '',2f12.4)')xm,ym
write(2,'(''rms x,y= '',2f12.4)')xrms,yrms
write(2,'(''mean abs x,y= '',2f12.4)')xmabs,ymabs
C find half-widths:
write(2,'(''xmin,max= '',2f12.4)')xmin,xmax
call sortxy(nrays,dxd)
dx=(xmax-xmin)/real(ns)
do is=1,ns
hist(is)=0
enddo
do i=1,nrays
nhist=(dxd(i)-xmin)/dx
nhist=nhist+1
c =position in histogram
if(nhist.lt.0)nhist=0
if(nhist.gt.ns)nhist=ns
hist(nhist)=hist(nhist)+1
C -histogram incremented at position nhist
enddo
cc write(2,'(''xf= '',5f12.4)')(dxd(i),i=1,nrays)
write(2,'(''histogram='',10i3)')(hist(i),i=1,ns)
C and same for y:
write(2,'(''ymin,max= '',2f12.4)')ymin,ymax
call sortxy(nrays,dyd)
dy=(ymax-ymin)/real(ns)
do is=1,ns
hist(is)=0
enddo
do i=1,nrays
nhist=(dyd(i)-ymin)/dy
nhist=nhist+1
if(nhist.lt.0)nhist=0
if(nhist.gt.ns)nhist=ns
hist(nhist)=hist(nhist)+1
enddo
cc write(2,'(''yf= '',5f12.4)')(dyd(i),i=1,nrays)
write(2,'(''histogram='',10i3)')(hist(i),i=1,ns)
c
c convolute with a uniform disc
nd=10
c =number of strips into which the uniform disc is divided
ndh=nd/2
nd=2*ndh
c -must be even
write(2,'(''disc, diam, x,yrms, x,yabs='')')
do idisc=1,ndisc
c (ndisc=number of uniform discs)
rdisc=disc_max*0.5*real(idisc)/real(ndisc)
c =radius of present disc
ddisc=2.*rdisc
do id=1,ndh
c id labels strip
h_inner=rdisc*real(id-1)/real(ndh)
h_outer=rdisc*real(id)/real(ndh)
c =distances of boundaries of strips from the centre
theta_inner=acos(h_inner/rdisc)
if(id.ne.ndh)then
theta_outer=acos(h_outer/rdisc)
else
theta_outer=0.
c (to avoid a possible numerical problem)
endif
s_inner=rdisc*sin(theta_inner)
s_outer=rdisc*sin(theta_outer)
c =extent from axis
area_inner=theta_inner/ppi-h_inner*s_inner/(ppi*rdisc**2)
area_outer=theta_outer/ppi-h_outer*s_outer/(ppi*rdisc**2)
c =areas lying outside lines, in units of pi*r**2
area_strip(id)=area_inner-area_outer
area_strip(id+ndh)=area_strip(id)
c =area of strip, normalised
d_strip(id)=0.5*(h_inner+h_outer)
d_strip(id+ndh)=-d_strip(id)
c =mean distance from origin
enddo
cc write(2,'(i3,2f12.6)')(i,area_strip(i),d_strip(i),i=1,nd)
d_strip_min=d_strip(nd)
d_strip_max=d_strip(ndh)
c then find rms etc:
c xm,ym are same as before
xmin=1.e10
xmax=-1.e10
ymin=1.e10
do i=1,nrays
xmin=min(xmin,(dxd(i)+d_strip_min))
xmax=max(xmax,(dxd(i)+d_strip_max))
ymin=min(ymin,(dyd(i)+d_strip_min))
c (all y values are negative)
enddo
if(double)then
xmin=0.
endif
ymax=0.
xrms=0.
yrms=0.
xmabs=0.
ymabs=0.
c =mean absolute deviations
do i=1,nrays
do id=1,nd
c -ie cycle over strips
xrms=xrms+area_strip(id)*(d_strip(id)+dxd(i)-xm)**2
yrms=yrms+area_strip(id)*(d_strip(id)+dyd(i)-ym)**2
xmabs=xmabs+area_strip(id)*abs(d_strip(id)+dxd(i)-xm)
ymabs=ymabs+area_strip(id)*abs(d_strip(id)+dyd(i)-ym)
enddo
enddo
xrms=sqrt(xrms/real(nrays))
yrms=sqrt(yrms/real(nrays))
xmabs=xmabs/real(nrays)
ymabs=ymabs/real(nrays)
write(2,'(5f8.4)')ddisc,xrms,yrms,xmabs,ymabs
c quick test of rms of uniform disc:
u=0.
do id=1,nd
u=u+area_strip(id)*d_strip(id)**2
cc write(2,'(''area,d= '',2f12.4)')area_strip(id),d_strip(id)
enddo
u=sqrt(u)
cc write(2,'(1pe9.2)')u
enddo
stop
end
C
subroutine sortxy(n,a)
c from num recipes
dimension a(500)
do j=2,n
at=a(j)
do i=j-1,1,-1
if(a(i).le.at)goto10
a(i+1)=a(i)
enddo
i=0
10 a(i+1)=at
enddo
return
end