#include <math.h>
#include <complex>
#include <iostream>
#include <fstream>
#include "globals.h"
//#include "cufft.h"
using namespace std;

#define pi 3.14159

typedef complex<float> scComplex;

int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
    complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
int bessjyv(double v,double x,double &vm,double *jv,double *yv,
    double *djv,double *dyv);

complex<double> Jl_neg(complex<double> x)
{
	//this function computes the bessel function of the first kind Jl(x) for l = -0.5
	return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}

double Jl_neg(double x)
{
	//this function computes the bessel function of the first kind Jl(x) for l = -0.5
	return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
}

double Yl_neg(double x)
{
	//this function computes the bessel function of the second kind Yl(x) for l = -0.5;
	return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
}

void computeB(complex<float>* B, float radius, complex<double> refIndex, float lambda, int Nl)
{
	double k = (2*pi)/lambda;
	int b = 2;

	//allocate space for the real bessel functions
	double* jv = (double*)malloc(sizeof(double)*(Nl+b));
	double* yv = (double*)malloc(sizeof(double)*(Nl+b));
	double* jvp = (double*)malloc(sizeof(double)*(Nl+b));
	double* yvp = (double*)malloc(sizeof(double)*(Nl+b));

	//allocate space for the complex bessel functions
	complex<double>* cjv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
	complex<double>* cyv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
	complex<double>* cjvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
	complex<double>* cyvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));

	double kr = k*radius;
	complex<double> knr = k*refIndex*(double)radius;
	complex<double> n = refIndex;

	//compute the bessel functions for k*r
	double vm;// = Nl - 1;
	bessjyv((Nl)+0.5, kr, vm, jv, yv, jvp, yvp);
	//cout<<"Nl: "<<Nl<<"  vm: "<<vm<<endl;
	//printf("Nl: %f, vm: %f\n", (float)Nl, (float)vm);

	//compute the bessel functions for k*n*r
	cbessjyva((Nl)+0.5, knr, vm, cjv, cyv, cjvp, cyvp);

	//scale factor for spherical bessel functions
	double scale_kr = sqrt(pi/(2.0*kr));
	complex<double> scale_knr = sqrt(pi/(2.0*knr));

	complex<double> numer, denom;
	double j_kr;
	double y_kr;
	complex<double> j_knr;
	complex<double> j_d_knr;
	double j_d_kr;
	complex<double> h_kr;
	complex<double> h_d_kr;
	complex<double> h_neg;
	complex<double> h_pos;

	//cout<<"B coefficients:"<<endl;
	for(int l=0; l<Nl; l++)
	{
		//compute the spherical bessel functions
		j_kr = jv[l] * scale_kr;
		y_kr = yv[l] * scale_kr;
		j_knr = cjv[l] * scale_knr;

		//compute the Hankel function
		h_kr = complex<double>(j_kr, y_kr);

		//compute the derivatives
		if(l == 0)
		{
			//spherical bessel functions for l=0
			j_d_kr = scale_kr * (Jl_neg(kr) - (jv[l] + kr*jv[l+1])/kr )/2.0;
			j_d_knr = scale_knr * ( Jl_neg(knr) - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
			h_neg = complex<double>(scale_kr*Jl_neg(kr), scale_kr*Yl_neg(kr));
			h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
			h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
		}
		else
		{
			//spherical bessel functions
			j_d_kr = scale_kr * (jv[l-1] - (jv[l] + kr*jv[l+1])/kr )/2.0;
			j_d_knr = scale_knr * ( cjv[l-1] - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
			h_neg = complex<double>(scale_kr*jv[l-1], scale_kr*yv[l-1]);
			h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
			h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
		}

		numer = j_kr*j_d_knr*n - j_knr*j_d_kr;
		denom = j_knr*h_d_kr - h_kr*j_d_knr*n;
		complex<double> temp =  numer/denom;

		B[l] = scComplex(temp.real(), temp.imag());
		//cout<<B[l]<<endl;
	}

	free(jv);
	free(yv);
	free(jvp);
	free(yvp);
	free(cjv);
	free(cyv);
	free(cjvp);
	free(cyvp);
}

void Legendre(float* P, float x, int Nl)
{
	//computes the legendre polynomials from orders 0 to Nl-1
	P[0] = 1;
	if(Nl == 1) return;
	P[1] = x;
	for(int l = 2; l < Nl; l++)
	{
		P[l] = ((2*l - 1)*x*P[l-1] - (l - 1)*P[l-2])/l;
	}	

}

complex<float> integrateUi(float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
{
	/*This function integrates the incident field of magnitude M in the far zone
	in order to evaluate the field at the central pixel of a detector.
	cNAi = condenser inner angle
	cNAo = condenser outer angle
	oNAi = objective inner angle
	oNAo = objective outer angle
	M = field magnitude*/

	float alphaIn = max(cAngleI, oAngleI);
	float alphaOut = min(cAngleO,oAngleO);

	complex<float> Ui;
	if(alphaIn > alphaOut)
		Ui = complex<float>(0.0, 0.0);
	else
		Ui = complex<float>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);

	return Ui;

}

void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO)
{
	/*This function computes the condenser integral in order to build the field of incident light
	alpha = list of Nl floating point values representing the condenser alpha as a function of l
	Nl = number of orders in the incident field
	cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/

	//compute the Legendre polynomials for the condenser aperature
	float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAo, cos(cAngleO), Nl+1);
	float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAi, cos(cAngleI), Nl+1);

	for(int l=0; l<Nl; l++)
	{
		//integration term
		if(l == 0)
			alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
		else
			alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];

		alpha[l] *= 2 * pi;
	}

}

complex<float> integrateUs(float r, float lambda, complex<float> eta, 
						   float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
{
	/*This function integrates the incident field of magnitude M in the far zone
	in order to evaluate the field at the central pixel of a detector.
	r = sphere radius
	lambda = wavelength
	eta = index of refraction
	cNAi = condenser inner NA
	cNAo = condenser outer NA
	oNAi = objective inner NA
	oNAo = objective outer NA
	M = field magnitude*/

	//compute the required number of orders
	float k = 2*pi/lambda;
	int Nl = ceil( k + 4 * exp(log(k*r)/3) + 3 );

	//compute the material coefficients B
	complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
	//compute the Legendre polynomials for the condenser and objective aperatures
	float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAo, cos(cAngleO), Nl+1);
	float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAi, cos(cAngleI), Nl+1);

	float* PoNAo = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PoNAo, cos(oAngleO), Nl+1);
	float* PoNAi = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PoNAi, cos(oAngleI), Nl+1);

	//store the index of refraction;
	complex<double> IR(eta.real(), eta.imag());
	
	//compute the scattering coefficients
	computeB(B, r, IR, lambda, Nl);

	//aperature terms for the condenser (alpha) and objective (beta)
	float alpha;
	float beta;
	float c;
	complex<float> Us(0.0, 0.0);

	for(int l=0; l<Nl; l++)
	{
		//integration term
		if(l == 0)
		{
			alpha = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
			beta = -PoNAo[l+1] + PoNAo[0] + PoNAi[l+1] - PoNAi[0];
		}
		else
		{
			alpha = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
			beta = -PoNAo[l+1] + PoNAo[l-1] + PoNAi[l+1] - PoNAi[l-1];
		}
		c = (2*pi)/(2.0 * l + 1.0);
		Us += c * alpha * beta * B[l] * M;
		
			
	}
	free(PcNAo);
	free(PcNAi);
	free(PoNAo);
	free(PoNAi);
	free(B);

	return Us;

}

void pointSpectrum()
{
	PD.StartTimer(SIMULATE_SPECTRUM);
	//clear the previous spectrum
	SimSpectrum.clear();

	float dNu = 2.0f;
	float lambda;
	
	//compute the angles based on NA
	float cAngleI = asin(cNAi);
	float cAngleO = asin(cNAo);
	float oAngleI = asin(oNAi);
	float oAngleO = asin(oNAo);

	//implement a reflection-mode system if necessary
	if(opticsMode == ReflectionOpticsType){
		
		//set the condenser to match the objective
		cAngleI = oAngleI;
		cAngleO = oAngleO;

		//invert the objective
		oAngleO = pi - cAngleI;
		oAngleI = pi - cAngleO;
	}

	//integrate the incident field at the detector position
	complex<float> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
	float I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
	I0 *= scaleI0;

	

	float I;
	SpecPair temp;
	float nu;
	complex<float> eta;
	complex<float> Us, U;

	float vecLen = 0.0;
	for(int i=0; i<EtaK.size(); i++)
	{
		nu = EtaK[i].nu;
		lambda = 10000.0f/nu;
		if(applyMaterial)
			eta = complex<float>(EtaN[i].A, EtaK[i].A);
		else
			eta = complex<float>(baseIR, 0.0);


		//integrate the scattered field at the detector position
		Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
		U = Us + Ui;
		float I = U.real() * U.real() + U.imag() * U.imag();
		
		temp.nu = nu;

		//set the spectrum value based on the current display type
		if(dispSimType == AbsorbanceSpecType)
			temp.A = -log10(I/I0);
		else
			temp.A = I;

		if(dispNormalize)
			vecLen += temp.A * temp.A;
		
		SimSpectrum.push_back(temp);		
	}
	vecLen = sqrt(vecLen);

	if(dispNormalize)
		for(int i=0; i<SimSpectrum.size(); i++)
			SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;

	PD.EndTimer(SIMULATE_SPECTRUM);
}

/*
complex<float> sampleUs(complex<float>* B, float* Alpha, int Nl, float r, 
						   float cAngleI, float cAngleO, float theta, float M = 2*pi)
{
	/*This function takes a point sample of the scattered field in the far zone
	in order to evaluate the field at the central pixel of a detector.
	r = sphere radius
	lambda = wavelength
	eta = index of refraction
	cNAi = condenser inner NA
	cNAo = condenser outer NA
	theta = angle of the sample
	M = field magnitude*/

/*
	//compute the material coefficients B
	//complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);

	//compute the Legendre polynomials for theta (at the objective)
	float* Ptheta = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(Ptheta, cos(theta), Nl+1);

	//complex<double> IR(eta.real(), eta.imag());

	//aperature terms for the condenser (alpha) and objective (beta)
	float beta;
	float c;
	complex<float> Us(0.0, 0.0);

	for(int l=0; l<Nl; l++)
	{

		complex<float> numer(0.0, -(l*pi)/2.0);
		Us += B[l] * exp(numer) * Ptheta[l] * Alpha[l] * pow(complex<float>(0.0, 1.0), l);
		
			
	}
	//printf("Ptheta: %f\n", Ptheta[Nl-1]);

	return Us;

}*/


/*
void detectorSpectrum(int numSamples)
{
	//integrate across the objective aperature and calculate the resulting intensity on a detector
	PD.StartTimer(SIMULATE_SPECTRUM);
	//clear the previous spectrum
	SimSpectrum.clear();

	float dNu = 2.0f;
	float lambda;
	
	//compute the angles based on NA
	float cAngleI = asin(cNAi);
	float cAngleO = asin(cNAo);
	float oAngleI = asin(oNAi);
	float oAngleO = asin(oNAo);

	//implement a reflection-mode system if necessary
	if(opticsMode == ReflectionOpticsType){
		
		//set the condenser to match the objective
		cAngleI = oAngleI;
		cAngleO = oAngleO;

		//invert the objective
		oAngleO = pi - cAngleI;
		oAngleI = pi - cAngleO;
	}

	//compute Nl (maximum order of the spectrum)
	//****************************************************************************
	float maxNu = EtaK.back().nu;
	float maxLambda = 10000.0f/maxNu;
	float k = 2*pi/maxLambda;
	int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );
	int nLambda = EtaK.size();

	//compute alpha (condenser integral)
	//****************************************************************************
	//compute the Legendre polynomials for the condenser aperature
	float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAo, cos(cAngleO), Nl+1);
	float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
	Legendre(PcNAi, cos(cAngleI), Nl+1);

	//allocate space for the alpha array
	float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
	computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
	
	for(int l=0; l<Nl; l++)
	{
		//integration term
		if(l == 0)
			alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
		else
			alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];

		alpha[l] *= 2 * pi;
	}

	//compute the information based on wavelength and sphere

	//evaluate the incident field intensity
	float I0 = 0.0;
	float theta;
	float dTheta = (oAngleO - oAngleI)/numSamples;
	complex<float> Ui;

	Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
	I0 = Ui.real()*2*pi;

	float I;
	SpecPair temp;
	float nu;
	complex<float> eta;
	complex<float> Us, U;
	
	

	float vecLen = 0.0;
	for(int i=0; i<EtaK.size(); i++)
	{
		//compute information based on wavelength and material
		nu = EtaK[i].nu;
		lambda = 10000.0f/nu;
		if(applyMaterial)
			eta = complex<float>(EtaN[i].A, EtaK[i].A);
		else
			eta = complex<float>(baseIR, 0.0);

		//allocate memory for the scattering coefficients
		complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);		
		
		complex<double> IR(eta.real(), eta.imag());
		computeB(B, radius, IR, lambda, Nl);


		//integrate the scattered field at the detector position
		I = 0.0;
		
		for(int iTheta = 0; iTheta < numSamples; iTheta++)
		{
			theta = oAngleI + iTheta * dTheta;
			Us = sampleUs(B, alpha, Nl, radius, cAngleI, cAngleO, theta, 2*pi);
			

			//calculate the intensity and add
			if(theta >= cAngleI && theta <= cAngleO)
				U = Us + 2*(float)pi;
			else
				U = Us;

			I += (U.real() * U.real() + U.imag() * U.imag()) * sin(theta) * 2 * pi * dTheta;
			
		}
		
		temp.nu = nu;

		if(i == 0)
			printf("I: %f        I0: %f\n", I, I0);

		//set the spectrum value based on the current display type
		if(dispSimType == AbsorbanceSpecType)
			temp.A = -log10(I/I0);
		else
			temp.A = I;

		if(dispNormalize)
			vecLen += temp.A * temp.A;
		//temp.A = Us.real();
		SimSpectrum.push_back(temp);	

		free(B);
	}
	vecLen = sqrt(vecLen);

	if(dispNormalize)
		for(int i=0; i<SimSpectrum.size(); i++)
			SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;

	free(alpha);

	PD.EndTimer(SIMULATE_SPECTRUM);
}*/

void updateSpectrum(float* I, float I0, int n)
{
	SimSpectrum.clear();
	SpecPair temp;

	//update the displayed spectrum based on the computed intensity I
	for(int i=0; i<n; i++)
	{
		temp.nu = EtaK[i].nu;

		//set the spectrum value based on the current display type
		if(dispSimType == AbsorbanceSpecType)
			temp.A = -log10(I[i]/I0);
		else
			temp.A = I[i];

		SimSpectrum.push_back(temp);
	}
}

void computeCassegrainAngles(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO)
{
	//compute the angles based on NA
	cAngleI = asin(cNAi);
	cAngleO = asin(cNAo);
	oAngleI = asin(oNAi);
	oAngleO = asin(oNAo);

	//implement a reflection-mode system if necessary
	if(opticsMode == ReflectionOpticsType){
		
		//set the condenser to match the objective
		cAngleI = oAngleI;
		cAngleO = oAngleO;

		//invert the objective
		oAngleO = pi - cAngleI;
		oAngleI = pi - cAngleO;
	}


}

int computeNl()
{
	float maxNu = EtaK.back().nu;
	float maxLambda = 10000.0f/maxNu;
	float k = 2*pi/maxLambda;
	int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );

	return Nl;
}

void computeBArray(complex<float>* B, int Nl, int nLambda)
{
	float nu;
	complex<float> eta;
	float* Lambda = (float*)malloc(sizeof(float) * nLambda);

	//for each wavenumber nu
	for(int i=0; i<EtaK.size(); i++)
	{
		//compute information based on wavelength and material
		nu = EtaK[i].nu;
		Lambda[i] = 10000.0f/nu;
		if(applyMaterial)
			eta = complex<float>(EtaN[i].A, EtaK[i].A);
		else
			eta = complex<float>(baseIR, 0.0);

		//allocate memory for the scattering coefficients
		//complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);		

		complex<double> IR(eta.real(), eta.imag());
		computeB(&B[i * Nl], radius, IR, Lambda[i], Nl);
	}
}

void computeOpticalParameters(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO, float& I0, float* alpha, int Nl)
{
	computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO);	

	//evaluate the incident field intensity
	I0 = 0.0;
	float theta;
	complex<float> Ui;

	Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
	I0 = Ui.real()*2*pi;

	//compute alpha (condenser integral)
	computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
}

void gpuDetectorSpectrum(int numSamples)
{
	//integrate across the objective aperature and calculate the resulting intensity on a detector
	PD.StartTimer(SIMULATE_SPECTRUM);
	//clear the previous spectrum
	SimSpectrum.clear();	

	//compute Nl (maximum order of the spectrum)
	int Nl = computeNl();

	float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
	float cAngleI, cAngleO, oAngleI, oAngleO, I0;
	computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);

	//allocate space for a list of wavelengths
	int nLambda = EtaK.size();
	
	//allocate space for the 2D array (Nl x nu) of scattering coefficients
	complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
	computeBArray(B, Nl, nLambda);
	

	//allocate temporary space for the spectrum
	float* I = (float*)malloc(sizeof(float) * EtaK.size());

	//compute the spectrum on the GPU
	PD.StartTimer(SIMULATE_GPU);
	cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
	PD.EndTimer(SIMULATE_GPU);

	updateSpectrum(I, I0, nLambda);

	PD.EndTimer(SIMULATE_SPECTRUM);

}

void SimulateSpectrum()
{
	if(pointDetector)
		pointSpectrum();
	else
		gpuDetectorSpectrum(objectiveSamples);
		//detectorSpectrum(objectiveSamples);
}

float absorbanceDistortion(){

	//compute the mean of the spectrum
	float sumSim = 0.0;
	for(int i=0; i<SimSpectrum.size(); i++)
	{
		sumSim += SimSpectrum[i].A;
	}
	float meanSim = sumSim/SimSpectrum.size();

	//compute the distortion (MSE from the mean)
	float sumSE = 0.0;
	for(int i=0; i<SimSpectrum.size(); i++)
	{
		sumSE += pow(SimSpectrum[i].A - meanSim, 2);
	}
	float MSE = sumSE / SimSpectrum.size();

	return MSE;
}

float intensityDistortion(){

	//compute the magnitude of the spectrum
	float sumSim = 0.0;
	for(int i=0; i<SimSpectrum.size(); i++)
	{
		sumSim += SimSpectrum[i].A * SimSpectrum[i].A;
	}
	float magSim = sqrt(sumSim);

	//compute the distortion (MSE from the mean)
	float sumSE = 0.0;
	for(int i=0; i<SimSpectrum.size(); i++)
	{
		sumSE += (SimSpectrum[i].A/magSim) * (1.0/SimSpectrum.size());
	}
	float MSE = sumSE;

	return MSE;
}

void MinimizeDistortion(){
	ofstream outFile("distortion.txt");

	//set the parameters for the distortion simulation
	float step = 0.001;

	oNAi = 0.2;
	oNAo = 0.5;

	//compute the optical parameters
	//compute Nl (maximum order of the spectrum)
	int Nl = computeNl();

	float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
	float cAngleI, cAngleO, oAngleI, oAngleO, I0;
	
	//allocate space for a list of wavelengths
	int nLambda = EtaK.size();

	//allocate temporary space for the spectrum
	float* I = (float*)malloc(sizeof(float) * EtaK.size());

	//calculate the material parameters
	//allocate space for the 2D array (Nl x nu) of scattering coefficients
	complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
	computeBArray(B, Nl, nLambda);



	float D;
	float e = 0.001;
	for(float i=0.0; i<=oNAo-step; i+=step)
	{
		
		for(float o=oNAi+step; o<=1.0; o+=step)
		{
			
			
			//set the current optical parameters
			cNAi = i;
			cNAo = o;

			//compute the optical parameters
			computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);

			//simulate the spectrum
			cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);
			updateSpectrum(I, I0, nLambda);

			if(dispSimType == AbsorbanceSpecType)
			{
				if(i + e >= o || i + e >= oNAo || oNAi + e >= o || oNAi + e >= oNAo)
					D = 0.0;
				else
					D = absorbanceDistortion();
			}
			else
			{
				if(i >= o || oNAi >= oNAo)
					D=0;
				else
					D = intensityDistortion();
			}
			outFile<<D<<"       ";
		}
		outFile<<endl;
		cout<<i<<endl;
	}
	outFile.close();
}