realslong
begin
integer a line, b line, c line, d line,
e line, f line, g prime line,
h line
real array spectral line(1:8)
integer current surfaces
real array test case(1:4, 1:4)
integer curvature radius, index of refraction, dispersion,
edge thickness
real clear aperture, aberr lspher, aberr osc, aberr lchrom,
max lspher, max osc, max lchrom, radius of curvature,
object distance, ray height, axis slope angle,
from index, to index
integer number of iterations, iteration
integer paraxial, marginal ray, paraxial ray
integer od field, sa field
real array od sa(1:2, 1:2)
real od cline, od fline
real fn cot(real x)
result = 1.0 / (sin(x) / cos(x))
end {cot}
real fn arc sin(real x)
result = 2 * arc tan1(x / (1 + sqrt(1 - x^2)))
end {arc sin}
routine out big real(integer channel, real x)
real fraction
fraction = x
end {out big real}
routine transit surface
real iang, rang, iang sin, rang sin,
old axis slope angle, sagitta
if paraxial = paraxial ray then start
if radius of curvature # 0 then start
if object distance = 0 then start
axis slope angle = 0
iang sin = ray height / radius of curvature
finish else c
iang sin = ((object distance - c
radius of curvature) / radius of curvature) * c
axis slope angle
rang sin = (from index / to index) * iang sin
old axis slope angle = axis slope angle
axis slope angle = axis slope angle + c
iang sin - rang sin
if object distance # 0 then c
ray height = object distance * old axis slope angle
object distance = ray height / axis slope angle
else
object distance = object distance * (to index / from index)
axis slope angle = axis slope angle * (from index / to index)
finish
else
if radius of curvature # 0 then start
if object distance = 0 then start
axis slope angle = 0
iang sin = ray height / radius of curvature
finish else c
iang sin = ((object distance - c
radius of curvature) / radius of curvature) * c
sin(axis slope angle)
iang = arc sin(iang sin)
rang sin = (from index / to index) * iang sin
old axis slope angle = axis slope angle
axis slope angle = axis slope angle + c
iang - arcsin(rang sin)
! for some reason imp2021 likes this but not the original!:
sagitta = 2 * radius of curvature * (sin((old axis slope angle + iang) / 2) ^ 2)
! sagitta = sin((old axis slope angle + iang) / 2)
! sagitta = 2 * radius of curvature * (sagitta ^ 2)
!%comment sagitta = sin((old axis slope angle + iang) / 2)
!%comment sagitta = 2 * radius of curvature * (sagitta ^ 2)
object distance = ((radius of curvature * c
sin(old axis slope angle + iang)) * c
cot(axis slope angle)) + sagitta
else
rang = -arcsin((from index / to index) * c
sin(axis slope angle))
object distance = object distance * ((to index * c
cos(-rang)) / (from index * c
cos(axis slope angle)))
axis slope angle = -rang
finish
finish
end {transit surface}
routine trace line(integer line, real ray h)
integer i
object distance = 0
ray height = ray h
from index = 1
for i = 1, 1, current surfaces cycle
radius of curvature = test case(i, curvature radius)
to index = test case(i, index of refraction)
if to index > 1 then c
to index = to index + ((spectral line(d line) - c
spectral line(line)) / c
(spectral line(c line) - spectral line(f line))) * c
((test case(i, index of refraction) - 1.0) / c
test case(i, dispersion))
transit surface
from index = to index
if i < current surfaces then c
object distance = object distance - c
test case(i, edge thickness)
repeat
end {trace line}
comment END GLOBAL DECLARATONS
comment Spectral lines
a line = 1; b line = 2; c line = 3; d line = 4
e line = 5; f line = 6; g prime line = 7
h line = 8
spectral line(a line) = 7621.0
spectral line(b line) = 6869.955
spectral line(c line) = 6562.8160
spectral line(d line) = 5895.944
spectral line(e line) = 5269.557
spectral line(f line) = 4861.344
spectral line(g prime line) = 4340.477
spectral line(h line) = 3968.494
comment Initialise the test case array
current surfaces = 4
curvature radius = 1; index of refraction = 2; dispersion = 3
edge thickness = 4
test case(1, curvature radius) = 27.05
test case(1, index of refraction) = 1.5137
test case(1, dispersion) = 63.6
test case(1, edge thickness) = 0.52
test case(2, curvature radius) = -16.68
test case(2, index of refraction) = 1.0
test case(2, dispersion) = 0.0
test case(2, edge thickness) = 0.138
test case(3, curvature radius) = -16.68
test case(3, index of refraction) = 1.6164
test case(3, dispersion) = 36.7
test case(3, edge thickness) = 0.38
test case(4, curvature radius) = -78.1
test case(4, index of refraction) = 1.0
test case(4, dispersion) = 0.0
test case(4, edge thickness) = 0.0
marginal ray = 1; paraxial ray = 2
od field = 1; sa field = 2
comment number of iterations = 40263820;
number of iterations = 1
clear aperture = 4
for iteration = 1, 1, number of iterations cycle
for paraxial = marginal ray, 1, paraxial ray cycle
trace line(d line, clear aperture / 2)
od sa(paraxial, od field) = object distance
od sa(paraxial, sa field) = axis slope angle
repeat
comment Trace marginal ray in C
paraxial = marginal ray
trace line(c line, clear aperture / 2)
od cline = object distance
comment Trace marginal ray in F
trace line(f line, clear aperture / 2)
od fline = object distance
aberr lspher = od sa(paraxial ray, od field) - c
od sa(marginal ray, od field)
aberr osc = 1 - ((od sa(paraxial ray, od field) * c
od sa(paraxial ray, sa field)) / c
(sin(od sa(marginal ray, sa field)) * c
od sa(marginal ray, od field)))
aberr lchrom = od fline - od cline
max lspher = 0.0000926 / sin(od sa(marginal ray, sa field)) ^ 2
max lchrom = max lspher
max osc = 0.0025
repeat
comment Print the analysis of the ray trace
print string(" Marginal ray ")
outbigreal(1, od sa(marginal ray, od field))
print string(" ")
print(od sa(marginal ray, sa field), 5, 11); newline
print string(" Paraxial ray ")
outbigreal(1, od sa(paraxial ray, od field))
print string(" ")
print(od sa(paraxial ray, sa field), 5, 11); newline
print string("Longitudinal spherical aberration: ")
print(aberr lspher, 5, 11); newline
print string(" (Maximum permissible): ")
print(max lspher, 5, 11); newline
print string("Offense against sine condition (coma): ")
print(aberr osc, 5, 11); newline
print string(" (Maximum permissible): ")
print(max osc, 5, 11); newline
print string("Axial chromatic aberration: ")
print(aberr lchrom, 5, 11); newline
print string(" (Maximum permissible): ")
print(max lchrom, 5, 11); newline
end of program