jc@saber.UUCP (John Cincotta) (07/18/85)
*** REPLACE THIS LINE WITH YOUR MESSAGE ***
here is an astronomical raytrace program it is a changed version of
the program that was in sky and tell a few years back it has been
changed from basic to c it has also been modifyed for doing multiple
passes through the system at angle angle/2 and on axis
any help figering out what some of the variables are and
how this program could be modifyed for a complete raytrace would
be appreaciated
the first file is input to do a 17.5" f5 neutonian
1degree off axis .5degrees and on axis
the second file is the output from the program for above
the third is the c program
8.75 0 1.0
-.005714285 -87.5 -1 1 3
-1. 0 0 0 0
0 0 -1 -1 0
-.005714285 -87.5 -1 1 3
-1. 0 0 0 0
0 0 -1 -1 0
-.005714285 -87.5 -1 1 3
-1. 0 0 0 0
0 0 -1 -1 0
-.005714285 -87.5 -1 1 3
-1. 0 0 0 0
0 0 -1 -1 0
-.005714285 -87.5 -1 1 3
-1. 0 0 0 0
0 0 -1 -1 0
0 1.527318, 1.527318, 0.763601, 0.763601, 0.000000
1 1.527214, 1.527480, 0.763582, 0.763649, 0.000000
2 1.527166, 1.527699, 0.763592, 0.763725, 0.000000
3 1.527176, 1.527976, 0.763630, 0.763830, 0.000000
4 1.527243, 1.528310, 0.763697, 0.763963, 0.000000
5 1.527368, 1.528701, 0.763792, 0.764126, 0.000000
6 1.527549, 1.529150, 0.763917, 0.764317, 0.000000
7 1.527788, 1.529657, 0.764070, 0.764536, 0.000000
8 1.528084, 1.530221, 0.764251, 0.764785, 0.000000
9 1.528438, 1.530843, 0.764461, 0.765062, 0.000000
10 1.528849, 1.531523, 0.764700, 0.765368, 0.000000
11 1.529317, 1.532260, 0.764968, 0.765703, 0.000000
12 1.529842, 1.533056, 0.765264, 0.766067, 0.000000
13 1.530425, 1.533910, 0.765589, 0.766459, 0.000001
14 1.531065, 1.534821, 0.765943, 0.766881, 0.000001
15 1.531763, 1.535791, 0.766326, 0.767332, 0.000001
16 1.532518, 1.536820, 0.766737, 0.767812, 0.000001
17 1.533331, 1.537907, 0.767178, 0.768321, 0.000001
18 1.534201, 1.539052, 0.767647, 0.768859, 0.000001
19 1.535129, 1.540256, 0.768146, 0.769426, 0.000001
20 1.536115, 1.541519, 0.768673, 0.770023, 0.000001
1.531242 1.531242 0.765562 0.765562 0.000000
0 0.003924, 0.003924, 0.001961, 0.001961, 0.000000
1 0.004029, 0.003762, 0.001980, 0.001913, 0.000000
2 0.004076, 0.003543, 0.001970, 0.001837, 0.000000
3 0.004066, 0.003266, 0.001932, 0.001732, 0.000000
4 0.003999, 0.002932, 0.001865, 0.001599, 0.000000
5 0.003874, 0.002541, 0.001770, 0.001436, 0.000000
6 0.003693, 0.002092, 0.001645, 0.001246, 0.000000
7 0.003454, 0.001585, 0.001493, 0.001026, 0.000000
8 0.003158, 0.001021, 0.001311, 0.000777, 0.000000
9 0.002804, 0.000399, 0.001101, 0.000500, 0.000000
10 0.002394, -0.000281, 0.000862, 0.000194, 0.000000
11 0.001926, -0.001018, 0.000594, -0.000141, -0.000000
12 0.001400, -0.001814, 0.000298, -0.000505, -0.000000
13 0.000817, -0.002667, -0.000027, -0.000897, -0.000000
14 0.000177, -0.003579, -0.000381, -0.001319, -0.000000
15 -0.000520, -0.004549, -0.000764, -0.001770, -0.000000
16 -0.001276, -0.005578, -0.001175, -0.002250, -0.000000
17 -0.002088, -0.006664, -0.001616, -0.002759, -0.000000
18 -0.002959, -0.007810, -0.002085, -0.003297, -0.000000
19 -0.003887, -0.009014, -0.002584, -0.003864, -0.000000
20 -0.004873, -0.010277, -0.003111, -0.004461, -0.000000
#include <math.h>
double l[20][8]; /*the surface parameters */
/* first parameter of the aray
0 is object at infinity
1 is entry point
2 is first surface
3 ...next surfaces
second parameter of the aray
0 curvature
1 thickness
2 index of refraction
3 4 5 6 7 aspherisity parameters
*/
double average[6]; /*this is where the average for an offset is asembeled*/
double printer[6]; /*this is used to get the numbers to print the last aray */
double y[22]; /*ray height at each surface */
double output[21][5]; /*aray of ray heights for last surface */
double yy; /*current height*/
double d; /* */
double l1; /*cosines???? */
double l2; /* */
double m1; /*cosines???? */
double m2; /* */
double m3; /* */
double g1; /* */
double g2; /* */
double g3; /* */
double g4; /* */
double yy1; /* */
double z; /* */
double r; /* */
double r1; /* */
double d; /* */
double s; /* */
double aa; /*input angle */
double ab; /*entry point ray height */
double ad; /*entry point to first surface distance */
double an; /*tan of ray angle */
double d1; /* */
int k[20]; /*surface code*/
int k5; /*used to tell no more points ?????*/
int l5; /*last surface number*/
int q; /*count different ray heights*/
int j; /*current surface number*/
int verbose = 0 ; /*do you want any printout except last */
main(argc, argv)
int argc;
char **argv;
{
double t1,t2,t3,t4; /*temp work variables */
int i; /*general counter */
int s1; /*surface number used on input */
int pass; /*pass counter 1 to 5 */
if(argc>1)verbose = 1; /*if there were args they want printout */
if(verbose)printf ("raytrace pgm\n");
if(verbose)printf("enter etry point ray ht,entry pt -s3, input angle and k5\n");
scanf("%f %f %f %d",&ab,&ad,&aa,&k5);
an = tan(aa/57.29577951);
/*make 5 passes +- input angle,+- .5 input angle,on axsis */
for (pass =1;pass < 6;pass++){ /* different obliquetys*/
if(verbose)printf ("pass %d\n",pass);
for (i=0;i<20;i++){ /*zero out arays before entering data*/
for(j=0;j<8;j++){
l[i][j] = 0.0;
}
k[i] = 0;
}
if(pass < 3){ /*angel requested */
an = tan(aa/57.29577951);
}else if(pass < 5){ /*angel/2 requested */
an = tan(aa/(2*57.29577951));
}else{ /*on axis */
an = 0;
}
l[1][1] = ad;
l[0][2] = 1.0;
l[1][2] = 1.0;
if (k5 == 0){
for (i=2;i<20;i++){ /*get l[i][0-2],k,s1 */
if(verbose)printf("enter cv,th,index,skode,sufno\n");
scanf("%f %f %f %d %d",&l[i][0],&l[i][1],&l[i][2],&k[i],&s1);
if (k[i] > 0){ /* get l[1][3-7] */
if(verbose)printf("enter conic aspherisity 4,6,8,10\n");
scanf("%f %f %f %f %f",&l[i][3],&l[i][4],&l[i][5],&l[i][6],&l[i][7]);
if(l[i][3] > 0.0){
l[i][3] = (l[i][3] - 1.0);
}
}
if(s1 == 0){
l5 = i;
break;
}
}
for (i=0;i<=l5;i++){
l[i][3] = (l[i][3]+1.0);
}
}
/*this ends the inputing of data and begins the traceing of rays for
a given obliquity */
for (q=0;q<21;q++){ /*differnt ray heights*/
y[1] = q*ab/20.0; /*ray height*/
l[1][1] = ad; /*entry point to first surface*/
k[l5] = -1; /* makes the focal plane flat*/
k[1] = -1; /* makes the entry point flat*/
l[0][1] = 1.0e20;
y[0] = 0.0;
yy = y[1];
z = 0.0;
m2 = 1.0;
t1 = an; /* tan of input angle */
m1 = 1/sqrt(1+t1*t1);
l1 = t1*m1;
g1 = m1;
for(j = 2;j <=l5;j++){ /*different surfaces*/
d = (l[(j-1)][1]-z)/m1;
yy1 = yy + l1*d;
if (k[j] < 0){ /* this handels flat surfaces skode = -1 */
yy = yy1;
z = 0.0 ;
l2 = z;
m2 =1.0;
g3 = m1;
}else { /* this handels simple curved surfaces skode = 0 */
/* and does same before calling higher for aspheric surfaces skode = 1*/
t1 = l[j][0]*(1.0+(l[j][3]-1)*m1*m1);
t2 = m1 -l[j][0]*yy1*l1;
t3 = l[j][0]*yy1*yy1;
t4 = t2*t2-t1*t3;
if(t4<=0.0){ /*bail out if sqrt of -number*/
printf("failure on curved surface %d t4 = %f\n",j,t4);
exit(0);
}
g3 = sqrt(t4);
if(t2<0.0){
g3 = -g3;
}
d1 = t3/(t2+g3);
yy = yy1+l1*d1;
z = m1*d1;
if (k[j] > 0){
higher(); /*go here for surfaces beyond sphere*/
}
if (k[j] == 0){
l2 = -l[j][0]*yy;
m2 = 1-l[j][3]*l[j][0]*z;
}
}
t1= yy*yy; /*transfer ray from surface ??*/
y[j] = yy;
if (j <= l5){
g2 = (l2*l2)+(m2*m2);
r1= l[j-1][2]/l[j][2];
if(r1>=0.0){
r = r1;
t1 = r*g3/g2;
t2 = (r*r-1)/g2;
t3 = t1*t1-t2;
if (t3 <=0.0){ /*bail out if sqrt of -number*/
printf("trace failed surface no%d,%f\n",j,t3);
exit(0);
}
g1 = sqrt(t3);
if (t1<0.0){
g1 = -g1;
}
g4 =g1 -t1;
}else{
g1 = -g3 / g2;
r = 1;
g4 = g1 +g1;
}
l1 = r * l1 + (g4 * l2); /*recalculate cosines ???*/
m1 = r * m1 +(g4 * m2); /*recalculate cosines ???*/
}
} /*end of different surfaces*/
/* this is where we print out the report */
if(verbose){
printf("curvature,thickness,index,y\n");
for (i = 2;i <= l5;i++){
printf("%d,%10.6f,%10.6f,%10.6f,%10.6f\n",i,l[i][0],l[i][1],l[i][2],y[i]);
}
}
output[q][(pass-1)] = y[l5];
} /*differnt ray heights*/
ab = -ab; /* other side of the lense */
} /* different obliquetys*/
for(i=0;i<6;i++){
average[i] = 0.;
printer[i] = 0.;
}
for(i=0;i<21;i++){
average[0] += output[i][0];
average[1] += output[i][1];
average[2] += output[i][2];
average[3] += output[i][3];
average[4] += output[i][4];
printf("%d %10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n",i,output[i][0],output[i][1],output[i][2],output[i][3],output[i][4]);
}
average[0] = average[0]/21.;
average[1] = average[1]/21.;
average[2] = average[2]/21.;
average[3] = average[3]/21.;
average[4] = average[4]/21.;
printer[0] = (average[0]+average[1])/2;
average[0] = average[1] = printer[0];
printer[0] = (average[2]+average[3])/2;
average[2] = average[3] = printer[0];
printf("\n%10.6f %10.6f %10.6f %10.6f %10.6f\n\n",average[0],average[1],average[2],average[3],average[4]);
for(i=0;i<21;i++){
printer[0] =( average[0] - output[i][0]);
printer[1] =( average[1] - output[i][1]);
printer[2] =( average[2] - output[i][2]);
printer[3] =( average[3] - output[i][3]);
printer[4] =( average[4] - output[i][4]);
printf("%d %10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n",i,printer[0],printer[1],printer[2],printer[3],printer[4]);
}
}
higher()
{ /* this handels higher order curved surfaces skode = 1 */
/*it does it by sucessive approximation */
int i; /*iteration counter */
double t1,t2,t3,t4,t5,t6; /*temp work variables */
double f; /*partial accumulator for aspherisity calculations ??*/
m2 = 1.0;
for(i= 0 ;i<1000;i++){
t1 = yy*yy;
t2 = 1-l[j][3]*l[j][0]*l[j][0]*t1;
if(t2 <=0.0){ /*bail out if sqrt of -number*/
printf("aspheric test failed surface%d, t2 =%f\n",j,t2);
exit();
}
t3 = sqrt(t2);
t4 = l[j][0]/(1.0+t3);
t5 = l[j][0]/t3;
s = z-((((l[j][7]*t1+l[j][6])*t1+l[j][5])*t1+l[j][4])*t1+t4)*t1;
f = (((10*l[j][7]*t1+8*l[j][6])*t1+6*l[j][5])*t1+4*l[j][4])*t1;
f = f+t5;
l2 = -yy*f;
g3 = l2*l1+m1;
d1 = d1-s/g3;
yy = yy1+l1*d1;
z = m1*d1;
t6 = abs(s);
if(t6 <= .0000001){
return(0);
}
}
printf("aspheric test failed in 1000 passes %f\n",t6);
exit();
}