#include <perms.h>
// Edinburgh IMP77 Compiler - Version 8.4
//
// 1 %realslong
// 2
// 3 %begin
int main(int argc, char **argv) { /*%begin*/
_imp_initialise(argc, argv);
// 4 %integer a line, b line, c line, d line,
int ALINE;
int BLINE;
int CLINE;
int DLINE;
int ELINE;
int FLINE;
int GPRIMELINE;
int HLINE;
// 5+ e line, f line, g prime line,
// 6+ h line
// 7 %real %array spectral line(1:8)
double SPECTRALLINE[8];
struct _imp_dv SPECTRALLINE_dv = {SPECTRALLINE, sizeof(double) /*, dims, <bounds> */};
// 8 %integer current surfaces
int CURRENTSURFACES;
// 9 %real %array test case(1:4, 1:4)
double TESTCASE[4][4];
struct _imp_dv TESTCASE_dv = {TESTCASE, sizeof(double) /*, dims, <bounds> */};
// 10 %integer curvature radius, index of refraction, dispersion,
int CURVATURERADIUS;
int INDEXOFREFRACTION;
int DISPERSION;
int EDGETHICKNESS;
// 11+ edge thickness
// 12
// 13 %real clear aperture, aberr lspher, aberr osc, aberr lchrom,
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;
// 14+ max lspher, max osc, max lchrom, radius of curvature,
// 15+ object distance, ray height, axis slope angle,
// 16+ from index, to index
// 17 %integer number of iterations, iteration
int NUMBEROFITERATIONS;
int ITERATION;
// 18 %integer paraxial, marginal ray, paraxial ray
int PARAXIAL;
int MARGINALRAY;
int PARAXIALRAY;
// 19 %integer od field, sa field
int ODFIELD;
int SAFIELD;
// 20 %real %array od sa(1:2, 1:2)
double ODSA[2][2];
struct _imp_dv ODSA_dv = {ODSA, sizeof(double) /*, dims, <bounds> */};
// 21 %real od cline, od fline
double ODCLINE;
double ODFLINE;
// 22 %real %fn cot(%real x)
double COT(double X) {
// 23 %result = 1.0 / (sin(x) / cos(x))
return (double)(1.0 / ((double)(_imp_SIN(X) / _imp_COS(X))));
// 24 %end {cot}
} /* end of COT */
// 25 %real %fn arc sin(%real x)
double ARCSIN(double X) {
// 26 %result = 2 * arc tan1(x / (1 + sqrt(1 - x^2)))
return 2 * _imp_ARCTAN1((double)(X / (1 + _imp_SQRT(1 - _imp_REXP(X, 2)))));
// 27 %end {arc sin}
} /* end of ARCSIN */
// 28 %routine out big real(%integer channel, %real x)
void OUTBIGREAL(int CHANNEL, double X) {
// 29 %real fraction
double FRACTION;
// 30 fraction = x
FRACTION = X;
// 31 %end {out big real}
return;
} /* end of OUTBIGREAL */
// ?CHANNEL unused
// 32 %routine transit surface
void TRANSITSURFACE(void) {
// 33 %real iang, rang, iang sin, rang sin,
double IANG;
double RANG;
double IANGSIN;
double RANGSIN;
double OLDAXISSLOPEANGLE;
double SAGITTA;
// 34+ old axis slope angle, sagitta
// 35 %if paraxial = paraxial ray %then %start
if (PARAXIAL != PARAXIALRAY) goto L_0001;
// 36 %if radius of curvature # 0 %then %start
if (RADIUSOFCURVATURE == 0) goto L_0002;
// 37 %if object distance = 0 %then %start
if (OBJECTDISTANCE != 0) goto L_0003;
// 38 axis slope angle = 0
AXISSLOPEANGLE = 0;
// 39 iang sin = ray height / radius of curvature
IANGSIN = (double)(RAYHEIGHT / RADIUSOFCURVATURE);
// 40 %finish %else %c
// 41+ iang sin = ((object distance - %c
goto L_0004;
L_0003:;
IANGSIN = (double)((OBJECTDISTANCE - RADIUSOFCURVATURE) / RADIUSOFCURVATURE) * AXISSLOPEANGLE;
L_0004:;
// 42+ radius of curvature) / radius of curvature) * %c
// 43+ axis slope angle
// 44 rang sin = (from index / to index) * iang sin
RANGSIN = (double)(FROMINDEX / TOINDEX) * IANGSIN;
// 45 old axis slope angle = axis slope angle
OLDAXISSLOPEANGLE = AXISSLOPEANGLE;
// 46 axis slope angle = axis slope angle + %c
AXISSLOPEANGLE = AXISSLOPEANGLE + IANGSIN - RANGSIN;
// 47+ iang sin - rang sin
// 48 %if object distance # 0 %then %c
if (OBJECTDISTANCE == 0) goto L_0005;
RAYHEIGHT = OBJECTDISTANCE * OLDAXISSLOPEANGLE;
L_0005:;
// 49+ ray height = object distance * old axis slope angle
// 50 object distance = ray height / axis slope angle
OBJECTDISTANCE = (double)(RAYHEIGHT / AXISSLOPEANGLE);
// 51 %else
goto L_0006;
L_0002:;
// 52 object distance = object distance * (to index / from index)
OBJECTDISTANCE = OBJECTDISTANCE * ((double)(TOINDEX / FROMINDEX));
// 53 axis slope angle = axis slope angle * (from index / to index)
AXISSLOPEANGLE = AXISSLOPEANGLE * ((double)(FROMINDEX / TOINDEX));
// 54 %finish
L_0006:;
// 55 %else
goto L_0007;
L_0001:;
// 56 %if radius of curvature # 0 %then %start
if (RADIUSOFCURVATURE == 0) goto L_0008;
// 57 %if object distance = 0 %then %start
if (OBJECTDISTANCE != 0) goto L_0009;
// 58 axis slope angle = 0
AXISSLOPEANGLE = 0;
// 59 iang sin = ray height / radius of curvature
IANGSIN = (double)(RAYHEIGHT / RADIUSOFCURVATURE);
// 60 %finish %else %c
// 61+ iang sin = ((object distance - %c
goto L_000A;
L_0009:;
IANGSIN = (double)((OBJECTDISTANCE - RADIUSOFCURVATURE) / RADIUSOFCURVATURE) * _imp_SIN(AXISSLOPEANGLE);
L_000A:;
// 62+ radius of curvature) / radius of curvature) * %c
// 63+ sin(axis slope angle)
// 64 iang = arc sin(iang sin)
IANG = ARCSIN(IANGSIN);
// 65 rang sin = (from index / to index) * iang sin
RANGSIN = (double)(FROMINDEX / TOINDEX) * IANGSIN;
// 66 old axis slope angle = axis slope angle
OLDAXISSLOPEANGLE = AXISSLOPEANGLE;
// 67 axis slope angle = axis slope angle + %c
AXISSLOPEANGLE = AXISSLOPEANGLE + IANG - ARCSIN(RANGSIN);
// 68+ iang - arcsin(rang sin)
// 69
// 70 ! for some reason imp2021 likes this but not the original!:
// 71 sagitta = 2 * radius of curvature * (sin((old axis slope angle + iang) / 2) ^ 2)
SAGITTA = 2 * RADIUSOFCURVATURE * (_imp_REXP(_imp_SIN((double)((OLDAXISSLOPEANGLE + IANG) / 2)), 2));
// 72
// 73 ! sagitta = sin((old axis slope angle + iang) / 2)
// 74 ! sagitta = 2 * radius of curvature * (sagitta ^ 2)
// 75 !%comment sagitta = sin((old axis slope angle + iang) / 2)
// 76 !%comment sagitta = 2 * radius of curvature * (sagitta ^ 2)
// 77
// 78
// 79 object distance = ((radius of curvature * %c
OBJECTDISTANCE = RADIUSOFCURVATURE * _imp_SIN(OLDAXISSLOPEANGLE + IANG) * COT(AXISSLOPEANGLE) + SAGITTA;
// 80+ sin(old axis slope angle + iang)) * %c
// 81+ cot(axis slope angle)) + sagitta
// 82 %else
goto L_000B;
L_0008:;
// 83 rang = -arcsin((from index / to index) * %c
RANG = -ARCSIN((double)(FROMINDEX / TOINDEX) * _imp_SIN(AXISSLOPEANGLE));
// 84+ sin(axis slope angle))
// 85 object distance = object distance * ((to index * %c
OBJECTDISTANCE = OBJECTDISTANCE * ((double)((TOINDEX * _imp_COS(-RANG)) / (FROMINDEX * _imp_COS(AXISSLOPEANGLE))));
// 86+ cos(-rang)) / (from index * %c
// 87+ cos(axis slope angle)))
// 88 axis slope angle = -rang
AXISSLOPEANGLE = -RANG;
// 89 %finish
L_000B:;
// 90 %finish
L_0007:;
// 91 %end {transit surface}
return;
} /* end of TRANSITSURFACE */
// 92 %routine trace line(%integer line, %real ray h)
void TRACELINE(int LINE, double RAYH) {
// 93 %integer i
int I;
// 94 object distance = 0
OBJECTDISTANCE = 0;
// 95 ray height = ray h
RAYHEIGHT = RAYH;
// 96 from index = 1
FROMINDEX = 1;
// 97 %for i = 1, 1, current surfaces %cycle
// %for I = 1, 1, CURRENTSURFACES %cycle
{
int _final = CURRENTSURFACES;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle I = 1, 1, CURRENTSURFACES");
I = 0;
L_000C:
if (I == _final) goto L_000D;
I += 1;
// 98 radius of curvature = test case(i, curvature radius)
RADIUSOFCURVATURE = TESTCASE[I - 1][CURVATURERADIUS - 1];
// 99 to index = test case(i, index of refraction)
TOINDEX = TESTCASE[I - 1][INDEXOFREFRACTION - 1];
// 100 %if to index > 1 %then %c
if (TOINDEX <= 1) goto L_000F;
TOINDEX =
TOINDEX + (double)((SPECTRALLINE[DLINE - 1] - SPECTRALLINE[LINE - 1]) /
(SPECTRALLINE[CLINE - 1] - SPECTRALLINE[FLINE - 1])) *
((double)((TESTCASE[I - 1][INDEXOFREFRACTION - 1] - 1.0) / TESTCASE[I - 1][DISPERSION - 1]));
L_000F:;
// 101+ to index = to index + ((spectral line(d line) - %c
// 102+ spectral line(line)) / %c
// 103+ (spectral line(c line) - spectral line(f line))) * %c
// 104+ ((test case(i, index of refraction) - 1.0) / %c
// 105+ test case(i, dispersion))
// 106 transit surface
TRANSITSURFACE();
// 107 from index = to index
FROMINDEX = TOINDEX;
// 108 %if i < current surfaces %then %c
if (I >= CURRENTSURFACES) goto L_0010;
OBJECTDISTANCE = OBJECTDISTANCE - TESTCASE[I - 1][EDGETHICKNESS - 1];
L_0010:;
// 109+ object distance = object distance - %c
// 110+ test case(i, edge thickness)
// 111 %repeat
goto L_000C;
}
L_000D:;
// 112 %end {trace line}
return;
} /* end of TRACELINE */
// 113 %comment END GLOBAL DECLARATONS
// 114 %comment Spectral lines
// 115 a line = 1; b line = 2; c line = 3; d line = 4
ALINE = 1;
BLINE = 2;
CLINE = 3;
DLINE = 4;
// 116 e line = 5; f line = 6; g prime line = 7
ELINE = 5;
FLINE = 6;
GPRIMELINE = 7;
// 117 h line = 8
HLINE = 8;
// 118 spectral line(a line) = 7621.0
SPECTRALLINE[ALINE - 1] = 7621.0;
// 119 spectral line(b line) = 6869.955
SPECTRALLINE[BLINE - 1] = 6869.955;
// 120 spectral line(c line) = 6562.8160
SPECTRALLINE[CLINE - 1] = 6562.8160;
// 121 spectral line(d line) = 5895.944
SPECTRALLINE[DLINE - 1] = 5895.944;
// 122 spectral line(e line) = 5269.557
SPECTRALLINE[ELINE - 1] = 5269.557;
// 123 spectral line(f line) = 4861.344
SPECTRALLINE[FLINE - 1] = 4861.344;
// 124 spectral line(g prime line) = 4340.477
SPECTRALLINE[GPRIMELINE - 1] = 4340.477;
// 125 spectral line(h line) = 3968.494
SPECTRALLINE[HLINE - 1] = 3968.494;
// 126 %comment Initialise the test case array
// 127 current surfaces = 4
CURRENTSURFACES = 4;
// 128 curvature radius = 1; index of refraction = 2; dispersion = 3
CURVATURERADIUS = 1;
INDEXOFREFRACTION = 2;
DISPERSION = 3;
// 129 edge thickness = 4
EDGETHICKNESS = 4;
// 130 test case(1, curvature radius) = 27.05
TESTCASE[0][CURVATURERADIUS - 1] = 27.05;
// 131 test case(1, index of refraction) = 1.5137
TESTCASE[0][INDEXOFREFRACTION - 1] = 1.5137;
// 132 test case(1, dispersion) = 63.6
TESTCASE[0][DISPERSION - 1] = 63.6;
// 133 test case(1, edge thickness) = 0.52
TESTCASE[0][EDGETHICKNESS - 1] = 0.52;
// 134 test case(2, curvature radius) = -16.68
TESTCASE[1][CURVATURERADIUS - 1] = -16.68;
// 135 test case(2, index of refraction) = 1.0
TESTCASE[1][INDEXOFREFRACTION - 1] = 1.0;
// 136 test case(2, dispersion) = 0.0
TESTCASE[1][DISPERSION - 1] = 0.0;
// 137 test case(2, edge thickness) = 0.138
TESTCASE[1][EDGETHICKNESS - 1] = 0.138;
// 138 test case(3, curvature radius) = -16.68
TESTCASE[2][CURVATURERADIUS - 1] = -16.68;
// 139 test case(3, index of refraction) = 1.6164
TESTCASE[2][INDEXOFREFRACTION - 1] = 1.6164;
// 140 test case(3, dispersion) = 36.7
TESTCASE[2][DISPERSION - 1] = 36.7;
// 141 test case(3, edge thickness) = 0.38
TESTCASE[2][EDGETHICKNESS - 1] = 0.38;
// 142 test case(4, curvature radius) = -78.1
TESTCASE[3][CURVATURERADIUS - 1] = -78.1;
// 143 test case(4, index of refraction) = 1.0
TESTCASE[3][INDEXOFREFRACTION - 1] = 1.0;
// 144 test case(4, dispersion) = 0.0
TESTCASE[3][DISPERSION - 1] = 0.0;
// 145 test case(4, edge thickness) = 0.0
TESTCASE[3][EDGETHICKNESS - 1] = 0.0;
// 146 marginal ray = 1; paraxial ray = 2
MARGINALRAY = 1;
PARAXIALRAY = 2;
// 147 od field = 1; sa field = 2
ODFIELD = 1;
SAFIELD = 2;
// 148 %comment number of iterations = 40263820;
// 149 number of iterations = 1
NUMBEROFITERATIONS = 1;
// 150 clear aperture = 4
CLEARAPERTURE = 4;
// 151 %for iteration = 1, 1, number of iterations %cycle
// %for ITERATION = 1, 1, NUMBEROFITERATIONS %cycle
{
int _final = NUMBEROFITERATIONS;
if (invalid_for_loop(1, 1, _final)) _imp_signal(5, 1, 0, "Illegal cycle ITERATION = 1, 1, NUMBEROFITERATIONS");
ITERATION = 0;
L_0011:
if (ITERATION == _final) goto L_0012;
ITERATION += 1;
// 152 %for paraxial = marginal ray, 1, paraxial ray %cycle
// %for PARAXIAL = MARGINALRAY, 1, PARAXIALRAY %cycle
{
int _initial = MARGINALRAY, _final = PARAXIALRAY;
if (invalid_for_loop(_initial, 1, _final))
_imp_signal(5, 1, 0, "Illegal cycle PARAXIAL = MARGINALRAY, 1, PARAXIALRAY");
PARAXIAL = _initial - (1);
L_0013:
if (PARAXIAL == _final) goto L_0014;
PARAXIAL += 1;
// 153 trace line(d line, clear aperture / 2)
TRACELINE(DLINE, (double)(CLEARAPERTURE / 2));
// 154 od sa(paraxial, od field) = object distance
ODSA[PARAXIAL - 1][ODFIELD - 1] = OBJECTDISTANCE;
// 155 od sa(paraxial, sa field) = axis slope angle
ODSA[PARAXIAL - 1][SAFIELD - 1] = AXISSLOPEANGLE;
// 156 %repeat
goto L_0013;
}
L_0014:;
// 157 %comment Trace marginal ray in C
// 158 paraxial = marginal ray
PARAXIAL = MARGINALRAY;
// 159 trace line(c line, clear aperture / 2)
TRACELINE(CLINE, (double)(CLEARAPERTURE / 2));
// 160 od cline = object distance
ODCLINE = OBJECTDISTANCE;
// 161 %comment Trace marginal ray in F
// 162 trace line(f line, clear aperture / 2)
TRACELINE(FLINE, (double)(CLEARAPERTURE / 2));
// 163 od fline = object distance
ODFLINE = OBJECTDISTANCE;
// 164 aberr lspher = od sa(paraxial ray, od field) - %c
ABERRLSPHER = ODSA[PARAXIALRAY - 1][ODFIELD - 1] - ODSA[MARGINALRAY - 1][ODFIELD - 1];
// 165+ od sa(marginal ray, od field)
// 166 aberr osc = 1 - ((od sa(paraxial ray, od field) * %c
ABERROSC = 1 - (double)((ODSA[PARAXIALRAY - 1][ODFIELD - 1] * ODSA[PARAXIALRAY - 1][SAFIELD - 1]) /
(_imp_SIN(ODSA[MARGINALRAY - 1][SAFIELD - 1]) * ODSA[MARGINALRAY - 1][ODFIELD - 1]));
// 167+ od sa(paraxial ray, sa field)) / %c
// 168+ (sin(od sa(marginal ray, sa field)) * %c
// 169+ od sa(marginal ray, od field)))
// 170 aberr lchrom = od fline - od cline
ABERRLCHROM = ODFLINE - ODCLINE;
// 171 max lspher = 0.0000926 / sin(od sa(marginal ray, sa field)) ^ 2
MAXLSPHER = (double)(0.0000926 / (_imp_REXP(_imp_SIN(ODSA[MARGINALRAY - 1][SAFIELD - 1]), 2)));
// 172 max lchrom = max lspher
MAXLCHROM = MAXLSPHER;
// 173 max osc = 0.0025
MAXOSC = 0.0025;
// 174 %repeat
goto L_0011;
}
L_0012:;
// 175 %comment Print the analysis of the ray trace
// 176 print string(" Marginal ray ")
_imp_PRINTSTRING(_imp_str_literal(" Marginal ray "));
// 177 outbigreal(1, od sa(marginal ray, od field))
OUTBIGREAL(1, ODSA[MARGINALRAY - 1][ODFIELD - 1]);
// 178 print string(" ")
_imp_PRINTSTRING(_imp_str_literal(" "));
// 179 print(od sa(marginal ray, sa field), 5, 11); newline
_imp_PRINT(ODSA[MARGINALRAY - 1][SAFIELD - 1], 5, 11);
_imp_NEWLINE();
// 180 print string(" Paraxial ray ")
_imp_PRINTSTRING(_imp_str_literal(" Paraxial ray "));
// 181 outbigreal(1, od sa(paraxial ray, od field))
OUTBIGREAL(1, ODSA[PARAXIALRAY - 1][ODFIELD - 1]);
// 182 print string(" ")
_imp_PRINTSTRING(_imp_str_literal(" "));
// 183 print(od sa(paraxial ray, sa field), 5, 11); newline
_imp_PRINT(ODSA[PARAXIALRAY - 1][SAFIELD - 1], 5, 11);
_imp_NEWLINE();
// 184 print string("Longitudinal spherical aberration: ")
_imp_PRINTSTRING(_imp_str_literal("Longitudinal spherical aberration: "));
// 185 print(aberr lspher, 5, 11); newline
_imp_PRINT(ABERRLSPHER, 5, 11);
_imp_NEWLINE();
// 186 print string(" (Maximum permissible): ")
_imp_PRINTSTRING(_imp_str_literal(" (Maximum permissible): "));
// 187 print(max lspher, 5, 11); newline
_imp_PRINT(MAXLSPHER, 5, 11);
_imp_NEWLINE();
// 188 print string("Offense against sine condition (coma): ")
_imp_PRINTSTRING(_imp_str_literal("Offense against sine condition (coma): "));
// 189 print(aberr osc, 5, 11); newline
_imp_PRINT(ABERROSC, 5, 11);
_imp_NEWLINE();
// 190 print string(" (Maximum permissible): ")
_imp_PRINTSTRING(_imp_str_literal(" (Maximum permissible): "));
// 191 print(max osc, 5, 11); newline
_imp_PRINT(MAXOSC, 5, 11);
_imp_NEWLINE();
// 192 print string("Axial chromatic aberration: ")
_imp_PRINTSTRING(_imp_str_literal("Axial chromatic aberration: "));
// 193 print(aberr lchrom, 5, 11); newline
_imp_PRINT(ABERRLCHROM, 5, 11);
_imp_NEWLINE();
// 194 print string(" (Maximum permissible): ")
_imp_PRINTSTRING(_imp_str_literal(" (Maximum permissible): "));
// 195 print(max lchrom, 5, 11); newline
_imp_PRINT(MAXLCHROM, 5, 11);
_imp_NEWLINE();
// 196 %end %of %program
_imp_stop(0);
return 0;
end_0:;
}