home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c221 / 5.ddi / MWHC.005 / C0 < prev    next >
Encoding:
Text File  |  1992-02-10  |  23.5 KB  |  753 lines
  1. /*
  2.  
  3.     This is a complete optical design raytracing algorithm,
  4.     stripped of its user interface and recast into portable C.  It
  5.     not only determines execution speed on an extremely floating
  6.     point (including trig function) intensive real-world
  7.     application, it checks accuracy on an algorithm that is
  8.     exquisitely sensitive to errors.  The performance of this
  9.     program is typically far more sensitive to changes in the
  10.     efficiency of the trigonometric library routines than the
  11.     average floating point program.
  12.  
  13.     The benchmark may be compiled in two modes.  If the symbol
  14.     INTRIG is defined, built-in trigonometric and square root
  15.     routines will be used for all calculations.  Timings made with
  16.         INTRIG defined reflect the machine's basic floating point
  17.     performance for the arithmetic operators.  If INTRIG is not
  18.     defined, the system library <math.h> functions are used.  Results
  19.         with INTRIG not defined reflect the system's library performance
  20.     and/or floating point hardware support for trig functions and
  21.     square root.  Results with INTRIG defined are a good guide to
  22.     general floating point performance, while results with INTRIG
  23.     undefined indicate the performance of an application which is
  24.     math function intensive.
  25.  
  26.     Special note regarding errors in accuracy: this program has
  27.     generated numbers identical to the last digit it formats and
  28.     checks on the following machines, floating point architectures,
  29.     and languages:
  30.  
  31.     IBM PC / XT / AT  Lattice C IEEE 64 bit, 80 bit temporaries
  32.             High C    same, in line 80x87 code
  33.                           BASICA    "Double precision"
  34.             Quick BASIC IEEE double precision, software routines
  35.  
  36.     Sun 3      C       IEEE 64 bit, 80 bit temporaries,
  37.                  in-line 68881 code, in-line FPA code.
  38.  
  39.         MicroVAX II       C         Vax "G" format floating point
  40.  
  41.     Macintosh Plus   MPW C     SANE floating point, IEEE 64 bit format
  42.                  implemented in ROM.
  43.  
  44.     Inaccuracies reported by this program should be taken VERY
  45.     SERIOUSLY INDEED, as the program has been demonstrated to be
  46.     invariant under changes in floating point format, as long as
  47.     the format is a recognised double precision format.  If you
  48.     encounter errors, please remember that they are just as likely
  49.     to be in the floating point editing library or the
  50.     trigonometric libraries as in the low level operator code.
  51.  
  52.     The benchmark assumes that results are basically reliable, and
  53.     only tests the last result computed against the reference.  If
  54.         you're running on a suspect system you can compile this program
  55.     with ACCURACY defined.   This will generate a version which
  56.     executes as an infinite loop, performing the ray trace and
  57.     checking the results on every pass.  All incorrect results will
  58.     be reported.
  59.  
  60.     Representative timings are given below.  All have been normalised
  61.     as if run for 1000 iterations.
  62.  
  63.   Time in seconds        Computer, Compiler, and notes
  64.  Normal      INTRIG
  65.  
  66.  3466.00    4031.00   Commodore 128, 2 Mhz 8510 with software floating
  67.           point.   Abacus Software/Data-Becker Super-C 128,
  68.           version 3.00, run in fast (2 Mhz) mode.  Note:
  69.           the results generated by this system differed
  70.           from the reference results in the 8th to 10th
  71.           decimal place.
  72.  
  73.  3290.00     IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
  74.                         Run with the "/d" switch, software floating point.
  75.  
  76.  2131.50     IBM PC/AT 6 Mhz, Lattice C version 2.14, small model.
  77.           This version of Lattice compiles subroutine
  78.           calls which either do software floating point
  79.           or use the 80x87.  The machine on which I ran
  80.           this had an 80287, but the results were so bad
  81.           I wonder if it was being used.
  82.  
  83.  1598.00     Macintosh Plus, MPW C, SANE Software floating point.
  84.  
  85.   404.00     IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
  86.           Software floating point.
  87.  
  88.   165.15     IBM PC/AT 6 Mhz, Metaware High C version 1.3, small
  89.           model.   This was compiled to call subroutines for
  90.           floating point, and the machine contained an 80287
  91.           which was used by the subroutines.
  92.  
  93.   143.20     Macintosh II, MPW C, SANE calls.  I was unable to
  94.           determine whether SANE was using the 68881 chip or
  95.           not.
  96.  
  97.   121.80     Sun 3/160 16 Mhz, Sun C.  Compiled with -fsoft switch
  98.           which executes floating point in software.
  99.  
  100.    78.78     IBM RT PC (Model 6150).  IBM AIX 1.0 C compiler
  101.           with -O switch.
  102.  
  103.    66.36     206.35   IBM PC/AT 6Mhz, Metaware High C version 1.3, small
  104.           model.   This was compiled with in-line code for the
  105.           80287 math coprocessor.  Trig functions still call
  106.           library routines.
  107.  
  108.    17.18     Apollo DN-3000, 12 Mhz 68020 with 68881, compiled
  109.           with in-line code for the 68881 coprocessor.
  110.           According to Apollo, the library routines are chosen
  111.           at runtime based on coprocessor presence.  Since the
  112.           coprocessor was present, the library is supposed to
  113.           use in-line floating point code.
  114.  
  115.    15.55      27.56   VAXstation II GPX.  Compiled and executed under
  116.           VAX/VMS C.
  117.  
  118.    15.14      37.93   Macintosh II, Unix system V.  Green Hills 68020
  119.           Unix compiler with in-line code for the 68881
  120.           coprocessor (-O -ZI switches).
  121.  
  122.    12.69     Sun 3/160 16 Mhz, Sun C.  Compiled with -fswitch,
  123.           which calls a subroutine to select the fastest
  124.           floating point processor.  This was using the 68881.
  125.  
  126.    11.74      26.73   Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  127.           Metaware High C version 1.3, compiled with in-line
  128.           for the math coprocessor (but not optimised for the
  129.           80386/80387).  Trig functions still call library
  130.           routines.
  131.  
  132.     8.43      30.49   Sun 3/160 16 Mhz, Sun C.  Compiled with -f68881,
  133.           generating in-line MC68881 instructions.  Trig
  134.           functions still call library routines.
  135.  
  136.     6.29      25.17   Sun 3/260 25 Mhz, Sun C.  Compiled with -f68881,
  137.           generating in-line MC68881 instructions.  Trig
  138.           functions still call library routines.
  139.  
  140.     6.00      23.00   Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  141.           Metaware High C version 1.3, compiled with in-line
  142.           for the math coprocessor (optimised for the
  143.           80386/80387).
  144.  
  145.     5.9       23.00   Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  146.           Metaware High C version 1.4, compiled with in-line
  147.           for the math coprocessor (optimised for the
  148.           80386/80387).
  149.  
  150.     5.25      ---       Intel 386/20, 16 Mhz 80386 with 16 Mhz 80387.
  151.           Metaware High C version 1.4, compiled with in-line
  152.           for the math coprocessor (optimised for the
  153.           80386/80387).
  154.  
  155. */
  156.  
  157. #include <stdio.h>
  158. #ifndef INTRIG
  159. #include <math.h>
  160. #endif
  161.  
  162. extern exit();
  163.  
  164. #define cot(x) (1.0 / tan(x))
  165.  
  166. #define TRUE  1
  167. #define FALSE 0
  168.  
  169. #define max_surfaces 10
  170.  
  171. /*  Local variables  */
  172.  
  173. static char tbfr[132];
  174.  
  175. static short current_surfaces;
  176. static short paraxial;
  177.  
  178. static double clear_aperture;
  179.  
  180. static double aberr_lspher;
  181. static double aberr_osc;
  182. static double aberr_lchrom;
  183.  
  184. static double max_lspher;
  185. static double max_osc;
  186. static double max_lchrom;
  187.  
  188. static double radius_of_curvature;
  189. static double object_distance;
  190. static double ray_height;
  191. static double axis_slope_angle;
  192. static double from_index;
  193. static double to_index;
  194.  
  195. static double spectral_line[9];
  196. static double s[max_surfaces][5];
  197. static double od_sa[2][2];
  198.  
  199. static char outarr[8][80];      /* Computed output of program goes here */
  200.  
  201. int itercount;           /* The iteration counter for the main loop
  202.                    in the program is made global so that
  203.                    the compiler should not be allowed to
  204.                    optimise out the loop over the ray
  205.                    tracing code. */
  206. int niter = 1000;        /* Iteration counter */
  207.  
  208. static char *refarr[] = {    /* Reference results.  These happen to
  209.                    be derived from a run on Microsoft
  210.                    Quick BASIC on the IBM PC/AT. */
  211.  
  212.         "   Marginal ray          47.09479120920   0.04178472683",
  213.         "   Paraxial ray          47.08372160249   0.04177864821",
  214.         "Longitudinal spherical aberration:        -0.01106960671",
  215.         "    (Maximum permissible):                 0.05306749907",
  216.         "Offense against sine condition (coma):     0.00008954761",
  217.         "    (Maximum permissible):                 0.00250000000",
  218.         "Axial chromatic aberration:                0.00448229032",
  219.         "    (Maximum permissible):                 0.05306749907"
  220. };
  221.  
  222. /* The test case used in this program is the design for a 4 inch
  223.    achromatic telescope objective used as the example in Wyld's
  224.    classic work on ray tracing by hand, given in Amateur Telescope
  225.    Making, Volume 3. */
  226.  
  227. static double testcase[4][4] = {
  228.     {27.05, 1.5137, 63.6, 0.52},
  229.     {-16.68, 1, 0, 0.138},
  230.     {-16.68, 1.6164, 36.7, 0.38},
  231.     {-78.1, 1, 0, 0}
  232. };
  233.  
  234. /*  Internal trig functions (used only if INTRIG is defined).  These
  235.     standard functions may be enabled to obtain timings that reflect
  236.     the machine's floating point performance rather than the speed of
  237.     its trig function evaluation.  */
  238.  
  239. #ifdef INTRIG
  240.  
  241. #define fabs(x)  ((x < 0.0) ? -x : x)
  242.  
  243. #define pic 3.1415926535897932
  244.  
  245. /*  Commonly used constants  */
  246.  
  247. static double pi = pic,
  248.     twopi =pic * 2.0,
  249.     piover4 = pic / 4.0,
  250.     fouroverpi = 4.0 / pic,
  251.     piover2 = pic / 2.0;
  252.  
  253. /*  Coefficients for ATAN evaluation  */
  254.  
  255. static double atanc[] = {
  256.     0.0,
  257.     0.4636476090008061165,
  258.     0.7853981633974483094,
  259.     0.98279372324732906714,
  260.     1.1071487177940905022,
  261.     1.1902899496825317322,
  262.     1.2490457723982544262,
  263.     1.2924966677897852673,
  264.     1.3258176636680324644
  265. };
  266.  
  267. /*  aint(x)    Return integer part of number.  Truncates towards 0  */
  268.  
  269. double aint(x)
  270. double x;
  271. {
  272.     long l;
  273.  
  274.     /*  Note that this routine cannot handle the full floating point
  275.         number range.  This function should be in the machine-dependent
  276.         floating point library!  */
  277.  
  278.     l = x;
  279.     if ((int)(-0.5) != 0  &&  l < 0 )
  280.        l++;
  281.     x = l;
  282.     return x;
  283. }
  284.  
  285. /*  sin(x)     Return sine, x in radians  */
  286.  
  287. static double sin(x)
  288. double x;
  289. {
  290.     int sign;
  291.     double y, r, z;
  292.  
  293.     x = (((sign= (x < 0.0)) != 0) ? -x: x);
  294.  
  295.     if (x > twopi)
  296.        x -= (aint(x / twopi) * twopi);
  297.  
  298.     if (x > pi) {
  299.        x -= pi;
  300.        sign = !sign;
  301.     }
  302.  
  303.     if (x > piover2)
  304.        x = pi - x;
  305.  
  306.     if (x < piover4) {
  307.        y = x * fouroverpi;
  308.        z = y * y;
  309.        r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z -
  310.           0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z -
  311.           0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z -
  312.           0.0807455121882807815) * z + 0.785398163397448310);
  313.     } else {
  314.        y = (piover2 - x) * fouroverpi;
  315.        z = y * y;
  316.        r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z -
  317.           0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z -
  318.           0.325991886926687550E-3) * z + 0.0158543442438154109) * z -
  319.           0.308425137534042452) * z + 1.0;
  320.     }
  321.     return sign ? -r : r;
  322. }
  323.  
  324. /*  cos(x)     Return cosine, x in radians, by identity  */
  325.  
  326. static double cos(x)
  327. double x;
  328. {
  329.     x = (x < 0.0) ? -x : x;
  330.     if (x > twopi)          /* Do range reduction here to limit */
  331.        x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2    */
  332.     return sin(x + piover2);
  333. }
  334.  
  335. /*  tan(x)     Return tangent, x in radians, by identity  */
  336.  
  337. static double tan(x)
  338. double x;
  339. {
  340.     return sin(x) / cos(x);
  341. }
  342.  
  343. /*  sqrt(x)    Return square root.  Initial guess, then Newton-
  344.          Raphson refinement  */
  345.  
  346. double sqrt(x)
  347. double x;
  348. {
  349.     double c, cl, y;
  350.     int n;
  351.  
  352.     if (x == 0.0)
  353.        return 0.0;
  354.  
  355.     if (x < 0.0) {
  356.        fprintf(stderr,
  357.               "\nGood work!  You tried to take the square root of %g",
  358.          x);
  359.        fprintf(stderr,
  360.               "\nunfortunately, that is too complex for me to handle.\n");
  361.        exit(1);
  362.     }
  363.  
  364.     y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x);
  365.  
  366.     c = (y - x / y) / 2.0;
  367.     cl = 0.0;
  368.     for (n = 50; c != cl && n--;) {
  369.        y = y - c;
  370.        cl = c;
  371.        c = (y - x / y) / 2.0;
  372.     }
  373.     return y;
  374. }
  375.  
  376. /*  atan(x)    Return arctangent in radians,
  377.          range -pi/2 to pi/2  */
  378.  
  379. static double atan(x)
  380. double x;
  381. {
  382.     int sign, l, y;
  383.     double a, b, z;
  384.  
  385.     x = (((sign = (x < 0.0)) != 0) ? -x : x);
  386.     l = 0;
  387.  
  388.     if (x >= 4.0) {
  389.        l = -1;
  390.        x = 1.0 / x;
  391.        y = 0;
  392.        goto atl;
  393.     } else {
  394.        if (x < 0.25) {
  395.           y = 0;
  396.           goto atl;
  397.        }
  398.     }
  399.  
  400.     y = aint(x / 0.5);
  401.     z = y * 0.5;
  402.     x = (x - z) / (x * z + 1);
  403.  
  404. atl:
  405.     z = x * x;
  406.     b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z +
  407.         1277025750.0) * z + 1550674125.0) * z + 654729075.0;
  408.     a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z +
  409.         1332431100.0) * z + 654729075.0;
  410.     a = (a / b) * x + atanc[y];
  411.     if (l)
  412.        a=piover2 - a;
  413.     return sign ? -a : a;
  414. }
  415.  
  416. /*  atan2(y,x)    Return arctangent in radians of y/x,
  417.          range -pi to pi  */
  418.  
  419. static double atan2(y, x)
  420. double y, x;
  421. {
  422.     double temp;
  423.  
  424.     if (x == 0.0) {
  425.        if (y == 0.0)   /*  Special case: atan2(0,0) = 0  */
  426.           return 0.0;
  427.        else if (y > 0)
  428.           return piover2;
  429.        else
  430.           return -piover2;
  431.     }
  432.     temp = atan(y / x);
  433.     if (x < 0.0) {
  434.        if (y >= 0.0)
  435.           temp += pic;
  436.        else
  437.           temp -= pic;
  438.     }
  439.     return temp;
  440. }
  441.  
  442. /*  asin(x)    Return arcsine in radians of x  */
  443.  
  444. static double asin(x)
  445. double x;
  446. {
  447.     if (fabs(x)>1.0) {
  448.        fprintf(stderr,
  449.               "\nInverse trig functions lose much of their gloss when");
  450.        fprintf(stderr,
  451.               "\ntheir arguments are greater than 1, such as the");
  452.        fprintf(stderr,
  453.               "\nvalue %g you passed.\n", x);
  454.        exit(1);
  455.     }
  456.     return atan2(x, sqrt(1 - x * x));
  457. }
  458. #endif
  459.  
  460. /*  Perform ray trace in specific spectral line  */
  461.  
  462. static trace_line(line, ray_h)
  463. int line;
  464. double ray_h;
  465. {
  466.     int i;
  467.  
  468.     object_distance = 0.0;
  469.     ray_height = ray_h;
  470.     from_index = 1.0;
  471.  
  472.     for (i = 1; i <= current_surfaces; i++) {
  473.        radius_of_curvature = s[i][1];
  474.        to_index = s[i][2];
  475.        if (to_index > 1.0)
  476.           to_index = to_index + ((spectral_line[4] -
  477.         spectral_line[line]) /
  478.         (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
  479.         s[i][3]);
  480.        transit_surface();
  481.        from_index = to_index;
  482.        if (i < current_surfaces)
  483.           object_distance = object_distance - s[i][4];
  484.     }
  485. }
  486.  
  487. /*        Calculate passage through surface
  488.  
  489.           If the variable PARAXIAL is true, the trace through the
  490.           surface will be done using the paraxial approximations.
  491.           Otherwise, the normal trigonometric trace will be done.
  492.  
  493.           This routine takes the following inputs:
  494.  
  495.           RADIUS_OF_CURVATURE     Radius of curvature of surface
  496.                   being crossed.  If 0, surface is
  497.                   plane.
  498.  
  499.           OBJECT_DISTANCE      Distance of object focus from
  500.                   lens vertex.  If 0, incoming
  501.                   rays are parallel and
  502.                   the following must be specified:
  503.  
  504.           RAY_HEIGHT        Height of ray from axis.  Only
  505.                   relevant if OBJECT.DISTANCE == 0
  506.  
  507.           AXIS_SLOPE_ANGLE        Angle incoming ray makes with axis
  508.                   at intercept
  509.  
  510.           FROM_INDEX        Refractive index of medium being left
  511.  
  512.           TO_INDEX          Refractive index of medium being
  513.                   entered.
  514.  
  515.           The outputs are the following variables:
  516.  
  517.           OBJECT_DISTANCE      Distance from vertex to object focus
  518.                   after refraction.
  519.  
  520.           AXIS_SLOPE_ANGLE        Angle incoming ray makes with axis
  521.                   at intercept after refraction.
  522.  
  523. */
  524.  
  525. static transit_surface() {
  526.     double iang,         /* Incidence angle */
  527.            rang,         /* Refraction angle */
  528.            iang_sin,     /* Incidence angle sin */
  529.            rang_sin,     /* Refraction angle sin */
  530.            old_axis_slope_angle, sagitta;
  531.  
  532.     if (paraxial) {
  533.        if (radius_of_curvature != 0.0) {
  534.           if (object_distance == 0.0) {
  535.         axis_slope_angle = 0.0;
  536.         iang_sin = ray_height / radius_of_curvature;
  537.           } else
  538.         iang_sin = ((object_distance -
  539.            radius_of_curvature) / radius_of_curvature) *
  540.            axis_slope_angle;
  541.  
  542.           rang_sin = (from_index / to_index) *
  543.         iang_sin;
  544.           old_axis_slope_angle = axis_slope_angle;
  545.           axis_slope_angle = axis_slope_angle +
  546.         iang_sin - rang_sin;
  547.           if (object_distance != 0.0)
  548.         ray_height = object_distance * old_axis_slope_angle;
  549.           object_distance = ray_height / axis_slope_angle;
  550.           return;
  551.        }
  552.        object_distance = object_distance * (to_index / from_index);
  553.        axis_slope_angle = axis_slope_angle * (from_index / to_index);
  554.        return;
  555.     }
  556.  
  557.     if (radius_of_curvature != 0.0) {
  558.        if (object_distance == 0.0) {
  559.           axis_slope_angle = 0.0;
  560.           iang_sin = ray_height / radius_of_curvature;
  561.        } else {
  562.           iang_sin = ((object_distance -
  563.         radius_of_curvature) / radius_of_curvature) *
  564.         sin(axis_slope_angle);
  565.        }
  566.        iang = asin(iang_sin);
  567.        rang_sin = (from_index / to_index) *
  568.           iang_sin;
  569.        old_axis_slope_angle = axis_slope_angle;
  570.        axis_slope_angle = axis_slope_angle +
  571.           iang - asin(rang_sin);
  572.        sagitta = sin((old_axis_slope_angle + iang) / 2.0);
  573.        sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
  574.        object_distance = ((radius_of_curvature * sin(
  575.           old_axis_slope_angle + iang)) *
  576.           cot(axis_slope_angle)) + sagitta;
  577.        return;
  578.     }
  579.  
  580.     rang = -asin((from_index / to_index) *
  581.        sin(axis_slope_angle));
  582.     object_distance = object_distance * ((to_index *
  583.        cos(-rang)) / (from_index *
  584.        cos(axis_slope_angle)));
  585.     axis_slope_angle = -rang;
  586. }
  587.  
  588. /*  Initialise when called the first time  */
  589.  
  590. main(argc, argv)
  591. int argc;
  592. char *argv[];
  593. {
  594.     int i, j, k, errors;
  595.     double od_fline, od_cline;
  596. #ifdef ACCURACY
  597.     long passes;
  598. #endif
  599.  
  600.     spectral_line[1] = 7621.0;   /* A */
  601.     spectral_line[2] = 6869.955;    /* B */
  602.     spectral_line[3] = 6562.816;    /* C */
  603.     spectral_line[4] = 5895.944;    /* D */
  604.     spectral_line[5] = 5269.557;    /* E */
  605.     spectral_line[6] = 4861.344;    /* F */
  606.         spectral_line[7] = 4340.477;     /* G'*/
  607.     spectral_line[8] = 3968.494;    /* H */
  608.  
  609.     /* Process the number of iterations argument, if one is supplied. */
  610.  
  611.     if (argc > 1) {
  612.        niter = atoi(argv[1]);
  613.            if (*argv[1] == '-' || niter < 1) {
  614.               printf("This is John Walker's floating point accuracy and\n");
  615.               printf("performance benchmark program.  You call it with\n");
  616.               printf("\nfbench <itercount>\n\n");
  617.               printf("where <itercount> is the number of iterations\n");
  618.               printf("to be executed.  Archival timings should be made\n");
  619.               printf("with the iteration count set so that roughly five\n");
  620.               printf("minutes of execution is timed.\n");
  621.           exit(0);
  622.        }
  623.     }
  624.  
  625.     /* Load test case into working array */
  626.  
  627.     clear_aperture = 4.0;
  628.     current_surfaces = 4;
  629.     for (i = 0; i < current_surfaces; i++)
  630.        for (j = 0; j < 4; j++)
  631.           s[i + 1][j + 1] = testcase[i][j];
  632.  
  633. #ifdef ACCURACY
  634.         printf("Beginning execution of floating point accuracy test...\n");
  635.     passes = 0;
  636. #else
  637.         printf("Ready to begin John Walker's floating point accuracy\n");
  638.         printf("and performance benchmark.  %d iterations will be made.\n\n",
  639.        niter);
  640.         printf("Press return to begin benchmark:");
  641.     gets(tbfr);
  642. #endif
  643.  
  644.     /* Perform ray trace the specified number of times. */
  645.  
  646. #ifdef ACCURACY
  647.     while (TRUE) {
  648.        passes++;
  649.        if ((passes % 100L) == 0) {
  650.               printf("Pass %ld.\n", passes);
  651.        }
  652. #else
  653.     for (itercount = 0; itercount < niter; itercount++) {
  654. #endif
  655.  
  656.        for (paraxial = 0; paraxial <= 1; paraxial++) {
  657.  
  658.           /* Do main trace in D light */
  659.  
  660.           trace_line(4, clear_aperture / 2.0);
  661.           od_sa[paraxial][0] = object_distance;
  662.           od_sa[paraxial][1] = axis_slope_angle;
  663.        }
  664.        paraxial = FALSE;
  665.  
  666.        /* Trace marginal ray in C */
  667.  
  668.        trace_line(3, clear_aperture / 2.0);
  669.        od_cline = object_distance;
  670.  
  671.        /* Trace marginal ray in F */
  672.  
  673.        trace_line(6, clear_aperture / 2.0);
  674.        od_fline = object_distance;
  675.  
  676.        aberr_lspher = od_sa[1][0] - od_sa[0][0];
  677.        aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
  678.           (sin(od_sa[0][1]) * od_sa[0][0]);
  679.        aberr_lchrom = od_fline - od_cline;
  680.        max_lspher = sin(od_sa[0][1]);
  681.  
  682.        /* D light */
  683.  
  684.        max_lspher = 0.0000926 / (max_lspher * max_lspher);
  685.        max_osc = 0.0025;
  686.        max_lchrom = max_lspher;
  687. #ifndef ACCURACY
  688.     }
  689.  
  690.         printf("Stop the timer:\007");
  691.     gets(tbfr);
  692. #endif
  693.  
  694.     /* Now evaluate the accuracy of the results from the last ray trace */
  695.  
  696.         sprintf(outarr[0], "%15s   %21.11f  %14.11f",
  697.            "Marginal ray", od_sa[0][0], od_sa[0][1]);
  698.         sprintf(outarr[1], "%15s   %21.11f  %14.11f",
  699.            "Paraxial ray", od_sa[1][0], od_sa[1][1]);
  700.     sprintf(outarr[2],
  701.            "Longitudinal spherical aberration:      %16.11f",
  702.        aberr_lspher);
  703.     sprintf(outarr[3],
  704.            "    (Maximum permissible):              %16.11f",
  705.        max_lspher);
  706.     sprintf(outarr[4],
  707.            "Offense against sine condition (coma):  %16.11f",
  708.        aberr_osc);
  709.     sprintf(outarr[5],
  710.            "    (Maximum permissible):              %16.11f",
  711.        max_osc);
  712.     sprintf(outarr[6],
  713.            "Axial chromatic aberration:             %16.11f",
  714.        aberr_lchrom);
  715.     sprintf(outarr[7],
  716.            "    (Maximum permissible):              %16.11f",
  717.        max_lchrom);
  718.  
  719.     /* Now compare the edited results with the master values from
  720.        reference executions of this program. */
  721.  
  722.     errors = 0;
  723.     for (i = 0; i < 8; i++) {
  724.        if (strcmp(outarr[i], refarr[i]) != 0) {
  725. #ifdef ACCURACY
  726.               printf("\nError in pass %ld for results on line %d...\n",
  727.             passes, i + 1);
  728. #else
  729.               printf("\nError in results on line %d...\n", i + 1);
  730. #endif
  731.               printf("Expected:  \"%s\"\n", refarr[i]);
  732.               printf("Received:  \"%s\"\n", outarr[i]);
  733.               printf("(Errors)    ");
  734.           k = strlen(refarr[i]);
  735.           for (j = 0; j < k; j++) {
  736.                  printf("%c", refarr[i][j] == outarr[i][j] ? ' ' : '^');
  737.         if (refarr[i][j] != outarr[i][j])
  738.            errors++;
  739.           }
  740.               printf("\n");
  741.        }
  742.     }
  743. #ifdef ACCURACY
  744.     }
  745. #else
  746.     if (errors > 0) {
  747.            printf("\n%d error%s in results.  This is VERY SERIOUS.\n",
  748.               errors, errors > 1 ? "s" : "");
  749.     } else
  750.            printf("\nNo errors in results.\n");
  751. #endif
  752. }
  753.