optics / imie

  #include "globals.h"

  #include <stdlib.h>

  #define PI 3.14159
  

  double CalculateError(double* E)

  {
  	//Calculate the error between the Reference Spectrum and the Simulated Spectrum

  	double sumE = 0.0;

  	int nVals = RefSpectrum[currentSpec].size();

  	double nu;

  	for(int i=0; i<nVals; i++)
  	{
  		nu = RefSpectrum[currentSpec][i].nu;
  		E[i] = RefSpectrum[currentSpec][i].A + nu * refSlope - SimSpectrum[i].A;
  		sumE += E[i]*E[i];
  	}
  
  	return sumE/nVals;
  }
  

  void EstimateK(double* E)

  {
  	int nVals = RefSpectrum[currentSpec].size();

  	double nuStart = RefSpectrum[currentSpec].front().nu;
  	double nuEnd = RefSpectrum[currentSpec].back().nu;

  	double r = radius/10000.0;
  	double nu;
  	double dNu = (nuEnd - nuStart)/(nVals-1);
  	double eScale;

  	for(int i=0; i<nVals; i++)
  	{
  		nu = nuStart + i*2;
  
  		eScale = 1/(8*PI*r*nu);
  		EtaK[i].A = EtaK[i].A + eScale * E[i];
  		if(EtaK[i].A < 0.0) EtaK[i].A = 0.0;
  	}
  }
  
  void EstimateMaterial()
  {
  	/*This function estimates the material properties of a sphere based on the
  	input spectrum RefSpectrum and the optical properties of the system.
  	1) The material properties are stored in EtaK and EtaN
  	2) The best fit is stored in SimSpectrum*/
  
  	//initialize the material index of refraction
  	EtaN.clear();	
  	EtaK.clear();
  
  	//insert the default material properties
  	SpecPair temp;

  	for(int s=0; s<(int)RefSpectrum[currentSpec].size(); s++)

  	{
  		//the real part of the IR is the user-specified baseline IR
  		temp.nu = RefSpectrum[currentSpec][s].nu;
  		temp.A = baseIR;
  		EtaN.push_back(temp);
  		//the imaginary part of the IR is zero absorbance
  		temp.A = 0.0f;
  		EtaK.push_back(temp);
  	}
  
  
  	//allocate space to store the list of error values

  	double* E = (double*)malloc(sizeof(double) * RefSpectrum[currentSpec].size());

  	//copy the absorbance values into a linear array

  	double* k = (double*)malloc(sizeof(double) * EtaK.size());
  	double* n = (double*)malloc(sizeof(double) * EtaN.size());

  
  	//iterate to solve for both n and k

  	double sumE = 99999.9;
  	int j=0;

  	//clear the console
  	system("cls");

  	while(sumE > minMSE && j < maxFitIter)

  	{	
  		//simulate a spectrum based on the current IR		
  		SimulateSpectrum();
  
  		//calculate the error term
  		sumE = CalculateError(E);
  
  		//estimate the new absorbance
  		EstimateK(E);
  
  		//use Kramers-Kronig to compute n
  		

  		for(unsigned int i=0; i<EtaK.size(); i++)

  			k[i] = EtaK[i].A;
  		cudaKramersKronig(n, k, EtaK.size(), EtaK.front().nu, EtaK.back().nu, baseIR);
  
  		//copy the real part of the index of refraction into the vector
  		EtaN.clear();

  		for(int i=0; i<(int)EtaK.size(); i++)

  		{
  			temp.nu =  EtaK[i].nu;
  			temp.A = n[i];
  			EtaN.push_back(temp);
  		}
  
  		cout<<"    E = "<<sumE<<endl;

  		j++;

  		//SaveSpectrum(n, nVals, "simNj.txt");
  		//SaveSpectrum(k, nVals, "simKj.txt");
  		//SaveSpectrum(simSpec, nVals, "simSpec.txt");
  		//exit(1);
  		
  	}
  
  	free(E);
  	free(k);
  	free(n);
  
  
  
  
  }