Skip to content

Commit c998070

Browse files
mohawk2mohawk2
committed
rewrite as PP function
1 parent 4ae296f commit c998070

File tree

1 file changed

+116
-140
lines changed

1 file changed

+116
-140
lines changed

‎src/perl/fbench-pdl.pl

Lines changed: 116 additions & 140 deletions
Original file line numberDiff line numberDiff line change
@@ -75,128 +75,113 @@
7575
# classic work on ray tracing by hand, given in Amateur Telescope
7676
# Making, Volume 3, p597.
7777

78-
my @testcase = (
78+
my $testcase = pdl(
7979
# radius of curvature, refractive index, Abbe number, axial distance
8080
[ 27.05, 1.5137, 63.6, 0.52 ], # crown
8181
[ -16.68, 1, 0, 0.138 ],
8282
[ -16.68, 1.6164, 36.7, 0.38 ], # flint
8383
[ -78.1, 1, 0, 0 ]
8484
);
8585

86-
# Perform ray trace in specific spectral line
87-
sub trace_line {
88-
my ($paraxial, $line, $ray_h) = @_;
89-
$object_distance .= 0;
90-
$ray_height .= $ray_h;
91-
$from_index .= 1;
92-
for (@s) {
93-
($radius_of_curvature, $to_index, my $abbe_number, my $axial_distance) = map PDL::Core::topdl($_), @$_;
94-
if ($to_index > 1) {
95-
$to_index += (
96-
($spectral_line->slice("(4)") - $spectral_line->slice("($line)")) /
97-
($spectral_line->slice("(3)") - $spectral_line->slice("(6)"))
98-
) *
99-
(($to_index - 1) / $abbe_number);
100-
}
101-
$paraxial ? transit_surface_paraxial() : transit_surface();
102-
$from_index .= $to_index;
103-
$object_distance -= $axial_distance;
86+
use Inline Pdlpp => <<'EOF';
87+
pp_def('trace_line',
88+
Pars => '
89+
surface(elt=4,s); spectral_lines(l); int paraxial(); indx line(); ray_height();
90+
[o] object_distance(s); [o] slope_angle(s);',
91+
GenericTypes => ['D'],
92+
Code => <<'EOCODE',
93+
/*
94+
Inputs:
95+
radius_of_curvature Radius of curvature of surface
96+
being crossed. If 0, surface is plane.
97+
object_distance Distance of object focus from
98+
lens vertex. If 0, incoming rays are parallel and
99+
the following must be specified:
100+
ray_height Height of ray from axis. Only
101+
relevant if $object_distance == 0
102+
axis_slope_angle Angle incoming ray makes with axis at intercept
103+
from_index Refractive index of medium being left
104+
to_index Refractive index of medium being entered.
105+
Outputs:
106+
object_distance Distance from vertex to object focus after refraction.
107+
axis_slope_angle Angle incoming ray makes with axis
108+
at intercept after refraction.
109+
*/
110+
$GENERIC() o_d = 0, s_a, r_h = $ray_height(), from_index = 1;
111+
loop(s) %{
112+
/* Perform ray trace in specific spectral line */
113+
$GENERIC() radius_of_curvature = $surface(elt=>0, s=>s),
114+
to_index = $surface(elt=>1, s=>s),
115+
abbe_number = $surface(elt=>2, s=>s),
116+
axial_distance = $surface(elt=>3, s=>s);
117+
PDL_Indx which_line = $line();
118+
if (to_index > 1) {
119+
to_index += (
120+
($spectral_lines(l=>4) - $spectral_lines(l=>which_line)) /
121+
($spectral_lines(l=>3) - $spectral_lines(l=>6))
122+
) *
123+
((to_index - 1) / abbe_number);
104124
}
105-
}
106-
107-
# Calculate passage through surface using the paraxial approximations.
108-
#
109-
# This subroutine takes the following global inputs:
110-
#
111-
# $radius_of_curvature Radius of curvature of surface
112-
# being crossed. If 0, surface is
113-
# plane.
114-
#
115-
# $object_distance Distance of object focus from
116-
# lens vertex. If 0, incoming
117-
# rays are parallel and
118-
# the following must be specified:
119-
#
120-
# $ray_height Height of ray from axis. Only
121-
# relevant if $object_distance == 0
122-
#
123-
# $axis_slope_angle Angle incoming ray makes with axis
124-
# at intercept
125-
#
126-
# $from_index Refractive index of medium being left
127-
#
128-
# $to_index Refractive index of medium being
129-
# entered.
130-
#
131-
# The outputs are the following global variables:
132-
#
133-
# $object_distance Distance from vertex to object focus
134-
# after refraction.
135-
#
136-
# $axis_slope_angle Angle incoming ray makes with axis
137-
# at intercept after refraction.
138-
sub transit_surface_paraxial {
139-
my ($iang_sin, # Incidence angle sin
140-
$rang_sin, # Refraction angle sin
141-
);
142-
if ($radius_of_curvature != 0) {
143-
if ($object_distance == 0) {
144-
$axis_slope_angle .= 0;
145-
$iang_sin = $ray_height / $radius_of_curvature;
125+
$GENERIC() iang_sin /* Incidence angle sin */,
126+
rang_sin /* Refraction angle sin */;
127+
if ($paraxial()) {
128+
if (radius_of_curvature != 0) {
129+
if (o_d == 0) {
130+
s_a = 0;
131+
iang_sin = r_h / radius_of_curvature;
132+
} else {
133+
iang_sin = ((o_d -
134+
radius_of_curvature) / radius_of_curvature) *
135+
s_a;
136+
}
137+
rang_sin = (from_index / to_index) *
138+
iang_sin;
139+
$GENERIC() old_axis_slope_angle = s_a;
140+
s_a += iang_sin - rang_sin;
141+
if (o_d != 0)
142+
r_h = o_d * old_axis_slope_angle;
143+
o_d = r_h / s_a;
146144
} else {
147-
$iang_sin = (($object_distance -
148-
$radius_of_curvature) / $radius_of_curvature) *
149-
$axis_slope_angle;
150-
}
151-
$rang_sin = ($from_index / $to_index) *
152-
$iang_sin;
153-
my $old_axis_slope_angle = $axis_slope_angle->copy;
154-
$axis_slope_angle .= $axis_slope_angle +
155-
$iang_sin - $rang_sin;
156-
if ($object_distance != 0) {
157-
$ray_height .= $object_distance * $old_axis_slope_angle;
145+
o_d *= to_index / from_index;
146+
s_a *= from_index / to_index;
158147
}
159-
$object_distance .= $ray_height / $axis_slope_angle;
160-
return;
161-
}
162-
$object_distance *= ($to_index / $from_index);
163-
$axis_slope_angle *= ($from_index / $to_index);
164-
}
165-
166-
# Calculate passage through surface using normal trigonometric trace
167-
# inputs and output same as paraxial version above
168-
sub transit_surface {
169-
my ($iang_sin, # Incidence angle sin
170-
$rang_sin, # Refraction angle sin
171-
);
172-
if ($radius_of_curvature != 0) {
173-
if ($object_distance == 0) {
174-
$axis_slope_angle .= 0;
175-
$iang_sin = $ray_height / $radius_of_curvature;
176-
} else {
177-
$iang_sin = (($object_distance -
178-
$radius_of_curvature) / $radius_of_curvature) *
179-
sin($axis_slope_angle);
180-
}
181-
my $iang = asin($iang_sin); # Incidence angle
182-
$rang_sin = ($from_index / $to_index) *
183-
$iang_sin;
184-
my $old_axis_slope_angle = $axis_slope_angle->copy;
185-
$axis_slope_angle += $iang - asin($rang_sin);
186-
my $sagitta = sin(($old_axis_slope_angle + $iang) / 2);
187-
$sagitta .= 2 * $radius_of_curvature * $sagitta * $sagitta;
188-
$object_distance .= (($radius_of_curvature * sin(
189-
$old_axis_slope_angle + $iang)) /
190-
tan($axis_slope_angle)) + $sagitta;
191-
return;
148+
} else {
149+
if (radius_of_curvature != 0) {
150+
if (o_d == 0) {
151+
s_a = 0;
152+
iang_sin = r_h / radius_of_curvature;
153+
} else {
154+
iang_sin = ((o_d -
155+
radius_of_curvature) / radius_of_curvature) *
156+
sin(s_a);
192157
}
193-
my $rang = -asin(($from_index / $to_index) *
194-
sin($axis_slope_angle)); # Refraction angle
195-
$object_distance *= (($to_index *
196-
cos(-$rang)) / ($from_index *
197-
cos($axis_slope_angle)));
198-
$axis_slope_angle .= -$rang;
199-
}
158+
$GENERIC() iang = asin(iang_sin); /* Incidence angle */
159+
rang_sin = (from_index / to_index) *
160+
iang_sin;
161+
$GENERIC() old_axis_slope_angle = s_a;
162+
s_a += iang - asin(rang_sin);
163+
$GENERIC() sagitta = sin((old_axis_slope_angle + iang) / 2);
164+
sagitta = 2 * radius_of_curvature * sagitta * sagitta;
165+
o_d = ((radius_of_curvature * sin(
166+
old_axis_slope_angle + iang)) /
167+
tan(s_a)) + sagitta;
168+
} else {
169+
$GENERIC() rang = -asin((from_index / to_index) *
170+
sin(s_a)); /* Refraction angle */
171+
o_d *= ((to_index *
172+
cos(-rang)) / (from_index *
173+
cos(s_a)));
174+
s_a = -rang;
175+
}
176+
}
177+
from_index = to_index;
178+
o_d -= axial_distance;
179+
$object_distance() = o_d;
180+
$slope_angle() = s_a;
181+
%}
182+
EOCODE
183+
);
184+
EOF
200185

201186
# Process the number of iterations argument, if one is supplied.
202187

@@ -221,7 +206,6 @@ sub transit_surface {
221206
# Load test case into working array
222207

223208
my $clear_aperture = 4;
224-
@s = map [@$_], @testcase; # one-level copy
225209

226210
my $nik = $niter / 1000;
227211
print << "EOD";
@@ -241,34 +225,26 @@ sub transit_surface {
241225
my ($od_fline, $od_cline);
242226

243227
for (my $itercount = 0; $itercount < $niter; $itercount++) {
244-
245-
for my $paraxial (0, 1) {
246-
# Do main trace in D light
247-
trace_line($paraxial, 4, $clear_aperture / 2);
248-
$od_sa[$paraxial] = [$object_distance->copy, $axis_slope_angle->copy];
249-
}
250-
251-
# Trace marginal ray in C
252-
253-
trace_line(0, 3, $clear_aperture / 2);
254-
$od_cline = $object_distance->copy;
255-
256-
# Trace marginal ray in F
257-
258-
trace_line(0, 6, $clear_aperture / 2);
259-
$od_fline = $object_distance->copy;
260-
261-
$aberr_lspher = $od_sa[1][0] - $od_sa[0][0];
262-
$aberr_osc = 1 - ($od_sa[1][0] * $od_sa[1][1]) /
263-
(sin($od_sa[0][1]) * $od_sa[0][0]);
264-
$aberr_lchrom = $od_fline - $od_cline;
265-
$max_lspher = sin($od_sa[0][1]);
266-
267-
# D light
268-
269-
$max_lspher = 0.0000926 / ($max_lspher * $max_lspher);
270-
$max_osc = 0.0025;
271-
$max_lchrom = $max_lspher;
228+
my @inputs = (
229+
$testcase,
230+
$spectral_line,
231+
pdl(0, 1, 0, 0), # paraxial
232+
pdl(4, 4, 3, 6), # spectral line - main trace in D light, marginal in C,F
233+
pdl($clear_aperture / 2), # ray height, threads so no need to repeat
234+
);
235+
my ($od, $sa) = PDL::trace_line(@inputs);
236+
my $pdl_od_sa = pdl($od, $sa)->slice("(3)")->transpose; # slice as only last col is of interest
237+
@od_sa = @{ $pdl_od_sa->slice(",0:1")->unpdl };
238+
($od_cline, $od_fline) = @{ $pdl_od_sa->slice("(0),2:3")->unpdl };
239+
$aberr_lspher = $od_sa[1][0] - $od_sa[0][0];
240+
$aberr_osc = 1 - ($od_sa[1][0] * $od_sa[1][1]) /
241+
(sin($od_sa[0][1]) * $od_sa[0][0]);
242+
$aberr_lchrom = $od_fline - $od_cline;
243+
$max_lspher = sin($od_sa[0][1]);
244+
# D light
245+
$max_lspher = 0.0000926 / ($max_lspher * $max_lspher);
246+
$max_osc = 0.0025;
247+
$max_lchrom = $max_lspher;
272248
}
273249

274250
my $interval = tv_interval(\@t);

0 commit comments

Comments
 (0)