// This small program does some raytracing.  It tests Valgrind's handling of
// FP operations.  It apparently does a lot of trigonometry operations.

// Licensing: This program is closely based on the one of the same name from
// http://www.fourmilab.ch/.  The front page of that site says:
//
//   "Except for a few clearly-marked exceptions, all the material on this
//   site is in the public domain and may be used in any manner without
//   permission, restriction, attribution, or compensation."

/* This program can be used in two ways.  If INTRIG is undefined, sin,
   cos, tan, etc, will be used as supplied by <math.h>.  If it is
   defined, then the program calculates all this stuff from first
   principles (so to speak) and does not use the libc facilities.  For
   benchmarking purposes it seems better to avoid the libc stuff, so
   that the inner loops (sin, sqrt) present a workload independent of
   libc implementations on different platforms.  Hence: */

#define INTRIG 1


/*

        John Walker's Floating Point Benchmark, derived from...

	Marinchip Interactive Lens Design System

				     John Walker   December 1980

	By John Walker
	   http://www.fourmilab.ch/

	This  program may be used, distributed, and modified freely as
	long as the origin information is preserved.

	This  is  a  complete  optical	design	raytracing  algorithm,
	stripped of its user interface and recast into portable C.  It
	not only determines execution speed on an  extremely  floating
	point	(including   trig   function)	intensive   real-world
	application, it  checks  accuracy  on  an  algorithm  that  is
	exquisitely  sensitive	to  errors.   The  performance of this
	program is typically far more  sensitive  to  changes  in  the
	efficiency  of	the  trigonometric  library  routines than the
	average floating point program.

	The benchmark may be compiled in two  modes.   If  the	symbol
	INTRIG	is  defined,  built-in	trigonometric  and square root
	routines will be used for all calculations.  Timings made with
        INTRIG  defined  reflect  the  machine's  basic floating point
	performance for the arithmetic operators.  If  INTRIG  is  not
	defined,  the  system  library	<math.h>  functions  are used.
        Results with INTRIG not defined reflect the  system's  library
	performance  and/or  floating  point hardware support for trig
	functions and square root.  Results with INTRIG defined are  a
	good  guide  to  general  floating  point  performance,  while
	results with INTRIG undefined indicate the performance	of  an
	application which is math function intensive.

	Special  note  regarding  errors in accuracy: this program has
	generated numbers identical to the last digit it  formats  and
	checks on the following machines, floating point
	architectures, and languages:

	Marinchip 9900	  QBASIC    IBM 370 double-precision (REAL * 8) format

	IBM PC / XT / AT  Lattice C IEEE 64 bit, 80 bit temporaries
			  High C    same, in line 80x87 code
                          BASICA    "Double precision"
			  Quick BASIC IEEE double precision, software routines

	Sun 3		  C	    IEEE 64 bit, 80 bit temporaries,
				    in-line 68881 code, in-line FPA code.

        MicroVAX II       C         Vax "G" format floating point

	Macintosh Plus	  MPW C     SANE floating point, IEEE 64 bit format
				    implemented in ROM.

	Inaccuracies  reported	by  this  program should be taken VERY
	SERIOUSLY INDEED, as the program has been demonstrated	to  be
	invariant  under  changes in floating point format, as long as
	the format is a recognised double precision  format.   If  you
	encounter errors, please remember that they are just as likely
	to  be	in  the  floating  point  editing   library   or   the
	trigonometric  libraries  as  in  the low level operator code.

	The benchmark assumes that results are basically reliable, and
	only tests the last result computed against the reference.  If
        you're running on  a  suspect  system  you  can  compile  this
	program  with  ACCURACY defined.  This will generate a version
	which executes as an infinite loop, performing the  ray  trace
	and checking the results on every pass.  All incorrect results
	will be reported.

	Representative	timings  are  given  below.   All  have   been
	normalised as if run for 1000 iterations.

  Time in seconds		   Computer, Compiler, and notes
 Normal      INTRIG

 3466.00    4031.00	Commodore 128, 2 Mhz 8510 with software floating
			point.	Abacus Software/Data-Becker Super-C 128,
			version 3.00, run in fast (2 Mhz) mode.  Note:
			the results generated by this system differed
			from the reference results in the 8th to 10th
			decimal place.

 3290.00		IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
                        Run with the "/d" switch, software floating point.

 2131.50		IBM PC/AT 6 Mhz, Lattice C version 2.14, small model.
			This version of Lattice compiles subroutine
			calls which either do software floating point
			or use the 80x87.  The machine on which I ran
			this had an 80287, but the results were so bad
			I wonder if it was being used.

 1598.00		Macintosh Plus, MPW C, SANE Software floating point.

 1582.13		Marinchip 9900 2 Mhz, QBASIC compiler with software
			floating point.  This was a QBASIC version of the
			program which contained the identical algorithm.

  404.00		IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
			Software floating point.

  165.15		IBM PC/AT 6 Mhz, Metaware High C version 1.3, small
			model.	This was compiled to call subroutines for
			floating point, and the machine contained an 80287
			which was used by the subroutines.

  143.20		Macintosh II, MPW C, SANE calls.  I was unable to
			determine whether SANE was using the 68881 chip or
			not.

  121.80		Sun 3/160 16 Mhz, Sun C.  Compiled with -fsoft switch
			which executes floating point in software.

   78.78     110.11	IBM RT PC (Model 6150).  IBM AIX 1.0 C compiler
			with -O switch.

   75.2      254.0	Microsoft Quick C 1.0, in-line 8087 instructions,
			compiled with 80286 optimisation on.  (Switches
			were -Ol -FPi87-G2 -AS).  Small memory model.

   69.50		IBM PC/AT 6Mhz, Borland Turbo BASIC 1.0.  Compiled
                        in "8087 required" mode to generate in-line
			code for the math coprocessor.

   66.96		IBM PC/AT 6Mhz, Microsoft QuickBASIC 4.0.  This
			release of QuickBASIC compiles code for the
			80287 math coprocessor.

   66.36     206.35	IBM PC/AT 6Mhz, Metaware High C version 1.3, small
			model.	This was compiled with in-line code for the
			80287 math coprocessor.  Trig functions still call
			library routines.

   63.07     220.43	IBM PC/AT, 6Mhz, Borland Turbo C, in-line 8087 code,
			small model, word alignment, no stack checking,
			8086 code mode.

   17.18		Apollo DN-3000, 12 Mhz 68020 with 68881, compiled
			with in-line code for the 68881 coprocessor.
			According to Apollo, the library routines are chosen
			at runtime based on coprocessor presence.  Since the
			coprocessor was present, the library is supposed to
			use in-line floating point code.

   15.55      27.56	VAXstation II GPX.  Compiled and executed under
			VAX/VMS C.

   15.14      37.93	Macintosh II, Unix system V.  Green Hills 68020
			Unix compiler with in-line code for the 68881
			coprocessor (-O -ZI switches).

   12.69		Sun 3/160 16 Mhz, Sun C.  Compiled with -fswitch,
			which calls a subroutine to select the fastest
			floating point processor.  This was using the 68881.

   11.74      26.73	Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
			Metaware High C version 1.3, compiled with in-line
			for the math coprocessor (but not optimised for the
			80386/80387).  Trig functions still call library
			routines.

    8.43      30.49	Sun 3/160 16 Mhz, Sun C.  Compiled with -f68881,
			generating in-line MC68881 instructions.  Trig
			functions still call library routines.

    6.29      25.17	Sun 3/260 25 Mhz, Sun C.  Compiled with -f68881,
			generating in-line MC68881 instructions.  Trig
			functions still call library routines.

    4.57		Sun 3/260 25 Mhz, Sun FORTRAN 77.  Compiled with
			-O -f68881, generating in-line MC68881 instructions.
			Trig functions are compiled in-line.  This used
			the FORTRAN 77 version of the program, FBFORT77.F.

    4.00      14.20	Sun386i/25 Mhz model 250, Sun C compiler.

    4.00      14.00	Sun386i/25 Mhz model 250, Metaware C.

    3.10      12.00	Compaq 386/387 25 Mhz running SCO Xenix 2.
			Compiled with Metaware HighC 386, optimized
			for 386.

    3.00      12.00	Compaq 386/387 25MHZ optimized for 386/387.

    2.96       5.17	Sun 4/260, Sparc RISC processor.  Sun C,
			compiled with the -O2 switch for global
			optimisation.

    2.47		COMPAQ 486/25, secondary cache disabled, High C,
			486/387, inline f.p., small memory model.

    2.20       3.40	Data General Motorola 88000, 16 Mhz, Gnu C.

    1.56		COMPAQ 486/25, 128K secondary cache, High C, 486/387,
			inline f.p., small memory model.

    0.66       1.50	DEC Pmax, Mips processor.

    0.63       0.91	Sun SparcStation 2, Sun C (SunOS 4.1.1) with
                        -O4 optimisation and "/usr/lib/libm.il" inline
			floating point.

    0.60       1.07	Intel 860 RISC processor, 33 Mhz, Greenhills
			C compiler.

    0.40       0.90	Dec 3MAX, MIPS 3000 processor, -O4.

    0.31       0.90	IBM RS/6000, -O.

    0.1129     0.2119	Dell Dimension XPS P133c, Pentium 133 MHz,
			Windows 95, Microsoft Visual C 5.0.

    0.0883     0.2166	Silicon Graphics Indigo, MIPS R4400,
                        175 Mhz, "-O3".

    0.0351     0.0561	Dell Dimension XPS R100, Pentium II 400 MHz,
			Windows 98, Microsoft Visual C 5.0.

    0.0312     0.0542	Sun Ultra 2, UltraSPARC V9, 300 MHz, Solaris
			2.5.1.
			
    0.00862    0.01074  Dell Inspiron 9100, Pentium 4, 3.4 GHz, gcc -O3.

*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifndef INTRIG
#include <math.h>
#endif

#define cot(x) (1.0 / tan(x))

#define TRUE  1
#define FALSE 0

#define max_surfaces 10

/*  Local variables  */

/* static char tbfr[132]; */

static short current_surfaces;
static short paraxial;

static double clear_aperture;

static double aberr_lspher;
static double aberr_osc;
static double aberr_lchrom;

static double max_lspher;
static double max_osc;
static double max_lchrom;

static double radius_of_curvature;
static double object_distance;
static double ray_height;
static double axis_slope_angle;
static double from_index;
static double to_index;

static double spectral_line[9];
static double s[max_surfaces][5];
static double od_sa[2][2];

static char outarr[8][80];	   /* Computed output of program goes here */

int itercount;			   /* The iteration counter for the main loop
				      in the program is made global so that
				      the compiler should not be allowed to
				      optimise out the loop over the ray
				      tracing code. */

#ifndef ITERATIONS
#define ITERATIONS /*1000*/ /*500000*/ 80000
#endif
int niter = ITERATIONS; 	   /* Iteration counter */

static char *refarr[] = {	   /* Reference results.  These happen to
				      be derived from a run on Microsoft 
				      Quick BASIC on the IBM PC/AT. */

        "   Marginal ray          47.09479120920   0.04178472683",
        "   Paraxial ray          47.08372160249   0.04177864821",
        "Longitudinal spherical aberration:        -0.01106960671",
        "    (Maximum permissible):                 0.05306749907",
        "Offense against sine condition (coma):     0.00008954761",
        "    (Maximum permissible):                 0.00250000000",
        "Axial chromatic aberration:                0.00448229032",
        "    (Maximum permissible):                 0.05306749907"
};

/* The	test  case  used  in  this program is the  design for a 4 inch
   achromatic telescope  objective  used  as  the  example  in  Wyld's
   classic  work  on  ray  tracing by hand, given in Amateur Telescope
   Making, Volume 3.  */

static double testcase[4][4] = {
	{27.05, 1.5137, 63.6, 0.52},
	{-16.68, 1, 0, 0.138},
	{-16.68, 1.6164, 36.7, 0.38},
	{-78.1, 1, 0, 0}
};

/*  Internal trig functions (used only if INTRIG is  defined).	 These
    standard  functions  may be enabled to obtain timings that reflect
    the machine's floating point performance rather than the speed  of
    its trig function evaluation.  */

#ifdef INTRIG

/*  The following definitions should keep you from getting intro trouble
    with compilers which don't let you redefine intrinsic functions.  */

#define sin I_sin
#define cos I_cos
#define tan I_tan
#define sqrt I_sqrt
#define atan I_atan
#define atan2 I_atan2
#define asin I_asin

#define fabs(x)  ((x < 0.0) ? -x : x)

#define pic 3.1415926535897932

/*  Commonly used constants  */

static double pi = pic,
	twopi =pic * 2.0,
	piover4 = pic / 4.0,
	fouroverpi = 4.0 / pic,
	piover2 = pic / 2.0;

/*  Coefficients for ATAN evaluation  */

static double atanc[] = {
	0.0,
	0.4636476090008061165,
	0.7853981633974483094,
	0.98279372324732906714,
	1.1071487177940905022,
	1.1902899496825317322,
	1.2490457723982544262,
	1.2924966677897852673,
	1.3258176636680324644
};

/*  aint(x)	  Return integer part of number.  Truncates towards 0	 */

double aint(x)
double x;
{
	long l;

	/*  Note that this routine cannot handle the full floating point
	    number range.  This function should be in the machine-dependent
	    floating point library!  */

	l = x;
	if ((int)(-0.5) != 0  &&  l < 0 )
	   l++;
	x = l;
	return x;
}

/*  sin(x)	  Return sine, x in radians  */

static double sin(x)
double x;
{
	int sign;
	double y, r, z;

	x = (((sign= (x < 0.0)) != 0) ? -x: x);

	if (x > twopi)
	   x -= (aint(x / twopi) * twopi);

	if (x > pi) {
	   x -= pi;
	   sign = !sign;
	}

	if (x > piover2)
	   x = pi - x;

	if (x < piover4) {
	   y = x * fouroverpi;
	   z = y * y;
	   r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z -
	      0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z -
	      0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z -
	      0.0807455121882807815) * z + 0.785398163397448310);
	} else {
	   y = (piover2 - x) * fouroverpi;
	   z = y * y;
	   r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z -
	      0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z -
	      0.325991886926687550E-3) * z + 0.0158543442438154109) * z -
	      0.308425137534042452) * z + 1.0;
	}
	return sign ? -r : r;
}

/*  cos(x)	  Return cosine, x in radians, by identity  */

static double cos(x)
double x;
{
	x = (x < 0.0) ? -x : x;
	if (x > twopi)		      /* Do range reduction here to limit */
	   x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2    */
	return sin(x + piover2);
}

/*  tan(x)	  Return tangent, x in radians, by identity  */

static double tan(x)
double x;
{
	return sin(x) / cos(x);
}

/*  sqrt(x)	  Return square root.  Initial guess, then Newton-
		  Raphson refinement  */

double sqrt(x)
double x;
{
	double c, cl, y;
	int n;

	if (x == 0.0)
	   return 0.0;

	if (x < 0.0) {
	   fprintf(stderr,
              "\nGood work!  You tried to take the square root of %g",
	     x);
	   fprintf(stderr,
              "\nunfortunately, that is too complex for me to handle.\n");
	   exit(1);
	}

	y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x);

	c = (y - x / y) / 2.0;
	cl = 0.0;
	for (n = 50; c != cl && n--;) {
	   y = y - c;
	   cl = c;
	   c = (y - x / y) / 2.0;
	}
	return y;
}

/*  atan(x)	  Return arctangent in radians,
		  range -pi/2 to pi/2  */

static double atan(x)
double x;
{
	int sign, l, y;
	double a, b, z;

	x = (((sign = (x < 0.0)) != 0) ? -x : x);
	l = 0;

	if (x >= 4.0) {
	   l = -1;
	   x = 1.0 / x;
	   y = 0;
	   goto atl;
	} else {
	   if (x < 0.25) {
	      y = 0;
	      goto atl;
	   }
	}

	y = aint(x / 0.5);
	z = y * 0.5;
	x = (x - z) / (x * z + 1);

atl:
	z = x * x;
	b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z +
	    1277025750.0) * z + 1550674125.0) * z + 654729075.0;
	a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z +
	    1332431100.0) * z + 654729075.0;
	a = (a / b) * x + atanc[y];
	if (l)
	   a=piover2 - a;
	return sign ? -a : a;
}

/*  atan2(y,x)	  Return arctangent in radians of y/x,
		  range -pi to pi  */

static double atan2(y, x)
double y, x;
{
	double temp;

	if (x == 0.0) {
	   if (y == 0.0)   /*  Special case: atan2(0,0) = 0  */
	      return 0.0;
	   else if (y > 0)
	      return piover2;
	   else
	      return -piover2;
	}
	temp = atan(y / x);
	if (x < 0.0) {
	   if (y >= 0.0)
	      temp += pic;
	   else
	      temp -= pic;
	}
	return temp;
}

/*  asin(x)	  Return arcsine in radians of x  */

static double asin(x)
double x;
{
	if (fabs(x)>1.0) {
	   fprintf(stderr,
              "\nInverse trig functions lose much of their gloss when");
	   fprintf(stderr,
              "\ntheir arguments are greater than 1, such as the");
	   fprintf(stderr,
              "\nvalue %g you passed.\n", x);
	   exit(1);
	}
	return atan2(x, sqrt(1 - x * x));
}
#endif

/*	      Calculate passage through surface

	      If  the variable PARAXIAL is true, the trace through the
	      surface will be done using the paraxial  approximations.
	      Otherwise,  the normal trigonometric trace will be done.

	      This routine takes the following 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 variables:

	      OBJECT_DISTANCE		  Distance from vertex to object focus
					  after refraction.

	      AXIS_SLOPE_ANGLE		  Angle incoming ray makes with axis
					  at intercept after refraction.

*/

static void transit_surface() {
	double iang,		   /* Incidence angle */
	       rang,		   /* Refraction angle */
	       iang_sin,	   /* Incidence angle sin */
	       rang_sin,	   /* Refraction angle sin */
	       old_axis_slope_angle, sagitta;

	if (paraxial) {
	   if (radius_of_curvature != 0.0) {
	      if (object_distance == 0.0) {
		 axis_slope_angle = 0.0;
		 iang_sin = ray_height / radius_of_curvature;
	      } else
		 iang_sin = ((object_distance -
		    radius_of_curvature) / radius_of_curvature) *
		    axis_slope_angle;

	      rang_sin = (from_index / to_index) *
		 iang_sin;
	      old_axis_slope_angle = axis_slope_angle;
	      axis_slope_angle = axis_slope_angle +
		 iang_sin - rang_sin;
	      if (object_distance != 0.0)
		 ray_height = object_distance * old_axis_slope_angle;
	      object_distance = ray_height / axis_slope_angle;
	      return;
	   }
	   object_distance = object_distance * (to_index / from_index);
	   axis_slope_angle = axis_slope_angle * (from_index / to_index);
	   return;
	}

	if (radius_of_curvature != 0.0) {
	   if (object_distance == 0.0) {
	      axis_slope_angle = 0.0;
	      iang_sin = ray_height / radius_of_curvature;
	   } else {
	      iang_sin = ((object_distance -
		 radius_of_curvature) / radius_of_curvature) *
		 sin(axis_slope_angle);
	   }
	   iang = asin(iang_sin);
	   rang_sin = (from_index / to_index) *
	      iang_sin;
	   old_axis_slope_angle = axis_slope_angle;
	   axis_slope_angle = axis_slope_angle +
	      iang - asin(rang_sin);
	   sagitta = sin((old_axis_slope_angle + iang) / 2.0);
	   sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
	   object_distance = ((radius_of_curvature * sin(
	      old_axis_slope_angle + iang)) *
	      cot(axis_slope_angle)) + sagitta;
	   return;
	}

	rang = -asin((from_index / to_index) *
	   sin(axis_slope_angle));
	object_distance = object_distance * ((to_index *
	   cos(-rang)) / (from_index *
	   cos(axis_slope_angle)));
	axis_slope_angle = -rang;
}

/*  Perform ray trace in specific spectral line  */

static void trace_line(line, ray_h)
int line;
double ray_h;
{
	int i;

	object_distance = 0.0;
	ray_height = ray_h;
	from_index = 1.0;

	for (i = 1; i <= current_surfaces; i++) {
	   radius_of_curvature = s[i][1];
	   to_index = s[i][2];
	   if (to_index > 1.0)
	      to_index = to_index + ((spectral_line[4] -
		 spectral_line[line]) /
		 (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
		 s[i][3]);
	   transit_surface();
	   from_index = to_index;
	   if (i < current_surfaces)
	      object_distance = object_distance - s[i][4];
	}
}

/*  Initialise when called the first time  */

int main(argc, argv)
int argc;
char *argv[];
{
	int i, j, k, errors;
	double od_fline, od_cline;
#ifdef ACCURACY
	long passes;
#endif

	spectral_line[1] = 7621.0;	 /* A */
	spectral_line[2] = 6869.955;	 /* B */
	spectral_line[3] = 6562.816;	 /* C */
	spectral_line[4] = 5895.944;	 /* D */
	spectral_line[5] = 5269.557;	 /* E */
	spectral_line[6] = 4861.344;	 /* F */
        spectral_line[7] = 4340.477;     /* G'*/
	spectral_line[8] = 3968.494;	 /* H */

	/* Process the number of iterations argument, if one is supplied. */

	if (argc > 1) {
	   niter = atoi(argv[1]);
           if (*argv[1] == '-' || niter < 1) {
              printf("This is John Walker's floating point accuracy and\n");
              printf("performance benchmark program.  You call it with\n");
              printf("\nfbench <itercount>\n\n");
              printf("where <itercount> is the number of iterations\n");
              printf("to be executed.  Archival timings should be made\n");
              printf("with the iteration count set so that roughly five\n");
              printf("minutes of execution is timed.\n");
	      exit(0);
	   }
	}

	/* Load test case into working array */

	clear_aperture = 4.0;
	current_surfaces = 4;
	for (i = 0; i < current_surfaces; i++)
	   for (j = 0; j < 4; j++)
	      s[i + 1][j + 1] = testcase[i][j];

#ifdef ACCURACY
        printf("Beginning execution of floating point accuracy test...\n");
	passes = 0;
#else
        printf("Ready to begin John Walker's floating point accuracy\n");
        printf("and performance benchmark.  %d iterations will be made.\n\n",
	   niter);

        printf("\nMeasured run time in seconds should be divided by %.f\n", niter / 1000.0);
        printf("to normalise for reporting results.  For archival results,\n");
        printf("adjust iteration count so the benchmark runs about five minutes.\n\n");

        //printf("Press return to begin benchmark:");
	//gets(tbfr);
#endif

	/* Perform ray trace the specified number of times. */

#ifdef ACCURACY
	while (TRUE) {
	   passes++;
	   if ((passes % 100L) == 0) {
              printf("Pass %ld.\n", passes);
	   }
#else
	for (itercount = 0; itercount < niter; itercount++) {
#endif

	   for (paraxial = 0; paraxial <= 1; paraxial++) {

	      /* Do main trace in D light */

	      trace_line(4, clear_aperture / 2.0);
	      od_sa[paraxial][0] = object_distance;
	      od_sa[paraxial][1] = axis_slope_angle;
	   }
	   paraxial = FALSE;

	   /* Trace marginal ray in C */

	   trace_line(3, clear_aperture / 2.0);
	   od_cline = object_distance;

	   /* Trace marginal ray in F */

	   trace_line(6, clear_aperture / 2.0);
	   od_fline = object_distance;

	   aberr_lspher = od_sa[1][0] - od_sa[0][0];
	   aberr_osc = 1.0 - (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;
#ifndef ACCURACY
	}

        //printf("Stop the timer:\007");
	//gets(tbfr);
#endif

	/* Now evaluate the accuracy of the results from the last ray trace */

        sprintf(outarr[0], "%15s   %21.11f  %14.11f",
           "Marginal ray", od_sa[0][0], od_sa[0][1]);
        sprintf(outarr[1], "%15s   %21.11f  %14.11f",
           "Paraxial ray", od_sa[1][0], od_sa[1][1]);
	sprintf(outarr[2],
           "Longitudinal spherical aberration:      %16.11f",
	   aberr_lspher);
	sprintf(outarr[3],
           "    (Maximum permissible):              %16.11f",
	   max_lspher);
	sprintf(outarr[4],
           "Offense against sine condition (coma):  %16.11f",
	   aberr_osc);
	sprintf(outarr[5],
           "    (Maximum permissible):              %16.11f",
	   max_osc);
	sprintf(outarr[6],
           "Axial chromatic aberration:             %16.11f",
	   aberr_lchrom);
	sprintf(outarr[7],
           "    (Maximum permissible):              %16.11f",
	   max_lchrom);

	/* Now compare the edited results with the master values from
	   reference executions of this program. */

	errors = 0;
	for (i = 0; i < 8; i++) {
	   if (strcmp(outarr[i], refarr[i]) != 0) {
#ifdef ACCURACY
              printf("\nError in pass %ld for results on line %d...\n",
		     passes, i + 1);
#else
              printf("\nError in results on line %d...\n", i + 1);
#endif
              printf("Expected:  \"%s\"\n", refarr[i]);
              printf("Received:  \"%s\"\n", outarr[i]);
              printf("(Errors)    ");
	      k = strlen(refarr[i]);
	      for (j = 0; j < k; j++) {
                 printf("%c", refarr[i][j] == outarr[i][j] ? ' ' : '^');
		 if (refarr[i][j] != outarr[i][j])
		    errors++;
	      }
              printf("\n");
	   }
	}
#ifdef ACCURACY
	}
#else
	if (errors > 0) {
           printf("\n%d error%s in results.  This is VERY SERIOUS.\n",
              errors, errors > 1 ? "s" : "");
	} else
           printf("\nNo errors in results.\n");
#endif
	return 0;
}