optics / imie

#include <fstream>
#include <iostream>
using namespace std;
#include "interactivemie.h"
#include <QtGui/QApplication>
#include <QGraphicsScene>
#include <QGraphicsView>
#include <QGraphicsPixmapItem>
#include <QLayout>
#include "qtSpectrumDisplay.h"
//#include "qwtSpectrumDisplay.h"
#include "globals.h"
#include "rtsGUIConsole.h"
#include "PerformanceData.h"
#include <complex>
#include <sstream>
//#include <direct.h>


PerformanceData PD;

qtSpectrumDisplay* gpSpectrumDisplay;
//qwtSpectrumDisplay* SpectrumDisplay;

QGraphicsScene* distortionScene = NULL;
QGraphicsView* distortionWindow = NULL;
QGraphicsPixmapItem* pixmapItem = NULL;

vector<vector<SpecPair> > RefSpectrum;
vector<SpecPair> SimSpectrum;
vector<SpecPair> EtaK;
vector<SpecPair> EtaN;

//source spectrum
bool useSourceSpectrum = false;
vector<SpecPair> SourceSpectrum;

//resample the source based on the material samples (EtaK and EtaN)
vector<SpecPair> SourceResampled;

int currentSpec = 0;

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

double aMin = 0;
double aMax = 1;

double nMag = 1.0;
double kMax = 1.0;

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 = -1;

//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();
	//SpectrumDisplay->replot();
}

void LoadMaterial(string fileName, string materialName)
{
	//open the file
	ifstream inFile(fileName.c_str());
	if(!inFile)
	{
		cout<<"Error loading material: "<<fileName<<endl;
		return;
	}

	Material newMaterial;
	newMaterial.validN = false;
	newMaterial.validK = false;
	//read the header
	string units;
	getline(inFile, units, '\t');
	char c = 0;
	while(c != '\n')
	{
		inFile.get(c);
		if(c == 'n') newMaterial.validN = true;
		if(c == 'k') newMaterial.validK = true;
	}

	//read the entire refractive index (both real and imaginary)
	float nu;
	float n;
	float k;

	while(inFile>>nu)
	{
		n = 1.0;
		k = 0.0;
		if(newMaterial.validN)
			inFile>>n;
		if(newMaterial.validK)
			inFile>>k;
		//ignore the rest of the line
		inFile.ignore();

		newMaterial.nu.push_back(nu);
		newMaterial.eta.push_back(complex<double>(n, k));
	}

	//iterate through the material to compute the necessary bounds for display
	float minN = 9999;
	float maxN = -9999;
	float sumN = 0;
	float maxK = 0;
	for(int i = 0; i<newMaterial.nu.size(); i++)
	{
		n = newMaterial.eta[i].real();
		k = newMaterial.eta[i].imag();

		if(n < minN) minN = n;
		if(n > maxN) maxN = n;
		if(k > maxK) maxK = k;
		sumN += n;
	}
	float meanN = sumN / newMaterial.nu.size();

	nMag = max(maxN - meanN, meanN - minN);
	kMax = maxK;
	baseIR = meanN;

	//set the name of the material object
	newMaterial.name = materialName;

	//add it to the material list
	MaterialList.push_back(newMaterial);
	currentMaterial = MaterialList.size() - 1;
}

void LoadSource(string fileNameSource)
{
	SourceSpectrum = LoadSpectrum(fileNameSource);
}

void ResampleSource()
{
	//clear the current resampled spectrum
	SourceResampled.clear();

	//get the number of source and material samples
	int nMatSamples = EtaK.size();
	int nSourceSamples = SourceSpectrum.size();

	float nu, I;
	for(int i=0; i<nMatSamples; i++)
	{
		SpecPair newSample;
		newSample.nu = EtaK[i].nu;

		//iterate through the SourceSpectrum to find the bounding wavelengths
		if(newSample.nu < SourceSpectrum[0].nu)
			newSample.A = 0.0;
		else if(newSample.nu > SourceSpectrum[nSourceSamples-1].nu)
			newSample.A = 0.0;
		else
		{
			int k=1;
			while(SourceSpectrum[k].nu < newSample.nu)
				k++;

			//interpolate
			float a = (newSample.nu - SourceSpectrum[k-1].nu)/(SourceSpectrum[k].nu - SourceSpectrum[k-1].nu);
			newSample.A = a * SourceSpectrum[k].A + (1.0 - a) * SourceSpectrum[k-1].A;
		}

		//insert the new spectral point into the resampled spectrum
		SourceResampled.push_back(newSample);
		//cout<<newSample.nu<<"   "<<newSample.A<<endl;
	}

}

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();
	if(currentMaterial == -1) return;

	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;

	//resample the source spectrum
	if(SourceSpectrum.size() != 0)
		ResampleSource();

}



int main(int argc, char *argv[])
{


	//load the default project file (any previous optical settings)
	LoadState();


	//load the default materials
	LoadMaterial("etaToluene.txt", "Toluene");
	LoadMaterial("kPMMA.txt", "PMMA");
	LoadMaterial("kPolyethylene.txt", "Polyethylene");
	LoadMaterial("kPTFE.txt", "Teflon");


	//LoadMaterial("eta_TolueneK.txt", "eta_TolueneN.txt", "Toluene");
	//LoadMaterial("kPMMA.txt", "PMMA");
	//LoadMaterial("eta_polystyreneK.txt", "Polystyrene");
	//LoadMaterial("../../../../data/materials/rtsSU8_k.txt", "../../../../data/materials/rtsSU8_n.txt", "SU8");
	SetMaterial();


	//load a mid-infrared source
	LoadSource("source_midIR.txt");

	ResampleSource();

	//compute the analytical solution for the Mie scattered spectrum
	SimulateSpectrum();

	QApplication a(argc, argv);

	//SpectrumDisplay = new qwtSpectrumDisplay();

	InteractiveMie w;
	w.show() ;
	w.move(0, 0);
	QSize uiFrame = w.frameSize() + QSize(10, 10);
	cout<<"Frame: "<<uiFrame.width()<<endl;
	QSize uiNoFrame = w.size();
	cout<<"No Frame: "<<uiNoFrame.width()<<endl;
	int frameHeight = uiFrame.height() - uiNoFrame.height();

	//activate a console for output
	RedirectIOToConsole(0, uiFrame.height(), uiFrame.width(), 400);
	printf("Frame height: %d\n", frameHeight);

	//set the size and position of the spectrum window
	int visWinSize = uiFrame.height()/2 - frameHeight;
	//create the far field window
	gpSpectrumDisplay = new qtSpectrumDisplay();
	gpSpectrumDisplay->resize(visWinSize*2, visWinSize);
	gpSpectrumDisplay->move(uiFrame.width(), 0);
	gpSpectrumDisplay->show();



	//distortion dialog box
	distortionDialog = new qtDistortionDialog();
	distortionDialog->move(0, 0);

	//display the distortion map
	distortionScene = new QGraphicsScene();
	distortionWindow = new QGraphicsView(distortionScene);
	distortionWindow->move(uiFrame.width(), visWinSize);

	//refresh the UI
	w.refreshUI();

	return a.exec();
}