Blame view

EstimateMaterial.cpp 2.77 KB
da3d4e0e   dmayerich   Initial commit.
1
2
3
  #include "globals.h"
  #define PI 3.14159
  
52a5fe9d   dmayerich   Added double supp...
4
  double CalculateError(double* E)
da3d4e0e   dmayerich   Initial commit.
5
6
  {
  	//Calculate the error between the Reference Spectrum and the Simulated Spectrum
52a5fe9d   dmayerich   Added double supp...
7
  	double sumE = 0.0;
da3d4e0e   dmayerich   Initial commit.
8
  	int nVals = RefSpectrum[currentSpec].size();
52a5fe9d   dmayerich   Added double supp...
9
  	double nu;
da3d4e0e   dmayerich   Initial commit.
10
11
12
13
14
15
16
17
18
19
  	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;
  }
  
52a5fe9d   dmayerich   Added double supp...
20
  void EstimateK(double* E)
da3d4e0e   dmayerich   Initial commit.
21
22
  {
  	int nVals = RefSpectrum[currentSpec].size();
52a5fe9d   dmayerich   Added double supp...
23
24
  	double nuStart = RefSpectrum[currentSpec].front().nu;
  	double nuEnd = RefSpectrum[currentSpec].back().nu;
da3d4e0e   dmayerich   Initial commit.
25
  
52a5fe9d   dmayerich   Added double supp...
26
27
28
29
  	double r = radius/10000.0;
  	double nu;
  	double dNu = (nuEnd - nuStart)/(nVals-1);
  	double eScale;
da3d4e0e   dmayerich   Initial commit.
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
  	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;
52a5fe9d   dmayerich   Added double supp...
53
  	for(int s=0; s<(int)RefSpectrum[currentSpec].size(); s++)
da3d4e0e   dmayerich   Initial commit.
54
55
56
57
58
59
60
61
62
63
64
65
  	{
  		//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
52a5fe9d   dmayerich   Added double supp...
66
  	double* E = (double*)malloc(sizeof(double) * RefSpectrum[currentSpec].size());
da3d4e0e   dmayerich   Initial commit.
67
  	//copy the absorbance values into a linear array
52a5fe9d   dmayerich   Added double supp...
68
69
  	double* k = (double*)malloc(sizeof(double) * EtaK.size());
  	double* n = (double*)malloc(sizeof(double) * EtaN.size());
da3d4e0e   dmayerich   Initial commit.
70
71
  
  	//iterate to solve for both n and k
52a5fe9d   dmayerich   Added double supp...
72
73
  	double sumE = 99999.9;
  	int j=0;
da3d4e0e   dmayerich   Initial commit.
74
75
  	//clear the console
  	system("cls");
52a5fe9d   dmayerich   Added double supp...
76
  	while(sumE > minMSE && j < maxFitIter)
da3d4e0e   dmayerich   Initial commit.
77
78
79
80
81
82
83
84
85
86
87
88
  	{	
  		//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
  		
52a5fe9d   dmayerich   Added double supp...
89
  		for(unsigned int i=0; i<EtaK.size(); i++)
da3d4e0e   dmayerich   Initial commit.
90
91
92
93
94
  			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();
52a5fe9d   dmayerich   Added double supp...
95
  		for(int i=0; i<(int)EtaK.size(); i++)
da3d4e0e   dmayerich   Initial commit.
96
97
98
99
100
101
102
  		{
  			temp.nu =  EtaK[i].nu;
  			temp.A = n[i];
  			EtaN.push_back(temp);
  		}
  
  		cout<<"    E = "<<sumE<<endl;
52a5fe9d   dmayerich   Added double supp...
103
  		j++;
da3d4e0e   dmayerich   Initial commit.
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
  		//SaveSpectrum(n, nVals, "simNj.txt");
  		//SaveSpectrum(k, nVals, "simKj.txt");
  		//SaveSpectrum(simSpec, nVals, "simSpec.txt");
  		//exit(1);
  		
  	}
  
  	free(E);
  	free(k);
  	free(n);
  
  
  
  
  }