/* Plots an F-P transmission curve
 *
 * Fraser Clarke (fclarke@ast.cam.ac.uk)
 * Jan 2000 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cpgplot.h>
#include <ada.h>

#define COS(x) (cos(M_PI*(x)/180.0))
#define SIN(x) (sin(M_PI*(x)/180.0))
#define SQR(x) ((x)*(x))

double FPtrans(double lambda, double theta, double R, double T, double gap, double ri) {
	double trans;
	double t1,t2,t3;
	
	t1 = SQR(T)/SQR(1-R);
	t2 = 4.0*R/SQR(1-R);
	t3 = 2.0*M_PI*gap*ri*COS(theta)/lambda;

	trans = t1/(1.0 + t2*SQR(sin(t3)));
	
	return trans;
}

int main() {

	char tmp[255];
	double R,T,sep,refindex,theta;	// Input parameters
	double sw,cw,fsr,fwhm,step,peak;
	int order, i;
	Spectrum *profile;
		
	/* Get some parameters from the User */
	
	printf("(Intensity) Relectivity > ");
	fgets(tmp,255,stdin);
	R=atof(tmp);
	
	printf("(Intensity) Transmission > ");
	fgets(tmp,255,stdin);
	T=atof(tmp);
	
	printf("Plate Seperation (microns) > ");
	fgets(tmp,255,stdin);
	sep=atof(tmp)/1.0E6;
	
	printf("Refractive Index > ");
	fgets(tmp,255,stdin);
	refindex=atof(tmp);
	
//	printf("Central Wavelength > ");
//	fgets(tmp,255,stdin);
//	cw=atof(tmp);
	
	printf("Angle of incidence (degrees) > ");
	fgets(tmp,255,stdin);
	theta=atof(tmp);
	
	printf("Order > ");
	fgets(tmp,255,stdin);
	order=atoi(tmp);
	
	/* Calculate the central wavelength and FSR */
	
	cw = 2.0*sep*refindex*cos(M_PI*theta/180.0)/(double)order;
	fsr = cw/(double)order;
	fwhm = cw*(1.0-R)/((double)order*M_PI*sqrt(R));
	
	printf("\nCentral Wavelength : %.1f Angstroms\n",cw*1.0E10);
	printf("Free Spectral Range : %.1f Angstroms\n",fsr*1.0E10);
	printf("Full Width Half Maximum : %.2f Angstroms\n",fwhm*1.0E10);

	/* Open a PGPLOT device */
	
	cpgopen("?");

	profile = specalloc(1000);

	/* Now we will calculate the profile over the one FSR */
		
	sw = cw - fsr/2.0;
	step = fsr/1000.0;
	
	// Calculate the peak transimission so we can normalise the curve.
	peak = SQR(T)/SQR(1-R);
	
	for(i=0;i<1000;i++) {
		profile->wave[i]=sw+(double)i*step;
		profile->flux[i]=FPtrans(profile->wave[i],theta,R,T,sep,refindex)/peak;
		
		// Change the wavelength to Angstroms
		
		profile->wave[i]*=1.0e10;
	}
	
	plotspec(profile);
	cpglab("Wavelength (Angstroms)","Relative Transmission","Transmission curve of F-P over one FSR");
	
	/* Now we will calculate the profile over a smaller range between +/- 2HWHM */

	sw = cw - fwhm;
	step = 2.0*fwhm/1000.0;
	
	for(i=0;i<1000;i++) {
		profile->wave[i]=sw+(double)i*step;
		profile->flux[i]=FPtrans(profile->wave[i],theta,R,T,sep,refindex)/peak;
		
		// Change the wavelength to Angstroms
		
		profile->wave[i]*=1.0e10;	
	}
	
	plotspec(profile);
	cpglab("Wavelength (Angstroms)","Relative Transmission","Transmission curve of F-P over 2x FWHM");
	
	cpgclos();
	
	/* Bye */
	
	freespec(profile);
	
	return 1;
}