optics / imie

Commit 8ffb8373b71fd5d6e8327b372b67c50fb8df43f9

Authored by dmayerich
1 parent 4198e3be

Improved materials and the Kramers-Kronig solver.

Showing 28 changed files with 7997 additions and 7851 deletions   Show diff stats
Inline Side-by-side
CMakeLists.txt
  Show/Hide comments   View file @8ffb837
... ... @@ -25,6 +25,9 @@ find_package(GLUT REQUIRED)
25 25 #find GLEW
26 26 find_package(GLEW REQUIRED)
27 27  
  28 +#find Qwt
  29 +find_package(Qwt REQUIRED)
  30 +
28 31 #add Qt OpenGL stuff
29 32 set(QT_USE_QTOPENGL TRUE)
30 33  
... ... @@ -53,9 +56,10 @@ file(GLOB SRC_CU "*.cu")
53 56 #set up copying data files
54 57 configure_file(kPMMA.txt ${CMAKE_CURRENT_BINARY_DIR}/kPMMA.txt @ONLY)
55 58 configure_file(eta_polystyreneK.txt ${CMAKE_CURRENT_BINARY_DIR}/eta_polystyreneK.txt @ONLY)
56   -configure_file(eta_TolueneK.txt ${CMAKE_CURRENT_BINARY_DIR}/eta_TolueneK.txt @ONLY)
57   -configure_file(eta_TolueneN.txt ${CMAKE_CURRENT_BINARY_DIR}/eta_TolueneN.txt @ONLY)
  59 +configure_file(etaToluene.txt ${CMAKE_CURRENT_BINARY_DIR}/etaToluene.txt @ONLY)
58 60 configure_file(source_midIR.txt ${CMAKE_CURRENT_BINARY_DIR}/source_midIR.txt @ONLY)
  61 +configure_file(kPolyethylene.txt ${CMAKE_CURRENT_BINARY_DIR}/kPolyethylene.txt @ONLY)
  62 +configure_file(kPTFE.txt ${CMAKE_CURRENT_BINARY_DIR}/kPTFE.txt @ONLY)
59 63  
60 64 #determine which source files have to be moc'd
61 65 Qt4_wrap_cpp(UI_MOC ${SRC_H})
... ... @@ -72,7 +76,4 @@ source_group(QtUI FILES ${SRC_UI})
72 76 cuda_add_executable(IMie ${SRC_CPP} ${SRC_H} ${UI_H} ${UI_MOC} ${ALL_RCC} ${SRC_CU})
73 77  
74 78 #set the link libraries
75   -target_link_libraries(IMie ${QT_LIBRARIES} ${QT_QTOPENGL_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ${GLEW_LIBRARY})
76   -
77   -
78   -
  79 +target_link_libraries(IMie ${QT_LIBRARIES} ${QT_QTOPENGL_LIBRARY} ${OPENGL_gl_LIBRARY} ${OPENGL_glu_LIBRARY} ${GLEW_LIBRARY} ${QWT_LIBRARY})
79 80 \ No newline at end of file
... ...
EstimateMaterial.cpp
  Show/Hide comments   View file @8ffb837
1   -#include "globals.h"
2   -#include <stdlib.h>
3   -#define PI 3.14159
4   -
5   -double CalculateError(double* E)
6   -{
7   - //Calculate the error between the Reference Spectrum and the Simulated Spectrum
8   - double sumE = 0.0;
9   - int nVals = RefSpectrum[currentSpec].size();
10   - double nu;
11   - for(int i=0; i<nVals; i++)
12   - {
13   - nu = RefSpectrum[currentSpec][i].nu;
14   - E[i] = RefSpectrum[currentSpec][i].A + nu * refSlope - SimSpectrum[i].A;
15   - sumE += E[i]*E[i];
16   - }
17   -
18   - return sumE/nVals;
19   -}
20   -
21   -void EstimateK(double* E)
22   -{
23   - int nVals = RefSpectrum[currentSpec].size();
24   - double nuStart = RefSpectrum[currentSpec].front().nu;
25   - double nuEnd = RefSpectrum[currentSpec].back().nu;
26   -
27   - double r = radius/10000.0;
28   - double nu;
29   - double dNu = (nuEnd - nuStart)/(nVals-1);
30   - double eScale;
31   - for(int i=0; i<nVals; i++)
32   - {
33   - nu = nuStart + i*2;
34   -
35   - eScale = 1/(8*PI*r*nu);
36   - EtaK[i].A = EtaK[i].A + eScale * E[i];
37   - if(EtaK[i].A < 0.0) EtaK[i].A = 0.0;
38   - }
39   -}
40   -
41   -void EstimateMaterial()
42   -{
43   - /*This function estimates the material properties of a sphere based on the
44   - input spectrum RefSpectrum and the optical properties of the system.
45   - 1) The material properties are stored in EtaK and EtaN
46   - 2) The best fit is stored in SimSpectrum*/
47   -
48   - //initialize the material index of refraction
49   - EtaN.clear();
50   - EtaK.clear();
51   -
52   - //insert the default material properties
53   - SpecPair temp;
54   - for(int s=0; s<(int)RefSpectrum[currentSpec].size(); s++)
55   - {
56   - //the real part of the IR is the user-specified baseline IR
57   - temp.nu = RefSpectrum[currentSpec][s].nu;
58   - temp.A = baseIR;
59   - EtaN.push_back(temp);
60   - //the imaginary part of the IR is zero absorbance
61   - temp.A = 0.0f;
62   - EtaK.push_back(temp);
63   - }
64   -
65   -
66   - //allocate space to store the list of error values
67   - double* E = (double*)malloc(sizeof(double) * RefSpectrum[currentSpec].size());
68   - //copy the absorbance values into a linear array
69   - double* k = (double*)malloc(sizeof(double) * EtaK.size());
70   - double* n = (double*)malloc(sizeof(double) * EtaN.size());
71   -
72   - //iterate to solve for both n and k
73   - double sumE = 99999.9;
74   - int j=0;
75   - //clear the console
76   - system("cls");
77   - while(sumE > minMSE && j < maxFitIter)
78   - {
79   - //simulate a spectrum based on the current IR
80   - SimulateSpectrum();
81   -
82   - //calculate the error term
83   - sumE = CalculateError(E);
84   -
85   - //estimate the new absorbance
86   - EstimateK(E);
87   -
88   - //use Kramers-Kronig to compute n
89   -
90   - for(unsigned int i=0; i<EtaK.size(); i++)
91   - k[i] = EtaK[i].A;
92   - cudaKramersKronig(n, k, EtaK.size(), EtaK.front().nu, EtaK.back().nu, baseIR);
93   -
94   - //copy the real part of the index of refraction into the vector
95   - EtaN.clear();
96   - for(int i=0; i<(int)EtaK.size(); i++)
97   - {
98   - temp.nu = EtaK[i].nu;
99   - temp.A = n[i];
100   - EtaN.push_back(temp);
101   - }
102   -
103   - cout<<" E = "<<sumE<<endl;
104   - j++;
105   - //SaveSpectrum(n, nVals, "simNj.txt");
106   - //SaveSpectrum(k, nVals, "simKj.txt");
107   - //SaveSpectrum(simSpec, nVals, "simSpec.txt");
108   - //exit(1);
109   -
110   - }
111   -
112   - free(E);
113   - free(k);
114   - free(n);
115   -
116   -
117   -
118   -
119   -}
120 1 \ No newline at end of file
  2 +#include "globals.h"
  3 +#include <stdlib.h>
  4 +#define PI 3.14159
  5 +
  6 +double CalculateError(double* E)
  7 +{
  8 + //Calculate the error between the Reference Spectrum and the Simulated Spectrum
  9 + double sumE = 0.0;
  10 + int nVals = RefSpectrum[currentSpec].size();
  11 + double nu;
  12 + for(int i=0; i<nVals; i++)
  13 + {
  14 + nu = RefSpectrum[currentSpec][i].nu;
  15 + E[i] = RefSpectrum[currentSpec][i].A + nu * refSlope - SimSpectrum[i].A;
  16 + sumE += E[i]*E[i];
  17 + }
  18 +
  19 + return sumE/nVals;
  20 +}
  21 +
  22 +void EstimateK(double* E)
  23 +{
  24 + int nVals = RefSpectrum[currentSpec].size();
  25 + double nuStart = RefSpectrum[currentSpec].front().nu;
  26 + double nuEnd = RefSpectrum[currentSpec].back().nu;
  27 +
  28 + double r = radius/10000.0;
  29 + double nu;
  30 + double dNu = (nuEnd - nuStart)/(nVals-1);
  31 + double eScale;
  32 + for(int i=0; i<nVals; i++)
  33 + {
  34 + nu = nuStart + i*2;
  35 +
  36 + eScale = 1/(8*PI*r*nu);
  37 + EtaK[i].A = EtaK[i].A + eScale * E[i];
  38 + if(EtaK[i].A < 0.0) EtaK[i].A = 0.0;
  39 + }
  40 +}
  41 +
  42 +void EstimateMaterial()
  43 +{
  44 + /*This function estimates the material properties of a sphere based on the
  45 + input spectrum RefSpectrum and the optical properties of the system.
  46 + 1) The material properties are stored in EtaK and EtaN
  47 + 2) The best fit is stored in SimSpectrum*/
  48 +
  49 + //initialize the material index of refraction
  50 + EtaN.clear();
  51 + EtaK.clear();
  52 +
  53 + //insert the default material properties
  54 + SpecPair temp;
  55 + for(int s=0; s<(int)RefSpectrum[currentSpec].size(); s++)
  56 + {
  57 + //the real part of the IR is the user-specified baseline IR
  58 + temp.nu = RefSpectrum[currentSpec][s].nu;
  59 + temp.A = baseIR;
  60 + EtaN.push_back(temp);
  61 + //the imaginary part of the IR is zero absorbance
  62 + temp.A = 0.0f;
  63 + EtaK.push_back(temp);
  64 + }
  65 +
  66 +
  67 + //allocate space to store the list of error values
  68 + double* E = (double*)malloc(sizeof(double) * RefSpectrum[currentSpec].size());
  69 + //copy the absorbance values into a linear array
  70 + double* k = (double*)malloc(sizeof(double) * EtaK.size());
  71 + double* n = (double*)malloc(sizeof(double) * EtaN.size());
  72 +
  73 + //iterate to solve for both n and k
  74 + double sumE = 99999.9;
  75 + int j=0;
  76 + //clear the console
  77 + system("cls");
  78 + while(sumE > minMSE && j < maxFitIter)
  79 + {
  80 + //simulate a spectrum based on the current IR
  81 + SimulateSpectrum();
  82 +
  83 + //calculate the error term
  84 + sumE = CalculateError(E);
  85 +
  86 + //estimate the new absorbance
  87 + EstimateK(E);
  88 +
  89 + //use Kramers-Kronig to compute n
  90 +
  91 + for(unsigned int i=0; i<EtaK.size(); i++)
  92 + k[i] = EtaK[i].A;
  93 + cudaKramersKronig(n, k, EtaK.size(), EtaK.front().nu, EtaK.back().nu, baseIR);
  94 +
  95 + //copy the real part of the index of refraction into the vector
  96 + EtaN.clear();
  97 + for(int i=0; i<(int)EtaK.size(); i++)
  98 + {
  99 + temp.nu = EtaK[i].nu;
  100 + temp.A = n[i];
  101 + EtaN.push_back(temp);
  102 + }
  103 +
  104 + cout<<" E = "<<sumE<<endl;
  105 + j++;
  106 + //SaveSpectrum(n, nVals, "simNj.txt");
  107 + //SaveSpectrum(k, nVals, "simKj.txt");
  108 + //SaveSpectrum(simSpec, nVals, "simSpec.txt");
  109 + //exit(1);
  110 +
  111 + }
  112 +
  113 + free(E);
  114 + free(k);
  115 + free(n);
  116 +
  117 +
  118 +
  119 +
  120 +}
... ...
FileIO.cpp
  Show/Hide comments   View file @8ffb837
... ... @@ -6,191 +6,198 @@ using namespace std;
6 6  
7 7 vector<SpecPair> LoadSpectrum(string filename)
8 8 {
9   - //load a spectrum from a file and resample to 2wn intervals
10   -
11   - //create the spectrum
12   - vector<SpecPair> S;
13   -
14   - //open the file
15   - ifstream inFile(filename.c_str());
16   -
17   - //read all elements of the file
18   - SpecPair temp;
19   - while(!inFile.eof()){
20   - inFile>>temp.nu;
21   - inFile>>temp.A;
22   - S.push_back(temp);
23   - }
24   -
25   - //compute the minimum and maximum input wavenumbers
26   - double inMin = S.front().nu;
27   - double inMax = S.back().nu;
28   -
29   - int nuMin = (int)ceil(inMin);
30   - int nuMax = (int)floor(inMax);
31   -
32   - //make sure both are either even or odd
33   - if(nuMin % 2 != nuMax % 2)
34   - nuMax--;
35   -
36   - //compute the number of values in the resampled spectrum
37   - int nVals = (nuMax - nuMin)/2 + 1;
38   -
39   - //allocate space for the spectrum
40   - vector<SpecPair> outSpec;
41   -
42   - double nu, highVal, lowVal, a;
43   - int j=1;
44   - for(int i=0; i<nVals; i++)
45   - {
46   - nu = nuMin + i * 2;
47   - temp.nu = nu;
48   -
49   - //handle the boundary cases
50   - if(nu < inMin || nu > inMax)
51   - temp.A = 0.0;
52   - else
53   - {
54   - //move to the correct position in the input array
55   - while(j < (int)S.size()-1 && S[j].nu <= nu)
56   - j++;
57   -
58   -
59   - lowVal = S[j-1].nu;
60   - highVal = S[j].nu;
61   - a = (nu - lowVal)/(highVal - lowVal);
62   - temp.A = S[j-1].A * (1.0 - a) + S[j].A * a;
63   - }
64   - outSpec.push_back(temp);
65   - }
66   -
67   -
68   -
69   - return outSpec;
  9 + //load a spectrum from a file and resample to 2wn intervals
  10 +
  11 + //create the spectrum
  12 + vector<SpecPair> S;
  13 +
  14 + //open the file
  15 + ifstream inFile(filename.c_str());
  16 + if(!inFile)
  17 + {
  18 + cout<<"Error loading spectrum: "<<inFile<<endl;
  19 + return S;
  20 + }
  21 +
  22 + //read all elements of the file
  23 + SpecPair temp;
  24 + while(!inFile.eof()) {
  25 + inFile>>temp.nu;
  26 + inFile>>temp.A;
  27 + S.push_back(temp);
  28 + }
  29 +
  30 + //compute the minimum and maximum input wavenumbers
  31 + double inMin = S.front().nu;
  32 + double inMax = S.back().nu;
  33 +
  34 + int nuMin = (int)ceil(inMin);
  35 + int nuMax = (int)floor(inMax);
  36 +
  37 + //make sure both are either even or odd
  38 + if(nuMin % 2 != nuMax % 2)
  39 + nuMax--;
  40 +
  41 + //compute the number of values in the resampled spectrum
  42 + int nVals = (nuMax - nuMin)/2 + 1;
  43 +
  44 + //allocate space for the spectrum
  45 + vector<SpecPair> outSpec;
  46 +
  47 + double nu, highVal, lowVal, a;
  48 + int j=1;
  49 + for(int i=0; i<nVals; i++)
  50 + {
  51 + nu = nuMin + i * 2;
  52 + temp.nu = nu;
  53 +
  54 + //handle the boundary cases
  55 + if(nu < inMin || nu > inMax)
  56 + temp.A = 0.0;
  57 + else
  58 + {
  59 + //move to the correct position in the input array
  60 + while(j < (int)S.size()-1 && S[j].nu <= nu)
  61 + j++;
  62 +
  63 +
  64 + lowVal = S[j-1].nu;
  65 + highVal = S[j].nu;
  66 + a = (nu - lowVal)/(highVal - lowVal);
  67 + temp.A = S[j-1].A * (1.0 - a) + S[j].A * a;
  68 + }
  69 + outSpec.push_back(temp);
  70 + }
  71 +
  72 +
  73 +
  74 + return outSpec;
70 75 }
71 76  
72 77 vector<SpecPair> SetReferenceSpectrum(char* text)
73 78 {
74   - stringstream inString(text);
  79 + stringstream inString(text);
75 80  
76   - //create the spectrum
77   - vector<SpecPair> S;
  81 + //create the spectrum
  82 + vector<SpecPair> S;
78 83  
79   - SpecPair temp;
80   - while(!inString.eof()){
81   - inString>>temp.nu;
82   - inString>>temp.A;
83   - S.push_back(temp);
84   - }
  84 + SpecPair temp;
  85 + while(!inString.eof()) {
  86 + inString>>temp.nu;
  87 + inString>>temp.A;
  88 + S.push_back(temp);
  89 + }
85 90  
86   - return S;
  91 + return S;
87 92 }
88 93  
89   -void SaveK(string fileName)
  94 +/*void SaveK(string fileName)
90 95 {
91   - ofstream outFile(fileName.c_str());
92   - for(unsigned int i=0; i<EtaK.size(); i++)
93   - {
94   - outFile<<EtaK[i].nu<<" ";
95   - outFile<<EtaK[i].A<<endl;
96   - }
97   - outFile.close();
98   -}
99   -
100   -void SaveN(string fileName)
  96 + ofstream outFile(fileName.c_str());
  97 + for(unsigned int i=0; i<EtaK.size(); i++)
  98 + {
  99 + outFile<<EtaK[i].nu<<" ";
  100 + outFile<<EtaK[i].A<<endl;
  101 + }
  102 + outFile.close();
  103 +}*/
  104 +
  105 +void SaveMaterial(string fileName)
101 106 {
102   - ofstream outFile(fileName.c_str());
103   - for(unsigned int i=0; i<EtaN.size(); i++)
104   - {
105   - outFile<<EtaN[i].nu<<" ";
106   - outFile<<EtaN[i].A<<endl;
107   - }
108   - outFile.close();
  107 + ofstream outFile(fileName.c_str());
  108 + outFile<<"nu"<<'\t'<<"n"<<'\t'<<"k"<<endl;
  109 + for(unsigned int i=0; i<EtaN.size(); i++)
  110 + {
  111 + outFile<<EtaN[i].nu<<'\t';
  112 + outFile<<EtaN[i].A<<'\t';
  113 + outFile<<EtaK[i].A<<endl;
  114 + }
  115 + outFile.close();
109 116 }
110 117  
111 118 void SaveSimulation(string fileName)
112 119 {
113   - ofstream outFile(fileName.c_str());
114   - outFile.precision(10);
115   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
116   - {
117   - outFile<<SimSpectrum[i].nu<<" ";
118   - outFile<<SimSpectrum[i].A<<endl;
119   - }
120   - outFile.close();
  120 + ofstream outFile(fileName.c_str());
  121 + outFile.precision(10);
  122 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
  123 + {
  124 + outFile<<SimSpectrum[i].nu<<" ";
  125 + outFile<<SimSpectrum[i].A<<endl;
  126 + }
  127 + outFile.close();
121 128 }
122 129  
123 130 void SaveState()
124 131 {
125   - ofstream outFile("main.prj");
126   - //Window Parameters
127   - outFile<<nuMin<<endl;
128   - outFile<<nuMax<<endl;
129   - outFile<<aMin<<endl;
130   - outFile<<aMax<<endl;
131   - outFile<<dNu<<endl;
132   -
133   - //material parameters
134   - outFile<<radius<<endl;
135   - outFile<<baseIR<<endl;
136   - outFile<<cA<<endl;
137   -
138   - //optical parameters
139   - outFile<<cNAi<<endl;
140   - outFile<<cNAo<<endl;
141   - outFile<<oNAi<<endl;
142   - outFile<<oNAo<<endl;
143   -
144   - outFile.close();
  132 + ofstream outFile("main.prj");
  133 + //Window Parameters
  134 + outFile<<nuMin<<endl;
  135 + outFile<<nuMax<<endl;
  136 + outFile<<aMin<<endl;
  137 + outFile<<aMax<<endl;
  138 + outFile<<dNu<<endl;
  139 +
  140 + //material parameters
  141 + outFile<<radius<<endl;
  142 + outFile<<baseIR<<endl;
  143 + outFile<<cA<<endl;
  144 +
  145 + //optical parameters
  146 + outFile<<cNAi<<endl;
  147 + outFile<<cNAo<<endl;
  148 + outFile<<oNAi<<endl;
  149 + outFile<<oNAo<<endl;
  150 +
  151 + outFile.close();
145 152  
146 153 }
147 154  
148 155 void LoadState()
149 156 {
150   - ifstream inFile("main.prj");
151   - //Window Parameters
152   - inFile>>nuMin;
153   - inFile>>nuMax;
154   - inFile>>aMin;
155   - inFile>>aMax;
156   - inFile>>dNu;
157   -
158   - //material parameters
159   - inFile>>radius;
160   - inFile>>baseIR;
161   - inFile>>cA;
162   -
163   - //optical parameters
164   - inFile>>cNAi;
165   - inFile>>cNAo;
166   - inFile>>oNAi;
167   - inFile>>oNAo;
168   -
169   - inFile.close();
  157 + ifstream inFile("main.prj");
  158 + //Window Parameters
  159 + inFile>>nuMin;
  160 + inFile>>nuMax;
  161 + inFile>>aMin;
  162 + inFile>>aMax;
  163 + inFile>>dNu;
  164 +
  165 + //material parameters
  166 + inFile>>radius;
  167 + inFile>>baseIR;
  168 + inFile>>cA;
  169 +
  170 + //optical parameters
  171 + inFile>>cNAi;
  172 + inFile>>cNAo;
  173 + inFile>>oNAi;
  174 + inFile>>oNAo;
  175 +
  176 + inFile.close();
170 177  
171 178 }
172 179  
173 180 void SetDefaults()
174 181 {
175 182  
176   - nuMin = 800;
177   - nuMax = 4000;
178   - dNu = 2;
  183 + nuMin = 800;
  184 + nuMax = 4000;
  185 + dNu = 2;
179 186  
180   - aMin = 0;
181   - aMax = 1;
  187 + aMin = 0;
  188 + aMax = 1;
182 189  
183 190  
184   - //material parameters
185   - radius = 4.0f;
186   - baseIR = 1.49f;
187   - cA = 1.0;
188   - vector<SpecPair> KMaterial;
189   - vector<SpecPair> NMaterial;
  191 + //material parameters
  192 + radius = 4.0f;
  193 + baseIR = 1.49f;
  194 + cA = 1.0;
  195 + vector<SpecPair> KMaterial;
  196 + vector<SpecPair> NMaterial;
190 197  
191   - //optical parameters
192   - cNAi = 0.0;
193   - cNAo = 0.6;
194   - oNAi = 0.0;
195   - oNAo = 0.6;
196   -}
197 198 \ No newline at end of file
  199 + //optical parameters
  200 + cNAi = 0.0;
  201 + cNAo = 0.6;
  202 + oNAi = 0.0;
  203 + oNAo = 0.6;
  204 +}
... ...
GAMMA.cpp
  Show/Hide comments   View file @8ffb837
1   -// gamma.cpp -- computation of gamma function.
2   -// Algorithms and coefficient values from "Computation of Special
3   -// Functions", Zhang and Jin, John Wiley and Sons, 1996.
4   -//
5   -// (C) 2003, C. Bond. All rights reserved.
6   -//
7   -// Returns gamma function of argument 'x'.
8   -//
9   -// NOTE: Returns 1e308 if argument is a negative integer or 0,
10   -// or if argument exceeds 171.
11   -//
12   -#define _USE_MATH_DEFINES
13   -#include <math.h>
14   -double gamma(double x)
15   -{
16   - int i,k,m;
17   - double ga,gr,r,z;
18   -
19   - static double g[] = {
20   - 1.0,
21   - 0.5772156649015329,
22   - -0.6558780715202538,
23   - -0.420026350340952e-1,
24   - 0.1665386113822915,
25   - -0.421977345555443e-1,
26   - -0.9621971527877e-2,
27   - 0.7218943246663e-2,
28   - -0.11651675918591e-2,
29   - -0.2152416741149e-3,
30   - 0.1280502823882e-3,
31   - -0.201348547807e-4,
32   - -0.12504934821e-5,
33   - 0.1133027232e-5,
34   - -0.2056338417e-6,
35   - 0.6116095e-8,
36   - 0.50020075e-8,
37   - -0.11812746e-8,
38   - 0.1043427e-9,
39   - 0.77823e-11,
40   - -0.36968e-11,
41   - 0.51e-12,
42   - -0.206e-13,
43   - -0.54e-14,
44   - 0.14e-14};
45   -
46   - if (x > 171.0) return 1e308; // This value is an overflow flag.
47   - if (x == (int)x) {
48   - if (x > 0.0) {
49   - ga = 1.0; // use factorial
50   - for (i=2;i<x;i++) {
51   - ga *= i;
52   - }
53   - }
54   - else
55   - ga = 1e308;
56   - }
57   - else {
58   - if (fabs(x) > 1.0) {
59   - z = fabs(x);
60   - m = (int)z;
61   - r = 1.0;
62   - for (k=1;k<=m;k++) {
63   - r *= (z-k);
64   - }
65   - z -= m;
66   - }
67   - else
68   - z = x;
69   - gr = g[24];
70   - for (k=23;k>=0;k--) {
71   - gr = gr*z+g[k];
72   - }
73   - ga = 1.0/(gr*z);
74   - if (fabs(x) > 1.0) {
75   - ga *= r;
76   - if (x < 0.0) {
77   - ga = -M_PI/(x*ga*sin(M_PI*x));
78   - }
79   - }
80   - }
81   - return ga;
82   -}
  1 +// gamma.cpp -- computation of gamma function.
  2 +// Algorithms and coefficient values from "Computation of Special
  3 +// Functions", Zhang and Jin, John Wiley and Sons, 1996.
  4 +//
  5 +// (C) 2003, C. Bond. All rights reserved.
  6 +//
  7 +// Returns gamma function of argument 'x'.
  8 +//
  9 +// NOTE: Returns 1e308 if argument is a negative integer or 0,
  10 +// or if argument exceeds 171.
  11 +//
  12 +#define _USE_MATH_DEFINES
  13 +#include <math.h>
  14 +double gamma(double x)
  15 +{
  16 + int i,k,m;
  17 + double ga,gr,r,z;
  18 +
  19 + static double g[] = {
  20 + 1.0,
  21 + 0.5772156649015329,
  22 + -0.6558780715202538,
  23 + -0.420026350340952e-1,
  24 + 0.1665386113822915,
  25 + -0.421977345555443e-1,
  26 + -0.9621971527877e-2,
  27 + 0.7218943246663e-2,
  28 + -0.11651675918591e-2,
  29 + -0.2152416741149e-3,
  30 + 0.1280502823882e-3,
  31 + -0.201348547807e-4,
  32 + -0.12504934821e-5,
  33 + 0.1133027232e-5,
  34 + -0.2056338417e-6,
  35 + 0.6116095e-8,
  36 + 0.50020075e-8,
  37 + -0.11812746e-8,
  38 + 0.1043427e-9,
  39 + 0.77823e-11,
  40 + -0.36968e-11,
  41 + 0.51e-12,
  42 + -0.206e-13,
  43 + -0.54e-14,
  44 + 0.14e-14
  45 + };
  46 +
  47 + if (x > 171.0) return 1e308; // This value is an overflow flag.
  48 + if (x == (int)x) {
  49 + if (x > 0.0) {
  50 + ga = 1.0; // use factorial
  51 + for (i=2; i<x; i++) {
  52 + ga *= i;
  53 + }
  54 + }
  55 + else
  56 + ga = 1e308;
  57 + }
  58 + else {
  59 + if (fabs(x) > 1.0) {
  60 + z = fabs(x);
  61 + m = (int)z;
  62 + r = 1.0;
  63 + for (k=1; k<=m; k++) {
  64 + r *= (z-k);
  65 + }
  66 + z -= m;
  67 + }
  68 + else
  69 + z = x;
  70 + gr = g[24];
  71 + for (k=23; k>=0; k--) {
  72 + gr = gr*z+g[k];
  73 + }
  74 + ga = 1.0/(gr*z);
  75 + if (fabs(x) > 1.0) {
  76 + ga *= r;
  77 + if (x < 0.0) {
  78 + ga = -M_PI/(x*ga*sin(M_PI*x));
  79 + }
  80 + }
  81 + }
  82 + return ga;
  83 +}
... ...
PerformanceData.h
  Show/Hide comments   View file @8ffb837
1   -// add the following to a cpp file:
2   -// PerformanceData PD;
3   -
4   -
5   -#pragma once
6   -#include <ostream>
7   -using namespace std;
8   -
9   -enum PerformanceDataType
10   -{
11   - PD_DISPLAY=0,
12   - PD_SPS,
13   - PD_UNUSED0,
14   -
15   - //my stuff
16   - SIMULATE_SPECTRUM,
17   - SIMULATE_GPU,
18   - KRAMERS_KRONIG,
19   -
20   -
21   -
22   - //end my stuff
23   - PERFORMANCE_DATA_TYPE_COUNT
24   -};
25   -
26   -static char PDTypeNames[][255] = {
27   - "Display ",
28   - "Simulation Total ",
29   - " ----------------- ",
30   - //my stuff
31   - "Simulate Spectrum ",
32   - " GPU Portion ",
33   - "Kramers-Kronig ",
34   -
35   - //end my stuff
36   -
37   -};
38   -#ifdef WIN32
39   -#include <stdio.h>
40   -#include <windows.h>
41   -#include <float.h>
42   -
43   -#include <iostream>
44   -#include <iomanip>
45   -
46   -//-------------------------------------------------------------------------------
47   -
48   -class PerformanceData
49   -{
50   -public:
51   - PerformanceData() { ClearAll(); QueryPerformanceFrequency(&cps); }
52   - ~PerformanceData(){}
53   -
54   - void ClearAll()
55   - {
56   - for ( int i=0; i<PERFORMANCE_DATA_TYPE_COUNT; i++ ) {
57   - for ( int j=0; j<256; j++ ) times[i][j] = 0;
58   - pos[i] = 0;
59   - minTime[i] = 0xFFFFFFFF;
60   - maxTime[i] = 0;
61   - totalTime[i] = 0;
62   - dataReady[i] = false;
63   - }
64   - }
65   -
66   - void StartTimer( int type ) { QueryPerformanceCounter( &startTime[type] );}
67   - void EndTimer( int type ) {
68   - LARGE_INTEGER endTime;
69   - QueryPerformanceCounter( &endTime );
70   - double t = (double)(endTime.QuadPart - startTime[type].QuadPart);
71   - //unsigned int t = GetTickCount() - startTime[type];
72   - if ( t < minTime[type] ) minTime[type] = t;
73   - if ( t > maxTime[type] ) maxTime[type] = t;
74   - totalTime[type] -= times[type][ pos[type] ];
75   - times[type][ pos[type] ] = t;
76   - totalTime[type] += t;
77   - pos[type]++;
78   - if ( pos[type] == 0 ) dataReady[type] = true;
79   - }
80   -
81   - void PrintResult( ostream &os,int i=PERFORMANCE_DATA_TYPE_COUNT)
82   - {
83   - os.setf(ios::fixed);
84   - if ((i<PERFORMANCE_DATA_TYPE_COUNT)&&(i>=0)){
85   - double a = GetAvrgTime(i);
86   - if ( a )
87   - os<< PDTypeNames[i]<<" : avrg="<<setw(8)<<setprecision(3)<<a<<"\tmin="<<setw(8)<<setprecision(3)<< GetMinTime(i) <<"\tmax="<<setw(8)<<setprecision(3)<< GetMaxTime(i) <<endl ;
88   - else
89   - os<< PDTypeNames[i]<<" : avrg= -----\tmin= -----\tmax= -----"<<endl;
90   - }
91   - }
92   -
93   - void PrintResults( ostream &os)
94   - {
95   - for ( int i=0; i<PERFORMANCE_DATA_TYPE_COUNT; i++ )
96   - PrintResult(os,i);
97   - }
98   -
99   - double GetLastTime( int type ) { return times[type][pos[type]]; }
100   - double GetAvrgTime( int type ) { double a = 1000.0 * totalTime[type] / (float)cps.QuadPart / ( (dataReady[type]) ? 256.0 : (double)pos[type] ); return (_finite(a))? a:0; }
101   - double GetMinTime( int type ) { return 1000.0 * minTime[type] / (float)cps.LowPart; }
102   - double GetMaxTime( int type ) { return 1000.0 * maxTime[type] / (float)cps.LowPart; }
103   -
104   -private:
105   - double times[PERFORMANCE_DATA_TYPE_COUNT][256];
106   - unsigned char pos[PERFORMANCE_DATA_TYPE_COUNT];
107   - LARGE_INTEGER startTime[PERFORMANCE_DATA_TYPE_COUNT];
108   - double minTime[ PERFORMANCE_DATA_TYPE_COUNT ];
109   - double maxTime[ PERFORMANCE_DATA_TYPE_COUNT ];
110   - double totalTime[ PERFORMANCE_DATA_TYPE_COUNT ];
111   - bool dataReady[ PERFORMANCE_DATA_TYPE_COUNT ];
112   - LARGE_INTEGER cps;
113   -};
114   -
115   -//-------------------------------------------------------------------------------
116   -#else
117   -
118   -class PerformanceData{
119   -public:
120   - PerformanceData() {;};
121   - ~PerformanceData(){;};
122   - void ClearAll(){;};
123   - void StartTimer( int type ) {;};
124   - void EndTimer( int type ) {;};
125   - void PrintResults( ostream &os){;};
126   - void PrintResult( ostream &os, int i=PERFORMANCE_DATA_TYPE_COUNT){;};
127   - double GetLastTime( int type ) { return 0.0; };
128   - double GetAvrgTime( int type ) { return 0.0; };
129   - double GetMinTime( int type ) { return 0.0; };
130   - double GetMaxTime( int type ) { return 0.0; };
131   -};
132   -
133   -#endif
134   -//-------------------------------------------------------------------------------
135   -
136   -extern PerformanceData PD;
137   -
138   -//-------------------------------------------------------------------------------
  1 +// add the following to a cpp file:
  2 +// PerformanceData PD;
  3 +
  4 +
  5 +#pragma once
  6 +#include <ostream>
  7 +using namespace std;
  8 +
  9 +enum PerformanceDataType
  10 +{
  11 + PD_DISPLAY=0,
  12 + PD_SPS,
  13 + PD_UNUSED0,
  14 +
  15 + //my stuff
  16 + SIMULATE_SPECTRUM,
  17 + SIMULATE_GPU,
  18 + KRAMERS_KRONIG,
  19 +
  20 +
  21 +
  22 + //end my stuff
  23 + PERFORMANCE_DATA_TYPE_COUNT
  24 +};
  25 +
  26 +static char PDTypeNames[][255] = {
  27 + "Display ",
  28 + "Simulation Total ",
  29 + " ----------------- ",
  30 + //my stuff
  31 + "Simulate Spectrum ",
  32 + " GPU Portion ",
  33 + "Kramers-Kronig ",
  34 +
  35 + //end my stuff
  36 +
  37 +};
  38 +#ifdef WIN32
  39 +#include <stdio.h>
  40 +#include <windows.h>
  41 +#include <float.h>
  42 +
  43 +#include <iostream>
  44 +#include <iomanip>
  45 +
  46 +//-------------------------------------------------------------------------------
  47 +
  48 +class PerformanceData
  49 +{
  50 +public:
  51 + PerformanceData() {
  52 + ClearAll();
  53 + QueryPerformanceFrequency(&cps);
  54 + }
  55 + ~PerformanceData() {}
  56 +
  57 + void ClearAll()
  58 + {
  59 + for ( int i=0; i<PERFORMANCE_DATA_TYPE_COUNT; i++ ) {
  60 + for ( int j=0; j<256; j++ ) times[i][j] = 0;
  61 + pos[i] = 0;
  62 + minTime[i] = 0xFFFFFFFF;
  63 + maxTime[i] = 0;
  64 + totalTime[i] = 0;
  65 + dataReady[i] = false;
  66 + }
  67 + }
  68 +
  69 + void StartTimer( int type ) {
  70 + QueryPerformanceCounter( &startTime[type] );
  71 + }
  72 + void EndTimer( int type ) {
  73 + LARGE_INTEGER endTime;
  74 + QueryPerformanceCounter( &endTime );
  75 + double t = (double)(endTime.QuadPart - startTime[type].QuadPart);
  76 + //unsigned int t = GetTickCount() - startTime[type];
  77 + if ( t < minTime[type] ) minTime[type] = t;
  78 + if ( t > maxTime[type] ) maxTime[type] = t;
  79 + totalTime[type] -= times[type][ pos[type] ];
  80 + times[type][ pos[type] ] = t;
  81 + totalTime[type] += t;
  82 + pos[type]++;
  83 + if ( pos[type] == 0 ) dataReady[type] = true;
  84 + }
  85 +
  86 + void PrintResult( ostream &os,int i=PERFORMANCE_DATA_TYPE_COUNT)
  87 + {
  88 + os.setf(ios::fixed);
  89 + if ((i<PERFORMANCE_DATA_TYPE_COUNT)&&(i>=0)) {
  90 + double a = GetAvrgTime(i);
  91 + if ( a )
  92 + os<< PDTypeNames[i]<<" : avrg="<<setw(8)<<setprecision(3)<<a<<"\tmin="<<setw(8)<<setprecision(3)<< GetMinTime(i) <<"\tmax="<<setw(8)<<setprecision(3)<< GetMaxTime(i) <<endl ;
  93 + else
  94 + os<< PDTypeNames[i]<<" : avrg= -----\tmin= -----\tmax= -----"<<endl;
  95 + }
  96 + }
  97 +
  98 + void PrintResults( ostream &os)
  99 + {
  100 + for ( int i=0; i<PERFORMANCE_DATA_TYPE_COUNT; i++ )
  101 + PrintResult(os,i);
  102 + }
  103 +
  104 + double GetLastTime( int type ) {
  105 + return times[type][pos[type]];
  106 + }
  107 + double GetAvrgTime( int type ) {
  108 + double a = 1000.0 * totalTime[type] / (float)cps.QuadPart / ( (dataReady[type]) ? 256.0 : (double)pos[type] );
  109 + return (_finite(a))? a:0;
  110 + }
  111 + double GetMinTime( int type ) {
  112 + return 1000.0 * minTime[type] / (float)cps.LowPart;
  113 + }
  114 + double GetMaxTime( int type ) {
  115 + return 1000.0 * maxTime[type] / (float)cps.LowPart;
  116 + }
  117 +
  118 +private:
  119 + double times[PERFORMANCE_DATA_TYPE_COUNT][256];
  120 + unsigned char pos[PERFORMANCE_DATA_TYPE_COUNT];
  121 + LARGE_INTEGER startTime[PERFORMANCE_DATA_TYPE_COUNT];
  122 + double minTime[ PERFORMANCE_DATA_TYPE_COUNT ];
  123 + double maxTime[ PERFORMANCE_DATA_TYPE_COUNT ];
  124 + double totalTime[ PERFORMANCE_DATA_TYPE_COUNT ];
  125 + bool dataReady[ PERFORMANCE_DATA_TYPE_COUNT ];
  126 + LARGE_INTEGER cps;
  127 +};
  128 +
  129 +//-------------------------------------------------------------------------------
  130 +#else
  131 +
  132 +class PerformanceData {
  133 +public:
  134 + PerformanceData() {
  135 + ;
  136 + };
  137 + ~PerformanceData() {
  138 + ;
  139 + };
  140 + void ClearAll() {
  141 + ;
  142 + };
  143 + void StartTimer( int type ) {
  144 + ;
  145 + };
  146 + void EndTimer( int type ) {
  147 + ;
  148 + };
  149 + void PrintResults( ostream &os) {
  150 + ;
  151 + };
  152 + void PrintResult( ostream &os, int i=PERFORMANCE_DATA_TYPE_COUNT) {
  153 + ;
  154 + };
  155 + double GetLastTime( int type ) {
  156 + return 0.0;
  157 + };
  158 + double GetAvrgTime( int type ) {
  159 + return 0.0;
  160 + };
  161 + double GetMinTime( int type ) {
  162 + return 0.0;
  163 + };
  164 + double GetMaxTime( int type ) {
  165 + return 0.0;
  166 + };
  167 +};
  168 +
  169 +#endif
  170 +//-------------------------------------------------------------------------------
  171 +
  172 +extern PerformanceData PD;
  173 +
  174 +//-------------------------------------------------------------------------------
... ...
SimulateSpectrum.cpp
  Show/Hide comments   View file @8ffb837
1   -#include <math.h>
2   -#include <complex>
3   -#include <iostream>
4   -#include <fstream>
5   -#include "globals.h"
6   -#include <QProgressDialog>
7   -#include <stdlib.h>
8   -//#include "cufft.h"
9   -using namespace std;
10   -
11   -#define pi 3.14159
12   -
13   -typedef complex<double> scComplex;
14   -
15   -extern int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
16   - complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
17   -extern int bessjyv(double v,double x,double &vm,double *jv,double *yv,
18   - double *djv,double *dyv);
19   -
20   -complex<double> Jl_neg(complex<double> x)
21   -{
22   - //this function computes the bessel function of the first kind Jl(x) for l = -0.5
23   - return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
24   -}
25   -
26   -double Jl_neg(double x)
27   -{
28   - //this function computes the bessel function of the first kind Jl(x) for l = -0.5
29   - return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
30   -}
31   -
32   -double Yl_neg(double x)
33   -{
34   - //this function computes the bessel function of the second kind Yl(x) for l = -0.5;
35   - return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
36   -}
37   -
38   -void computeB(complex<double>* B, double radius, complex<double> refIndex, double lambda, int Nl)
39   -{
40   - double k = (2*pi)/lambda;
41   - int b = 2;
42   -
43   - //allocate space for the real bessel functions
44   - double* jv = (double*)malloc(sizeof(double)*(Nl+b));
45   - double* yv = (double*)malloc(sizeof(double)*(Nl+b));
46   - double* jvp = (double*)malloc(sizeof(double)*(Nl+b));
47   - double* yvp = (double*)malloc(sizeof(double)*(Nl+b));
48   -
49   - //allocate space for the complex bessel functions
50   - complex<double>* cjv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
51   - complex<double>* cyv = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
52   - complex<double>* cjvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
53   - complex<double>* cyvp = (complex<double>*)malloc(sizeof(complex<double>)*(Nl+b));
54   -
55   - double kr = k*radius;
56   - complex<double> knr = k*refIndex*(double)radius;
57   - complex<double> n = refIndex;
58   -
59   - //compute the bessel functions for k*r
60   - double vm;// = Nl - 1;
61   - bessjyv((Nl)+0.5, kr, vm, jv, yv, jvp, yvp);
62   - //cout<<"Nl: "<<Nl<<" vm: "<<vm<<endl;
63   - //printf("Nl: %f, vm: %f\n", (float)Nl, (float)vm);
64   -
65   - //compute the bessel functions for k*n*r
66   - cbessjyva((Nl)+0.5, knr, vm, cjv, cyv, cjvp, cyvp);
67   -
68   - //scale factor for spherical bessel functions
69   - double scale_kr = sqrt(pi/(2.0*kr));
70   - complex<double> scale_knr = sqrt(pi/(2.0*knr));
71   -
72   - complex<double> numer, denom;
73   - double j_kr;
74   - double y_kr;
75   - complex<double> j_knr;
76   - complex<double> j_d_knr;
77   - double j_d_kr;
78   - complex<double> h_kr;
79   - complex<double> h_d_kr;
80   - complex<double> h_neg;
81   - complex<double> h_pos;
82   -
83   - //cout<<"B coefficients:"<<endl;
84   - for(int l=0; l<Nl; l++)
85   - {
86   - //compute the spherical bessel functions
87   - j_kr = jv[l] * scale_kr;
88   - y_kr = yv[l] * scale_kr;
89   - j_knr = cjv[l] * scale_knr;
90   -
91   - //compute the Hankel function
92   - h_kr = complex<double>(j_kr, y_kr);
93   -
94   - //compute the derivatives
95   - if(l == 0)
96   - {
97   - //spherical bessel functions for l=0
98   - j_d_kr = scale_kr * (Jl_neg(kr) - (jv[l] + kr*jv[l+1])/kr )/2.0;
99   - j_d_knr = scale_knr * ( Jl_neg(knr) - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
100   - h_neg = complex<double>(scale_kr*Jl_neg(kr), scale_kr*Yl_neg(kr));
101   - h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
102   - h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
103   - }
104   - else
105   - {
106   - //spherical bessel functions
107   - j_d_kr = scale_kr * (jv[l-1] - (jv[l] + kr*jv[l+1])/kr )/2.0;
108   - j_d_knr = scale_knr * ( cjv[l-1] - (cjv[l] + knr*cjv[l+1])/knr )/2.0;
109   - h_neg = complex<double>(scale_kr*jv[l-1], scale_kr*yv[l-1]);
110   - h_pos = complex<double>(scale_kr*jv[l+1], scale_kr*yv[l+1]);
111   - h_d_kr = (h_neg - (h_kr + kr*h_pos)/kr)/2.0;
112   - }
113   -
114   - numer = j_kr*j_d_knr*n - j_knr*j_d_kr;
115   - denom = j_knr*h_d_kr - h_kr*j_d_knr*n;
116   - B[l] = numer/denom;
117   -
118   - //B[l] = scComplex(temp.real(), temp.imag());
119   - //cout<<B[l]<<endl;
120   - }
121   -
122   - free(jv);
123   - free(yv);
124   - free(jvp);
125   - free(yvp);
126   - free(cjv);
127   - free(cyv);
128   - free(cjvp);
129   - free(cyvp);
130   -}
131   -
132   -void Legendre(double* P, double x, int Nl)
133   -{
134   - //computes the legendre polynomials from orders 0 to Nl-1
135   - P[0] = 1;
136   - if(Nl == 1) return;
137   - P[1] = x;
138   - for(int l = 2; l < Nl; l++)
139   - {
140   - P[l] = ((2*l - 1)*x*P[l-1] - (l - 1)*P[l-2])/l;
141   - }
142   -
143   -}
144   -
145   -complex<double> integrateUi(double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
146   -{
147   - /*This function integrates the incident field of magnitude M in the far zone
148   - in order to evaluate the field at the central pixel of a detector.
149   - cNAi = condenser inner angle
150   - cNAo = condenser outer angle
151   - oNAi = objective inner angle
152   - oNAo = objective outer angle
153   - M = field magnitude*/
154   -
155   - double alphaIn = max(cAngleI, oAngleI);
156   - double alphaOut = min(cAngleO,oAngleO);
157   -
158   - complex<double> Ui;
159   - if(alphaIn > alphaOut)
160   - Ui = complex<double>(0.0, 0.0);
161   - else
162   - Ui = complex<double>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);
163   -
164   - return Ui;
165   -
166   -}
167   -
168   -void computeCondenserAlpha(double* alpha, int Nl, double cAngleI, double cAngleO)
169   -{
170   - /*This function computes the condenser integral in order to build the field of incident light
171   - alpha = list of Nl floating point values representing the condenser alpha as a function of l
172   - Nl = number of orders in the incident field
173   - cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/
174   -
175   - //compute the Legendre polynomials for the condenser aperature
176   - double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
177   - Legendre(PcNAo, cos(cAngleO), Nl+1);
178   - double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
179   - Legendre(PcNAi, cos(cAngleI), Nl+1);
180   -
181   - for(int l=0; l<Nl; l++)
182   - {
183   - //integration term
184   - if(l == 0)
185   - alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
186   - else
187   - alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
188   -
189   - alpha[l] *= 2 * pi;
190   - }
191   -
192   -}
193   -
194   -complex<double> integrateUs(double r, double lambda, complex<double> eta,
195   - double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
196   -{
197   - /*This function integrates the incident field of magnitude M in the far zone
198   - in order to evaluate the field at the central pixel of a detector.
199   - r = sphere radius
200   - lambda = wavelength
201   - eta = index of refraction
202   - cNAi = condenser inner NA
203   - cNAo = condenser outer NA
204   - oNAi = objective inner NA
205   - oNAo = objective outer NA
206   - M = field magnitude*/
207   -
208   - //compute the required number of orders
209   - double k = 2*pi/lambda;
210   - int Nl = (int)ceil( k + 4 * exp(log(k*r)/3) + 3 );
211   -
212   - //compute the material coefficients B
213   - complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>)*Nl);
214   - //compute the Legendre polynomials for the condenser and objective aperatures
215   - double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
216   - Legendre(PcNAo, cos(cAngleO), Nl+1);
217   - double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
218   - Legendre(PcNAi, cos(cAngleI), Nl+1);
219   -
220   - double* PoNAo = (double*)malloc(sizeof(double)*(Nl+1));
221   - Legendre(PoNAo, cos(oAngleO), Nl+1);
222   - double* PoNAi = (double*)malloc(sizeof(double)*(Nl+1));
223   - Legendre(PoNAi, cos(oAngleI), Nl+1);
224   -
225   - //store the index of refraction;
226   - complex<double> IR(eta.real(), eta.imag());
227   -
228   - //compute the scattering coefficients
229   - computeB(B, r, IR, lambda, Nl);
230   -
231   - //aperature terms for the condenser (alpha) and objective (beta)
232   - double alpha;
233   - double beta;
234   - double c;
235   - complex<double> Us(0.0, 0.0);
236   -
237   - for(int l=0; l<Nl; l++)
238   - {
239   - //integration term
240   - if(l == 0)
241   - {
242   - alpha = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
243   - beta = -PoNAo[l+1] + PoNAo[0] + PoNAi[l+1] - PoNAi[0];
244   - }
245   - else
246   - {
247   - alpha = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
248   - beta = -PoNAo[l+1] + PoNAo[l-1] + PoNAi[l+1] - PoNAi[l-1];
249   - }
250   - c = (2*pi)/(2.0 * l + 1.0);
251   - Us += c * alpha * beta * B[l] * M;
252   -
253   -
254   - }
255   - free(PcNAo);
256   - free(PcNAi);
257   - free(PoNAo);
258   - free(PoNAi);
259   - free(B);
260   -
261   - return Us;
262   -
263   -}
264   -
265   -void pointSpectrum()
266   -{
267   - PD.StartTimer(SIMULATE_SPECTRUM);
268   - //clear the previous spectrum
269   - SimSpectrum.clear();
270   -
271   - double dNu = 2.0f;
272   - double lambda;
273   -
274   - //compute the angles based on NA
275   - double cAngleI = asin(cNAi);
276   - double cAngleO = asin(cNAo);
277   - double oAngleI = asin(oNAi);
278   - double oAngleO = asin(oNAo);
279   -
280   - //implement a reflection-mode system if necessary
281   - if(opticsMode == ReflectionOpticsType){
282   -
283   - //set the condenser to match the objective
284   - cAngleI = oAngleI;
285   - cAngleO = oAngleO;
286   -
287   - //invert the objective
288   - oAngleO = pi - cAngleI;
289   - oAngleI = pi - cAngleO;
290   - }
291   -
292   - //integrate the incident field at the detector position
293   - complex<double> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
294   - double I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
295   - I0 *= scaleI0;
296   -
297   -
298   -
299   - //double I;
300   - SpecPair temp;
301   - double nu;
302   - complex<double> eta;
303   - complex<double> Us, U;
304   -
305   - double vecLen = 0.0;
306   - for(unsigned int i=0; i<EtaK.size(); i++)
307   - {
308   - nu = EtaK[i].nu;
309   - lambda = 10000.0f/nu;
310   - if(applyMaterial)
311   - eta = complex<double>(EtaN[i].A, EtaK[i].A);
312   - else
313   - eta = complex<double>(baseIR, 0.0);
314   -
315   -
316   - //integrate the scattered field at the detector position
317   - Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
318   - U = Us + Ui;
319   - double I = U.real() * U.real() + U.imag() * U.imag();
320   -
321   - temp.nu = nu;
322   -
323   - //set the spectrum value based on the current display type
324   - if(dispSimType == AbsorbanceSpecType)
325   - temp.A = -log10(I/I0);
326   - else
327   - temp.A = I;
328   -
329   - if(dispNormalize)
330   - vecLen += temp.A * temp.A;
331   -
332   - SimSpectrum.push_back(temp);
333   - }
334   - vecLen = sqrt(vecLen);
335   -
336   - if(dispNormalize)
337   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
338   - SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;
339   -
340   - PD.EndTimer(SIMULATE_SPECTRUM);
341   -}
342   -
343   -void updateSpectrum(double* I, double I0, int n)
344   -{
345   - SimSpectrum.clear();
346   - SpecPair temp;
347   -
348   - //update the displayed spectrum based on the computed intensity I
349   - for(int i=0; i<n; i++)
350   - {
351   - temp.nu = EtaK[i].nu;
352   -
353   - //set the spectrum value based on the current display type
354   - if(dispSimType == AbsorbanceSpecType)
355   - {
356   - temp.A = -log10(I[i]/I0);
357   - //cout<<temp.nu<<" "<<I[i]<<endl;
358   - }
359   - else
360   - {
361   - temp.A = I[i];
362   - if(useSourceSpectrum)
363   - temp.A *= SourceResampled[i].A;
364   - }
365   -
366   - SimSpectrum.push_back(temp);
367   - }
368   -}
369   -
370   -void computeCassegrainAngles(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO)
371   -{
372   - //compute the angles based on NA
373   - cAngleI = asin(cNAi);
374   - cAngleO = asin(cNAo);
375   - oAngleI = asin(oNAi);
376   - oAngleO = asin(oNAo);
377   -
378   - //implement a reflection-mode system if necessary
379   - if(opticsMode == ReflectionOpticsType){
380   -
381   - //set the condenser to match the objective
382   - cAngleI = oAngleI;
383   - cAngleO = oAngleO;
384   -
385   - //invert the objective
386   - oAngleO = pi - cAngleI;
387   - oAngleI = pi - cAngleO;
388   - }
389   -
390   -
391   -}
392   -
393   -int computeNl()
394   -{
395   - double maxNu = EtaK.back().nu;
396   - double maxLambda = 10000.0f/maxNu;
397   - double k = 2*pi/maxLambda;
398   - int Nl = (int)ceil( k + 4 * exp(log(k*radius)/3) + 3 );
399   -
400   - return Nl;
401   -}
402   -
403   -void computeBArray(complex<double>* B, int Nl, int nLambda)
404   -{
405   - double nu;
406   - complex<double> eta;
407   - double* Lambda = (double*)malloc(sizeof(double) * nLambda);
408   -
409   - //for each wavenumber nu
410   - for(unsigned int i=0; i<EtaK.size(); i++)
411   - {
412   - //compute information based on wavelength and material
413   - nu = EtaK[i].nu;
414   - Lambda[i] = 10000.0f/nu;
415   - if(applyMaterial)
416   - eta = complex<double>(EtaN[i].A, EtaK[i].A);
417   - else
418   - eta = complex<double>(baseIR, 0.0);
419   -
420   - //allocate memory for the scattering coefficients
421   - //complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
422   -
423   - complex<double> IR(eta.real(), eta.imag());
424   - computeB(&B[i * Nl], radius, IR, Lambda[i], Nl);
425   - }
426   -}
427   -
428   -void computeOpticalParameters(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO, double& I0, double* alpha, int Nl)
429   -{
430   - computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO);
431   -
432   - //evaluate the incident field intensity
433   - I0 = 0.0;
434   - complex<double> Ui;
435   -
436   - Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
437   - I0 = Ui.real()*2*pi;
438   -
439   - //compute alpha (condenser integral)
440   - computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
441   -}
442   -
443   -void gpuDetectorSpectrum(int numSamples)
444   -{
445   - //integrate across the objective aperature and calculate the resulting intensity on a detector
446   - PD.StartTimer(SIMULATE_SPECTRUM);
447   - //clear the previous spectrum
448   - SimSpectrum.clear();
449   -
450   - //compute Nl (maximum order of the spectrum)
451   - int Nl = computeNl();
452   -
453   - double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
454   - double cAngleI, cAngleO, oAngleI, oAngleO, I0;
455   - computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
456   -
457   - //allocate space for a list of wavelengths
458   - int nLambda = EtaK.size();
459   -
460   - //allocate space for the 2D array (Nl x nu) of scattering coefficients
461   - complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
462   - computeBArray(B, Nl, nLambda);
463   -
464   -
465   -
466   - //allocate temporary space for the spectrum
467   - double* I = (double*)malloc(sizeof(double) * EtaK.size());
468   -
469   - //compute the spectrum on the GPU
470   - PD.StartTimer(SIMULATE_GPU);
471   - cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
472   - PD.EndTimer(SIMULATE_GPU);
473   -
474   - updateSpectrum(I, I0, nLambda);
475   -
476   - PD.EndTimer(SIMULATE_SPECTRUM);
477   -
478   -}
479   -
480   -void SimulateSpectrum()
481   -{
482   - if(pointDetector)
483   - pointSpectrum();
484   - else
485   - gpuDetectorSpectrum(objectiveSamples);
486   - //detectorSpectrum(objectiveSamples);
487   -}
488   -
489   -double absorbanceDistortion(){
490   -
491   - //compute the mean of the spectrum
492   - double sumSim = 0.0;
493   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
494   - {
495   - sumSim += SimSpectrum[i].A;
496   - }
497   - double meanSim = sumSim/SimSpectrum.size();
498   -
499   - //compute the distortion (MSE from the mean)
500   - double sumSE = 0.0;
501   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
502   - {
503   - sumSE += pow(SimSpectrum[i].A - meanSim, 2);
504   - }
505   - double MSE = sumSE / SimSpectrum.size();
506   -
507   - return MSE;
508   -}
509   -
510   -double intensityDistortion(){
511   -
512   - //compute the mean intensity of the spectrum
513   - double sumSim = 0.0;
514   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
515   - {
516   - sumSim += pow(SimSpectrum[i].A, 2);
517   - }
518   - double magSim = sqrt(sumSim);
519   -
520   - //compute the distortion (MSE from the mean)
521   - double proj = 0.0;
522   - for(unsigned int i=0; i<SimSpectrum.size(); i++)
523   - {
524   - proj += (SimSpectrum[i].A/magSim) * (1.0/SimSpectrum.size());
525   - }
526   - double error = proj;
527   -
528   - return error;
529   -}
530   -
531   -void DistortionMap(float* distortionMap, int nSteps){
532   - ofstream outFile("distortion.txt");
533   -
534   - //set the parameters for the distortion simulation
535   - double range = 0.4;
536   - double step = (range)/(nSteps-1);
537   -
538   - oNAi = 0.2;
539   - oNAo = 0.5;
540   -
541   - double startNAi = 0.0;
542   - double startNAo = 0.3;
543   -
544   - //compute the optical parameters
545   - //compute Nl (maximum order of the spectrum)
546   - int Nl = computeNl();
547   -
548   - double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
549   - double cAngleI, cAngleO, oAngleI, oAngleO, I0;
550   -
551   - //allocate space for a list of wavelengths
552   - int nLambda = EtaK.size();
553   -
554   - //allocate temporary space for the spectrum
555   - double* I = (double*)malloc(sizeof(double) * EtaK.size());
556   -
557   - //calculate the material parameters
558   - //allocate space for the 2D array (Nl x nu) of scattering coefficients
559   - complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
560   - computeBArray(B, Nl, nLambda);
561   -
562   - QProgressDialog progress("Computing distortion map...", "Stop", 0, nSteps * nSteps);
563   - progress.setWindowModality(Qt::WindowModal);
564   -
565   - double D;
566   - double e = 0.001;
567   - int i, o;
568   - for(i=0; i<nSteps; i++)
569   - {
570   - for(o=0; o<nSteps; o++)
571   - {
572   - //update the progress bar and check for an exit
573   - progress.setValue(i * nSteps + o);
574   - if (progress.wasCanceled())
575   - break;
576   -
577   - //set the current optical parameters
578   - cNAi = startNAi + i * step;
579   - cNAo = startNAo + o * step;
580   - //cout<<cNAi<<" "<<cNAo<<endl;
581   -
582   - //set the current optical parameters
583   - //cNAi = i;
584   - //cNAo = o;
585   -
586   - //compute the optical parameters
587   - computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
588   -
589   - //simulate the spectrum
590   - cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);
591   - updateSpectrum(I, I0, nLambda);
592   -
593   - if(dispSimType == AbsorbanceSpecType)
594   - {
595   - if(cNAi >= cNAo || cNAi >= oNAo || oNAi >= cNAo || oNAi >= oNAo)
596   - D = -1.0;
597   - else
598   - D = absorbanceDistortion();
599   - }
600   - else
601   - {
602   - if(cNAi >= cNAo || oNAi >= oNAo)
603   - D = -1.0;
604   - else
605   - D = intensityDistortion();
606   - }
607   - distortionMap[o * nSteps + i] = D;
608   - outFile<<D<<" ";
609   - }
610   - outFile<<endl;
611   - //cout<<i<<endl;
612   - }
613   -
614   - progress.setValue(nSteps * nSteps);
615   -
616   - outFile.close();
617   -}
  1 +#include <math.h>
  2 +#include <complex>
  3 +#include <iostream>
  4 +#include <fstream>
  5 +#include "globals.h"
  6 +#include <QProgressDialog>
  7 +#include <stdlib.h>
  8 +//#include "cufft.h"
  9 +using namespace std;
  10 +
  11 +#define pi 3.14159
  12 +
  13 +typedef complex<double> scComplex;
  14 +
  15 +extern int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
  16 + complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
  17 +extern int bessjyv(double v,double x,double &vm,double *jv,double *yv,
  18 + double *djv,double *dyv);
  19 +
  20 +complex<double> Jl_neg(complex<double> x)
  21 +{
  22 + //this function computes the bessel function of the first kind Jl(x) for l = -0.5
  23 + return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
  24 +}
  25 +
  26 +double Jl_neg(double x)
  27 +{
  28 + //this function computes the bessel function of the first kind Jl(x) for l = -0.5
  29 + return ( sqrt(2.0/pi) * cos(x) )/sqrt(x);
  30 +}
  31 +
  32 +double Yl_neg(double x)
  33 +{
  34 + //this function computes the bessel function of the second kind Yl(x) for l = -0.5;
  35 + return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
  36 +}
  37 +
  38 +void computeB(complex<double>* B, double radius, complex<double> refIndex, double lambda, int Nl)
  39 +{
  40 + double k = (2*pi)/lambda;
  41 + int b = 2;
  42 +
  43 + //allocate space for the real bessel functions
  44 + double* jv = (double*)malloc(sizeof(double)*(Nl+b));
  45 + double* yv = (double*)malloc(sizeof(double)*(Nl+b));
  46 + double* jvp = (double*)malloc(sizeof(double)*(Nl+b));