Commit 52a5fe9d0778ae4c18411ad1dcee2785e451e78c

Authored by dmayerich
1 parent 212273d5

Added double support.

Fixed all known errors and warnings, with the exception of:
tmpxft_000013d0_00000000-4_cudaMain.ptx, line 86; warning : Double is not supported. Demoting to float

I'm not sure what is causing that, but it doesn't seem to be introducing any mathematical errors.  If anybody has an idea, please let me know!
BESSIK.CPP
... ... @@ -348,7 +348,7 @@ int bessikv(double v,double x,double &vm,double *iv,double *kv,
348 348 double *ivp,double *kvp)
349 349 {
350 350 double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum;
351   - double f,f0,f1,f2,ct,cs,wa,gan,ww,w0,v0n;
  351 + double f,f1,f2,ct,cs,wa,gan,ww,w0,v0n;
352 352 double r1,r2,bk0,bk1,bk2,a2,cb;
353 353 int n,k,kz,m;
354 354  
... ...
BESSJY.CPP
... ... @@ -274,7 +274,7 @@ int msta1(double x,int mp)
274 274 n1 = n0+5;
275 275 f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp;
276 276 for (i=0;i<20;i++) {
277   - nn = n1-(n1-n0)/(1.0-f0/f1);
  277 + nn = (int)(n1-(n1-n0)/(1.0-f0/f1));
278 278 f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp;
279 279 if (abs(nn-n1) < 1) break;
280 280 n0 = n1;
... ... @@ -305,7 +305,7 @@ int msta2(double x,int n,int mp)
305 305 n1 = n0+5;
306 306 f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj;
307 307 for (i=0;i<20;i++) {
308   - nn = n1-(n1-n0)/(1.0-f0/f1);
  308 + nn = (int)(n1-(n1-n0)/(1.0-f0/f1));
309 309 f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj;
310 310 if (abs(nn-n1) < 1) break;
311 311 n0 = n1;
... ... @@ -404,7 +404,7 @@ int bessjyna(int n,double x,int &amp;nm,double *jn,double *yn,
404 404 int bessjynb(int n,double x,int &nm,double *jn,double *yn,
405 405 double *jnp,double *ynp)
406 406 {
407   - double t1,t2,f,f0,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
  407 + double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv;
408 408 double ec,bs,byk,p0,p1,q0,q1;
409 409 static double a[] = {
410 410 -0.7031250000000000e-1,
... ... @@ -536,7 +536,7 @@ int bessjyv(double v,double x,double &amp;vm,double *jv,double *yv,
536 536 {
537 537 double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2;
538 538 double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq;
539   - double b,ec,w0,w1,bjy0,bjy1,bju0,bju1,pv0,pv1,byvk;
  539 + double b,ec,w0,w1,bju0,bju1,pv0,pv1,byvk;
540 540 int j,k,l,m,n,kz;
541 541  
542 542 x2 = x*x;
... ...
CBESSJY.CPP
... ... @@ -189,7 +189,7 @@ int cbessjyna(int n,complex&lt;double&gt; z,int &amp;nm,complex&lt;double&gt; *cj,
189 189 complex<double> cs,cg0,cg1,cyk,cyl1,cyl2,cylk,cp11,cp12,cp21,cp22;
190 190 complex<double> ch0,ch1,ch2;
191 191 double a0,yak,ya1,ya0,wa;
192   - int m,k,kz,lb,lb0;
  192 + int m,k,lb,lb0;
193 193  
194 194 if (n < 0) return 1;
195 195 a0 = abs(z);
... ...
EstimateMaterial.cpp
1 1 #include "globals.h"
2 2 #define PI 3.14159
3 3  
4   -float CalculateError(float* E)
  4 +double CalculateError(double* E)
5 5 {
6 6 //Calculate the error between the Reference Spectrum and the Simulated Spectrum
7   - float sumE = 0.0;
  7 + double sumE = 0.0;
8 8 int nVals = RefSpectrum[currentSpec].size();
9   - float nu;
  9 + double nu;
10 10 for(int i=0; i<nVals; i++)
11 11 {
12 12 nu = RefSpectrum[currentSpec][i].nu;
... ... @@ -17,16 +17,16 @@ float CalculateError(float* E)
17 17 return sumE/nVals;
18 18 }
19 19  
20   -void EstimateK(float* E)
  20 +void EstimateK(double* E)
21 21 {
22 22 int nVals = RefSpectrum[currentSpec].size();
23   - float nuStart = RefSpectrum[currentSpec].front().nu;
24   - float nuEnd = RefSpectrum[currentSpec].back().nu;
  23 + double nuStart = RefSpectrum[currentSpec].front().nu;
  24 + double nuEnd = RefSpectrum[currentSpec].back().nu;
25 25  
26   - float r = radius/10000.0;
27   - float nu;
28   - float dNu = (nuEnd - nuStart)/(nVals-1);
29   - float eScale;
  26 + double r = radius/10000.0;
  27 + double nu;
  28 + double dNu = (nuEnd - nuStart)/(nVals-1);
  29 + double eScale;
30 30 for(int i=0; i<nVals; i++)
31 31 {
32 32 nu = nuStart + i*2;
... ... @@ -50,7 +50,7 @@ void EstimateMaterial()
50 50  
51 51 //insert the default material properties
52 52 SpecPair temp;
53   - for(int s=0; s<RefSpectrum[currentSpec].size(); s++)
  53 + for(int s=0; s<(int)RefSpectrum[currentSpec].size(); s++)
54 54 {
55 55 //the real part of the IR is the user-specified baseline IR
56 56 temp.nu = RefSpectrum[currentSpec][s].nu;
... ... @@ -63,17 +63,17 @@ void EstimateMaterial()
63 63  
64 64  
65 65 //allocate space to store the list of error values
66   - float* E = (float*)malloc(sizeof(float) * RefSpectrum[currentSpec].size());
  66 + double* E = (double*)malloc(sizeof(double) * RefSpectrum[currentSpec].size());
67 67 //copy the absorbance values into a linear array
68   - float* k = (float*)malloc(sizeof(float) * EtaK.size());
69   - float* n = (float*)malloc(sizeof(float) * EtaN.size());
  68 + double* k = (double*)malloc(sizeof(double) * EtaK.size());
  69 + double* n = (double*)malloc(sizeof(double) * EtaN.size());
70 70  
71 71 //iterate to solve for both n and k
72   - float sumE = 99999.9;
73   - int i=0;
  72 + double sumE = 99999.9;
  73 + int j=0;
74 74 //clear the console
75 75 system("cls");
76   - while(sumE > minMSE && i < maxFitIter)
  76 + while(sumE > minMSE && j < maxFitIter)
77 77 {
78 78 //simulate a spectrum based on the current IR
79 79 SimulateSpectrum();
... ... @@ -86,13 +86,13 @@ void EstimateMaterial()
86 86  
87 87 //use Kramers-Kronig to compute n
88 88  
89   - for(int i=0; i<EtaK.size(); i++)
  89 + for(unsigned int i=0; i<EtaK.size(); i++)
90 90 k[i] = EtaK[i].A;
91 91 cudaKramersKronig(n, k, EtaK.size(), EtaK.front().nu, EtaK.back().nu, baseIR);
92 92  
93 93 //copy the real part of the index of refraction into the vector
94 94 EtaN.clear();
95   - for(int i=0; i<EtaK.size(); i++)
  95 + for(int i=0; i<(int)EtaK.size(); i++)
96 96 {
97 97 temp.nu = EtaK[i].nu;
98 98 temp.A = n[i];
... ... @@ -100,7 +100,7 @@ void EstimateMaterial()
100 100 }
101 101  
102 102 cout<<" E = "<<sumE<<endl;
103   - i++;
  103 + j++;
104 104 //SaveSpectrum(n, nVals, "simNj.txt");
105 105 //SaveSpectrum(k, nVals, "simKj.txt");
106 106 //SaveSpectrum(simSpec, nVals, "simSpec.txt");
... ...
FileIO.cpp
... ... @@ -23,11 +23,11 @@ vector&lt;SpecPair&gt; LoadSpectrum(string filename)
23 23 }
24 24  
25 25 //compute the minimum and maximum input wavenumbers
26   - float inMin = S.front().nu;
27   - float inMax = S.back().nu;
  26 + double inMin = S.front().nu;
  27 + double inMax = S.back().nu;
28 28  
29   - int nuMin = ceil(inMin);
30   - int nuMax = floor(inMax);
  29 + int nuMin = (int)ceil(inMin);
  30 + int nuMax = (int)floor(inMax);
31 31  
32 32 //make sure both are either even or odd
33 33 if(nuMin % 2 != nuMax % 2)
... ... @@ -39,7 +39,7 @@ vector&lt;SpecPair&gt; LoadSpectrum(string filename)
39 39 //allocate space for the spectrum
40 40 vector<SpecPair> outSpec;
41 41  
42   - float nu, highVal, lowVal, a;
  42 + double nu, highVal, lowVal, a;
43 43 int j=1;
44 44 for(int i=0; i<nVals; i++)
45 45 {
... ... @@ -52,7 +52,7 @@ vector&lt;SpecPair&gt; LoadSpectrum(string filename)
52 52 else
53 53 {
54 54 //move to the correct position in the input array
55   - while(j < S.size()-1 && S[j].nu <= nu)
  55 + while(j < (int)S.size()-1 && S[j].nu <= nu)
56 56 j++;
57 57  
58 58 lowVal = S[j-1].nu;
... ... @@ -88,7 +88,7 @@ vector&lt;SpecPair&gt; SetReferenceSpectrum(char* text)
88 88 void SaveK(string fileName)
89 89 {
90 90 ofstream outFile(fileName.c_str());
91   - for(int i=0; i<EtaK.size(); i++)
  91 + for(unsigned int i=0; i<EtaK.size(); i++)
92 92 {
93 93 outFile<<EtaK[i].nu<<" ";
94 94 outFile<<EtaK[i].A<<endl;
... ... @@ -99,7 +99,7 @@ void SaveK(string fileName)
99 99 void SaveN(string fileName)
100 100 {
101 101 ofstream outFile(fileName.c_str());
102   - for(int i=0; i<EtaN.size(); i++)
  102 + for(unsigned int i=0; i<EtaN.size(); i++)
103 103 {
104 104 outFile<<EtaN[i].nu<<" ";
105 105 outFile<<EtaN[i].A<<endl;
... ... @@ -110,7 +110,7 @@ void SaveN(string fileName)
110 110 void SaveSimulation(string fileName)
111 111 {
112 112 ofstream outFile(fileName.c_str());
113   - for(int i=0; i<SimSpectrum.size(); i++)
  113 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
114 114 {
115 115 outFile<<SimSpectrum[i].nu<<" ";
116 116 outFile<<SimSpectrum[i].A<<endl;
... ...
PerformanceData.h
... ... @@ -76,7 +76,7 @@ public:
76 76 void EndTimer( int type ) {
77 77 LARGE_INTEGER endTime;
78 78 QueryPerformanceCounter( &endTime );
79   - double t = endTime.QuadPart - startTime[type].QuadPart;
  79 + double t = (double)(endTime.QuadPart - startTime[type].QuadPart);
80 80 //unsigned int t = GetTickCount() - startTime[type];
81 81 if ( t < minTime[type] ) minTime[type] = t;
82 82 if ( t > maxTime[type] ) maxTime[type] = t;
... ...
  1 +I've just added support for double-precision values. This allows simulation of extremely small (nano-scale) particles, however a GPU with a compute capability of at least 1.3 is required.
  2 +
1 3 TrueEyes requires the following libraries:
2 4  
3 5 GLUT (OpenGL Utility Toolkit) http://www.opengl.org/resources/libraries/glut/
... ...
SimulateSpectrum.cpp
... ... @@ -8,7 +8,7 @@ using namespace std;
8 8  
9 9 #define pi 3.14159
10 10  
11   -typedef complex<float> scComplex;
  11 +typedef complex<double> scComplex;
12 12  
13 13 int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
14 14 complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
... ... @@ -33,7 +33,7 @@ double Yl_neg(double x)
33 33 return ( sqrt(2.0/pi) * sin(x) )/sqrt(x);
34 34 }
35 35  
36   -void computeB(complex<float>* B, float radius, complex<double> refIndex, float lambda, int Nl)
  36 +void computeB(complex<double>* B, double radius, complex<double> refIndex, double lambda, int Nl)
37 37 {
38 38 double k = (2*pi)/lambda;
39 39 int b = 2;
... ... @@ -111,9 +111,9 @@ void computeB(complex&lt;float&gt;* B, float radius, complex&lt;double&gt; refIndex, float l
111 111  
112 112 numer = j_kr*j_d_knr*n - j_knr*j_d_kr;
113 113 denom = j_knr*h_d_kr - h_kr*j_d_knr*n;
114   - complex<double> temp = numer/denom;
  114 + B[l] = numer/denom;
115 115  
116   - B[l] = scComplex(temp.real(), temp.imag());
  116 + //B[l] = scComplex(temp.real(), temp.imag());
117 117 //cout<<B[l]<<endl;
118 118 }
119 119  
... ... @@ -127,7 +127,7 @@ void computeB(complex&lt;float&gt;* B, float radius, complex&lt;double&gt; refIndex, float l
127 127 free(cyvp);
128 128 }
129 129  
130   -void Legendre(float* P, float x, int Nl)
  130 +void Legendre(double* P, double x, int Nl)
131 131 {
132 132 //computes the legendre polynomials from orders 0 to Nl-1
133 133 P[0] = 1;
... ... @@ -140,7 +140,7 @@ void Legendre(float* P, float x, int Nl)
140 140  
141 141 }
142 142  
143   -complex<float> integrateUi(float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
  143 +complex<double> integrateUi(double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
144 144 {
145 145 /*This function integrates the incident field of magnitude M in the far zone
146 146 in order to evaluate the field at the central pixel of a detector.
... ... @@ -150,20 +150,20 @@ complex&lt;float&gt; integrateUi(float cAngleI, float cAngleO, float oAngleI, float oA
150 150 oNAo = objective outer angle
151 151 M = field magnitude*/
152 152  
153   - float alphaIn = max(cAngleI, oAngleI);
154   - float alphaOut = min(cAngleO,oAngleO);
  153 + double alphaIn = max(cAngleI, oAngleI);
  154 + double alphaOut = min(cAngleO,oAngleO);
155 155  
156   - complex<float> Ui;
  156 + complex<double> Ui;
157 157 if(alphaIn > alphaOut)
158   - Ui = complex<float>(0.0, 0.0);
  158 + Ui = complex<double>(0.0, 0.0);
159 159 else
160   - Ui = complex<float>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);
  160 + Ui = complex<double>(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f);
161 161  
162 162 return Ui;
163 163  
164 164 }
165 165  
166   -void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO)
  166 +void computeCondenserAlpha(double* alpha, int Nl, double cAngleI, double cAngleO)
167 167 {
168 168 /*This function computes the condenser integral in order to build the field of incident light
169 169 alpha = list of Nl floating point values representing the condenser alpha as a function of l
... ... @@ -171,9 +171,9 @@ void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO)
171 171 cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/
172 172  
173 173 //compute the Legendre polynomials for the condenser aperature
174   - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
  174 + double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
175 175 Legendre(PcNAo, cos(cAngleO), Nl+1);
176   - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
  176 + double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
177 177 Legendre(PcNAi, cos(cAngleI), Nl+1);
178 178  
179 179 for(int l=0; l<Nl; l++)
... ... @@ -189,8 +189,8 @@ void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO)
189 189  
190 190 }
191 191  
192   -complex<float> integrateUs(float r, float lambda, complex<float> eta,
193   - float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi)
  192 +complex<double> integrateUs(double r, double lambda, complex<double> eta,
  193 + double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi)
194 194 {
195 195 /*This function integrates the incident field of magnitude M in the far zone
196 196 in order to evaluate the field at the central pixel of a detector.
... ... @@ -204,20 +204,20 @@ complex&lt;float&gt; integrateUs(float r, float lambda, complex&lt;float&gt; eta,
204 204 M = field magnitude*/
205 205  
206 206 //compute the required number of orders
207   - float k = 2*pi/lambda;
208   - int Nl = ceil( k + 4 * exp(log(k*r)/3) + 3 );
  207 + double k = 2*pi/lambda;
  208 + int Nl = (int)ceil( k + 4 * exp(log(k*r)/3) + 3 );
209 209  
210 210 //compute the material coefficients B
211   - complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
  211 + complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>)*Nl);
212 212 //compute the Legendre polynomials for the condenser and objective aperatures
213   - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
  213 + double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1));
214 214 Legendre(PcNAo, cos(cAngleO), Nl+1);
215   - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
  215 + double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1));
216 216 Legendre(PcNAi, cos(cAngleI), Nl+1);
217 217  
218   - float* PoNAo = (float*)malloc(sizeof(float)*(Nl+1));
  218 + double* PoNAo = (double*)malloc(sizeof(double)*(Nl+1));
219 219 Legendre(PoNAo, cos(oAngleO), Nl+1);
220   - float* PoNAi = (float*)malloc(sizeof(float)*(Nl+1));
  220 + double* PoNAi = (double*)malloc(sizeof(double)*(Nl+1));
221 221 Legendre(PoNAi, cos(oAngleI), Nl+1);
222 222  
223 223 //store the index of refraction;
... ... @@ -227,10 +227,10 @@ complex&lt;float&gt; integrateUs(float r, float lambda, complex&lt;float&gt; eta,
227 227 computeB(B, r, IR, lambda, Nl);
228 228  
229 229 //aperature terms for the condenser (alpha) and objective (beta)
230   - float alpha;
231   - float beta;
232   - float c;
233   - complex<float> Us(0.0, 0.0);
  230 + double alpha;
  231 + double beta;
  232 + double c;
  233 + complex<double> Us(0.0, 0.0);
234 234  
235 235 for(int l=0; l<Nl; l++)
236 236 {
... ... @@ -266,14 +266,14 @@ void pointSpectrum()
266 266 //clear the previous spectrum
267 267 SimSpectrum.clear();
268 268  
269   - float dNu = 2.0f;
270   - float lambda;
  269 + double dNu = 2.0f;
  270 + double lambda;
271 271  
272 272 //compute the angles based on NA
273   - float cAngleI = asin(cNAi);
274   - float cAngleO = asin(cNAo);
275   - float oAngleI = asin(oNAi);
276   - float oAngleO = asin(oNAo);
  273 + double cAngleI = asin(cNAi);
  274 + double cAngleO = asin(cNAo);
  275 + double oAngleI = asin(oNAi);
  276 + double oAngleO = asin(oNAo);
277 277  
278 278 //implement a reflection-mode system if necessary
279 279 if(opticsMode == ReflectionOpticsType){
... ... @@ -288,33 +288,33 @@ void pointSpectrum()
288 288 }
289 289  
290 290 //integrate the incident field at the detector position
291   - complex<float> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
292   - float I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
  291 + complex<double> Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
  292 + double I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag();
293 293 I0 *= scaleI0;
294 294  
295 295  
296 296  
297   - float I;
  297 + //double I;
298 298 SpecPair temp;
299   - float nu;
300   - complex<float> eta;
301   - complex<float> Us, U;
  299 + double nu;
  300 + complex<double> eta;
  301 + complex<double> Us, U;
302 302  
303   - float vecLen = 0.0;
304   - for(int i=0; i<EtaK.size(); i++)
  303 + double vecLen = 0.0;
  304 + for(unsigned int i=0; i<EtaK.size(); i++)
305 305 {
306 306 nu = EtaK[i].nu;
307 307 lambda = 10000.0f/nu;
308 308 if(applyMaterial)
309   - eta = complex<float>(EtaN[i].A, EtaK[i].A);
  309 + eta = complex<double>(EtaN[i].A, EtaK[i].A);
310 310 else
311   - eta = complex<float>(baseIR, 0.0);
  311 + eta = complex<double>(baseIR, 0.0);
312 312  
313 313  
314 314 //integrate the scattered field at the detector position
315 315 Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
316 316 U = Us + Ui;
317   - float I = U.real() * U.real() + U.imag() * U.imag();
  317 + double I = U.real() * U.real() + U.imag() * U.imag();
318 318  
319 319 temp.nu = nu;
320 320  
... ... @@ -332,202 +332,13 @@ void pointSpectrum()
332 332 vecLen = sqrt(vecLen);
333 333  
334 334 if(dispNormalize)
335   - for(int i=0; i<SimSpectrum.size(); i++)
  335 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
336 336 SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;
337 337  
338 338 PD.EndTimer(SIMULATE_SPECTRUM);
339 339 }
340 340  
341   -/*
342   -complex<float> sampleUs(complex<float>* B, float* Alpha, int Nl, float r,
343   - float cAngleI, float cAngleO, float theta, float M = 2*pi)
344   -{
345   - /*This function takes a point sample of the scattered field in the far zone
346   - in order to evaluate the field at the central pixel of a detector.
347   - r = sphere radius
348   - lambda = wavelength
349   - eta = index of refraction
350   - cNAi = condenser inner NA
351   - cNAo = condenser outer NA
352   - theta = angle of the sample
353   - M = field magnitude*/
354   -
355   -/*
356   - //compute the material coefficients B
357   - //complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
358   -
359   - //compute the Legendre polynomials for theta (at the objective)
360   - float* Ptheta = (float*)malloc(sizeof(float)*(Nl+1));
361   - Legendre(Ptheta, cos(theta), Nl+1);
362   -
363   - //complex<double> IR(eta.real(), eta.imag());
364   -
365   - //aperature terms for the condenser (alpha) and objective (beta)
366   - float beta;
367   - float c;
368   - complex<float> Us(0.0, 0.0);
369   -
370   - for(int l=0; l<Nl; l++)
371   - {
372   -
373   - complex<float> numer(0.0, -(l*pi)/2.0);
374   - Us += B[l] * exp(numer) * Ptheta[l] * Alpha[l] * pow(complex<float>(0.0, 1.0), l);
375   -
376   -
377   - }
378   - //printf("Ptheta: %f\n", Ptheta[Nl-1]);
379   -
380   - return Us;
381   -
382   -}*/
383   -
384   -
385   -/*
386   -void detectorSpectrum(int numSamples)
387   -{
388   - //integrate across the objective aperature and calculate the resulting intensity on a detector
389   - PD.StartTimer(SIMULATE_SPECTRUM);
390   - //clear the previous spectrum
391   - SimSpectrum.clear();
392   -
393   - float dNu = 2.0f;
394   - float lambda;
395   -
396   - //compute the angles based on NA
397   - float cAngleI = asin(cNAi);
398   - float cAngleO = asin(cNAo);
399   - float oAngleI = asin(oNAi);
400   - float oAngleO = asin(oNAo);
401   -
402   - //implement a reflection-mode system if necessary
403   - if(opticsMode == ReflectionOpticsType){
404   -
405   - //set the condenser to match the objective
406   - cAngleI = oAngleI;
407   - cAngleO = oAngleO;
408   -
409   - //invert the objective
410   - oAngleO = pi - cAngleI;
411   - oAngleI = pi - cAngleO;
412   - }
413   -
414   - //compute Nl (maximum order of the spectrum)
415   - //****************************************************************************
416   - float maxNu = EtaK.back().nu;
417   - float maxLambda = 10000.0f/maxNu;
418   - float k = 2*pi/maxLambda;
419   - int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );
420   - int nLambda = EtaK.size();
421   -
422   - //compute alpha (condenser integral)
423   - //****************************************************************************
424   - //compute the Legendre polynomials for the condenser aperature
425   - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1));
426   - Legendre(PcNAo, cos(cAngleO), Nl+1);
427   - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1));
428   - Legendre(PcNAi, cos(cAngleI), Nl+1);
429   -
430   - //allocate space for the alpha array
431   - float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
432   - computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO);
433   -
434   - for(int l=0; l<Nl; l++)
435   - {
436   - //integration term
437   - if(l == 0)
438   - alpha[l] = -PcNAo[l+1] + PcNAo[0] + PcNAi[l+1] - PcNAi[0];
439   - else
440   - alpha[l] = -PcNAo[l+1] + PcNAo[l-1] + PcNAi[l+1] - PcNAi[l-1];
441   -
442   - alpha[l] *= 2 * pi;
443   - }
444   -
445   - //compute the information based on wavelength and sphere
446   -
447   - //evaluate the incident field intensity
448   - float I0 = 0.0;
449   - float theta;
450   - float dTheta = (oAngleO - oAngleI)/numSamples;
451   - complex<float> Ui;
452   -
453   - Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
454   - I0 = Ui.real()*2*pi;
455   -
456   - float I;
457   - SpecPair temp;
458   - float nu;
459   - complex<float> eta;
460   - complex<float> Us, U;
461   -
462   -
463   -
464   - float vecLen = 0.0;
465   - for(int i=0; i<EtaK.size(); i++)
466   - {
467   - //compute information based on wavelength and material
468   - nu = EtaK[i].nu;
469   - lambda = 10000.0f/nu;
470   - if(applyMaterial)
471   - eta = complex<float>(EtaN[i].A, EtaK[i].A);
472   - else
473   - eta = complex<float>(baseIR, 0.0);
474   -
475   - //allocate memory for the scattering coefficients
476   - complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
477   -
478   - complex<double> IR(eta.real(), eta.imag());
479   - computeB(B, radius, IR, lambda, Nl);
480   -
481   -
482   - //integrate the scattered field at the detector position
483   - I = 0.0;
484   -
485   - for(int iTheta = 0; iTheta < numSamples; iTheta++)
486   - {
487   - theta = oAngleI + iTheta * dTheta;
488   - Us = sampleUs(B, alpha, Nl, radius, cAngleI, cAngleO, theta, 2*pi);
489   -
490   -
491   - //calculate the intensity and add
492   - if(theta >= cAngleI && theta <= cAngleO)
493   - U = Us + 2*(float)pi;
494   - else
495   - U = Us;
496   -
497   - I += (U.real() * U.real() + U.imag() * U.imag()) * sin(theta) * 2 * pi * dTheta;
498   -
499   - }
500   -
501   - temp.nu = nu;
502   -
503   - if(i == 0)
504   - printf("I: %f I0: %f\n", I, I0);
505   -
506   - //set the spectrum value based on the current display type
507   - if(dispSimType == AbsorbanceSpecType)
508   - temp.A = -log10(I/I0);
509   - else
510   - temp.A = I;
511   -
512   - if(dispNormalize)
513   - vecLen += temp.A * temp.A;
514   - //temp.A = Us.real();
515   - SimSpectrum.push_back(temp);
516   -
517   - free(B);
518   - }
519   - vecLen = sqrt(vecLen);
520   -
521   - if(dispNormalize)
522   - for(int i=0; i<SimSpectrum.size(); i++)
523   - SimSpectrum[i].A = (SimSpectrum[i].A / vecLen) * dispNormFactor;
524   -
525   - free(alpha);
526   -
527   - PD.EndTimer(SIMULATE_SPECTRUM);
528   -}*/
529   -
530   -void updateSpectrum(float* I, float I0, int n)
  341 +void updateSpectrum(double* I, double I0, int n)
531 342 {
532 343 SimSpectrum.clear();
533 344 SpecPair temp;
... ... @@ -547,7 +358,7 @@ void updateSpectrum(float* I, float I0, int n)
547 358 }
548 359 }
549 360  
550   -void computeCassegrainAngles(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO)
  361 +void computeCassegrainAngles(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO)
551 362 {
552 363 //compute the angles based on NA
553 364 cAngleI = asin(cNAi);
... ... @@ -572,30 +383,30 @@ void computeCassegrainAngles(float&amp; cAngleI, float&amp; cAngleO, float&amp; oAngleI, flo
572 383  
573 384 int computeNl()
574 385 {
575   - float maxNu = EtaK.back().nu;
576   - float maxLambda = 10000.0f/maxNu;
577   - float k = 2*pi/maxLambda;
578   - int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 );
  386 + double maxNu = EtaK.back().nu;
  387 + double maxLambda = 10000.0f/maxNu;
  388 + double k = 2*pi/maxLambda;
  389 + int Nl = (int)ceil( k + 4 * exp(log(k*radius)/3) + 3 );
579 390  
580 391 return Nl;
581 392 }
582 393  
583   -void computeBArray(complex<float>* B, int Nl, int nLambda)
  394 +void computeBArray(complex<double>* B, int Nl, int nLambda)
584 395 {
585   - float nu;
586   - complex<float> eta;
587   - float* Lambda = (float*)malloc(sizeof(float) * nLambda);
  396 + double nu;
  397 + complex<double> eta;
  398 + double* Lambda = (double*)malloc(sizeof(double) * nLambda);
588 399  
589 400 //for each wavenumber nu
590   - for(int i=0; i<EtaK.size(); i++)
  401 + for(unsigned int i=0; i<EtaK.size(); i++)
591 402 {
592 403 //compute information based on wavelength and material
593 404 nu = EtaK[i].nu;
594 405 Lambda[i] = 10000.0f/nu;
595 406 if(applyMaterial)
596   - eta = complex<float>(EtaN[i].A, EtaK[i].A);
  407 + eta = complex<double>(EtaN[i].A, EtaK[i].A);
597 408 else
598   - eta = complex<float>(baseIR, 0.0);
  409 + eta = complex<double>(baseIR, 0.0);
599 410  
600 411 //allocate memory for the scattering coefficients
601 412 //complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>)*Nl);
... ... @@ -605,14 +416,13 @@ void computeBArray(complex&lt;float&gt;* B, int Nl, int nLambda)
605 416 }
606 417 }
607 418  
608   -void computeOpticalParameters(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO, float& I0, float* alpha, int Nl)
  419 +void computeOpticalParameters(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO, double& I0, double* alpha, int Nl)
609 420 {
610 421 computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO);
611 422  
612 423 //evaluate the incident field intensity
613 424 I0 = 0.0;
614   - float theta;
615   - complex<float> Ui;
  425 + complex<double> Ui;
616 426  
617 427 Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi);
618 428 I0 = Ui.real()*2*pi;
... ... @@ -631,24 +441,24 @@ void gpuDetectorSpectrum(int numSamples)
631 441 //compute Nl (maximum order of the spectrum)
632 442 int Nl = computeNl();
633 443  
634   - float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
635   - float cAngleI, cAngleO, oAngleI, oAngleO, I0;
  444 + double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
  445 + double cAngleI, cAngleO, oAngleI, oAngleO, I0;
636 446 computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
637 447  
638 448 //allocate space for a list of wavelengths
639 449 int nLambda = EtaK.size();
640 450  
641 451 //allocate space for the 2D array (Nl x nu) of scattering coefficients
642   - complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
  452 + complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
643 453 computeBArray(B, Nl, nLambda);
644 454  
645 455  
646 456 //allocate temporary space for the spectrum
647   - float* I = (float*)malloc(sizeof(float) * EtaK.size());
  457 + double* I = (double*)malloc(sizeof(double) * EtaK.size());
648 458  
649 459 //compute the spectrum on the GPU
650 460 PD.StartTimer(SIMULATE_GPU);
651   - cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
  461 + cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples);
652 462 PD.EndTimer(SIMULATE_GPU);
653 463  
654 464 updateSpectrum(I, I0, nLambda);
... ... @@ -666,44 +476,44 @@ void SimulateSpectrum()
666 476 //detectorSpectrum(objectiveSamples);
667 477 }
668 478  
669   -float absorbanceDistortion(){
  479 +double absorbanceDistortion(){
670 480  
671 481 //compute the mean of the spectrum
672   - float sumSim = 0.0;
673   - for(int i=0; i<SimSpectrum.size(); i++)
  482 + double sumSim = 0.0;
  483 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
674 484 {
675 485 sumSim += SimSpectrum[i].A;
676 486 }
677   - float meanSim = sumSim/SimSpectrum.size();
  487 + double meanSim = sumSim/SimSpectrum.size();
678 488  
679 489 //compute the distortion (MSE from the mean)
680   - float sumSE = 0.0;
681   - for(int i=0; i<SimSpectrum.size(); i++)
  490 + double sumSE = 0.0;
  491 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
682 492 {
683 493 sumSE += pow(SimSpectrum[i].A - meanSim, 2);
684 494 }
685   - float MSE = sumSE / SimSpectrum.size();
  495 + double MSE = sumSE / SimSpectrum.size();
686 496  
687 497 return MSE;
688 498 }
689 499  
690   -float intensityDistortion(){
  500 +double intensityDistortion(){
691 501  
692 502 //compute the magnitude of the spectrum
693   - float sumSim = 0.0;
694   - for(int i=0; i<SimSpectrum.size(); i++)
  503 + double sumSim = 0.0;
  504 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
695 505 {
696 506 sumSim += SimSpectrum[i].A * SimSpectrum[i].A;
697 507 }
698   - float magSim = sqrt(sumSim);
  508 + double magSim = sqrt(sumSim);
699 509  
700 510 //compute the distortion (MSE from the mean)
701   - float sumSE = 0.0;
702   - for(int i=0; i<SimSpectrum.size(); i++)
  511 + double sumSE = 0.0;
  512 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
703 513 {
704 514 sumSE += (SimSpectrum[i].A/magSim) * (1.0/SimSpectrum.size());
705 515 }
706   - float MSE = sumSE;
  516 + double MSE = sumSE;
707 517  
708 518 return MSE;
709 519 }
... ... @@ -712,7 +522,7 @@ void MinimizeDistortion(){
712 522 ofstream outFile("distortion.txt");
713 523  
714 524 //set the parameters for the distortion simulation
715   - float step = 0.001;
  525 + double step = 0.001;
716 526  
717 527 oNAi = 0.2;
718 528 oNAo = 0.5;
... ... @@ -721,28 +531,28 @@ void MinimizeDistortion(){
721 531 //compute Nl (maximum order of the spectrum)
722 532 int Nl = computeNl();
723 533  
724   - float* alpha = (float*)malloc(sizeof(float)*(Nl + 1));
725   - float cAngleI, cAngleO, oAngleI, oAngleO, I0;
  534 + double* alpha = (double*)malloc(sizeof(double)*(Nl + 1));
  535 + double cAngleI, cAngleO, oAngleI, oAngleO, I0;
726 536  
727 537 //allocate space for a list of wavelengths
728 538 int nLambda = EtaK.size();
729 539  
730 540 //allocate temporary space for the spectrum
731   - float* I = (float*)malloc(sizeof(float) * EtaK.size());
  541 + double* I = (double*)malloc(sizeof(double) * EtaK.size());
732 542  
733 543 //calculate the material parameters
734 544 //allocate space for the 2D array (Nl x nu) of scattering coefficients
735   - complex<float>* B = (complex<float>*)malloc(sizeof(complex<float>) * Nl * nLambda);
  545 + complex<double>* B = (complex<double>*)malloc(sizeof(complex<double>) * Nl * nLambda);
736 546 computeBArray(B, Nl, nLambda);
737 547  
738 548  
739 549  
740   - float D;
741   - float e = 0.001;
742   - for(float i=0.0; i<=oNAo-step; i+=step)
  550 + double D;
  551 + double e = 0.001;
  552 + for(double i=0.0; i<=oNAo-step; i+=step)
743 553 {
744 554  
745   - for(float o=oNAi+step; o<=1.0; o+=step)
  555 + for(double o=oNAi+step; o<=1.0; o+=step)
746 556 {
747 557  
748 558  
... ... @@ -754,7 +564,7 @@ void MinimizeDistortion(){
754 564 computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl);
755 565  
756 566 //simulate the spectrum
757   - cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);
  567 + cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples);
758 568 updateSpectrum(I, I0, nLambda);
759 569  
760 570 if(dispSimType == AbsorbanceSpecType)
... ...
cudaKK.h
1   -__device__ float g(float v0, float v1)
  1 +__device__ double g(double v0, double v1)
2 2 {
3 3 return (v0 + v1)*log(abs(v0+v1)) + (v0-v1)*log(abs(v0-v1));
4 4 }
5 5  
6   -__device__ float hfin(float v0, float v1, float dv)
  6 +__device__ double hfin(double v0, double v1, double dv)
7 7 {
8   - float e = 0.001;
9   - float t0 = g(v0+e, v1-dv)/dv;
10   - float t1 = 2*g(v0+e, v1)/dv;
11   - float t2 = g(v0+e, v1+dv)/dv;
  8 + double e = 0.001;
  9 + double t0 = g(v0+e, v1-dv)/dv;
  10 + double t1 = 2*g(v0+e, v1)/dv;
  11 + double t2 = g(v0+e, v1+dv)/dv;
12 12  
13 13 return -1.0/PI * (t0 - t1 + t2);
14 14 }
15 15  
16   -__global__ void devKramersKronig(float* gpuN, float* gpuK, int numVals, float nuStart, float nuEnd, float nOffset)
  16 +__global__ void devKramersKronig(double* gpuN, double* gpuK, int numVals, double nuStart, double nuEnd, double nOffset)
17 17 {
18 18 int i = blockIdx.x * blockDim.x + threadIdx.x;
19 19  
20 20 if(i >= numVals) return;
21   - float nuDelta = (nuEnd - nuStart)/(numVals - 1);
  21 + double nuDelta = (nuEnd - nuStart)/(numVals - 1);
22 22  
23   - float nu = nuStart + i*nuDelta;
24   - float n = 0.0;
25   - float jNu;
26   - float jK;
  23 + double nu = nuStart + i*nuDelta;
  24 + double n = 0.0;
  25 + double jNu;
  26 + double jK;
27 27 for(int j=1; j<numVals-1; j++)
28 28 {
29 29 jNu = nuStart + j*nuDelta;
... ... @@ -35,55 +35,55 @@ __global__ void devKramersKronig(float* gpuN, float* gpuK, int numVals, float nu
35 35  
36 36 }
37 37  
38   -void cudaKramersKronig(float* cpuN, float* cpuK, int nVals, float nuStart, float nuEnd, float nOffset)
  38 +void cudaKramersKronig(double* cpuN, double* cpuK, int nVals, double nuStart, double nuEnd, double nOffset)
39 39 {
40   - float* gpuK;
41   - HANDLE_ERROR(cudaMalloc(&gpuK, sizeof(float)*nVals));
42   - HANDLE_ERROR(cudaMemcpy(gpuK, cpuK, sizeof(float)*nVals, cudaMemcpyHostToDevice));
43   - float* gpuN;
44   - HANDLE_ERROR(cudaMalloc(&gpuN, sizeof(float)*nVals));
  40 + double* gpuK;
  41 + HANDLE_ERROR(cudaMalloc(&gpuK, sizeof(double)*nVals));
  42 + HANDLE_ERROR(cudaMemcpy(gpuK, cpuK, sizeof(double)*nVals, cudaMemcpyHostToDevice));
  43 + double* gpuN;
  44 + HANDLE_ERROR(cudaMalloc(&gpuN, sizeof(double)*nVals));
45 45  
46 46 dim3 block(BLOCK_SIZE*BLOCK_SIZE);
47 47 dim3 grid(nVals/block.x + 1);
48 48 devKramersKronig<<<grid, block>>>(gpuN, gpuK, nVals, nuStart, nuEnd, nOffset);
49 49  
50   - HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(float)*nVals, cudaMemcpyDeviceToHost));
  50 + HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(double)*nVals, cudaMemcpyDeviceToHost));
51 51  
52 52 //free resources
53 53 HANDLE_ERROR(cudaFree(gpuK));
54 54 HANDLE_ERROR(cudaFree(gpuN));
55 55 }
56 56  
57   -__global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl,
58   - int nSamples, float oThetaI, float oThetaO, float cThetaI, float cThetaO)
  57 +__global__ void devComputeSpectrum(double* I, double2* B, double* alpha, int Nl,
  58 + int nSamples, double oThetaI, double oThetaO, double cThetaI, double cThetaO)
59 59 {
60 60 int i = blockIdx.x * blockDim.x + threadIdx.x;
61 61  
62 62 //compute the delta-theta value
63   - float dTheta = (oThetaO - oThetaI)/nSamples;
  63 + double dTheta = (oThetaO - oThetaI)/nSamples;
64 64  
65 65 //allocate space for the Legendre polynomials
66   - float Ptheta[2];
67   -
68   - float cosTheta, theta;
69   - cuComplex Us;
70   - cuComplex UsSample;
71   - cuComplex U;
72   - cuComplex Ui;
73   - Ui.x = 2*PI;
74   - Ui.y = 0.0;
75   - cuComplex numer;
  66 + double Ptheta[2];
  67 +
  68 + double cosTheta, theta;
  69 + cuDoubleComplex Us;
  70 + cuDoubleComplex UsSample;
  71 + cuDoubleComplex U;
  72 + //cuComplex Ui;
  73 + //Ui.x = 2*PI;
  74 + //Ui.y = 0.0;
  75 + cuDoubleComplex numer;
76 76 numer.x = 0.0;
77   - cuComplex exp_numer;
78   - cuComplex iL;
79   - cuComplex imag;
  77 + cuDoubleComplex exp_numer;
  78 + cuDoubleComplex iL;
  79 + cuDoubleComplex imag;
80 80 imag.x = 0.0; imag.y = 1.0;
81   - float realFac;
82   - cuComplex complexFac;
83   - float PlTheta;
84   - float Isum = 0.0;
85   - float maxVal = 0;
86   - float val;
  81 + double realFac;
  82 + cuDoubleComplex complexFac;
  83 + double PlTheta;
  84 + double Isum = 0.0;
  85 + //float maxVal = 0;
  86 + //float val;
87 87 for(int iTheta = 0; iTheta < nSamples; iTheta++)
88 88 {
89 89 //calculate theta
... ... @@ -131,12 +131,7 @@ __global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl,
131 131 //increment the imaginary exponent i^l
132 132 iL = cMult(iL, imag);
133 133  
134   - //val = cMag(Us);
135   - //if(val > maxVal)
136   - // maxVal = val;
137   -
138   - //Us += B[l] * exp(numer) * Ptheta[l] * alpha * M * pow(complex<float>(0.0, 1.0), l);
139   -
  134 +
140 135 }
141 136  
142 137 //sum the scattered and incident fields
... ... @@ -150,20 +145,20 @@ __global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl,
150 145 I[i] = Isum;
151 146 }
152 147  
153   -void cudaComputeSpectrum(float* cpuI, float* cpuB, float* cpuAlpha,
154   - int Nl, int nLambda, float oThetaI, float oThetaO, float cThetaI, float cThetaO, int nSamples)
  148 +void cudaComputeSpectrum(double* cpuI, double* cpuB, double* cpuAlpha,
  149 + int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples)
155 150 {
156 151 //copy everything to the GPU
157   - float2* gpuB;
158   - HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(float2) * nLambda * Nl));
159   - HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(float2) * nLambda * Nl, cudaMemcpyHostToDevice));
  152 + double2* gpuB;
  153 + HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(double2) * nLambda * Nl));
  154 + HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(double2) * nLambda * Nl, cudaMemcpyHostToDevice));
160 155  
161   - float* gpuAlpha;
162   - HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(float) * Nl));
163   - HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(float) * Nl, cudaMemcpyHostToDevice));
  156 + double* gpuAlpha;
  157 + HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(double) * Nl));
  158 + HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(double) * Nl, cudaMemcpyHostToDevice));
164 159  
165   - float* gpuI;
166   - HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(float) * nLambda));
  160 + double* gpuI;
  161 + HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(double) * nLambda));
167 162  
168 163  
169 164 //call the kernel to compute the spectrum
... ... @@ -171,10 +166,10 @@ void cudaComputeSpectrum(float* cpuI, float* cpuB, float* cpuAlpha,
171 166 dim3 grid(nLambda/block.x + 1);
172 167  
173 168 //devComputeSpectrum
174   - devComputeSpectrum<<<grid, block>>>(gpuI, (float2*)gpuB, gpuAlpha, Nl,
  169 + devComputeSpectrum<<<grid, block>>>(gpuI, (double2*)gpuB, gpuAlpha, Nl,
175 170 nSamples, oThetaI, oThetaO, cThetaI, cThetaO);
176 171  
177   - HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(float) * nLambda, cudaMemcpyDeviceToHost));
  172 + HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(double) * nLambda, cudaMemcpyDeviceToHost));
178 173  
179 174 HANDLE_ERROR(cudaFree(gpuB));
180 175 HANDLE_ERROR(cudaFree(gpuAlpha));
... ...
cudaMain.cu
... ... @@ -5,45 +5,45 @@
5 5 #define PI 3.14159
6 6 #define BLOCK_SIZE 16
7 7  
8   -__device__ cuComplex cMult(cuComplex a, cuComplex b)
  8 +__device__ cuDoubleComplex cMult(cuDoubleComplex a, cuDoubleComplex b)
9 9 {
10   - cuComplex result;
  10 + cuDoubleComplex result;
11 11 result.x = a.x * b.x - a.y * b.y;
12 12 result.y = a.x * b.y + a.y * b.x;
13 13  
14 14 return result;
15 15 }
16 16  
17   -__device__ cuComplex cMult(cuComplex a, float b)
  17 +__device__ cuDoubleComplex cMult(cuDoubleComplex a, float b)
18 18 {
19   - cuComplex result;
  19 + cuDoubleComplex result;
20 20 result.x = a.x * b;
21 21 result.y = a.y * b;
22 22  
23 23 return result;
24 24 }
25 25  
26   -__device__ cuComplex cAdd(cuComplex a, cuComplex b)
  26 +__device__ cuDoubleComplex cAdd(cuDoubleComplex a, cuDoubleComplex b)
27 27 {
28   - cuComplex r;
  28 + cuDoubleComplex r;
29 29 r.x = a.x + b.x;
30 30 r.y = a.y + b.y;
31 31  
32 32 return r;
33 33 }
34 34  
35   -__device__ cuComplex cAdd(cuComplex a, float b)
  35 +__device__ cuDoubleComplex cAdd(cuDoubleComplex a, float b)
36 36 {
37   - cuComplex r;
  37 + cuDoubleComplex r;
38 38 r.x = a.x + b;
39 39 r.y = a.y;
40 40  
41 41 return r;
42 42 }
43 43  
44   -__device__ cuComplex cExp(cuComplex a)
  44 +__device__ cuDoubleComplex cExp(cuDoubleComplex a)
45 45 {
46   - cuComplex r;
  46 + cuDoubleComplex r;
47 47  
48 48 r.x = exp(a.x) * cos(a.y);
49 49 r.y = exp(a.x) * sin(a.y);
... ... @@ -51,9 +51,9 @@ __device__ cuComplex cExp(cuComplex a)
51 51 return r;
52 52 }
53 53  
54   -__device__ float cMag(cuComplex a)
  54 +__device__ double cMag(cuDoubleComplex a)
55 55 {
56   - float r = sqrt(a.x * a.x + a.y * a.y);
  56 + double r = sqrt(a.x * a.x + a.y * a.y);
57 57 return r;
58 58 }
59 59  
... ...
globals.h
... ... @@ -8,14 +8,16 @@
8 8 #include <complex>
9 9 using namespace std;
10 10  
  11 +typedef float rtsFloat;
  12 +
11 13 struct SpecPair{
12   - float nu;
13   - float A;
  14 + double nu;
  15 + double A;
14 16 };
15 17  
16 18 struct Material{
17   - vector<float> nu;
18   - vector<complex<float>> eta;
  19 + vector<double> nu;
  20 + vector<complex<double>> eta;
19 21 string name;
20 22 };
21 23  
... ... @@ -47,16 +49,16 @@ void FitDisplay();
47 49 //Update Functions
48 50 void UpdateDisplay();
49 51 void SimulateSpectrum();
50   -void cudaKramersKronig(float* cpuN, float* cpuK, int nVals, float nuStart, float nuEnd, float nOffset);
51   -void cudaComputeSpectrum(float* cpuI, float* cpuB, float* alpha,
52   - int Nl, int nLambda, float oThetaI, float oThetaO, float cThetaI, float cThetaO, int nSamples);
  52 +void cudaKramersKronig(double* cpuN, double* cpuK, int nVals, double nuStart, double nuEnd, double nOffset);
  53 +void cudaComputeSpectrum(double* cpuI, double* cpuB, double* alpha,
  54 + int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples);
53 55  
54 56 //Window Parameters
55   -extern float nuMin;
56   -extern float nuMax;
57   -extern float aMin;
58   -extern float aMax;
59   -extern float dNu;
  57 +extern double nuMin;
  58 +extern double nuMax;
  59 +extern double aMin;
  60 +extern double aMax;
  61 +extern double dNu;
60 62 extern bool dispRefSpec;
61 63 extern bool dispSimSpec;
62 64 extern bool dispSimK;
... ... @@ -65,16 +67,16 @@ extern bool dispSimN;
65 67 extern bool dispMatN;
66 68 extern SpecType dispSimType;
67 69 extern bool dispNormalize;
68   -extern float dispNormFactor;
  70 +extern double dispNormFactor;
69 71  
70 72  
71   -extern float dispScaleK;
72   -extern float dispScaleN;
  73 +extern double dispScaleK;
  74 +extern double dispScaleN;
73 75  
74 76 //material parameters
75   -extern float radius;
76   -extern float baseIR;
77   -extern float cA;
  77 +extern double radius;
  78 +extern double baseIR;
  79 +extern double cA;
78 80 extern vector<SpecPair> EtaK;
79 81 extern vector<SpecPair> EtaN;
80 82 extern bool applyMaterial;
... ... @@ -86,22 +88,22 @@ void ChangeAbsorbance();
86 88 void SetMaterial();
87 89  
88 90 //optical parameters
89   -extern float cNAi;
90   -extern float cNAo;
91   -extern float oNAi;
92   -extern float oNAo;
  91 +extern double cNAi;
  92 +extern double cNAo;
  93 +extern double oNAi;
  94 +extern double oNAo;
93 95 extern OpticsType opticsMode;
94 96 extern bool pointDetector;
95 97 extern int objectiveSamples;
96 98  
97 99 //fitting parameters
98   -extern float minMSE;
  100 +extern double minMSE;
99 101 extern int maxFitIter;
100 102 void EstimateMaterial();
101   -extern float scaleI0;
102   -extern float refSlope;
  103 +extern double scaleI0;
  104 +extern double refSlope;
103 105  
104   -float ComputeDistortion();
  106 +double ComputeDistortion();
105 107 void MinimizeDistortion();
106 108  
107 109  
... ...
interactivemie.h
... ... @@ -48,7 +48,7 @@ public:
48 48  
49 49 //material selection combo box
50 50 ui.cmbMaterial->clear();
51   - for(int i=0; i<MaterialList.size(); i++)
  51 + for(unsigned int i=0; i<MaterialList.size(); i++)
52 52 ui.cmbMaterial->addItem(MaterialList[i].name.c_str(), i);
53 53 ui.cmbMaterial->setCurrentIndex(currentMaterial);
54 54  
... ...
interactivemie.ui
... ... @@ -350,6 +350,9 @@
350 350 <height>22</height>
351 351 </rect>
352 352 </property>
  353 + <property name="decimals">
  354 + <number>6</number>
  355 + </property>
353 356 <property name="minimum">
354 357 <double>-9999.000000000000000</double>
355 358 </property>
... ... @@ -362,7 +365,7 @@
362 365 </widget>
363 366 <widget class="QSpinBox" name="spinNuMin">
364 367 <property name="enabled">
365   - <bool>false</bool>
  368 + <bool>true</bool>
366 369 </property>
367 370 <property name="geometry">
368 371 <rect>
... ... @@ -372,6 +375,9 @@
372 375 <height>22</height>
373 376 </rect>
374 377 </property>
  378 + <property name="minimum">
  379 + <number>1</number>
  380 + </property>
375 381 <property name="maximum">
376 382 <number>4000</number>
377 383 </property>
... ... @@ -404,6 +410,9 @@
404 410 <height>22</height>
405 411 </rect>
406 412 </property>
  413 + <property name="decimals">
  414 + <number>6</number>
  415 + </property>
407 416 <property name="minimum">
408 417 <double>-9999.000000000000000</double>
409 418 </property>
... ... @@ -419,7 +428,7 @@
419 428 </widget>
420 429 <widget class="QSpinBox" name="spinNuMax">
421 430 <property name="enabled">
422   - <bool>false</bool>
  431 + <bool>true</bool>
423 432 </property>
424 433 <property name="geometry">
425 434 <rect>
... ... @@ -429,8 +438,11 @@
429 438 <height>22</height>
430 439 </rect>
431 440 </property>
  441 + <property name="minimum">
  442 + <number>1</number>
  443 + </property>
432 444 <property name="maximum">
433   - <number>4000</number>
  445 + <number>100000</number>
434 446 </property>
435 447 <property name="singleStep">
436 448 <number>10</number>
... ... @@ -509,6 +521,9 @@
509 521 <height>22</height>
510 522 </rect>
511 523 </property>
  524 + <property name="decimals">
  525 + <number>6</number>
  526 + </property>
512 527 <property name="minimum">
513 528 <double>-99.989999999999995</double>
514 529 </property>
... ...
main.cpp
... ... @@ -19,15 +19,15 @@ vector&lt;SpecPair&gt; EtaK;
19 19 vector<SpecPair> EtaN;
20 20 int currentSpec = 0;
21 21  
22   -float nuMin = 800;
23   -float nuMax = 4000;
24   -float dNu = 2;
  22 +double nuMin = 800;
  23 +double nuMax = 4000;
  24 +double dNu = 2;
25 25  
26   -float aMin = 0;
27   -float aMax = 1;
  26 +double aMin = 0;
  27 +double aMax = 1;
28 28  
29   -float scaleI0 = 1.0;
30   -float refSlope = 0.0;
  29 +double scaleI0 = 1.0;
  30 +double refSlope = 0.0;
31 31  
32 32 bool dispRefSpec = true;
33 33 bool dispSimSpec = true;
... ... @@ -35,17 +35,17 @@ bool dispSimK = true;
35 35 bool dispMatK = true;
36 36 bool dispSimN = true;
37 37 bool dispMatN = true;
38   -float dispScaleK = 1.0;
39   -float dispScaleN = 1.0;
  38 +double dispScaleK = 1.0;
  39 +double dispScaleN = 1.0;
40 40 SpecType dispSimType = AbsorbanceSpecType;
41 41 bool dispNormalize = false;
42   -float dispNormFactor = 1.0;
  42 +double dispNormFactor = 1.0;
43 43  
44 44  
45 45 //material parameters
46   -float radius = 4.0f;
47   -float baseIR = 1.49f;
48   -float cA = 1.0;
  46 +double radius = 4.0f;
  47 +double baseIR = 1.49f;
  48 +double cA = 1.0;
49 49 //vector<SpecPair> KMaterial;
50 50 //vector<SpecPair> NMaterial;
51 51 bool applyMaterial = true;
... ... @@ -53,16 +53,16 @@ vector&lt;Material&gt; MaterialList;
53 53 int currentMaterial = 0;
54 54  
55 55 //optical parameters
56   -float cNAi = 0.0;
57   -float cNAo = 0.6;
58   -float oNAi = 0.0;
59   -float oNAo = 0.6;
  56 +double cNAi = 0.0;
  57 +double cNAo = 0.6;
  58 +double oNAi = 0.0;
  59 +double oNAo = 0.6;
60 60 OpticsType opticsMode = TransmissionOpticsType;
61 61 bool pointDetector = false;
62 62 int objectiveSamples = 200;
63 63  
64 64 //fitting parameters
65   -float minMSE = 0.00001;
  65 +double minMSE = 0.00001;
66 66 int maxFitIter = 20;
67 67  
68 68 void TempSimSpectrum()
... ... @@ -71,7 +71,7 @@ void TempSimSpectrum()
71 71 for(int i=800; i<4000; i++)
72 72 {
73 73 temp.nu = i;
74   - temp.A = sin((float)i/200);
  74 + temp.A = sin((double)i/200);
75 75 SimSpectrum.push_back(temp);
76 76 }
77 77 }
... ... @@ -94,11 +94,11 @@ void LoadMaterial(string fileNameK, string fileNameN, string materialName)
94 94 exit(1);
95 95 }
96 96  
97   - complex<float> eta;
98   - int j;
99   - for(int i=0; i<KMaterial.size(); i++){
  97 + complex<double> eta;
  98 + //int j;
  99 + for(unsigned int i=0; i<KMaterial.size(); i++){
100 100 newMaterial.nu.push_back(KMaterial[i].nu);
101   - eta = complex<float>(NMaterial[i].A, KMaterial[i].A);
  101 + eta = complex<double>(NMaterial[i].A, KMaterial[i].A);
102 102 newMaterial.eta.push_back(eta);
103 103 }
104 104 MaterialList.push_back(newMaterial);
... ... @@ -114,15 +114,15 @@ void LoadMaterial(string fileNameK, string materialName){
114 114  
115 115 //compute the real IR using Kramers Kronig
116 116 //copy the absorbance values into a linear array
117   - float* k = (float*)malloc(sizeof(float) * KMaterial.size());
118   - float* n = (float*)malloc(sizeof(float) * KMaterial.size());
119   - for(int i=0; i<KMaterial.size(); i++)
  117 + double* k = (double*)malloc(sizeof(double) * KMaterial.size());
  118 + double* n = (double*)malloc(sizeof(double) * KMaterial.size());
  119 + for(unsigned int i=0; i<KMaterial.size(); i++)
120 120 k[i] = KMaterial[i].A;
121 121  
122 122 //use Kramers Kronig to determine the real part of the index of refraction
123 123 cudaKramersKronig(n, k, KMaterial.size(), KMaterial[0].nu, KMaterial.back().nu, baseIR);
124 124 SpecPair temp;
125   - for(int i=0; i<KMaterial.size(); i++)
  125 + for(unsigned int i=0; i<KMaterial.size(); i++)
126 126 {
127 127 temp.nu = KMaterial[i].nu;
128 128 temp.A = n[i];
... ... @@ -132,11 +132,10 @@ void LoadMaterial(string fileNameK, string materialName){
132 132 //create the material
133 133 Material newMaterial;
134 134 newMaterial.name = materialName;
135   - complex<float> eta;
136   - int j;
137   - for(int i=0; i<KMaterial.size(); i++){
  135 + complex<double> eta;
  136 + for(unsigned int i=0; i<KMaterial.size(); i++){
138 137 newMaterial.nu.push_back(KMaterial[i].nu);
139   - eta = complex<float>(NMaterial[i].A, KMaterial[i].A);
  138 + eta = complex<double>(NMaterial[i].A, KMaterial[i].A);
140 139 newMaterial.eta.push_back(eta);
141 140 }
142 141  
... ... @@ -144,12 +143,12 @@ void LoadMaterial(string fileNameK, string materialName){
144 143 }
145 144  
146 145 void FitDisplay(){
147   - float minA = 99999.0;
148   - float maxA = -99999.0;
149   - float k, n;
  146 + double minA = 99999.0;
  147 + double maxA = -99999.0;
  148 + double k, n;
150 149  
151 150 if(dispSimSpec)
152   - for(int i=0; i<SimSpectrum.size(); i++)
  151 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
153 152 {
154 153 if(SimSpectrum[i].A < minA)
155 154 minA = SimSpectrum[i].A;
... ... @@ -158,7 +157,7 @@ void FitDisplay(){
158 157 }
159 158  
160 159 if(dispRefSpec && RefSpectrum.size() > 0)
161   - for(int i=0; i<RefSpectrum[currentSpec].size(); i++)
  160 + for(unsigned int i=0; i<RefSpectrum[currentSpec].size(); i++)
162 161 {
163 162 if(RefSpectrum[currentSpec][i].A < minA)
164 163 minA = RefSpectrum[currentSpec][i].A;
... ... @@ -166,7 +165,7 @@ void FitDisplay(){
166 165 maxA = RefSpectrum[currentSpec][i].A;
167 166 }
168 167 if(dispMatK)
169   - for(int i=0; i<EtaK.size(); i++)
  168 + for(unsigned int i=0; i<EtaK.size(); i++)
170 169 {
171 170 k = MaterialList[currentMaterial].eta[i].imag() * dispScaleK;
172 171 if(k < minA)
... ... @@ -175,7 +174,7 @@ void FitDisplay(){
175 174 maxA = k;
176 175 }
177 176 if(dispSimK)
178   - for(int i=0; i<EtaK.size(); i++)
  177 + for(unsigned int i=0; i<EtaK.size(); i++)
179 178 {
180 179 k = EtaK[i].A * dispScaleK;
181 180 if(k < minA)
... ... @@ -184,7 +183,7 @@ void FitDisplay(){
184 183 maxA = k;
185 184 }
186 185 if(dispMatN)
187   - for(int i=0; i<EtaN.size(); i++)
  186 + for(unsigned int i=0; i<EtaN.size(); i++)
188 187 {
189 188 n = (MaterialList[currentMaterial].eta[i].real() - baseIR) * dispScaleN;
190 189 if(n < minA)
... ... @@ -193,7 +192,7 @@ void FitDisplay(){
193 192 maxA = n;
194 193 }
195 194 if(dispSimN)
196   - for(int i=0; i<EtaN.size(); i++)
  195 + for(unsigned int i=0; i<EtaN.size(); i++)
197 196 {
198 197 n = (EtaN[i].A - baseIR) * dispScaleN;
199 198 if(n < minA)
... ... @@ -213,10 +212,10 @@ void ChangeAbsorbance(){
213 212  
214 213 //copy the absorbance values into a linear array
215 214 int nSamples = MaterialList[currentMaterial].eta.size();
216   - float startNu = MaterialList[currentMaterial].nu.front();
217   - float endNu = MaterialList[currentMaterial].nu.back();
218   - float* k = (float*)malloc(sizeof(float) * nSamples);
219   - float* n = (float*)malloc(sizeof(float) * nSamples);
  215 + double startNu = MaterialList[currentMaterial].nu.front();
  216 + double endNu = MaterialList[currentMaterial].nu.back();
  217 + double* k = (double*)malloc(sizeof(double) * nSamples);
  218 + double* n = (double*)malloc(sizeof(double) * nSamples);
220 219 for(int i=0; i<nSamples; i++)
221 220 k[i] = MaterialList[currentMaterial].eta[i].imag() * cA;
222 221  
... ... @@ -230,7 +229,7 @@ void ChangeAbsorbance(){
230 229 SpecPair temp;
231 230  
232 231 //load the imaginary IR from the absorbance data
233   - float nu;
  232 + double nu;
234 233 for(int i=0; i<nSamples; i++){
235 234 nu = MaterialList[currentMaterial].nu[i];
236 235 if(nu >= nuMin && nu <= nuMax){
... ... @@ -253,7 +252,7 @@ void SetMaterial()
253 252 EtaN.clear();
254 253  
255 254 int nSamples = MaterialList[currentMaterial].eta.size();
256   - float nu;
  255 + double nu;
257 256 SpecPair temp;
258 257 for(int i=0; i<nSamples; i++){
259 258 nu = MaterialList[currentMaterial].nu[i];
... ...
qtSpectrumDisplay.cpp
... ... @@ -104,7 +104,7 @@
104 104 {
105 105 glColor3f(1.0, 1.0, 1.0);
106 106 glBegin(GL_LINE_STRIP);
107   - for(int i=0; i<SimSpectrum.size(); i++)
  107 + for(unsigned int i=0; i<SimSpectrum.size(); i++)
108 108 glVertex2f(SimSpectrum[i].nu, SimSpectrum[i].A);
109 109 glEnd();
110 110 }
... ... @@ -113,7 +113,7 @@
113 113 glColor3f(0.5, 0.5, 0.5);
114 114 glBegin(GL_LINE_STRIP);
115 115 float nu;
116   - for(int i=0; i<RefSpectrum[currentSpec].size(); i++)
  116 + for(unsigned int i=0; i<RefSpectrum[currentSpec].size(); i++)
117 117 {
118 118 nu = RefSpectrum[currentSpec][i].nu;
119 119 glVertex2f(nu, RefSpectrum[currentSpec][i].A + nu * refSlope);
... ... @@ -139,7 +139,7 @@
139 139 {
140 140 glColor3f(1.0, 1.0, 0.0);
141 141 glBegin(GL_LINE_STRIP);
142   - for(int i=0; i<EtaK.size(); i++){
  142 + for(unsigned int i=0; i<EtaK.size(); i++){
143 143 glVertex2f(EtaK[i].nu, EtaK[i].A * dispScaleK);
144 144 }
145 145 glEnd();
... ... @@ -160,7 +160,7 @@
160 160 {
161 161 glColor3f(0.0, 1.0, 1.0);
162 162 glBegin(GL_LINE_STRIP);
163   - for(int i=0; i<EtaN.size(); i++)
  163 + for(unsigned int i=0; i<EtaN.size(); i++)
164 164 glVertex2f(EtaN[i].nu, (EtaN[i].A - baseIR) * dispScaleN);
165 165 glEnd();
166 166 }
... ...