// 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; | |

} |