#include "impstrings.h" #define unless(cond) if (!(cond))// reals long // %realslong double SQRT(double NUM) { // %externalrealfn sqrt(%real num) return 0.0; // %result = 0.0 } // %end {sqrt} double ARCTAN(double ANGLE) { // %externalrealfn arctan(%real angle) return 0.0; // %result = 0.0 } // %end {arc tan} double SIN(double ANGLE) { // %externalrealfn sin(%real angle) return 0.0; // %result = 0.0 } // %end {sin} double COS(double ANGLE) { // %externalrealfn cos(%real angle) return 0.0; // %result = 0.0 } // %end {cos} void PRINTFL(double R, int P, int D) { // %externalroutine printfl(%real r, %integer p, d) } // %end int main(int _imp_argc, char **_imp_argv) // %begin { int ALINE; int BLINE; int CLINE; int DLINE; int ELINE; int FLINE; int GPRIMELINE; int HLINE; // h line double SPECTRALLINE[(8) - (1) + 1]; // %real %array spectral line(1:8) int CURRENTSURFACES; // %integer current surfaces double TESTCASE[(4) - (1) + 1][(4) - (1) + 1]; // %real %array test case(1:4, 1:4) int CURVATURERADIUS; int INDEXOFREFRACTION; int DISPERSION; int EDGETHICKNESS; double CLEARAPERTURE; double ABERRLSPHER; double ABERROSC; double ABERRLCHROM; double MAXLSPHER; double MAXOSC; double MAXLCHROM; double RADIUSOFCURVATURE; double OBJECTDISTANCE; double RAYHEIGHT; double AXISSLOPEANGLE; double FROMINDEX; double TOINDEX; int NUMBEROFITERATIONS; int ITERATION; // %integer number of iterations, iteration int PARAXIAL; int MARGINALRAY; int PARAXIALRAY; // %integer paraxial, marginal ray, paraxial ray int ODFIELD; int SAFIELD; // %integer od field, sa field double ODSA[(2) - (1) + 1][(2) - (1) + 1]; // %real %array od sa(1:2, 1:2) double ODCLINE; double ODFLINE; // %real od cline, od fline auto double COT(double X) { // %real %fn cot(%real x) return 1.0 / (SIN / COS); // %result = 1.0 / (sin(x) / cos(x)) } // %end {cot} auto double ARCSIN(double X) { // %real %fn arc sin(%real x) return 2 * ARCTAN; // %result = 2 * arc tan(x / (1 + sqrt(1 - x^2))) } // %end {arc sin} auto void OUTBIGREAL(int CHANNEL, double X) { // %routine out big real(%integer channel, %real x) double FRACTION; // %real fraction FRACTION = X; // fraction = x } // %end {out big real} auto void TRANSITSURFACE(void) { // %routine transit surface double IANG; double RANG; double IANGSIN; double RANGSIN; double OLDAXISSLOPEANGLE; double SAGITTA; if (PARAXIAL == PARAXIALRAY) { // %if paraxial = paraxial ray %then %start if (RADIUSOFCURVATURE != 0) { // %if radius of curvature # 0 %then %start if (OBJECTDISTANCE == 0) { // %if object distance = 0 %then %start AXISSLOPEANGLE = 0; // axis slope angle = 0 IANGSIN = RAYHEIGHT / RADIUSOFCURVATURE; // iang sin = ray height / radius of curvature } else IANGSIN = ((OBJECTDISTANCE - RADIUSOFCURVATURE) / RADIUSOFCURVATURE) * AXISSLOPEANGLE; // axis slope angle RANGSIN = (FROMINDEX / TOINDEX) * IANGSIN; // rang sin = (from index / to index) * iang sin OLDAXISSLOPEANGLE = AXISSLOPEANGLE; // old axis slope angle = axis slope angle AXISSLOPEANGLE = AXISSLOPEANGLE + IANGSIN - RANGSIN; // iang sin - rang sin if (OBJECTDISTANCE != 0) RAYHEIGHT = OBJECTDISTANCE * OLDAXISSLOPEANGLE; ; // ray height = object distance * old axis slope angle OBJECTDISTANCE = RAYHEIGHT / AXISSLOPEANGLE; // object distance = ray height / axis slope angle } else { // %else OBJECTDISTANCE = OBJECTDISTANCE * (TOINDEX / FROMINDEX); // object distance = object distance * (to index / from index) AXISSLOPEANGLE = AXISSLOPEANGLE * (FROMINDEX / TOINDEX); // axis slope angle = axis slope angle * (from index / to index) } // %finish } else { // %else if (RADIUSOFCURVATURE != 0) { // %if radius of curvature # 0 %then %start if (OBJECTDISTANCE == 0) { // %if object distance = 0 %then %start AXISSLOPEANGLE = 0; // axis slope angle = 0 IANGSIN = RAYHEIGHT / RADIUSOFCURVATURE; // iang sin = ray height / radius of curvature } else IANGSIN = ((OBJECTDISTANCE - RADIUSOFCURVATURE) / RADIUSOFCURVATURE) * SIN; // sin(axis slope angle) IANG = ARCSIN; // iang = arc sin(iang sin) RANGSIN = (FROMINDEX / TOINDEX) * IANGSIN; // rang sin = (from index / to index) * iang sin OLDAXISSLOPEANGLE = AXISSLOPEANGLE; // old axis slope angle = axis slope angle AXISSLOPEANGLE = AXISSLOPEANGLE + IANG - ARCSIN; // iang - arcsin(rang sin) SAGITTA = 2 * RADIUSOFCURVATURE * (SIN rexp() 2); // sagitta = 2 * radius of curvature * (sin((old axis slope angle + iang) / 2) ^ 2) OBJECTDISTANCE = ((RADIUSOFCURVATURE * SIN) * COT) + SAGITTA; // cot(axis slope angle)) + sagitta } else { // %else RANG = -ARCSIN; // sin(axis slope angle)) OBJECTDISTANCE = OBJECTDISTANCE * ((TOINDEX * COS) / (FROMINDEX * COS)); // cos(axis slope angle))) AXISSLOPEANGLE = -RANG; // axis slope angle = -rang } // %finish } // %finish } // %end {transit surface} auto void TRACELINE(int LINE, double RAYH) { // %routine trace line(%integer line, %real ray h) int I; // %integer i OBJECTDISTANCE = 0; // object distance = 0 RAYHEIGHT = RAYH; // ray height = ray h FROMINDEX = 1; // from index = 1 for (I = 1; I <= CURRENTSURFACES; I += 1) { // %for i = 1, 1, current surfaces %cycle RADIUSOFCURVATURE = TESTCASE; // radius of curvature = test case(i, curvature radius) TOINDEX = TESTCASE; // to index = test case(i, index of refraction) if (TOINDEX > 1) TOINDEX = TOINDEX + ((SPECTRALLINE - SPECTRALLINE) / % c(SPECTRALLINE - SPECTRALLINE)) * % c((TESTCASE - 1.0) / TESTCASE); ; // test case(i, dispersion)) TRANSITSURFACE; // transit surface FROMINDEX = TOINDEX; // from index = to index if (I < CURRENTSURFACES) OBJECTDISTANCE = OBJECTDISTANCE - TESTCASE; ; // test case(i, edge thickness) } // %repeat } // %end {trace line} // END GLOBAL DECLARATONS // %comment END GLOBAL DECLARATONS // Spectral lines // %comment Spectral lines ALINE = 1; BLINE = 2; CLINE = 3; DLINE = 4; // a line = 1; b line = 2; c line = 3; d line = 4 ELINE = 5; FLINE = 6; GPRIMELINE = 7; // e line = 5; f line = 6; g prime line = 7 HLINE = 8; // h line = 8 SPECTRALLINE(ALINE) = 7621.0; // spectral line(a line) = 7621.0 SPECTRALLINE(BLINE) = 6869.955; // spectral line(b line) = 6869.955 SPECTRALLINE(CLINE) = 6562.8160; // spectral line(c line) = 6562.8160 SPECTRALLINE(DLINE) = 5895.944; // spectral line(d line) = 5895.944 SPECTRALLINE(ELINE) = 5269.557; // spectral line(e line) = 5269.557 SPECTRALLINE(FLINE) = 4861.344; // spectral line(f line) = 4861.344 SPECTRALLINE(GPRIMELINE) = 4340.477; // spectral line(g prime line) = 4340.477 SPECTRALLINE(HLINE) = 3968.494; // spectral line(h line) = 3968.494 // Initialise the test case array // %comment Initialise the test case array CURRENTSURFACES = 4; // current surfaces = 4 CURVATURERADIUS = 1; INDEXOFREFRACTION = 2; DISPERSION = 3; // curvature radius = 1; index of refraction = 2; dispersion = 3 EDGETHICKNESS = 4; // edge thickness = 4 TESTCASE(1, CURVATURERADIUS) = 27.05; // test case(1, curvature radius) = 27.05 TESTCASE(1, INDEXOFREFRACTION) = 1.5137; // test case(1, index of refraction) = 1.5137 TESTCASE(1, DISPERSION) = 63.6; // test case(1, dispersion) = 63.6 TESTCASE(1, EDGETHICKNESS) = 0.52; // test case(1, edge thickness) = 0.52 TESTCASE(2, CURVATURERADIUS) = -16.68; // test case(2, curvature radius) = -16.68 TESTCASE(2, INDEXOFREFRACTION) = 1.0; // test case(2, index of refraction) = 1.0 TESTCASE(2, DISPERSION) = 0.0; // test case(2, dispersion) = 0.0 TESTCASE(2, EDGETHICKNESS) = 0.138; // test case(2, edge thickness) = 0.138 TESTCASE(3, CURVATURERADIUS) = -16.68; // test case(3, curvature radius) = -16.68 TESTCASE(3, INDEXOFREFRACTION) = 1.6164; // test case(3, index of refraction) = 1.6164 TESTCASE(3, DISPERSION) = 36.7; // test case(3, dispersion) = 36.7 TESTCASE(3, EDGETHICKNESS) = 0.38; // test case(3, edge thickness) = 0.38 TESTCASE(4, CURVATURERADIUS) = -78.1; // test case(4, curvature radius) = -78.1 TESTCASE(4, INDEXOFREFRACTION) = 1.0; // test case(4, index of refraction) = 1.0 TESTCASE(4, DISPERSION) = 0.0; // test case(4, dispersion) = 0.0 TESTCASE(4, EDGETHICKNESS) = 0.0; // test case(4, edge thickness) = 0.0 MARGINALRAY = 1; PARAXIALRAY = 2; // marginal ray = 1; paraxial ray = 2 ODFIELD = 1; SAFIELD = 2; // od field = 1; sa field = 2 // number of iterations = 40263820; // %comment number of iterations = 40263820; NUMBEROFITERATIONS = 1; // number of iterations = 1 CLEARAPERTURE = 4; // clear aperture = 4 for (ITERATION = 1; ITERATION <= NUMBEROFITERATIONS; ITERATION += 1) { // %for iteration = 1, 1, number of iterations %cycle for (PARAXIAL = MARGINALRAY; PARAXIAL <= PARAXIALRAY; PARAXIAL += 1) { // %for paraxial = marginal ray, 1, paraxial ray %cycle TRACELINE(DLINE, CLEARAPERTURE / 2); // trace line(d line, clear aperture / 2) ODSA(PARAXIAL, ODFIELD) = OBJECTDISTANCE; // od sa(paraxial, od field) = object distance ODSA(PARAXIAL, SAFIELD) = AXISSLOPEANGLE; // od sa(paraxial, sa field) = axis slope angle } // %repeat // Trace marginal ray in C // %comment Trace marginal ray in C PARAXIAL = MARGINALRAY; // paraxial = marginal ray TRACELINE(CLINE, CLEARAPERTURE / 2); // trace line(c line, clear aperture / 2) ODCLINE = OBJECTDISTANCE; // od cline = object distance // Trace marginal ray in F // %comment Trace marginal ray in F TRACELINE(FLINE, CLEARAPERTURE / 2); // trace line(f line, clear aperture / 2) ODFLINE = OBJECTDISTANCE; // od fline = object distance ABERRLSPHER = ODSA - ODSA; // od sa(marginal ray, od field) // ! // ! WARNING: the %c is making it through to 'c', presumably as white space? Find and fix. // ! ABERROSC = 1 - ((ODSA * ODSA) / % c(SIN * ODSA)); // od sa(marginal ray, od field))) ABERRLCHROM = ODFLINE - ODCLINE; // aberr lchrom = od fline - od cline MAXLSPHER = 0.0000926 / SIN rexp() 2; // max lspher = 0.0000926 / sin(od sa(marginal ray, sa field)) ^ 2 MAXLCHROM = MAXLSPHER; // max lchrom = max lspher MAXOSC = 0.0025; // max osc = 0.0025 } // %repeat // Print the analysis of the ray trace // %comment Print the analysis of the ray trace PRINTSTRING(" Marginal ray "); // print string(" Marginal ray ") OUTBIGREAL(1, ODSA); // outbigreal(1, od sa(marginal ray, od field)) PRINTSTRING(" "); // print string(" ") PRINTFL(ODSA, 1, 1); NEWLINE; // printfl(od sa(marginal ray, sa field), 1, 1); newline PRINTSTRING(" Paraxial ray "); // print string(" Paraxial ray ") OUTBIGREAL(1, ODSA); // outbigreal(1, od sa(paraxial ray, od field)) PRINTSTRING(" "); // print string(" ") PRINTFL(ODSA, 1, 1); NEWLINE; // printfl(od sa(paraxial ray, sa field), 1, 1); newline PRINTSTRING("Longitudinal spherical aberration: "); // print string("Longitudinal spherical aberration: ") PRINTFL(ABERRLSPHER, 1, 1); NEWLINE; // printfl(aberr lspher, 1, 1); newline PRINTSTRING(" (Maximum permissible): "); // print string(" (Maximum permissible): ") PRINTFL(MAXLSPHER, 1, 1); NEWLINE; // printfl(max lspher, 1, 1); newline PRINTSTRING("Offense against sine condition (coma): "); // print string("Offense against sine condition (coma): ") PRINTFL(ABERROSC, 1, 1); NEWLINE; // printfl(aberr osc, 1, 1); newline PRINTSTRING(" (Maximum permissible): "); // print string(" (Maximum permissible): ") PRINTFL(MAXOSC, 1, 1); NEWLINE; // printfl(max osc, 1, 1); newline PRINTSTRING("Axial chromatic aberration: "); // print string("Axial chromatic aberration: ") PRINTFL(ABERRLCHROM, 1, 1); NEWLINE; // printfl(aberr lchrom, 1, 1); newline PRINTSTRING(" (Maximum permissible): "); // print string(" (Maximum permissible): ") PRINTFL(MAXLCHROM, 1, 1); NEWLINE; // printfl(max lchrom, 1, 1); newline } // %end %of %program