/* SCHCORR is Schmidt corrector simulation software by Dominic-Luc Webb. 

Last updated 15 Nov 2009

Written for Unix/linux/BSD. Note the lines below (Start/End editing) where you 
edit the constants according to your optical configuration. They are all in 
the one place. Variable names are same as those in Rutten & vanVenrooij "Telescope 
Optics: Evaluation and Design", section 21.11. Compile with gcc using:

bash# gcc -o schcorr schcorr.c -lm

Data output is schcorr.dat. Schcorr also outputs a gnuplot file called schcorr.gnu. 
Run "gnuplot schcorr.gnu" to get schcorr.ps and schcorr.pbm graphics output. The ps 
file prints directly to a postscipt printer or can be viewed by ghostscript, etc 
and the pbm file can be openned by UNIX graphics software like XV if you want to 
view or convert format.

Schcorr prints values of all constants including corrector profile (terms A,B and C) 
to screem by default. That same output can instead be saved to file by redirect:

bash# schcorr > output_filename

*/

#include <stdio.h>
#include <math.h>

FILE *fp, *fp2;

main()
{


float i;
float T, Dc, h_n, r, n, h_p, p_h, p_h_temp, f, a, b, c, x, z, A, B, C, Oslo_RADIUS;

if((fp=fopen("schcorr.dat", "wt"))==NULL) {
printf("can't open file\n");
return 1;
}

if((fp2=fopen("schcorr.gnu", "wt"))==NULL) {
printf("can't open file\n");
return 1;
}


/* Start editing - Comstants below need to be edited for specific design, T is thickness */
T = 6.000000;
Dc = 250;
h_n = 0.866*Dc/2;
r = -961;  
n = 1.5211;
/* End editing - no more editing after this */ 


h_p = h_n / r;

p_h_temp = 1 - pow(h_p,2);
p_h = 2 * h_p * pow(p_h_temp,0.5);

f = h_n / p_h;

a = r * (1 - 0.5 / pow(p_h_temp,0.5));
b = r - a;
c = r * pow(p_h_temp,0.5);
x = (r+a-b-c)/(1-n);

printf("\nh_n = %f, ph = %f, f = %f\na = %f, b = %f\nc = %f, x = %f\n", h_n,p_h,f,a,b,c,x);

A = (f - r/2) / (f *r*(n-1));
B = (3*x-2*pow(h_n,2)*A) / pow(h_n,4);
C = (pow(h_n,2)*A-2*x) / pow(h_n,6);
Oslo_RADIUS = 0.5/A;

printf("A = %1.6e, B = %1.6e, C = %1.6e\n", A, B, C);
printf("Oslo RADIUS term (0.5 / A) = %1.6e\n\n", Oslo_RADIUS);
printf("Values for Oslo - full curve on first/front surface:\n");
printf("RADIUS = %1.6e, B (4th order) = %1.6e, C (6th order) = %1.6e\n\n", Oslo_RADIUS*-1, B*-1, C*-1);

printf("Values for Oslo - full curve on second/back surface:\n");
printf("RADIUS = %1.6e, B (4th order) = %1.6e, C (6th order) = %1.6e\n\n", Oslo_RADIUS, B, C);

printf("(NEEDS MORE TESTING) Values for Oslo - half curve on each side:\n");
printf("FRONT: RADIUS = %1.6e, B (4th order) = %1.6e, C (6th order) = %1.6e\n", Oslo_RADIUS*-2, B*-0.5, C*-0.5);
printf("BACK: RADIUS = %1.6e, B (4th order) = %1.6e, C (6th order) = %1.6e\n\n", Oslo_RADIUS*2, B*0.5, C*0.5);
printf("Oslo HINT (as of Oslo version 6.4.3):\n  RADIUS is in Lens > Surface Data Spreadsheet\n");
printf("  B (4th) and C (6th) order terms are in Lens > Surface Data Spreadsheet > SPECIAL > Polynomial Asphere > Standard Asphere\n\n");

for (i = 0; i < (Dc / 2);) {
i=i+0.1;
z = T + (A*pow(i,2) + B*pow(i,4) + C*pow(i,6));
fprintf(fp, "%2.6f %2.6f\n", i, z);
}

z = T + (A*pow(h_n,2) + B*pow(h_n,4) + C*pow(h_n,6));

fprintf(fp2, "set title 'Dominics Schmidt corrector simulator - thickness as function of radius'\n");
fprintf(fp2, "set label '* %3.2f' at %f, %f\n", h_n, h_n, T*0.99);
fprintf(fp2, "set xlabel 'radius (mm)'\n");
fprintf(fp2, "set ylabel 'thickness (mm)'\n");
fprintf(fp2, "set data style lines\n");
fprintf(fp2, "plot [0:%f] [%f:%f] 'schcorr.dat'\n", Dc/2, z, T+0.0001);
fprintf(fp2, "set terminal postscript\n");
fprintf(fp2, "set output 'schcorr.ps'\n");
fprintf(fp2, "replot\n");
fprintf(fp2, "set terminal pbm\n");
fprintf(fp2, "set output 'schcorr.pbm'\n");
fprintf(fp2, "replot\n");



fclose(fp);
fclose(fp2);
return 0;
}