1ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// This small program does some raytracing. It tests Valgrind's handling of 2ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// FP operations. It apparently does a lot of trigonometry operations. 3ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 4ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// Licensing: This program is closely based on the one of the same name from 5ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// http://www.fourmilab.ch/. The front page of that site says: 6ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// 7ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// "Except for a few clearly-marked exceptions, all the material on this 8ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// site is in the public domain and may be used in any manner without 9ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown// permission, restriction, attribution, or compensation." 10ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 11ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* This program can be used in two ways. If INTRIG is undefined, sin, 12ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cos, tan, etc, will be used as supplied by <math.h>. If it is 13ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown defined, then the program calculates all this stuff from first 14ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown principles (so to speak) and does not use the libc facilities. For 15ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown benchmarking purposes it seems better to avoid the libc stuff, so 16ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown that the inner loops (sin, sqrt) present a workload independent of 17ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown libc implementations on different platforms. Hence: */ 18ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 19ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define INTRIG 1 20ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 21ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 22ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* 23ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 24ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown John Walker's Floating Point Benchmark, derived from... 25ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 26ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Marinchip Interactive Lens Design System 27ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 28ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown John Walker December 1980 29ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 30ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown By John Walker 31ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown http://www.fourmilab.ch/ 32ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 33ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown This program may be used, distributed, and modified freely as 34ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown long as the origin information is preserved. 35ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 36ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown This is a complete optical design raytracing algorithm, 37ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown stripped of its user interface and recast into portable C. It 38ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown not only determines execution speed on an extremely floating 39ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown point (including trig function) intensive real-world 40ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown application, it checks accuracy on an algorithm that is 41ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown exquisitely sensitive to errors. The performance of this 42ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown program is typically far more sensitive to changes in the 43ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown efficiency of the trigonometric library routines than the 44ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown average floating point program. 45ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 46ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown The benchmark may be compiled in two modes. If the symbol 47ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown INTRIG is defined, built-in trigonometric and square root 48ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown routines will be used for all calculations. Timings made with 49ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown INTRIG defined reflect the machine's basic floating point 50ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown performance for the arithmetic operators. If INTRIG is not 51ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown defined, the system library <math.h> functions are used. 52ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Results with INTRIG not defined reflect the system's library 53ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown performance and/or floating point hardware support for trig 54ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown functions and square root. Results with INTRIG defined are a 55ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown good guide to general floating point performance, while 56ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown results with INTRIG undefined indicate the performance of an 57ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown application which is math function intensive. 58ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 59ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Special note regarding errors in accuracy: this program has 60ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown generated numbers identical to the last digit it formats and 61ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown checks on the following machines, floating point 62ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown architectures, and languages: 63ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 64ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Marinchip 9900 QBASIC IBM 370 double-precision (REAL * 8) format 65ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 66ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown IBM PC / XT / AT Lattice C IEEE 64 bit, 80 bit temporaries 67ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown High C same, in line 80x87 code 68ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown BASICA "Double precision" 69ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Quick BASIC IEEE double precision, software routines 70ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 71ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Sun 3 C IEEE 64 bit, 80 bit temporaries, 72ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown in-line 68881 code, in-line FPA code. 73ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 74ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown MicroVAX II C Vax "G" format floating point 75ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 76ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Macintosh Plus MPW C SANE floating point, IEEE 64 bit format 77ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown implemented in ROM. 78ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 79ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Inaccuracies reported by this program should be taken VERY 80ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown SERIOUSLY INDEED, as the program has been demonstrated to be 81ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown invariant under changes in floating point format, as long as 82ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the format is a recognised double precision format. If you 83ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown encounter errors, please remember that they are just as likely 84ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown to be in the floating point editing library or the 85ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown trigonometric libraries as in the low level operator code. 86ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 87ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown The benchmark assumes that results are basically reliable, and 88ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown only tests the last result computed against the reference. If 89ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown you're running on a suspect system you can compile this 90ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown program with ACCURACY defined. This will generate a version 91ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown which executes as an infinite loop, performing the ray trace 92ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown and checking the results on every pass. All incorrect results 93ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown will be reported. 94ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 95ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Representative timings are given below. All have been 96ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown normalised as if run for 1000 iterations. 97ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 98ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Time in seconds Computer, Compiler, and notes 99ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Normal INTRIG 100ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 101ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 3466.00 4031.00 Commodore 128, 2 Mhz 8510 with software floating 102ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown point. Abacus Software/Data-Becker Super-C 128, 103ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown version 3.00, run in fast (2 Mhz) mode. Note: 104ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the results generated by this system differed 105ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown from the reference results in the 8th to 10th 106ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown decimal place. 107ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 108ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 3290.00 IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00. 109ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Run with the "/d" switch, software floating point. 110ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 111ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 2131.50 IBM PC/AT 6 Mhz, Lattice C version 2.14, small model. 112ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown This version of Lattice compiles subroutine 113ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown calls which either do software floating point 114ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown or use the 80x87. The machine on which I ran 115ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown this had an 80287, but the results were so bad 116ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown I wonder if it was being used. 117ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 118ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1598.00 Macintosh Plus, MPW C, SANE Software floating point. 119ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 120ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1582.13 Marinchip 9900 2 Mhz, QBASIC compiler with software 121ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown floating point. This was a QBASIC version of the 122ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown program which contained the identical algorithm. 123ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 124ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 404.00 IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0. 125ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Software floating point. 126ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 127ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 165.15 IBM PC/AT 6 Mhz, Metaware High C version 1.3, small 128ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown model. This was compiled to call subroutines for 129ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown floating point, and the machine contained an 80287 130ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown which was used by the subroutines. 131ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 132ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 143.20 Macintosh II, MPW C, SANE calls. I was unable to 133ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown determine whether SANE was using the 68881 chip or 134ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown not. 135ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 136ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 121.80 Sun 3/160 16 Mhz, Sun C. Compiled with -fsoft switch 137ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown which executes floating point in software. 138ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 139ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 78.78 110.11 IBM RT PC (Model 6150). IBM AIX 1.0 C compiler 140ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown with -O switch. 141ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 142ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 75.2 254.0 Microsoft Quick C 1.0, in-line 8087 instructions, 143ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown compiled with 80286 optimisation on. (Switches 144ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown were -Ol -FPi87-G2 -AS). Small memory model. 145ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 146ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 69.50 IBM PC/AT 6Mhz, Borland Turbo BASIC 1.0. Compiled 147ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown in "8087 required" mode to generate in-line 148ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown code for the math coprocessor. 149ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 150ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 66.96 IBM PC/AT 6Mhz, Microsoft QuickBASIC 4.0. This 151ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown release of QuickBASIC compiles code for the 152ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 80287 math coprocessor. 153ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 154ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 66.36 206.35 IBM PC/AT 6Mhz, Metaware High C version 1.3, small 155ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown model. This was compiled with in-line code for the 156ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 80287 math coprocessor. Trig functions still call 157ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown library routines. 158ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 159ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 63.07 220.43 IBM PC/AT, 6Mhz, Borland Turbo C, in-line 8087 code, 160ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown small model, word alignment, no stack checking, 161ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 8086 code mode. 162ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 163ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 17.18 Apollo DN-3000, 12 Mhz 68020 with 68881, compiled 164ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown with in-line code for the 68881 coprocessor. 165ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown According to Apollo, the library routines are chosen 166ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown at runtime based on coprocessor presence. Since the 167ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown coprocessor was present, the library is supposed to 168ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown use in-line floating point code. 169ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 170ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 15.55 27.56 VAXstation II GPX. Compiled and executed under 171ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown VAX/VMS C. 172ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 173ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 15.14 37.93 Macintosh II, Unix system V. Green Hills 68020 174ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Unix compiler with in-line code for the 68881 175ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown coprocessor (-O -ZI switches). 176ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 177ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 12.69 Sun 3/160 16 Mhz, Sun C. Compiled with -fswitch, 178ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown which calls a subroutine to select the fastest 179ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown floating point processor. This was using the 68881. 180ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 181ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 11.74 26.73 Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387. 182ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Metaware High C version 1.3, compiled with in-line 183ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for the math coprocessor (but not optimised for the 184ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 80386/80387). Trig functions still call library 185ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown routines. 186ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 187ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 8.43 30.49 Sun 3/160 16 Mhz, Sun C. Compiled with -f68881, 188ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown generating in-line MC68881 instructions. Trig 189ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown functions still call library routines. 190ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 191ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 6.29 25.17 Sun 3/260 25 Mhz, Sun C. Compiled with -f68881, 192ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown generating in-line MC68881 instructions. Trig 193ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown functions still call library routines. 194ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 195ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 4.57 Sun 3/260 25 Mhz, Sun FORTRAN 77. Compiled with 196ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown -O -f68881, generating in-line MC68881 instructions. 197ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Trig functions are compiled in-line. This used 198ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the FORTRAN 77 version of the program, FBFORT77.F. 199ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 200ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 4.00 14.20 Sun386i/25 Mhz model 250, Sun C compiler. 201ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 202ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 4.00 14.00 Sun386i/25 Mhz model 250, Metaware C. 203ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 204ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 3.10 12.00 Compaq 386/387 25 Mhz running SCO Xenix 2. 205ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Compiled with Metaware HighC 386, optimized 206ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for 386. 207ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 208ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 3.00 12.00 Compaq 386/387 25MHZ optimized for 386/387. 209ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 210ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 2.96 5.17 Sun 4/260, Sparc RISC processor. Sun C, 211ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown compiled with the -O2 switch for global 212ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown optimisation. 213ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 214ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 2.47 COMPAQ 486/25, secondary cache disabled, High C, 215ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 486/387, inline f.p., small memory model. 216ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 217ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 2.20 3.40 Data General Motorola 88000, 16 Mhz, Gnu C. 218ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 219ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.56 COMPAQ 486/25, 128K secondary cache, High C, 486/387, 220ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown inline f.p., small memory model. 221ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 222ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.66 1.50 DEC Pmax, Mips processor. 223ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 224ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.63 0.91 Sun SparcStation 2, Sun C (SunOS 4.1.1) with 225ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown -O4 optimisation and "/usr/lib/libm.il" inline 226ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown floating point. 227ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 228ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.60 1.07 Intel 860 RISC processor, 33 Mhz, Greenhills 229ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown C compiler. 230ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 231ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.40 0.90 Dec 3MAX, MIPS 3000 processor, -O4. 232ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 233ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.31 0.90 IBM RS/6000, -O. 234ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 235ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.1129 0.2119 Dell Dimension XPS P133c, Pentium 133 MHz, 236ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Windows 95, Microsoft Visual C 5.0. 237ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 238ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.0883 0.2166 Silicon Graphics Indigo�, MIPS R4400, 239ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 175 Mhz, "-O3". 240ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 241ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.0351 0.0561 Dell Dimension XPS R100, Pentium II 400 MHz, 242ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Windows 98, Microsoft Visual C 5.0. 243ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 244ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.0312 0.0542 Sun Ultra 2, UltraSPARC V9, 300 MHz, Solaris 245ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 2.5.1. 246ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 247ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.00862 0.01074 Dell Inspiron 9100, Pentium 4, 3.4 GHz, gcc -O3. 248ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 249ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown*/ 250ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 251ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 252ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#include <stdio.h> 253ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#include <stdlib.h> 254ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#include <string.h> 255ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifndef INTRIG 256ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#include <math.h> 257ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 258ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 259ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define cot(x) (1.0 / tan(x)) 260ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 261ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define TRUE 1 262ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define FALSE 0 263ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 264ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define max_surfaces 10 265ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 266ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Local variables */ 267ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 268ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* static char tbfr[132]; */ 269ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 270ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic short current_surfaces; 271ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic short paraxial; 272ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 273ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double clear_aperture; 274ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 275ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double aberr_lspher; 276ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double aberr_osc; 277ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double aberr_lchrom; 278ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 279ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double max_lspher; 280ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double max_osc; 281ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double max_lchrom; 282ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 283ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double radius_of_curvature; 284ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double object_distance; 285ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double ray_height; 286ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double axis_slope_angle; 287ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double from_index; 288ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double to_index; 289ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 290ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double spectral_line[9]; 291ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double s[max_surfaces][5]; 292ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double od_sa[2][2]; 293ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 294ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic char outarr[8][80]; /* Computed output of program goes here */ 295ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 296ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownint itercount; /* The iteration counter for the main loop 297ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown in the program is made global so that 298ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the compiler should not be allowed to 299ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown optimise out the loop over the ray 300ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown tracing code. */ 301ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 302ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifndef ITERATIONS 303ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define ITERATIONS /*1000*/ /*500000*/ 80000 304ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 305ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownint niter = ITERATIONS; /* Iteration counter */ 306ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 307ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic char *refarr[] = { /* Reference results. These happen to 308ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown be derived from a run on Microsoft 309ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Quick BASIC on the IBM PC/AT. */ 310ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 311ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " Marginal ray 47.09479120920 0.04178472683", 312ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " Paraxial ray 47.08372160249 0.04177864821", 313ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Longitudinal spherical aberration: -0.01106960671", 314ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): 0.05306749907", 315ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Offense against sine condition (coma): 0.00008954761", 316ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): 0.00250000000", 317ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Axial chromatic aberration: 0.00448229032", 318ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): 0.05306749907" 319ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown}; 320ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 321ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* The test case used in this program is the design for a 4 inch 322ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown achromatic telescope objective used as the example in Wyld's 323ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown classic work on ray tracing by hand, given in Amateur Telescope 324ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Making, Volume 3. */ 325ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 326ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double testcase[4][4] = { 327ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown {27.05, 1.5137, 63.6, 0.52}, 328ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown {-16.68, 1, 0, 0.138}, 329ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown {-16.68, 1.6164, 36.7, 0.38}, 330ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown {-78.1, 1, 0, 0} 331ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown}; 332ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 333ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Internal trig functions (used only if INTRIG is defined). These 334ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown standard functions may be enabled to obtain timings that reflect 335ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the machine's floating point performance rather than the speed of 336ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown its trig function evaluation. */ 337ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 338ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef INTRIG 339ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 340ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* The following definitions should keep you from getting intro trouble 341ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown with compilers which don't let you redefine intrinsic functions. */ 342ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 343ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define sin I_sin 344ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define cos I_cos 345ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define tan I_tan 346ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define sqrt I_sqrt 347ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define atan I_atan 348ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define atan2 I_atan2 349ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define asin I_asin 350ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 351ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define fabs(x) ((x < 0.0) ? -x : x) 352ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 353ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#define pic 3.1415926535897932 354ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 355ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Commonly used constants */ 356ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 357ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double pi = pic, 358ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown twopi =pic * 2.0, 359ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown piover4 = pic / 4.0, 360ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fouroverpi = 4.0 / pic, 361ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown piover2 = pic / 2.0; 362ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 363ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Coefficients for ATAN evaluation */ 364ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 365ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double atanc[] = { 366ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.0, 367ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.4636476090008061165, 368ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.7853981633974483094, 369ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.98279372324732906714, 370ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.1071487177940905022, 371ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.1902899496825317322, 372ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.2490457723982544262, 373ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.2924966677897852673, 374ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1.3258176636680324644 375ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown}; 376ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 377ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* aint(x) Return integer part of number. Truncates towards 0 */ 378ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 379ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble aint(x) 380ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 381ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 382ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown long l; 383ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 384ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Note that this routine cannot handle the full floating point 385ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown number range. This function should be in the machine-dependent 386ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown floating point library! */ 387ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 388ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown l = x; 389ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if ((int)(-0.5) != 0 && l < 0 ) 390ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown l++; 391ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = l; 392ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return x; 393ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 394ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 395ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* sin(x) Return sine, x in radians */ 396ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 397ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double sin(x) 398ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 399ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 400ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown int sign; 401ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double y, r, z; 402ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 403ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = (((sign= (x < 0.0)) != 0) ? -x: x); 404ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 405ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x > twopi) 406ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x -= (aint(x / twopi) * twopi); 407ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 408ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x > pi) { 409ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x -= pi; 410ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sign = !sign; 411ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 412ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 413ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x > piover2) 414ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = pi - x; 415ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 416ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x < piover4) { 417ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = x * fouroverpi; 418ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown z = y * y; 419ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z - 420ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z - 421ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z - 422ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.0807455121882807815) * z + 0.785398163397448310); 423ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } else { 424ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = (piover2 - x) * fouroverpi; 425ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown z = y * y; 426ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z - 427ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z - 428ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.325991886926687550E-3) * z + 0.0158543442438154109) * z - 429ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 0.308425137534042452) * z + 1.0; 430ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 431ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return sign ? -r : r; 432ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 433ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 434ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* cos(x) Return cosine, x in radians, by identity */ 435ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 436ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double cos(x) 437ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 438ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 439ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = (x < 0.0) ? -x : x; 440ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x > twopi) /* Do range reduction here to limit */ 441ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2 */ 442ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return sin(x + piover2); 443ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 444ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 445ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* tan(x) Return tangent, x in radians, by identity */ 446ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 447ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double tan(x) 448ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 449ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 450ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return sin(x) / cos(x); 451ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 452ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 453ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* sqrt(x) Return square root. Initial guess, then Newton- 454ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Raphson refinement */ 455ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 456ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble sqrt(x) 457ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 458ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 459ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double c, cl, y; 460ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown int n; 461ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 462ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x == 0.0) 463ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return 0.0; 464ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 465ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x < 0.0) { 466ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fprintf(stderr, 467ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "\nGood work! You tried to take the square root of %g", 468ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x); 469ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fprintf(stderr, 470ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "\nunfortunately, that is too complex for me to handle.\n"); 471ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown exit(1); 472ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 473ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 474ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x); 475ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 476ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown c = (y - x / y) / 2.0; 477ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cl = 0.0; 478ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (n = 50; c != cl && n--;) { 479ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = y - c; 480ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cl = c; 481ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown c = (y - x / y) / 2.0; 482ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 483ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return y; 484ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 485ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 486ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* atan(x) Return arctangent in radians, 487ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown range -pi/2 to pi/2 */ 488ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 489ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double atan(x) 490ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 491ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 492ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown int sign, l, y; 493ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double a, b, z; 494ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 495ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = (((sign = (x < 0.0)) != 0) ? -x : x); 496ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown l = 0; 497ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 498ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x >= 4.0) { 499ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown l = -1; 500ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = 1.0 / x; 501ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = 0; 502ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown goto atl; 503ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } else { 504ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x < 0.25) { 505ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = 0; 506ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown goto atl; 507ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 508ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 509ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 510ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown y = aint(x / 0.5); 511ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown z = y * 0.5; 512ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown x = (x - z) / (x * z + 1); 513ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 514ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownatl: 515ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown z = x * x; 516ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z + 517ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1277025750.0) * z + 1550674125.0) * z + 654729075.0; 518ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z + 519ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 1332431100.0) * z + 654729075.0; 520ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown a = (a / b) * x + atanc[y]; 521ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (l) 522ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown a=piover2 - a; 523ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return sign ? -a : a; 524ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 525ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 526ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* atan2(y,x) Return arctangent in radians of y/x, 527ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown range -pi to pi */ 528ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 529ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double atan2(y, x) 530ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble y, x; 531ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 532ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double temp; 533ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 534ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x == 0.0) { 535ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (y == 0.0) /* Special case: atan2(0,0) = 0 */ 536ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return 0.0; 537ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown else if (y > 0) 538ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return piover2; 539ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown else 540ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return -piover2; 541ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 542ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown temp = atan(y / x); 543ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (x < 0.0) { 544ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (y >= 0.0) 545ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown temp += pic; 546ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown else 547ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown temp -= pic; 548ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 549ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return temp; 550ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 551ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 552ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* asin(x) Return arcsine in radians of x */ 553ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 554ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic double asin(x) 555ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble x; 556ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 557ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (fabs(x)>1.0) { 558ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fprintf(stderr, 559ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "\nInverse trig functions lose much of their gloss when"); 560ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fprintf(stderr, 561ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "\ntheir arguments are greater than 1, such as the"); 562ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown fprintf(stderr, 563ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "\nvalue %g you passed.\n", x); 564ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown exit(1); 565ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 566ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return atan2(x, sqrt(1 - x * x)); 567ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 568ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 569ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 570ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Calculate passage through surface 571ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 572ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown If the variable PARAXIAL is true, the trace through the 573ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown surface will be done using the paraxial approximations. 574ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown Otherwise, the normal trigonometric trace will be done. 575ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 576ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown This routine takes the following inputs: 577ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 578ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown RADIUS_OF_CURVATURE Radius of curvature of surface 579ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown being crossed. If 0, surface is 580ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown plane. 581ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 582ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown OBJECT_DISTANCE Distance of object focus from 583ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown lens vertex. If 0, incoming 584ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rays are parallel and 585ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown the following must be specified: 586ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 587ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown RAY_HEIGHT Height of ray from axis. Only 588ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown relevant if OBJECT.DISTANCE == 0 589ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 590ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown AXIS_SLOPE_ANGLE Angle incoming ray makes with axis 591ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown at intercept 592ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 593ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown FROM_INDEX Refractive index of medium being left 594ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 595ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown TO_INDEX Refractive index of medium being 596ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown entered. 597ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 598ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown The outputs are the following variables: 599ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 600ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown OBJECT_DISTANCE Distance from vertex to object focus 601ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown after refraction. 602ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 603ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown AXIS_SLOPE_ANGLE Angle incoming ray makes with axis 604ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown at intercept after refraction. 605ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 606ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown*/ 607ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 608ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic void transit_surface() { 609ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double iang, /* Incidence angle */ 610ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rang, /* Refraction angle */ 611ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin, /* Incidence angle sin */ 612ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rang_sin, /* Refraction angle sin */ 613ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown old_axis_slope_angle, sagitta; 614ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 615ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (paraxial) { 616ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (radius_of_curvature != 0.0) { 617ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (object_distance == 0.0) { 618ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = 0.0; 619ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin = ray_height / radius_of_curvature; 620ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } else 621ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin = ((object_distance - 622ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown radius_of_curvature) / radius_of_curvature) * 623ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle; 624ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 625ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rang_sin = (from_index / to_index) * 626ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin; 627ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown old_axis_slope_angle = axis_slope_angle; 628ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = axis_slope_angle + 629ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin - rang_sin; 630ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (object_distance != 0.0) 631ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown ray_height = object_distance * old_axis_slope_angle; 632ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = ray_height / axis_slope_angle; 633ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return; 634ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 635ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = object_distance * (to_index / from_index); 636ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = axis_slope_angle * (from_index / to_index); 637ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return; 638ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 639ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 640ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (radius_of_curvature != 0.0) { 641ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (object_distance == 0.0) { 642ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = 0.0; 643ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin = ray_height / radius_of_curvature; 644ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } else { 645ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin = ((object_distance - 646ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown radius_of_curvature) / radius_of_curvature) * 647ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sin(axis_slope_angle); 648ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 649ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang = asin(iang_sin); 650ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rang_sin = (from_index / to_index) * 651ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang_sin; 652ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown old_axis_slope_angle = axis_slope_angle; 653ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = axis_slope_angle + 654ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown iang - asin(rang_sin); 655ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sagitta = sin((old_axis_slope_angle + iang) / 2.0); 656ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sagitta = 2.0 * radius_of_curvature*sagitta*sagitta; 657ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = ((radius_of_curvature * sin( 658ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown old_axis_slope_angle + iang)) * 659ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cot(axis_slope_angle)) + sagitta; 660ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return; 661ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 662ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 663ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown rang = -asin((from_index / to_index) * 664ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sin(axis_slope_angle)); 665ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = object_distance * ((to_index * 666ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cos(-rang)) / (from_index * 667ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown cos(axis_slope_angle))); 668ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown axis_slope_angle = -rang; 669ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 670ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 671ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Perform ray trace in specific spectral line */ 672ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 673ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownstatic void trace_line(line, ray_h) 674ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownint line; 675ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Browndouble ray_h; 676ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 677ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown int i; 678ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 679ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = 0.0; 680ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown ray_height = ray_h; 681ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown from_index = 1.0; 682ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 683ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (i = 1; i <= current_surfaces; i++) { 684ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown radius_of_curvature = s[i][1]; 685ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown to_index = s[i][2]; 686ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (to_index > 1.0) 687ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown to_index = to_index + ((spectral_line[4] - 688ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[line]) / 689ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) / 690ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown s[i][3]); 691ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown transit_surface(); 692ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown from_index = to_index; 693ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (i < current_surfaces) 694ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown object_distance = object_distance - s[i][4]; 695ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 696ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 697ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 698ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown/* Initialise when called the first time */ 699ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 700ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownint main(argc, argv) 701ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownint argc; 702ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brownchar *argv[]; 703ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown{ 704ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown int i, j, k, errors; 705ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown double od_fline, od_cline; 706ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef ACCURACY 707ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown long passes; 708ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 709ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 710ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[1] = 7621.0; /* A */ 711ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[2] = 6869.955; /* B */ 712ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[3] = 6562.816; /* C */ 713ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[4] = 5895.944; /* D */ 714ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[5] = 5269.557; /* E */ 715ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[6] = 4861.344; /* F */ 716ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[7] = 4340.477; /* G'*/ 717ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown spectral_line[8] = 3968.494; /* H */ 718ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 719ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Process the number of iterations argument, if one is supplied. */ 720ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 721ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (argc > 1) { 722ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown niter = atoi(argv[1]); 723ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (*argv[1] == '-' || niter < 1) { 724ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("This is John Walker's floating point accuracy and\n"); 725ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("performance benchmark program. You call it with\n"); 726ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\nfbench <itercount>\n\n"); 727ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("where <itercount> is the number of iterations\n"); 728ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("to be executed. Archival timings should be made\n"); 729ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("with the iteration count set so that roughly five\n"); 730ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("minutes of execution is timed.\n"); 731ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown exit(0); 732ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 733ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 734ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 735ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Load test case into working array */ 736ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 737ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown clear_aperture = 4.0; 738ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown current_surfaces = 4; 739ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (i = 0; i < current_surfaces; i++) 740ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (j = 0; j < 4; j++) 741ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown s[i + 1][j + 1] = testcase[i][j]; 742ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 743ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef ACCURACY 744ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("Beginning execution of floating point accuracy test...\n"); 745ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown passes = 0; 746ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#else 747ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("Ready to begin John Walker's floating point accuracy\n"); 748ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("and performance benchmark. %d iterations will be made.\n\n", 749ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown niter); 750ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 751ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\nMeasured run time in seconds should be divided by %.f\n", niter / 1000.0); 752ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("to normalise for reporting results. For archival results,\n"); 753ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("adjust iteration count so the benchmark runs about five minutes.\n\n"); 754ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 755ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown //printf("Press return to begin benchmark:"); 756ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown //gets(tbfr); 757ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 758ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 759ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Perform ray trace the specified number of times. */ 760ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 761ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef ACCURACY 762ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown while (TRUE) { 763ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown passes++; 764ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if ((passes % 100L) == 0) { 765ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("Pass %ld.\n", passes); 766ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 767ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#else 768ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (itercount = 0; itercount < niter; itercount++) { 769ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 770ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 771ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (paraxial = 0; paraxial <= 1; paraxial++) { 772ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 773ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Do main trace in D light */ 774ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 775ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown trace_line(4, clear_aperture / 2.0); 776ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown od_sa[paraxial][0] = object_distance; 777ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown od_sa[paraxial][1] = axis_slope_angle; 778ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 779ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown paraxial = FALSE; 780ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 781ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Trace marginal ray in C */ 782ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 783ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown trace_line(3, clear_aperture / 2.0); 784ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown od_cline = object_distance; 785ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 786ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Trace marginal ray in F */ 787ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 788ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown trace_line(6, clear_aperture / 2.0); 789ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown od_fline = object_distance; 790ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 791ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_lspher = od_sa[1][0] - od_sa[0][0]; 792ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) / 793ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown (sin(od_sa[0][1]) * od_sa[0][0]); 794ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_lchrom = od_fline - od_cline; 795ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_lspher = sin(od_sa[0][1]); 796ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 797ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* D light */ 798ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 799ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_lspher = 0.0000926 / (max_lspher * max_lspher); 800ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_osc = 0.0025; 801ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_lchrom = max_lspher; 802ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifndef ACCURACY 803ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 804ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 805ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown //printf("Stop the timer:\007"); 806ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown //gets(tbfr); 807ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 808ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 809ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Now evaluate the accuracy of the results from the last ray trace */ 810ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 811ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[0], "%15s %21.11f %14.11f", 812ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Marginal ray", od_sa[0][0], od_sa[0][1]); 813ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[1], "%15s %21.11f %14.11f", 814ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Paraxial ray", od_sa[1][0], od_sa[1][1]); 815ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[2], 816ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Longitudinal spherical aberration: %16.11f", 817ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_lspher); 818ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[3], 819ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): %16.11f", 820ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_lspher); 821ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[4], 822ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Offense against sine condition (coma): %16.11f", 823ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_osc); 824ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[5], 825ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): %16.11f", 826ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_osc); 827ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[6], 828ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown "Axial chromatic aberration: %16.11f", 829ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown aberr_lchrom); 830ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown sprintf(outarr[7], 831ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown " (Maximum permissible): %16.11f", 832ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown max_lchrom); 833ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 834ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown /* Now compare the edited results with the master values from 835ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown reference executions of this program. */ 836ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown 837ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown errors = 0; 838ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (i = 0; i < 8; i++) { 839ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (strcmp(outarr[i], refarr[i]) != 0) { 840ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef ACCURACY 841ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\nError in pass %ld for results on line %d...\n", 842ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown passes, i + 1); 843ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#else 844ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\nError in results on line %d...\n", i + 1); 845ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 846ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("Expected: \"%s\"\n", refarr[i]); 847ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("Received: \"%s\"\n", outarr[i]); 848ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("(Errors) "); 849ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown k = strlen(refarr[i]); 850ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown for (j = 0; j < k; j++) { 851ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("%c", refarr[i][j] == outarr[i][j] ? ' ' : '^'); 852ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (refarr[i][j] != outarr[i][j]) 853ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown errors++; 854ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 855ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\n"); 856ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 857ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 858ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#ifdef ACCURACY 859ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } 860ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#else 861ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown if (errors > 0) { 862ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\n%d error%s in results. This is VERY SERIOUS.\n", 863ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown errors, errors > 1 ? "s" : ""); 864ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown } else 865ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown printf("\nNo errors in results.\n"); 866ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown#endif 867ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown return 0; 868ed07e00d438c74b7a23c01bfffde77e3968305e4Jeff Brown} 869