mirror of
https://github.com/ioacademy-jikim/debugging
synced 2025-06-08 00:16:11 +00:00
869 lines
25 KiB
C
869 lines
25 KiB
C
// 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<67>, 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;
|
||
}
|