#include int _imp_mainep(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:8*/]; // %real %array spectral line(1:8) int Currentsurfaces; // %integer current surfaces double Testcase[4 /*1:4*/][4 /*1:4*/]; // %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:2*/][2 /*1:2*/]; // %real %array od sa(1:2, 1:2) double Odcline; double Odfline; // %real od cline, od fline double Cot(double X) { // %real %fn cot(%real x) return (1.0 / (Sin(X) / Cos(X))); // %result = 1.0 / (sin(x) / cos(x)) } // %end {cot} double Arcsin(double X) { // %real %fn arc sin(%real x) return (2 * Arctan1(X / (1 + Sqrt(1 - REXP(X, 2))))); // %result = 2 * arc tan1(x / (1 + sqrt(1 - x^2))) } // %end {arc sin} 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} 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) { // %if radius of curvature # 0 // %then %start if (!Objectdistance) { // %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; // %finish %else // iang sin = ((object distance - radius of curvature) / // radius of curvature) * 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; // axis slope angle = axis slope // angle + iang sin - rang sin if (Objectdistance) Rayheight = Objectdistance * Oldaxisslopeangle; // %if object distance # 0 %then 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) { // %if radius of curvature # 0 // %then %start if (!Objectdistance) { // %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(Axisslopeangle); // %finish %else iang sin = ((object // distance - radius of curvature) / // radius of curvature) * sin(axis // slope angle) Iang = Arcsin(Iangsin); // 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(Rangsin); // axis slope angle = axis slope angle + // iang - arcsin(rang sin) Sagitta = 2 * Radiusofcurvature * (REXP(Sin((Oldaxisslopeangle + Iang) / 2), 2)); // sagitta = 2 * radius of curvature * (sin((old // axis slope angle + iang) / 2) ^ 2) Objectdistance = ((Radiusofcurvature * Sin(Oldaxisslopeangle + Iang)) * Cot(Axisslopeangle)) + Sagitta; // object distance = ((radius of curvature * sin(old // axis slope angle + iang)) * cot(axis slope angle)) // + sagitta } else { // %else Rang = -Arcsin( (Fromindex / Toindex) * Sin(Axisslopeangle)); // rang = -arcsin((from index / to // index) * sin(axis slope angle)) Objectdistance = Objectdistance * ((Toindex * Cos(-Rang)) / (Fromindex * Cos(Axisslopeangle))); // object distance = object // distance * ((to index * // cos(-rang)) / (from index * // cos(axis slope angle))) Axisslopeangle = -Rang; // axis slope angle = -rang } // %finish } // %finish } // %end {transit surface} 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++) { // %for i = 1, 1, current surfaces %cycle Radiusofcurvature = Testcase[I] [Curvatureradius]; // radius of curvature = test // case(i, curvature radius) Toindex = Testcase[I] [Indexofrefraction]; // to index = test case(i, // index of refraction) if (Toindex > 1) Toindex += ((Spectralline[Dline] - Spectralline[Line]) / (Spectralline[Cline] - Spectralline[Fline])) * ((Testcase[I][Indexofrefraction] - 1.0) / Testcase[I][Dispersion]); // %if to index > 1 %then to index = // to index + ((spectral line(d // line) - spectral line(line)) / // (spectral line(c line) - spectral // line(f line))) * ((test case(i, // index of refraction) - 1.0) / // test case(i, dispersion)) Transitsurface(); // transit surface Fromindex = Toindex; // from index = to index if (I < Currentsurfaces) Objectdistance -= Testcase[I] [Edgethickness]; // %if i < current surfaces %then // object distance = object distance - // 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++) { // %for iteration = 1, 1, number of iterations // %cycle for (Paraxial = Marginalray; Paraxial <= Paraxialray; Paraxial++) { // %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[Paraxialray][Odfield] - Odsa[Marginalray] [Odfield]; // aberr lspher = od sa(paraxial ray, od field) - 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[Paraxialray][Odfield] * Odsa[Paraxialray][Safield]) / (Sin(Odsa[Marginalray][Safield]) * Odsa[Marginalray] [Odfield])); // aberr osc = 1 - ((od sa(paraxial ray, od // field) * od sa(paraxial ray, sa field)) / // (sin(od sa(marginal ray, sa field)) * od // sa(marginal ray, od field))) Aberrlchrom = Odfline - Odcline; // aberr lchrom = od fline - od cline Maxlspher = 0.0000926 / REXP(Sin(Odsa[Marginalray][Safield]), 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(_imp_str_literal( " Marginal ray ")); // print string(" Marginal ray ") Outbigreal( 1, Odsa[Marginalray] [Odfield]); // outbigreal(1, od sa(marginal ray, od field)) Printstring(_imp_str_literal(" ")); // print string(" ") Print(Odsa[Marginalray][Safield], 5, 11); Newline(); // print(od sa(marginal ray, sa field), 5, 11); newline Printstring(_imp_str_literal( " Paraxial ray ")); // print string(" Paraxial ray ") Outbigreal( 1, Odsa[Paraxialray] [Odfield]); // outbigreal(1, od sa(paraxial ray, od field)) Printstring(_imp_str_literal(" ")); // print string(" ") Print(Odsa[Paraxialray][Safield], 5, 11); Newline(); // print(od sa(paraxial ray, sa field), 5, 11); newline Printstring(_imp_str_literal( "Longitudinal spherical aberration: ")); // print // string("Longitudinal // spherical aberration: ") Print(Aberrlspher, 5, 11); Newline(); // print(aberr lspher, 5, 11); newline Printstring(_imp_str_literal( " (Maximum permissible): ")); // print string(" (Maximum // permissible): ") Print(Maxlspher, 5, 11); Newline(); // print(max lspher, 5, 11); newline Printstring(_imp_str_literal( "Offense against sine condition (coma): ")); // print // string("Offense // against sine // condition (coma): ") Print(Aberrosc, 5, 11); Newline(); // print(aberr osc, 5, 11); newline Printstring(_imp_str_literal( " (Maximum permissible): ")); // print string(" (Maximum // permissible): ") Print(Maxosc, 5, 11); Newline(); // print(max osc, 5, 11); newline Printstring(_imp_str_literal( "Axial chromatic aberration: ")); // print string("Axial chromatic // aberration: ") Print(Aberrlchrom, 5, 11); Newline(); // print(aberr lchrom, 5, 11); newline Printstring(_imp_str_literal( " (Maximum permissible): ")); // print string(" (Maximum // permissible): ") Print(Maxlchrom, 5, 11); Newline(); // print(max lchrom, 5, 11); newline exit(0); return (1); } // %end %of %program