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