//This program use the Hermite Interpolation to either
//	-plot relative intensity as a function of rho
//	-find the location of the first fringe
#include <stdio.h>
#include <math.h>
#include "numcip.h"

double **bes;

void plotairy(void);
void findfringe(void);
double fbes(double r);
double airy(double);
void initbessel(void);
double dbes1(int n);

main()
{
	int choice;
	//user io
	printf("Please choose one of the two following functions\n");
	printf("\t[0]Plot the relative intensity of the airy pattern\n");
	printf("\t[1]Find the location of the first fringe\n");
	scanf("%d",&choice);
	
	initbessel();

	if(choice)
		findfringe();
	else
		plotairy();

	return 0;
}

void plotairy(void)
//this will plot the airy pattern at 0.1 increments from 0 to 10.0
{
	double rho;
	double lastpt;
	int i;
	FILE *outfile;

	outfile=fopen("airyplot.dat","w");
	fprintf(outfile,"rho\tairy\n");

	//at rho=0, by l'hospital rule,
	//	the relative intensity has a value of 1.0
	fprintf(outfile,"%lg\t%lg\n",0.0,1.0);


	rho=0.1;
	for(i=1;i<99;i++)
		{
		fprintf(outfile,"%lg\t%lg\n",rho,airy(rho));
		rho+=0.1;
		}
	//plotting the last point for the graph
	lastpt=bes[10][1]/5.0;
	fprintf(outfile,"%lg\t%lg\n",10.0,lastpt*lastpt);

}

#define tol 1e-12
void findfringe(void)
//this will use the Bisection Method to find the first fringe
//	the first fringe is located at the first zero of J_1(x)
{
	FILE *outfile;
	double newrho,rho,oldrho;
	double mbesrho,besrho;

	outfile=fopen("airyzero.dat","w");
	fprintf(outfile,"The location of the first fringe is at\n");

	//from the graph of J_1(x) on page 91, the first zero is located
	//	close to rho=3.0.  We will use the following bracket for our root
	rho=5;
	oldrho=3;
	while(fabs(rho-oldrho)>tol)
		{
		besrho=fbes(rho);//bessel at rho
		//here we take bisection step
		newrho=(rho+oldrho)/2.0;
		mbesrho=fbes(newrho);//bessel at midpoint
		if(mbesrho*besrho<0.0)//check to see in which half
			oldrho=newrho;
		else
			rho=newrho;
		};
	fprintf(outfile,"\t%20.15lg\n",rho);
	printf("The fringe is at: %20.15lg\n",rho);
	scanf("%ld",&rho);
	
}

double fbes(double r)
//give interpolated value of bessel
{
	int where;
	double x1,x2;
	double y;

	where=(int)r;//this will give index to the left pt in table for bessel
	
	x1=r-where;//x-x1
	x2=r-(where+1);//x-x2
	
	//J_1 Interpolation by Hermite Method
	//first two terms with function
	y=(1.0+2.0*x1)*x2*x2*bes[where][1]+(1.0-2.0*x2)*x1*x1*bes[where+1][1];
	//last two terms with derivative of J_1
	y+=(x1*x2*(x2*dbes1(where)+x1*dbes1(where+1)));
	//the denominator is 1.0
	//in this example, x1-x2=1

	return y;
}

double airy(double r)
//evaluate airy value by calling fbes()
{
	double y;
	
	//airy function
	y=2.0*fbes(r)/r;
	y=y*y;

	return y;
}

double dbes1(int n)
//this evaluate the derivative of the bessel function
{
	return (bes[n][0]-bes[n][2])/2.0;
}

void initbessel(void)
//bes[n][]-x value
//bes[][m]-order number
{
	bes=dmatrix(0,10,0,2);//this allocates space for bes (numcip.c)

	bes[0][0]=1.0;
	bes[0][1]=0.0;
	bes[0][2]=0.0;

	bes[1][0]=0.7651976866;
	bes[1][1]=0.4400505857;
	bes[1][2]=0.1149034849;

	bes[2][0]=0.2238907791;
	bes[2][1]=0.5767248078;
	bes[2][2]=0.3528340286;

	bes[3][0]=-0.2600519549;
	bes[3][1]=0.3390589585;
	bes[3][2]=0.4860912606;
	
	bes[4][0]=-0.3971498099;
	bes[4][1]=-0.0660433280;
	bes[4][2]=0.3641281459;

	bes[5][0]=-0.1775967713;
	bes[5][1]=-0.3275791376;
	bes[5][2]=0.0465651163;

	bes[6][0]=0.1506452573;
	bes[6][1]=-0.2766838581;
	bes[6][2]=-0.2428732100;
	
	bes[7][0]=0.3000792705;
	bes[7][1]=-0.0046828235;
	bes[7][2]=-0.3014172201;

	bes[8][0]=0.1716508071;
	bes[8][1]=0.2346363469;
	bes[8][2]=-0.1129917204;

	bes[9][0]=-0.0903336112;
	bes[9][1]=0.2453117866;
	bes[9][2]=0.1448473415;

	bes[10][0]=-0.2459357645;
	bes[10][1]=0.0434727462;
	bes[10][2]=0.2546303137;
}