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;

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

float aMin = 0;
float aMax = 1;

float scaleI0 = 1.0;
float refSlope = 0.0;

bool dispRefSpec = true;
bool dispSimSpec = true;
bool dispSimK = true;
bool dispMatK = true;
bool dispSimN = true;
bool dispMatN = true;
float dispScaleK = 1.0;
float dispScaleN = 1.0;
SpecType dispSimType = AbsorbanceSpecType;
bool dispNormalize = false;
float dispNormFactor = 1.0;


//material parameters
float radius = 4.0f;
float baseIR = 1.49f;
float cA = 1.0;
//vector<SpecPair> KMaterial;
//vector<SpecPair> NMaterial;
bool applyMaterial = true;
vector<Material> MaterialList;
int currentMaterial = 0;

//optical parameters
float cNAi = 0.0;
float cNAo = 0.6;
float oNAi = 0.0;
float oNAo = 0.6;
OpticsType opticsMode = TransmissionOpticsType;
bool pointDetector = false;
int objectiveSamples = 200;

//fitting parameters
float minMSE = 0.00001;
int maxFitIter = 20;

void TempSimSpectrum()
{
	SpecPair temp;
	for(int i=800; i<4000; i++)
	{
		temp.nu = i;
		temp.A = sin((float)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<float> eta;
	int j;
	for(int i=0; i<KMaterial.size(); i++){
		newMaterial.nu.push_back(KMaterial[i].nu);
		eta = complex<float>(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
	float* k = (float*)malloc(sizeof(float) * KMaterial.size());
	float* n = (float*)malloc(sizeof(float) * KMaterial.size());
	for(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(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<float> eta;
	int j;
	for(int i=0; i<KMaterial.size(); i++){
		newMaterial.nu.push_back(KMaterial[i].nu);
		eta = complex<float>(NMaterial[i].A, KMaterial[i].A);
		newMaterial.eta.push_back(eta);
	}

	MaterialList.push_back(newMaterial);
}

void FitDisplay(){
	float minA = 99999.0;
	float maxA = -99999.0;
	float k, n;

	if(dispSimSpec)
		for(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(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(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(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(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(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();
	float startNu = MaterialList[currentMaterial].nu.front();
	float endNu = MaterialList[currentMaterial].nu.back();
	float* k = (float*)malloc(sizeof(float) * nSamples);
	float* n = (float*)malloc(sizeof(float) * 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
	float 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();
	float nu;
	SpecPair temp;
	for(int i=0; i<nSamples; i++){
		nu = MaterialList[currentMaterial].nu[i];
		if(nu >= nuMin && nu <= nuMax){
			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();
}