diff --git a/BESSIK.CPP b/BESSIK.CPP index af85a13..73dee21 100644 --- a/BESSIK.CPP +++ b/BESSIK.CPP @@ -348,7 +348,7 @@ int bessikv(double v,double x,double &vm,double *iv,double *kv, double *ivp,double *kvp) { double x2,v0,piv,vt,a1,v0p,gap,r,bi0,ca,sum; - double f,f0,f1,f2,ct,cs,wa,gan,ww,w0,v0n; + double f,f1,f2,ct,cs,wa,gan,ww,w0,v0n; double r1,r2,bk0,bk1,bk2,a2,cb; int n,k,kz,m; diff --git a/BESSJY.CPP b/BESSJY.CPP index 4e58505..f1c1c74 100644 --- a/BESSJY.CPP +++ b/BESSJY.CPP @@ -274,7 +274,7 @@ int msta1(double x,int mp) n1 = n0+5; f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-mp; for (i=0;i<20;i++) { - nn = n1-(n1-n0)/(1.0-f0/f1); + nn = (int)(n1-(n1-n0)/(1.0-f0/f1)); f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-mp; if (abs(nn-n1) < 1) break; n0 = n1; @@ -305,7 +305,7 @@ int msta2(double x,int n,int mp) n1 = n0+5; f1 = 0.5*log10(6.28*n1)-n1*log10(1.36*a0/n1)-obj; for (i=0;i<20;i++) { - nn = n1-(n1-n0)/(1.0-f0/f1); + nn = (int)(n1-(n1-n0)/(1.0-f0/f1)); f = 0.5*log10(6.28*nn)-nn*log10(1.36*a0/nn)-obj; if (abs(nn-n1) < 1) break; n0 = n1; @@ -404,7 +404,7 @@ int bessjyna(int n,double x,int &nm,double *jn,double *yn, int bessjynb(int n,double x,int &nm,double *jn,double *yn, double *jnp,double *ynp) { - double t1,t2,f,f0,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv; + double t1,t2,f,f1,f2,bj0,bj1,bjk,by0,by1,cu,s0,su,sv; double ec,bs,byk,p0,p1,q0,q1; static double a[] = { -0.7031250000000000e-1, @@ -536,7 +536,7 @@ int bessjyv(double v,double x,double &vm,double *jv,double *yv, { double v0,vl,vg,vv,a,a0,r,x2,bjv0,bjv1,bjvl,f,f0,f1,f2; double r0,r1,ck,cs,cs0,cs1,sk,qx,px,byv0,byv1,rp,xk,rq; - double b,ec,w0,w1,bjy0,bjy1,bju0,bju1,pv0,pv1,byvk; + double b,ec,w0,w1,bju0,bju1,pv0,pv1,byvk; int j,k,l,m,n,kz; x2 = x*x; diff --git a/CBESSJY.CPP b/CBESSJY.CPP index 103b878..645bc4c 100644 --- a/CBESSJY.CPP +++ b/CBESSJY.CPP @@ -189,7 +189,7 @@ int cbessjyna(int n,complex z,int &nm,complex *cj, complex cs,cg0,cg1,cyk,cyl1,cyl2,cylk,cp11,cp12,cp21,cp22; complex ch0,ch1,ch2; double a0,yak,ya1,ya0,wa; - int m,k,kz,lb,lb0; + int m,k,lb,lb0; if (n < 0) return 1; a0 = abs(z); diff --git a/EstimateMaterial.cpp b/EstimateMaterial.cpp index 5564972..5ddd887 100644 --- a/EstimateMaterial.cpp +++ b/EstimateMaterial.cpp @@ -1,12 +1,12 @@ #include "globals.h" #define PI 3.14159 -float CalculateError(float* E) +double CalculateError(double* E) { //Calculate the error between the Reference Spectrum and the Simulated Spectrum - float sumE = 0.0; + double sumE = 0.0; int nVals = RefSpectrum[currentSpec].size(); - float nu; + double nu; for(int i=0; i minMSE && i < maxFitIter) + while(sumE > minMSE && j < maxFitIter) { //simulate a spectrum based on the current IR SimulateSpectrum(); @@ -86,13 +86,13 @@ void EstimateMaterial() //use Kramers-Kronig to compute n - for(int i=0; i LoadSpectrum(string filename) } //compute the minimum and maximum input wavenumbers - float inMin = S.front().nu; - float inMax = S.back().nu; + double inMin = S.front().nu; + double inMax = S.back().nu; - int nuMin = ceil(inMin); - int nuMax = floor(inMax); + int nuMin = (int)ceil(inMin); + int nuMax = (int)floor(inMax); //make sure both are either even or odd if(nuMin % 2 != nuMax % 2) @@ -39,7 +39,7 @@ vector LoadSpectrum(string filename) //allocate space for the spectrum vector outSpec; - float nu, highVal, lowVal, a; + double nu, highVal, lowVal, a; int j=1; for(int i=0; i LoadSpectrum(string filename) else { //move to the correct position in the input array - while(j < S.size()-1 && S[j].nu <= nu) + while(j < (int)S.size()-1 && S[j].nu <= nu) j++; lowVal = S[j-1].nu; @@ -88,7 +88,7 @@ vector SetReferenceSpectrum(char* text) void SaveK(string fileName) { ofstream outFile(fileName.c_str()); - for(int i=0; i maxTime[type] ) maxTime[type] = t; diff --git a/README b/README index 9a8fc52..0069b94 100644 --- a/README +++ b/README @@ -1,3 +1,5 @@ +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. + TrueEyes requires the following libraries: GLUT (OpenGL Utility Toolkit) http://www.opengl.org/resources/libraries/glut/ diff --git a/SimulateSpectrum.cpp b/SimulateSpectrum.cpp index 991807b..6948811 100644 --- a/SimulateSpectrum.cpp +++ b/SimulateSpectrum.cpp @@ -8,7 +8,7 @@ using namespace std; #define pi 3.14159 -typedef complex scComplex; +typedef complex scComplex; int cbessjyva(double v,complex z,double &vm,complex*cjv, complex*cyv,complex*cjvp,complex*cyvp); @@ -33,7 +33,7 @@ double Yl_neg(double x) return ( sqrt(2.0/pi) * sin(x) )/sqrt(x); } -void computeB(complex* B, float radius, complex refIndex, float lambda, int Nl) +void computeB(complex* B, double radius, complex refIndex, double lambda, int Nl) { double k = (2*pi)/lambda; int b = 2; @@ -111,9 +111,9 @@ void computeB(complex* B, float radius, complex refIndex, float l numer = j_kr*j_d_knr*n - j_knr*j_d_kr; denom = j_knr*h_d_kr - h_kr*j_d_knr*n; - complex temp = numer/denom; + B[l] = numer/denom; - B[l] = scComplex(temp.real(), temp.imag()); + //B[l] = scComplex(temp.real(), temp.imag()); //cout<* B, float radius, complex refIndex, float l free(cyvp); } -void Legendre(float* P, float x, int Nl) +void Legendre(double* P, double x, int Nl) { //computes the legendre polynomials from orders 0 to Nl-1 P[0] = 1; @@ -140,7 +140,7 @@ void Legendre(float* P, float x, int Nl) } -complex integrateUi(float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi) +complex integrateUi(double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi) { /*This function integrates the incident field of magnitude M in the far zone in order to evaluate the field at the central pixel of a detector. @@ -150,20 +150,20 @@ complex integrateUi(float cAngleI, float cAngleO, float oAngleI, float oA oNAo = objective outer angle M = field magnitude*/ - float alphaIn = max(cAngleI, oAngleI); - float alphaOut = min(cAngleO,oAngleO); + double alphaIn = max(cAngleI, oAngleI); + double alphaOut = min(cAngleO,oAngleO); - complex Ui; + complex Ui; if(alphaIn > alphaOut) - Ui = complex(0.0, 0.0); + Ui = complex(0.0, 0.0); else - Ui = complex(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f); + Ui = complex(M * 2 * pi * (cos(alphaIn) - cos(alphaOut)), 0.0f); return Ui; } -void computeCondenserAlpha(float* alpha, int Nl, float cAngleI, float cAngleO) +void computeCondenserAlpha(double* alpha, int Nl, double cAngleI, double cAngleO) { /*This function computes the condenser integral in order to build the field of incident light 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) cAngleI, cAngleO = inner and outer condenser angles (inner and outer NA)*/ //compute the Legendre polynomials for the condenser aperature - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1)); + double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PcNAo, cos(cAngleO), Nl+1); - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1)); + double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PcNAi, cos(cAngleI), Nl+1); for(int l=0; l integrateUs(float r, float lambda, complex eta, - float cAngleI, float cAngleO, float oAngleI, float oAngleO, float M = 2*pi) +complex integrateUs(double r, double lambda, complex eta, + double cAngleI, double cAngleO, double oAngleI, double oAngleO, double M = 2*pi) { /*This function integrates the incident field of magnitude M in the far zone in order to evaluate the field at the central pixel of a detector. @@ -204,20 +204,20 @@ complex integrateUs(float r, float lambda, complex eta, M = field magnitude*/ //compute the required number of orders - float k = 2*pi/lambda; - int Nl = ceil( k + 4 * exp(log(k*r)/3) + 3 ); + double k = 2*pi/lambda; + int Nl = (int)ceil( k + 4 * exp(log(k*r)/3) + 3 ); //compute the material coefficients B - complex* B = (complex*)malloc(sizeof(complex)*Nl); + complex* B = (complex*)malloc(sizeof(complex)*Nl); //compute the Legendre polynomials for the condenser and objective aperatures - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1)); + double* PcNAo = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PcNAo, cos(cAngleO), Nl+1); - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1)); + double* PcNAi = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PcNAi, cos(cAngleI), Nl+1); - float* PoNAo = (float*)malloc(sizeof(float)*(Nl+1)); + double* PoNAo = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PoNAo, cos(oAngleO), Nl+1); - float* PoNAi = (float*)malloc(sizeof(float)*(Nl+1)); + double* PoNAi = (double*)malloc(sizeof(double)*(Nl+1)); Legendre(PoNAi, cos(oAngleI), Nl+1); //store the index of refraction; @@ -227,10 +227,10 @@ complex integrateUs(float r, float lambda, complex eta, computeB(B, r, IR, lambda, Nl); //aperature terms for the condenser (alpha) and objective (beta) - float alpha; - float beta; - float c; - complex Us(0.0, 0.0); + double alpha; + double beta; + double c; + complex Us(0.0, 0.0); for(int l=0; l Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi); - float I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag(); + complex Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi); + double I0 = Ui.real() * Ui.real() + Ui.imag() * Ui.imag(); I0 *= scaleI0; - float I; + //double I; SpecPair temp; - float nu; - complex eta; - complex Us, U; + double nu; + complex eta; + complex Us, U; - float vecLen = 0.0; - for(int i=0; i(EtaN[i].A, EtaK[i].A); + eta = complex(EtaN[i].A, EtaK[i].A); else - eta = complex(baseIR, 0.0); + eta = complex(baseIR, 0.0); //integrate the scattered field at the detector position Us = integrateUs(radius, lambda, eta, cAngleI, cAngleO, oAngleI, oAngleO, 2*pi); U = Us + Ui; - float I = U.real() * U.real() + U.imag() * U.imag(); + double I = U.real() * U.real() + U.imag() * U.imag(); temp.nu = nu; @@ -332,202 +332,13 @@ void pointSpectrum() vecLen = sqrt(vecLen); if(dispNormalize) - for(int i=0; i sampleUs(complex* B, float* Alpha, int Nl, float r, - float cAngleI, float cAngleO, float theta, float M = 2*pi) -{ - /*This function takes a point sample of the scattered field in the far zone - in order to evaluate the field at the central pixel of a detector. - r = sphere radius - lambda = wavelength - eta = index of refraction - cNAi = condenser inner NA - cNAo = condenser outer NA - theta = angle of the sample - M = field magnitude*/ - -/* - //compute the material coefficients B - //complex* B = (complex*)malloc(sizeof(complex)*Nl); - - //compute the Legendre polynomials for theta (at the objective) - float* Ptheta = (float*)malloc(sizeof(float)*(Nl+1)); - Legendre(Ptheta, cos(theta), Nl+1); - - //complex IR(eta.real(), eta.imag()); - - //aperature terms for the condenser (alpha) and objective (beta) - float beta; - float c; - complex Us(0.0, 0.0); - - for(int l=0; l numer(0.0, -(l*pi)/2.0); - Us += B[l] * exp(numer) * Ptheta[l] * Alpha[l] * pow(complex(0.0, 1.0), l); - - - } - //printf("Ptheta: %f\n", Ptheta[Nl-1]); - - return Us; - -}*/ - - -/* -void detectorSpectrum(int numSamples) -{ - //integrate across the objective aperature and calculate the resulting intensity on a detector - PD.StartTimer(SIMULATE_SPECTRUM); - //clear the previous spectrum - SimSpectrum.clear(); - - float dNu = 2.0f; - float lambda; - - //compute the angles based on NA - float cAngleI = asin(cNAi); - float cAngleO = asin(cNAo); - float oAngleI = asin(oNAi); - float oAngleO = asin(oNAo); - - //implement a reflection-mode system if necessary - if(opticsMode == ReflectionOpticsType){ - - //set the condenser to match the objective - cAngleI = oAngleI; - cAngleO = oAngleO; - - //invert the objective - oAngleO = pi - cAngleI; - oAngleI = pi - cAngleO; - } - - //compute Nl (maximum order of the spectrum) - //**************************************************************************** - float maxNu = EtaK.back().nu; - float maxLambda = 10000.0f/maxNu; - float k = 2*pi/maxLambda; - int Nl = ceil( k + 4 * exp(log(k*radius)/3) + 3 ); - int nLambda = EtaK.size(); - - //compute alpha (condenser integral) - //**************************************************************************** - //compute the Legendre polynomials for the condenser aperature - float* PcNAo = (float*)malloc(sizeof(float)*(Nl+1)); - Legendre(PcNAo, cos(cAngleO), Nl+1); - float* PcNAi = (float*)malloc(sizeof(float)*(Nl+1)); - Legendre(PcNAi, cos(cAngleI), Nl+1); - - //allocate space for the alpha array - float* alpha = (float*)malloc(sizeof(float)*(Nl + 1)); - computeCondenserAlpha(alpha, Nl, cAngleI, cAngleO); - - for(int l=0; l Ui; - - Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi); - I0 = Ui.real()*2*pi; - - float I; - SpecPair temp; - float nu; - complex eta; - complex Us, U; - - - - float vecLen = 0.0; - for(int i=0; i(EtaN[i].A, EtaK[i].A); - else - eta = complex(baseIR, 0.0); - - //allocate memory for the scattering coefficients - complex* B = (complex*)malloc(sizeof(complex)*Nl); - - complex IR(eta.real(), eta.imag()); - computeB(B, radius, IR, lambda, Nl); - - - //integrate the scattered field at the detector position - I = 0.0; - - for(int iTheta = 0; iTheta < numSamples; iTheta++) - { - theta = oAngleI + iTheta * dTheta; - Us = sampleUs(B, alpha, Nl, radius, cAngleI, cAngleO, theta, 2*pi); - - - //calculate the intensity and add - if(theta >= cAngleI && theta <= cAngleO) - U = Us + 2*(float)pi; - else - U = Us; - - I += (U.real() * U.real() + U.imag() * U.imag()) * sin(theta) * 2 * pi * dTheta; - - } - - temp.nu = nu; - - if(i == 0) - printf("I: %f I0: %f\n", I, I0); - - //set the spectrum value based on the current display type - if(dispSimType == AbsorbanceSpecType) - temp.A = -log10(I/I0); - else - temp.A = I; - - if(dispNormalize) - vecLen += temp.A * temp.A; - //temp.A = Us.real(); - SimSpectrum.push_back(temp); - - free(B); - } - vecLen = sqrt(vecLen); - - if(dispNormalize) - for(int i=0; i* B, int Nl, int nLambda) +void computeBArray(complex* B, int Nl, int nLambda) { - float nu; - complex eta; - float* Lambda = (float*)malloc(sizeof(float) * nLambda); + double nu; + complex eta; + double* Lambda = (double*)malloc(sizeof(double) * nLambda); //for each wavenumber nu - for(int i=0; i(EtaN[i].A, EtaK[i].A); + eta = complex(EtaN[i].A, EtaK[i].A); else - eta = complex(baseIR, 0.0); + eta = complex(baseIR, 0.0); //allocate memory for the scattering coefficients //complex* B = (complex*)malloc(sizeof(complex)*Nl); @@ -605,14 +416,13 @@ void computeBArray(complex* B, int Nl, int nLambda) } } -void computeOpticalParameters(float& cAngleI, float& cAngleO, float& oAngleI, float& oAngleO, float& I0, float* alpha, int Nl) +void computeOpticalParameters(double& cAngleI, double& cAngleO, double& oAngleI, double& oAngleO, double& I0, double* alpha, int Nl) { computeCassegrainAngles(cAngleI, cAngleO, oAngleI, oAngleO); //evaluate the incident field intensity I0 = 0.0; - float theta; - complex Ui; + complex Ui; Ui = integrateUi(cAngleI, cAngleO, oAngleI, oAngleO, 2*pi); I0 = Ui.real()*2*pi; @@ -631,24 +441,24 @@ void gpuDetectorSpectrum(int numSamples) //compute Nl (maximum order of the spectrum) int Nl = computeNl(); - float* alpha = (float*)malloc(sizeof(float)*(Nl + 1)); - float cAngleI, cAngleO, oAngleI, oAngleO, I0; + double* alpha = (double*)malloc(sizeof(double)*(Nl + 1)); + double cAngleI, cAngleO, oAngleI, oAngleO, I0; computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl); //allocate space for a list of wavelengths int nLambda = EtaK.size(); //allocate space for the 2D array (Nl x nu) of scattering coefficients - complex* B = (complex*)malloc(sizeof(complex) * Nl * nLambda); + complex* B = (complex*)malloc(sizeof(complex) * Nl * nLambda); computeBArray(B, Nl, nLambda); //allocate temporary space for the spectrum - float* I = (float*)malloc(sizeof(float) * EtaK.size()); + double* I = (double*)malloc(sizeof(double) * EtaK.size()); //compute the spectrum on the GPU PD.StartTimer(SIMULATE_GPU); - cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples); + cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, numSamples); PD.EndTimer(SIMULATE_GPU); updateSpectrum(I, I0, nLambda); @@ -666,44 +476,44 @@ void SimulateSpectrum() //detectorSpectrum(objectiveSamples); } -float absorbanceDistortion(){ +double absorbanceDistortion(){ //compute the mean of the spectrum - float sumSim = 0.0; - for(int i=0; i* B = (complex*)malloc(sizeof(complex) * Nl * nLambda); + complex* B = (complex*)malloc(sizeof(complex) * Nl * nLambda); computeBArray(B, Nl, nLambda); - float D; - float e = 0.001; - for(float i=0.0; i<=oNAo-step; i+=step) + double D; + double e = 0.001; + for(double i=0.0; i<=oNAo-step; i+=step) { - for(float o=oNAi+step; o<=1.0; o+=step) + for(double o=oNAi+step; o<=1.0; o+=step) { @@ -754,7 +564,7 @@ void MinimizeDistortion(){ computeOpticalParameters(cAngleI, cAngleO, oAngleI, oAngleO, I0, alpha, Nl); //simulate the spectrum - cudaComputeSpectrum(I, (float*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples); + cudaComputeSpectrum(I, (double*)B, alpha, Nl, nLambda, oAngleI, oAngleO, cAngleI, cAngleO, objectiveSamples); updateSpectrum(I, I0, nLambda); if(dispSimType == AbsorbanceSpecType) diff --git a/cudaKK.h b/cudaKK.h index 6459fc6..8efcab4 100644 --- a/cudaKK.h +++ b/cudaKK.h @@ -1,29 +1,29 @@ -__device__ float g(float v0, float v1) +__device__ double g(double v0, double v1) { return (v0 + v1)*log(abs(v0+v1)) + (v0-v1)*log(abs(v0-v1)); } -__device__ float hfin(float v0, float v1, float dv) +__device__ double hfin(double v0, double v1, double dv) { - float e = 0.001; - float t0 = g(v0+e, v1-dv)/dv; - float t1 = 2*g(v0+e, v1)/dv; - float t2 = g(v0+e, v1+dv)/dv; + double e = 0.001; + double t0 = g(v0+e, v1-dv)/dv; + double t1 = 2*g(v0+e, v1)/dv; + double t2 = g(v0+e, v1+dv)/dv; return -1.0/PI * (t0 - t1 + t2); } -__global__ void devKramersKronig(float* gpuN, float* gpuK, int numVals, float nuStart, float nuEnd, float nOffset) +__global__ void devKramersKronig(double* gpuN, double* gpuK, int numVals, double nuStart, double nuEnd, double nOffset) { int i = blockIdx.x * blockDim.x + threadIdx.x; if(i >= numVals) return; - float nuDelta = (nuEnd - nuStart)/(numVals - 1); + double nuDelta = (nuEnd - nuStart)/(numVals - 1); - float nu = nuStart + i*nuDelta; - float n = 0.0; - float jNu; - float jK; + double nu = nuStart + i*nuDelta; + double n = 0.0; + double jNu; + double jK; for(int j=1; j>>(gpuN, gpuK, nVals, nuStart, nuEnd, nOffset); - HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(float)*nVals, cudaMemcpyDeviceToHost)); + HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(double)*nVals, cudaMemcpyDeviceToHost)); //free resources HANDLE_ERROR(cudaFree(gpuK)); HANDLE_ERROR(cudaFree(gpuN)); } -__global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl, - int nSamples, float oThetaI, float oThetaO, float cThetaI, float cThetaO) +__global__ void devComputeSpectrum(double* I, double2* B, double* alpha, int Nl, + int nSamples, double oThetaI, double oThetaO, double cThetaI, double cThetaO) { int i = blockIdx.x * blockDim.x + threadIdx.x; //compute the delta-theta value - float dTheta = (oThetaO - oThetaI)/nSamples; + double dTheta = (oThetaO - oThetaI)/nSamples; //allocate space for the Legendre polynomials - float Ptheta[2]; - - float cosTheta, theta; - cuComplex Us; - cuComplex UsSample; - cuComplex U; - cuComplex Ui; - Ui.x = 2*PI; - Ui.y = 0.0; - cuComplex numer; + double Ptheta[2]; + + double cosTheta, theta; + cuDoubleComplex Us; + cuDoubleComplex UsSample; + cuDoubleComplex U; + //cuComplex Ui; + //Ui.x = 2*PI; + //Ui.y = 0.0; + cuDoubleComplex numer; numer.x = 0.0; - cuComplex exp_numer; - cuComplex iL; - cuComplex imag; + cuDoubleComplex exp_numer; + cuDoubleComplex iL; + cuDoubleComplex imag; imag.x = 0.0; imag.y = 1.0; - float realFac; - cuComplex complexFac; - float PlTheta; - float Isum = 0.0; - float maxVal = 0; - float val; + double realFac; + cuDoubleComplex complexFac; + double PlTheta; + double Isum = 0.0; + //float maxVal = 0; + //float val; for(int iTheta = 0; iTheta < nSamples; iTheta++) { //calculate theta @@ -131,12 +131,7 @@ __global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl, //increment the imaginary exponent i^l iL = cMult(iL, imag); - //val = cMag(Us); - //if(val > maxVal) - // maxVal = val; - - //Us += B[l] * exp(numer) * Ptheta[l] * alpha * M * pow(complex(0.0, 1.0), l); - + } //sum the scattered and incident fields @@ -150,20 +145,20 @@ __global__ void devComputeSpectrum(float* I, float2* B, float* alpha, int Nl, I[i] = Isum; } -void cudaComputeSpectrum(float* cpuI, float* cpuB, float* cpuAlpha, - int Nl, int nLambda, float oThetaI, float oThetaO, float cThetaI, float cThetaO, int nSamples) +void cudaComputeSpectrum(double* cpuI, double* cpuB, double* cpuAlpha, + int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples) { //copy everything to the GPU - float2* gpuB; - HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(float2) * nLambda * Nl)); - HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(float2) * nLambda * Nl, cudaMemcpyHostToDevice)); + double2* gpuB; + HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(double2) * nLambda * Nl)); + HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(double2) * nLambda * Nl, cudaMemcpyHostToDevice)); - float* gpuAlpha; - HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(float) * Nl)); - HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(float) * Nl, cudaMemcpyHostToDevice)); + double* gpuAlpha; + HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(double) * Nl)); + HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(double) * Nl, cudaMemcpyHostToDevice)); - float* gpuI; - HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(float) * nLambda)); + double* gpuI; + HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(double) * nLambda)); //call the kernel to compute the spectrum @@ -171,10 +166,10 @@ void cudaComputeSpectrum(float* cpuI, float* cpuB, float* cpuAlpha, dim3 grid(nLambda/block.x + 1); //devComputeSpectrum - devComputeSpectrum<<>>(gpuI, (float2*)gpuB, gpuAlpha, Nl, + devComputeSpectrum<<>>(gpuI, (double2*)gpuB, gpuAlpha, Nl, nSamples, oThetaI, oThetaO, cThetaI, cThetaO); - HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(float) * nLambda, cudaMemcpyDeviceToHost)); + HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(double) * nLambda, cudaMemcpyDeviceToHost)); HANDLE_ERROR(cudaFree(gpuB)); HANDLE_ERROR(cudaFree(gpuAlpha)); diff --git a/cudaMain.cu b/cudaMain.cu index 23225b1..9af99ce 100644 --- a/cudaMain.cu +++ b/cudaMain.cu @@ -5,45 +5,45 @@ #define PI 3.14159 #define BLOCK_SIZE 16 -__device__ cuComplex cMult(cuComplex a, cuComplex b) +__device__ cuDoubleComplex cMult(cuDoubleComplex a, cuDoubleComplex b) { - cuComplex result; + cuDoubleComplex result; result.x = a.x * b.x - a.y * b.y; result.y = a.x * b.y + a.y * b.x; return result; } -__device__ cuComplex cMult(cuComplex a, float b) +__device__ cuDoubleComplex cMult(cuDoubleComplex a, float b) { - cuComplex result; + cuDoubleComplex result; result.x = a.x * b; result.y = a.y * b; return result; } -__device__ cuComplex cAdd(cuComplex a, cuComplex b) +__device__ cuDoubleComplex cAdd(cuDoubleComplex a, cuDoubleComplex b) { - cuComplex r; + cuDoubleComplex r; r.x = a.x + b.x; r.y = a.y + b.y; return r; } -__device__ cuComplex cAdd(cuComplex a, float b) +__device__ cuDoubleComplex cAdd(cuDoubleComplex a, float b) { - cuComplex r; + cuDoubleComplex r; r.x = a.x + b; r.y = a.y; return r; } -__device__ cuComplex cExp(cuComplex a) +__device__ cuDoubleComplex cExp(cuDoubleComplex a) { - cuComplex r; + cuDoubleComplex r; r.x = exp(a.x) * cos(a.y); r.y = exp(a.x) * sin(a.y); @@ -51,9 +51,9 @@ __device__ cuComplex cExp(cuComplex a) return r; } -__device__ float cMag(cuComplex a) +__device__ double cMag(cuDoubleComplex a) { - float r = sqrt(a.x * a.x + a.y * a.y); + double r = sqrt(a.x * a.x + a.y * a.y); return r; } diff --git a/globals.h b/globals.h index e7a3a1a..57eb52f 100644 --- a/globals.h +++ b/globals.h @@ -8,14 +8,16 @@ #include using namespace std; +typedef float rtsFloat; + struct SpecPair{ - float nu; - float A; + double nu; + double A; }; struct Material{ - vector nu; - vector> eta; + vector nu; + vector> eta; string name; }; @@ -47,16 +49,16 @@ void FitDisplay(); //Update Functions void UpdateDisplay(); void SimulateSpectrum(); -void cudaKramersKronig(float* cpuN, float* cpuK, int nVals, float nuStart, float nuEnd, float nOffset); -void cudaComputeSpectrum(float* cpuI, float* cpuB, float* alpha, - int Nl, int nLambda, float oThetaI, float oThetaO, float cThetaI, float cThetaO, int nSamples); +void cudaKramersKronig(double* cpuN, double* cpuK, int nVals, double nuStart, double nuEnd, double nOffset); +void cudaComputeSpectrum(double* cpuI, double* cpuB, double* alpha, + int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples); //Window Parameters -extern float nuMin; -extern float nuMax; -extern float aMin; -extern float aMax; -extern float dNu; +extern double nuMin; +extern double nuMax; +extern double aMin; +extern double aMax; +extern double dNu; extern bool dispRefSpec; extern bool dispSimSpec; extern bool dispSimK; @@ -65,16 +67,16 @@ extern bool dispSimN; extern bool dispMatN; extern SpecType dispSimType; extern bool dispNormalize; -extern float dispNormFactor; +extern double dispNormFactor; -extern float dispScaleK; -extern float dispScaleN; +extern double dispScaleK; +extern double dispScaleN; //material parameters -extern float radius; -extern float baseIR; -extern float cA; +extern double radius; +extern double baseIR; +extern double cA; extern vector EtaK; extern vector EtaN; extern bool applyMaterial; @@ -86,22 +88,22 @@ void ChangeAbsorbance(); void SetMaterial(); //optical parameters -extern float cNAi; -extern float cNAo; -extern float oNAi; -extern float oNAo; +extern double cNAi; +extern double cNAo; +extern double oNAi; +extern double oNAo; extern OpticsType opticsMode; extern bool pointDetector; extern int objectiveSamples; //fitting parameters -extern float minMSE; +extern double minMSE; extern int maxFitIter; void EstimateMaterial(); -extern float scaleI0; -extern float refSlope; +extern double scaleI0; +extern double refSlope; -float ComputeDistortion(); +double ComputeDistortion(); void MinimizeDistortion(); diff --git a/interactivemie.h b/interactivemie.h index 0d06408..02eab7e 100644 --- a/interactivemie.h +++ b/interactivemie.h @@ -48,7 +48,7 @@ public: //material selection combo box ui.cmbMaterial->clear(); - for(int i=0; iaddItem(MaterialList[i].name.c_str(), i); ui.cmbMaterial->setCurrentIndex(currentMaterial); diff --git a/interactivemie.ui b/interactivemie.ui index 32037b8..3e89aa7 100644 --- a/interactivemie.ui +++ b/interactivemie.ui @@ -350,6 +350,9 @@ 22 + + 6 + -9999.000000000000000 @@ -362,7 +365,7 @@ - false + true @@ -372,6 +375,9 @@ 22 + + 1 + 4000 @@ -404,6 +410,9 @@ 22 + + 6 + -9999.000000000000000 @@ -419,7 +428,7 @@ - false + true @@ -429,8 +438,11 @@ 22 + + 1 + - 4000 + 100000 10 @@ -509,6 +521,9 @@ 22 + + 6 + -99.989999999999995 diff --git a/main.cpp b/main.cpp index 386a95b..ea331fe 100644 --- a/main.cpp +++ b/main.cpp @@ -19,15 +19,15 @@ vector EtaK; vector EtaN; int currentSpec = 0; -float nuMin = 800; -float nuMax = 4000; -float dNu = 2; +double nuMin = 800; +double nuMax = 4000; +double dNu = 2; -float aMin = 0; -float aMax = 1; +double aMin = 0; +double aMax = 1; -float scaleI0 = 1.0; -float refSlope = 0.0; +double scaleI0 = 1.0; +double refSlope = 0.0; bool dispRefSpec = true; bool dispSimSpec = true; @@ -35,17 +35,17 @@ bool dispSimK = true; bool dispMatK = true; bool dispSimN = true; bool dispMatN = true; -float dispScaleK = 1.0; -float dispScaleN = 1.0; +double dispScaleK = 1.0; +double dispScaleN = 1.0; SpecType dispSimType = AbsorbanceSpecType; bool dispNormalize = false; -float dispNormFactor = 1.0; +double dispNormFactor = 1.0; //material parameters -float radius = 4.0f; -float baseIR = 1.49f; -float cA = 1.0; +double radius = 4.0f; +double baseIR = 1.49f; +double cA = 1.0; //vector KMaterial; //vector NMaterial; bool applyMaterial = true; @@ -53,16 +53,16 @@ vector MaterialList; int currentMaterial = 0; //optical parameters -float cNAi = 0.0; -float cNAo = 0.6; -float oNAi = 0.0; -float oNAo = 0.6; +double cNAi = 0.0; +double cNAo = 0.6; +double oNAi = 0.0; +double oNAo = 0.6; OpticsType opticsMode = TransmissionOpticsType; bool pointDetector = false; int objectiveSamples = 200; //fitting parameters -float minMSE = 0.00001; +double minMSE = 0.00001; int maxFitIter = 20; void TempSimSpectrum() @@ -71,7 +71,7 @@ void TempSimSpectrum() for(int i=800; i<4000; i++) { temp.nu = i; - temp.A = sin((float)i/200); + temp.A = sin((double)i/200); SimSpectrum.push_back(temp); } } @@ -94,11 +94,11 @@ void LoadMaterial(string fileNameK, string fileNameN, string materialName) exit(1); } - complex eta; - int j; - for(int i=0; i eta; + //int j; + for(unsigned int i=0; i(NMaterial[i].A, KMaterial[i].A); + eta = complex(NMaterial[i].A, KMaterial[i].A); newMaterial.eta.push_back(eta); } MaterialList.push_back(newMaterial); @@ -114,15 +114,15 @@ void LoadMaterial(string fileNameK, string materialName){ //compute the real IR using Kramers Kronig //copy the absorbance values into a linear array - float* k = (float*)malloc(sizeof(float) * KMaterial.size()); - float* n = (float*)malloc(sizeof(float) * KMaterial.size()); - for(int i=0; i eta; - int j; - for(int i=0; i eta; + for(unsigned int i=0; i(NMaterial[i].A, KMaterial[i].A); + eta = complex(NMaterial[i].A, KMaterial[i].A); newMaterial.eta.push_back(eta); } @@ -144,12 +143,12 @@ void LoadMaterial(string fileNameK, string materialName){ } void FitDisplay(){ - float minA = 99999.0; - float maxA = -99999.0; - float k, n; + double minA = 99999.0; + double maxA = -99999.0; + double k, n; if(dispSimSpec) - for(int i=0; i 0) - for(int i=0; i= nuMin && nu <= nuMax){ @@ -253,7 +252,7 @@ void SetMaterial() EtaN.clear(); int nSamples = MaterialList[currentMaterial].eta.size(); - float nu; + double nu; SpecPair temp; for(int i=0; i