From c9980708797da27511a0de2859998190f54d719c Mon Sep 17 00:00:00 2001 From: Ed J Date: Thu, 30 Sep 2021 18:35:15 +0100 Subject: [PATCH] rewrite as PP function --- src/perl/fbench-pdl.pl | 256 +++++++++++++++++++---------------------- 1 file changed, 116 insertions(+), 140 deletions(-) diff --git a/src/perl/fbench-pdl.pl b/src/perl/fbench-pdl.pl index b131357..92ced61 100644 --- a/src/perl/fbench-pdl.pl +++ b/src/perl/fbench-pdl.pl @@ -75,7 +75,7 @@ # classic work on ray tracing by hand, given in Amateur Telescope # Making, Volume 3, p597. -my @testcase = ( +my $testcase = pdl( # radius of curvature, refractive index, Abbe number, axial distance [ 27.05, 1.5137, 63.6, 0.52 ], # crown [ -16.68, 1, 0, 0.138 ], @@ -83,120 +83,105 @@ [ -78.1, 1, 0, 0 ] ); -# Perform ray trace in specific spectral line -sub trace_line { - my ($paraxial, $line, $ray_h) = @_; - $object_distance .= 0; - $ray_height .= $ray_h; - $from_index .= 1; - for (@s) { - ($radius_of_curvature, $to_index, my $abbe_number, my $axial_distance) = map PDL::Core::topdl($_), @$_; - if ($to_index > 1) { - $to_index += ( - ($spectral_line->slice("(4)") - $spectral_line->slice("($line)")) / - ($spectral_line->slice("(3)") - $spectral_line->slice("(6)")) - ) * - (($to_index - 1) / $abbe_number); - } - $paraxial ? transit_surface_paraxial() : transit_surface(); - $from_index .= $to_index; - $object_distance -= $axial_distance; +use Inline Pdlpp => <<'EOF'; +pp_def('trace_line', + Pars => ' + surface(elt=4,s); spectral_lines(l); int paraxial(); indx line(); ray_height(); + [o] object_distance(s); [o] slope_angle(s);', + GenericTypes => ['D'], + Code => <<'EOCODE', + /* + Inputs: + radius_of_curvature Radius of curvature of surface + being crossed. If 0, surface is plane. + object_distance Distance of object focus from + lens vertex. If 0, incoming rays are parallel and + the following must be specified: + ray_height Height of ray from axis. Only + relevant if $object_distance == 0 + axis_slope_angle Angle incoming ray makes with axis at intercept + from_index Refractive index of medium being left + to_index Refractive index of medium being entered. + Outputs: + object_distance Distance from vertex to object focus after refraction. + axis_slope_angle Angle incoming ray makes with axis + at intercept after refraction. + */ + $GENERIC() o_d = 0, s_a, r_h = $ray_height(), from_index = 1; + loop(s) %{ + /* Perform ray trace in specific spectral line */ + $GENERIC() radius_of_curvature = $surface(elt=>0, s=>s), + to_index = $surface(elt=>1, s=>s), + abbe_number = $surface(elt=>2, s=>s), + axial_distance = $surface(elt=>3, s=>s); + PDL_Indx which_line = $line(); + if (to_index > 1) { + to_index += ( + ($spectral_lines(l=>4) - $spectral_lines(l=>which_line)) / + ($spectral_lines(l=>3) - $spectral_lines(l=>6)) + ) * + ((to_index - 1) / abbe_number); } -} - -# Calculate passage through surface using the paraxial approximations. -# -# This subroutine takes the following global inputs: -# -# $radius_of_curvature Radius of curvature of surface -# being crossed. If 0, surface is -# plane. -# -# $object_distance Distance of object focus from -# lens vertex. If 0, incoming -# rays are parallel and -# the following must be specified: -# -# $ray_height Height of ray from axis. Only -# relevant if $object_distance == 0 -# -# $axis_slope_angle Angle incoming ray makes with axis -# at intercept -# -# $from_index Refractive index of medium being left -# -# $to_index Refractive index of medium being -# entered. -# -# The outputs are the following global variables: -# -# $object_distance Distance from vertex to object focus -# after refraction. -# -# $axis_slope_angle Angle incoming ray makes with axis -# at intercept after refraction. -sub transit_surface_paraxial { - my ($iang_sin, # Incidence angle sin - $rang_sin, # Refraction angle sin - ); - if ($radius_of_curvature != 0) { - if ($object_distance == 0) { - $axis_slope_angle .= 0; - $iang_sin = $ray_height / $radius_of_curvature; + $GENERIC() iang_sin /* Incidence angle sin */, + rang_sin /* Refraction angle sin */; + if ($paraxial()) { + if (radius_of_curvature != 0) { + if (o_d == 0) { + s_a = 0; + iang_sin = r_h / radius_of_curvature; + } else { + iang_sin = ((o_d - + radius_of_curvature) / radius_of_curvature) * + s_a; + } + rang_sin = (from_index / to_index) * + iang_sin; + $GENERIC() old_axis_slope_angle = s_a; + s_a += iang_sin - rang_sin; + if (o_d != 0) + r_h = o_d * old_axis_slope_angle; + o_d = r_h / s_a; } else { - $iang_sin = (($object_distance - - $radius_of_curvature) / $radius_of_curvature) * - $axis_slope_angle; - } - $rang_sin = ($from_index / $to_index) * - $iang_sin; - my $old_axis_slope_angle = $axis_slope_angle->copy; - $axis_slope_angle .= $axis_slope_angle + - $iang_sin - $rang_sin; - if ($object_distance != 0) { - $ray_height .= $object_distance * $old_axis_slope_angle; + o_d *= to_index / from_index; + s_a *= from_index / to_index; } - $object_distance .= $ray_height / $axis_slope_angle; - return; - } - $object_distance *= ($to_index / $from_index); - $axis_slope_angle *= ($from_index / $to_index); -} - -# Calculate passage through surface using normal trigonometric trace -# inputs and output same as paraxial version above -sub transit_surface { - my ($iang_sin, # Incidence angle sin - $rang_sin, # Refraction angle sin - ); - if ($radius_of_curvature != 0) { - if ($object_distance == 0) { - $axis_slope_angle .= 0; - $iang_sin = $ray_height / $radius_of_curvature; - } else { - $iang_sin = (($object_distance - - $radius_of_curvature) / $radius_of_curvature) * - sin($axis_slope_angle); - } - my $iang = asin($iang_sin); # Incidence angle - $rang_sin = ($from_index / $to_index) * - $iang_sin; - my $old_axis_slope_angle = $axis_slope_angle->copy; - $axis_slope_angle += $iang - asin($rang_sin); - my $sagitta = sin(($old_axis_slope_angle + $iang) / 2); - $sagitta .= 2 * $radius_of_curvature * $sagitta * $sagitta; - $object_distance .= (($radius_of_curvature * sin( - $old_axis_slope_angle + $iang)) / - tan($axis_slope_angle)) + $sagitta; - return; + } else { + if (radius_of_curvature != 0) { + if (o_d == 0) { + s_a = 0; + iang_sin = r_h / radius_of_curvature; + } else { + iang_sin = ((o_d - + radius_of_curvature) / radius_of_curvature) * + sin(s_a); } - my $rang = -asin(($from_index / $to_index) * - sin($axis_slope_angle)); # Refraction angle - $object_distance *= (($to_index * - cos(-$rang)) / ($from_index * - cos($axis_slope_angle))); - $axis_slope_angle .= -$rang; -} + $GENERIC() iang = asin(iang_sin); /* Incidence angle */ + rang_sin = (from_index / to_index) * + iang_sin; + $GENERIC() old_axis_slope_angle = s_a; + s_a += iang - asin(rang_sin); + $GENERIC() sagitta = sin((old_axis_slope_angle + iang) / 2); + sagitta = 2 * radius_of_curvature * sagitta * sagitta; + o_d = ((radius_of_curvature * sin( + old_axis_slope_angle + iang)) / + tan(s_a)) + sagitta; + } else { + $GENERIC() rang = -asin((from_index / to_index) * + sin(s_a)); /* Refraction angle */ + o_d *= ((to_index * + cos(-rang)) / (from_index * + cos(s_a))); + s_a = -rang; + } + } + from_index = to_index; + o_d -= axial_distance; + $object_distance() = o_d; + $slope_angle() = s_a; + %} +EOCODE +); +EOF # Process the number of iterations argument, if one is supplied. @@ -221,7 +206,6 @@ sub transit_surface { # Load test case into working array my $clear_aperture = 4; -@s = map [@$_], @testcase; # one-level copy my $nik = $niter / 1000; print << "EOD"; @@ -241,34 +225,26 @@ sub transit_surface { my ($od_fline, $od_cline); for (my $itercount = 0; $itercount < $niter; $itercount++) { - - for my $paraxial (0, 1) { - # Do main trace in D light - trace_line($paraxial, 4, $clear_aperture / 2); - $od_sa[$paraxial] = [$object_distance->copy, $axis_slope_angle->copy]; - } - - # Trace marginal ray in C - - trace_line(0, 3, $clear_aperture / 2); - $od_cline = $object_distance->copy; - - # Trace marginal ray in F - - trace_line(0, 6, $clear_aperture / 2); - $od_fline = $object_distance->copy; - - $aberr_lspher = $od_sa[1][0] - $od_sa[0][0]; - $aberr_osc = 1 - ($od_sa[1][0] * $od_sa[1][1]) / - (sin($od_sa[0][1]) * $od_sa[0][0]); - $aberr_lchrom = $od_fline - $od_cline; - $max_lspher = sin($od_sa[0][1]); - - # D light - - $max_lspher = 0.0000926 / ($max_lspher * $max_lspher); - $max_osc = 0.0025; - $max_lchrom = $max_lspher; + my @inputs = ( + $testcase, + $spectral_line, + pdl(0, 1, 0, 0), # paraxial + pdl(4, 4, 3, 6), # spectral line - main trace in D light, marginal in C,F + pdl($clear_aperture / 2), # ray height, threads so no need to repeat + ); + my ($od, $sa) = PDL::trace_line(@inputs); + my $pdl_od_sa = pdl($od, $sa)->slice("(3)")->transpose; # slice as only last col is of interest + @od_sa = @{ $pdl_od_sa->slice(",0:1")->unpdl }; + ($od_cline, $od_fline) = @{ $pdl_od_sa->slice("(0),2:3")->unpdl }; + $aberr_lspher = $od_sa[1][0] - $od_sa[0][0]; + $aberr_osc = 1 - ($od_sa[1][0] * $od_sa[1][1]) / + (sin($od_sa[0][1]) * $od_sa[0][0]); + $aberr_lchrom = $od_fline - $od_cline; + $max_lspher = sin($od_sa[0][1]); + # D light + $max_lspher = 0.0000926 / ($max_lspher * $max_lspher); + $max_osc = 0.0025; + $max_lchrom = $max_lspher; } my $interval = tv_interval(\@t);