optics / imie

  #include <fstream>
  using namespace std;
  #include "interactivemie.h"
  #include <QtGui/QApplication>
  #include "qtSpectrumDisplay.h"
  #include "globals.h"
  #include "rtsGUIConsole.h"
  #include "PerformanceData.h"
  #include <complex>

  //#include <direct.h>

  
  PerformanceData PD;
  
  qtSpectrumDisplay* gpSpectrumDisplay;
  

  vector<vector<SpecPair> > RefSpectrum;

  vector<SpecPair> SimSpectrum;
  vector<SpecPair> EtaK;
  vector<SpecPair> EtaN;
  int currentSpec = 0;
  

  double nuMin = 800;
  double nuMax = 4000;
  double dNu = 2;

  double aMin = 0;
  double aMax = 1;

  double scaleI0 = 1.0;
  double refSlope = 0.0;

  
  bool dispRefSpec = true;
  bool dispSimSpec = true;
  bool dispSimK = true;
  bool dispMatK = true;
  bool dispSimN = true;
  bool dispMatN = true;

  double dispScaleK = 1.0;
  double dispScaleN = 1.0;

  SpecType dispSimType = AbsorbanceSpecType;
  bool dispNormalize = false;

  double dispNormFactor = 1.0;

  
  
  //material parameters

  double radius = 4.0f;
  double baseIR = 1.49f;
  double cA = 1.0;

  //vector<SpecPair> KMaterial;
  //vector<SpecPair> NMaterial;
  bool applyMaterial = true;
  vector<Material> MaterialList;
  int currentMaterial = 0;
  
  //optical parameters

  double cNAi = 0.0;
  double cNAo = 0.6;
  double oNAi = 0.0;
  double oNAo = 0.6;

  OpticsType opticsMode = TransmissionOpticsType;
  bool pointDetector = false;
  int objectiveSamples = 200;
  
  //fitting parameters

  double minMSE = 0.00001;

  int maxFitIter = 20;
  
  void TempSimSpectrum()
  {
  	SpecPair temp;
  	for(int i=800; i<4000; i++)
  	{
  		temp.nu = i;

  		temp.A = sin((double)i/200);

  		SimSpectrum.push_back(temp);
  	}
  }
  
  void UpdateDisplay(){
  	gpSpectrumDisplay->updateGL();
  }
  
  void LoadMaterial(string fileNameK, string fileNameN, string materialName)
  {
  	Material newMaterial;
  	newMaterial.name = materialName;
  
  	vector<SpecPair> KMaterial = LoadSpectrum(fileNameK.c_str());
  	vector<SpecPair> NMaterial = LoadSpectrum(fileNameN.c_str());
  
  	//make sure that the sizes are the same
  	if(KMaterial.size() != NMaterial.size()){
  		cout<<"Error, material properties don't match."<<endl;
  		exit(1);
  	}
  

  	complex<double> eta;
  	//int j;
  	for(unsigned int i=0; i<KMaterial.size(); i++){

  		newMaterial.nu.push_back(KMaterial[i].nu);

  		eta = complex<double>(NMaterial[i].A, KMaterial[i].A);

  		newMaterial.eta.push_back(eta);
  	}
  	MaterialList.push_back(newMaterial);
  
  }
  
  void LoadMaterial(string fileNameK, string materialName){
  
  	//load the material absorbance
  	vector<SpecPair> KMaterial = LoadSpectrum(fileNameK.c_str());
  	vector<SpecPair> NMaterial;
  	//KMaterial = LoadSpectrum("eta_TolueneK.txt");
  
  	//compute the real IR using Kramers Kronig
  	//copy the absorbance values into a linear array

  	double* k = (double*)malloc(sizeof(double) * KMaterial.size());
  	double* n = (double*)malloc(sizeof(double) * KMaterial.size());
  	for(unsigned int i=0; i<KMaterial.size(); i++)

  		k[i] = KMaterial[i].A;
  
  	//use Kramers Kronig to determine the real part of the index of refraction
  	cudaKramersKronig(n, k, KMaterial.size(), KMaterial[0].nu, KMaterial.back().nu, baseIR);
  	SpecPair temp;

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

  	{
  		temp.nu =  KMaterial[i].nu;
  		temp.A = n[i];
  		NMaterial.push_back(temp);
  	}
  
  	//create the material
  	Material newMaterial;
  	newMaterial.name = materialName;

  	complex<double> eta;
  	for(unsigned int i=0; i<KMaterial.size(); i++){

  		newMaterial.nu.push_back(KMaterial[i].nu);

  		eta = complex<double>(NMaterial[i].A, KMaterial[i].A);

  		newMaterial.eta.push_back(eta);
  	}
  
  	MaterialList.push_back(newMaterial);
  }
  
  void FitDisplay(){

  	double minA = 99999.0;
  	double maxA = -99999.0;
  	double k, n;

  
  	if(dispSimSpec)

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

  		{
  			if(SimSpectrum[i].A < minA)
  				minA = SimSpectrum[i].A;
  			if(SimSpectrum[i].A > maxA)
  				maxA = SimSpectrum[i].A;
  		}
  
  	if(dispRefSpec && RefSpectrum.size() > 0)

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

  		{
  			if(RefSpectrum[currentSpec][i].A < minA)
  				minA = RefSpectrum[currentSpec][i].A;
  			if(RefSpectrum[currentSpec][i].A > maxA)
  				maxA = RefSpectrum[currentSpec][i].A;
  		}
  	if(dispMatK)

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

  		{
  			k = MaterialList[currentMaterial].eta[i].imag() * dispScaleK;
  			if(k < minA)
  				minA = k;
  			if(k > maxA)
  				maxA = k;
  		}
  	if(dispSimK)

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

  		{
  			k = EtaK[i].A * dispScaleK;
  			if(k < minA)
  				minA = k;
  			if(EtaK[i].A > maxA)
  				maxA = k;
  		}
  	if(dispMatN)

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

  		{
  			n = (MaterialList[currentMaterial].eta[i].real() - baseIR) * dispScaleN;
  			if(n < minA)
  				minA = n;
  			if(n > maxA)
  				maxA = n;
  		}
  	if(dispSimN)

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

  		{
  			n = (EtaN[i].A - baseIR) * dispScaleN;
  			if(n < minA)
  				minA = n;
  			if(n > maxA)
  				maxA = n;
  		}
  
  	aMin = minA;
  	aMax = maxA;
  	UpdateDisplay();
  }
  
  void ChangeAbsorbance(){
  
  	//compute the real part of the index of refraction
  	
  	//copy the absorbance values into a linear array
  	int nSamples = MaterialList[currentMaterial].eta.size();

  	double startNu = MaterialList[currentMaterial].nu.front();
  	double endNu = MaterialList[currentMaterial].nu.back();
  	double* k = (double*)malloc(sizeof(double) * nSamples);
  	double* n = (double*)malloc(sizeof(double) * nSamples);

  	for(int i=0; i<nSamples; i++)
  		k[i] = MaterialList[currentMaterial].eta[i].imag() * cA;
  
  	//NMaterial.clear();
  	EtaK.clear();
  	EtaN.clear();
  	//use Kramers Kronig to determine the real part of the index of refraction
  	cudaKramersKronig(n, k, nSamples, startNu, endNu, baseIR);
  
  	//copy the real part of the index of refraction into the vector
  	SpecPair temp;
  
  	//load the imaginary IR from the absorbance data

  	double nu;

  	for(int i=0; i<nSamples; i++){
  		nu = MaterialList[currentMaterial].nu[i];
  		if(nu >= nuMin && nu <= nuMax){
  			temp.nu = nu;
  			temp.A = k[i];
  			EtaK.push_back(temp);
  			//temp.A = NMaterial[i].A;
  			temp.A = n[i];
  			EtaN.push_back(temp);
  		}
  	}
  
  	free(k);
  	free(n);
  }
  
  void SetMaterial()
  {
  	EtaK.clear();
  	EtaN.clear();
  
  	int nSamples = MaterialList[currentMaterial].eta.size();

  	double nu;

  	SpecPair temp;

  
  	//initialize the current nuMin and nuMax values
  	nuMin = MaterialList[currentMaterial].nu[0];
  	nuMax = nuMin;

  	for(int i=0; i<nSamples; i++){
  		nu = MaterialList[currentMaterial].nu[i];

  		//if(nu >= nuMin && nu <= nuMax){
  
  		//update the min and max values for display
  		if(nu < nuMin) nuMin = nu;
  		if(nu > nuMax) nuMax = nu;
  

  			temp.nu = nu;
  			temp.A = MaterialList[currentMaterial].eta[i].imag();
  			EtaK.push_back(temp);
  			temp.A = MaterialList[currentMaterial].eta[i].real();
  			EtaN.push_back(temp);

  	}
  	cA = 1.0;
  
  }
  
  
  
  int main(int argc, char *argv[])
  {
  	
  
  	//load the default project file (any previous optical settings)
  	LoadState();
  	
  	//load the default materials
  	LoadMaterial("eta_TolueneK.txt", "eta_TolueneN.txt", "Toluene");
  	LoadMaterial("kPMMA.txt", "PMMA");
  	//LoadMaterial("../../../../data/materials/rtsSU8_k.txt", "../../../../data/materials/rtsSU8_n.txt", "SU8");
  	SetMaterial();
  
  	//compute the analytical solution for the Mie scattered spectrum
  	SimulateSpectrum();
  
  	QApplication a(argc, argv);
  	InteractiveMie w;
  	
  
  	w.show();
  	
  	
  	w.move(0, 0);
  	QRect frame = w.frameGeometry();
  	QRect inside = w.geometry();
  
  	//activate a console for output
  	RedirectIOToConsole(0, frame.height(), frame.width());
  
  	gpSpectrumDisplay = new qtSpectrumDisplay();
  	gpSpectrumDisplay->move(frame.width(), 0);
  	gpSpectrumDisplay->resize(2*inside.height(), inside.height());
  
  	gpSpectrumDisplay->show();
  
  	//refresh the UI
  	w.refreshUI();
  
  	return a.exec();
  }