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


#define DEG_IN_RAD    57.295779513
#define ROOT2         1.41421356237
#define A_IN_MM       1.0e7
#define MICR_IN_MM    1.0e3
#define GAMMA         31.   /* collimator-camera angle */

#define C_PIXEL       24.  /* Charlotte (Templeton) pixel size */
#define C_COLS        1024  /* Charlotte (Templeton) useful rows */  
#define W_PIXEL       15.
#define W_COLS        2048  
#define N_PIXEL       21.  /* Nellie */
#define N_COLS        2048
#define E_PIXEL       24.  /* Echelle */
#define E_COLS        1680  /* Echelle, rough vignetting limit */
#define SHORT_CAMF    85.  /* short camera focus, mm */
#define LONG_CAMF     200.
#define COLLF         429. /* collimator focus, mm */
#define TEL_SCALE     11.528  /* telescope scale, arcsec/mm at f7.5 */
#define GR1           600.   /* grating lines per mm */
#define BL1           8.63   /* blaze angle, degrees */
#define GR2           830.8
#define BL2           20.57  
#define GR3           1200.
#define BL3           17.45
#define PROJSLIT      1.5    /* 1.5 pixels is about the minimum slit size */

double tilt(lambda, lines, beta)
double lambda;  /* central wavelength */
double lines;   /* grating lines per mm */
double *beta;   /* angle toward camera from normal to grating */

/* Gives tilt of grating with respect to input beam 
to put wavelength alpha at angle gamma away from input beam --
gamma is the camera-collimator separation.  Angles are
given and returned in degrees. */
 
{
	double alpha;
	double cosgam, gamma, argu;
        double gr_ratio;  /* lambda over g, where g = spacing of grooves */
       
	gr_ratio = lambda * lines / A_IN_MM;

/*	printf("Each groove is %f wavelengths\n",(1./gr_ratio)); */
	
	gamma = GAMMA / DEG_IN_RAD; 
        cosgam = cos(gamma);
	

	argu = gr_ratio / (ROOT2 * pow((1.+cosgam),0.5));
	if(fabs(argu) <= 1.0) alpha = asin(argu);
		else {
			printf("No solution ... argument = %f\n",argu); 
			return(-99.);
		}
	alpha  =  DEG_IN_RAD * (alpha + gamma/2.);
	*beta = alpha - GAMMA; 
        return(alpha);
}


int get_setup(fcm, ln,  bl, pix, csz)

double *fcm, *ln, *bl, *pix, *csz;

{
	double fcam, lines, blaze, pixsize, chipsize;
	int port, grat, chip, ok;

        fcam = LONG_CAMF;
        lines = GR1; 
	blaze = BL1;
	pixsize = C_PIXEL;
	chipsize = C_COLS;

	printf("Type 1 for long port, 2 for short: ");
	scanf("%d",&port);
	if(port == 1) fcam = LONG_CAMF;
	else if (port == 2) fcam = SHORT_CAMF;
	else printf("Illegal choice, defaulting to long.\n");
        
	printf("Choose a grating:\n  1 -- 600 l/mm\n  2 -- 830.8 l/mm\n  3 -- 1200 l/mm\n");
	printf("Your choice: --> ");
	scanf("%d",&grat);
	if(grat == 1)       { lines = GR1; blaze = BL1; }
	else if (grat == 2) { lines = GR2; blaze = BL2; }
	else if (grat == 3) { lines = GR3; blaze = BL3; }
	else printf("Illegal choice, defaulting to %f\n",lines);

	printf("Choose a detector:\n  1 -- Templeton\n  2 -- Wilbur\n");
	printf("  3 -- Nellie\n  4 -- Echelle\n");
	printf("Your choice: --> ");
	scanf("%d",&chip);
	if(chip == 1) {
		pixsize = C_PIXEL;
		chipsize = C_COLS; 	
	}
	else if (chip == 2) {
		pixsize = W_PIXEL;
		chipsize = W_COLS; 	
	}
	else if (chip == 3) {
		pixsize = N_PIXEL;
		chipsize = N_COLS;
	}
	else if (chip == 4) {
		pixsize = E_PIXEL;
		chipsize = E_COLS;
	}
	else printf("Illegal choice, defaulting to Templetons's values.\n");

	printf("Camera focal length %5.1f mm\n",fcam);
	printf("Grating %7.2f lines / mm\n",lines);
	printf("Detector: %5.2f microns/pixel, %4.0f useful columns\n",
		pixsize,chipsize);

        *fcm = fcam;
	*ln  = lines;
	*bl = blaze;
	*pix = pixsize;
	*csz = chipsize;

}

double blazewln(lines, tilt, blaze)
double lines, tilt, blaze;

{
 	double wln;
	
	tilt = tilt / DEG_IN_RAD;  /* pass-by-value! */
	blaze = blaze / DEG_IN_RAD;

	wln = A_IN_MM * (sin(tilt) + sin(2. * blaze - tilt)) / lines;
	
	return(wln);
}
	
double dispersion(lines, fcam, beta, pixsize) 
double lines, fcam, beta, pixsize;

{
	double disp;
	
	disp = A_IN_MM * cos(beta / DEG_IN_RAD) / (lines * fcam);
	   /* in Angstroms per mm */
/*	printf("Dispersion %f A/mm\n",disp); */
	disp = disp * pixsize / MICR_IN_MM;
        return(disp);
}

double demagnification(fcam, fcoll, alpha, beta)
double fcam, fcoll, alpha, beta;
{
	double demag;
	
	demag = (fcam / fcoll) * cos(alpha / DEG_IN_RAD) /
				 cos(beta / DEG_IN_RAD);
	return(demag);
}

void ends(mid,disp,chipsize,shortw,longw)
double mid, disp, chipsize, *shortw, *longw;
{
	*shortw = mid - (chipsize/2.) * disp;
	*longw = mid + (chipsize/2.) * disp;
}
	
main() {
       
	double lambda = 5000.;
	double dlam = 100.;  /* for "numerical derivative" */
	double alpha;  /* grating angle */
	double alpha2;
	double beta;   /* angle between normal to grating and camera */
	double beta2;
	double trueblaze;  /* blaze wavelength at this incident angle */
	double disp, demag;
        double fcam, fcoll = COLLF, lines, blazeang, pixsize, chipsize;
	double match, slitscale;
	double shortw, longw;
	int more = 1;

	printf("\nModspec setup program -- John Thorstensen, April 1994.\n\n");

	printf("** NOTE ** Calculations are based on nominal parameters -- you should\n");
	printf("           allow a little room for error.\n\n");
while(more == 1) {  /* funny indentation, sorry */
	get_setup(&fcam, &lines, &blazeang, &pixsize, &chipsize);

	while(lambda > 0.) {
		printf("\nGive central wavelength, negative exits : ");
		scanf("%lf",&lambda);
		if(lambda <= 0.) break;
	        alpha = tilt(lambda,lines,&beta);
/*		printf("Alpha = %f\n",alpha);
                printf("Beta = %f\n",beta);  */
		disp = dispersion(lines, fcam, beta, pixsize); 
		trueblaze = blazewln(lines, alpha, blazeang);

		printf("\n");
		if(fcam == LONG_CAMF) 
		   printf("Port                        :    Long\n");
		else if (fcam == SHORT_CAMF)
		   printf("Port                        :    Short\n");
		printf("Grating : %7.2f lines/mm, blaze angle %5.2f degrees\n",
				lines, blazeang);
		if(pixsize == C_PIXEL && chipsize == C_COLS) 
		   printf("Detector                    :    Charlotte\n");
		else if (pixsize == W_PIXEL && chipsize == W_COLS) 
		   printf("Detector                    :    Wilbur\n");
		else if (pixsize == N_PIXEL && chipsize == N_COLS) 
		   printf("Detector                    :    Nellie [MAY VIGNETTE?]\n");
		else if (pixsize == E_PIXEL && chipsize == E_COLS) 
		   printf("Detector                    :    Echelle [vignettes]\n");
		else 
		   printf("Detector                    :    ???\n");

		printf("Central wavelength (A)      :   %5.0f\n",lambda);
		printf("Grating angle (degrees)     :    %5.2f\n",alpha);
		printf("Dispersion (A per pixel)    :   %5.2f\n",disp);
		demag = demagnification(fcam,fcoll,alpha,beta);
		printf("Slit demagnification factor :    %5.3f in dispersion direction\n",demag);
		match = TEL_SCALE * PROJSLIT * pixsize /(demag * MICR_IN_MM);
		printf("Matching (%4.1f pix) slit    :   %5.2f arcsec\n",
				PROJSLIT,match);
		slitscale = TEL_SCALE * pixsize * fcoll / (fcam * MICR_IN_MM);
		printf("scale along slit (arcsc/pix):    %5.3f\n",
			slitscale);
                ends(lambda,disp,chipsize,&shortw,&longw);
	        printf("Short end (Angstroms)       :   %5.0f\n",shortw);
		printf("Long end (Angstroms)        :   %5.0f\n",longw);
		printf("Effective blaze wln (A)     :   %5.0f\n",trueblaze);
		alpha2 = tilt(lambda + dlam, lines, &beta2);
		printf("Grating adjustment factor   :    %5.3f",
			(alpha2 - alpha),dlam);
		printf(" degrees per %4.0f Angstroms.\n\n",dlam);
	}
	printf("Type 1 to try other port, grating etc, or 0 to exit :");
	scanf("%d",&more);
	lambda = 5000.;  /* to allow looping */
}  /* closes do while more */	
}