From 396a5f1225210b5c52e5b6df9c9d96e7f63176e5 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Thu, 15 May 2014 12:16:59 -0500 Subject: [PATCH] added custom code for dealing with command-line arguments --- CMakeLists.txt | 21 +++++++++++++-------- dataTypes.h | 4 ++-- defaults.h | 36 +++++++++++++++++++++--------------- main.cpp | 40 ++++++++++++++++++++++++++++------------ nearfield.h | 4 ++-- nfScalarUf.cu | 58 +++++++++++++++++++++++++++++----------------------------- nfScalarUfLut.cu | 192 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------ nfScalarUpLut.cu | 276 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------ nfScalarUs.cu | 266 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------- options.h | 609 --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- planewave.h | 26 +++++++++++++------------- rts | 2 +- sphere.h | 3 ++- 13 files changed, 478 insertions(+), 1059 deletions(-) delete mode 100644 options.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 03f6820..2af58cf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,10 +12,10 @@ set(CMAKE_AUTOMOC ON) set(CMAKE_INCLUDE_CURRENT_DIR ON) #find BOOST -set(Boost_USE_STATIC_LIBS ON) -set(Boost_USE_MULTITHREADED ON) -set(Boost_USE_STATIC_RUNTIME OFF) -find_package( Boost 1.46.0 COMPONENTS program_options ) +#set(Boost_USE_STATIC_LIBS ON) +#set(Boost_USE_MULTITHREADED ON) +#set(Boost_USE_STATIC_RUNTIME OFF) +#find_package( Boost 1.46.0 COMPONENTS program_options ) #find the Qt5 find_package(Qt5Widgets REQUIRED) @@ -28,7 +28,7 @@ include_directories(${QT_INCLUDE_DIRECTORY}) find_package(CUDA) #ask the user for the RTS location -find_package(RTS REQUIRED) +#find_package(RTS REQUIRED) #set the include directories include_directories( @@ -37,8 +37,9 @@ include_directories( ${Qt5Core_INCLUDE_DIRS} ${Qt5Gui_INCLUDE_DIRS} # ${Qt5OpenGL_INCLUDE_DIRS} - ${RTS_INCLUDE_DIR} - ${Boost_INCLUDE_DIR} +# ${RTS_INCLUDE_DIR} +# ${Boost_INCLUDE_DIR} + ${CMAKE_CURRENT_SOURCE_DIR} ) #build position independent code for Qt (-fPIC) @@ -55,6 +56,9 @@ file(GLOB SRC_H "*.h") file(GLOB SRC_UI "*.ui") file(GLOB SRC_CU "*.cu") +#assign RTS source files +file(GLOB SRC_RTS "rts/source/*.cpp") + #determine which source files have to be moc'd Qt5_wrap_cpp(UI_MOC ${SRC_H}) Qt5_wrap_ui(UI_H ${SRC_UI}) @@ -70,6 +74,7 @@ cuda_add_executable(bimsim ${UI_H} ${SRC_UI} ${SRC_CU} + ${SRC_RTS} ) #specify which qt5 modules to use @@ -83,7 +88,7 @@ target_link_libraries(bimsim # ${Qt5OpenGL_LIBRARIES} ${CUDA_cufft_LIBRARY} ${CUDA_cublas_LIBRARY} - ${Boost_LIBRARIES} +# ${Boost_LIBRARIES} ) diff --git a/dataTypes.h b/dataTypes.h index f30bc88..a78cd14 100644 --- a/dataTypes.h +++ b/dataTypes.h @@ -14,10 +14,10 @@ typedef double ptype; #define BLOCK 256 #define SQRT_BLOCK 16 -#define PI 3.14159 +#define PI 3.14159f //a very small number -#define EPSILON 0.00001 +#define EPSILON 0.00001f //CUDA hybrid code - complex class should run on both the CPU and GPU diff --git a/defaults.h b/defaults.h index f0dd9d5..03c19c5 100644 --- a/defaults.h +++ b/defaults.h @@ -8,22 +8,26 @@ #define DEFAULT_SPHERE_A 1 //default near field parameters -#define DEFAULT_LAMBDA 1 -#define DEFAULT_AMPLITUDE 1 +#define DEFAULT_LAMBDA "1" +#define DEFAULT_AMPLITUDE "1" +#define DEFAULT_MATERIAL "1.4 0.05" #define DEFAULT_N 1.4 #define DEFAULT_K 0.5 -#define DEFAULT_FOCUS_X 0 -#define DEFAULT_FOCUS_Y 0 -#define DEFAULT_FOCUS_Z 0 +#define DEFAULT_FOCUS "0 0 0" +//#define DEFAULT_FOCUS_X "0" +//#define DEFAULT_FOCUS_Y "0" +//#define DEFAULT_FOCUS_Z "0" //#define DEFAULT_INCIDENT_ORDER 20 #define DEFAULT_STABILITY_PARM 1.4 //optics -#define DEFAULT_CONDENSER_MIN 0 -#define DEFAULT_CONDENSER_MAX 1 +//#define DEFAULT_CONDENSER_MIN "0.0" +//#define DEFAULT_CONDENSER_MAX "1.0" +#define DEFAULT_CONDENSER "0 1" -#define DEFAULT_OBJECTIVE_MIN 0 -#define DEFAULT_OBJECTIVE_MAX 1 +//#define DEFAULT_OBJECTIVE_MIN "0" +//#define DEFAULT_OBJECTIVE_MAX "1" +#define DEFAULT_OBJECTIVE "0 1" //incident light direction #define DEFAULT_K_THETA 0 @@ -35,15 +39,17 @@ #define DEFAULT_APPEND false //#define DEFAULT_OUTPUT_POINT fileoutStruct::imageObjective - +#define DEFAULT_PLANE_MIN "-5 0 -5" #define DEFAULT_PLANE_MIN_X -5 #define DEFAULT_PLANE_MIN_Y 0 #define DEFAULT_PLANE_MIN_Z -5 +#define DEFAULT_PLANE_MAX "5 0 5" #define DEFAULT_PLANE_MAX_X 5 #define DEFAULT_PLANE_MAX_Y 0 #define DEFAULT_PLANE_MAX_Z 5 +#define DEFAULT_PLANE_NORM "0 1 0" #define DEFAULT_PLANE_NORM_X 0 #define DEFAULT_PLANE_NORM_Y 1 #define DEFAULT_PLANE_NORM_Z 0 @@ -67,16 +73,16 @@ */ -#define DEFAULT_FIELD_ORDER 10 +#define DEFAULT_FIELD_ORDER "10" -#define DEFAULT_SAMPLES 400 +#define DEFAULT_SAMPLES "400" -#define DEFAULT_SLICE_RES 256 +#define DEFAULT_SLICE_RES "256" #define DEFAULT_SPHERE_THETA_R 1000 -#define DEFAULT_PADDING 1 -#define DEFAULT_SUPERSAMPLE 1 +#define DEFAULT_PADDING "1" +#define DEFAULT_SUPERSAMPLE "1" #define DEFAULT_INTENSITY_FILE "out_i.bmp" #define DEFAULT_TRANSMITTANCE_FILE "" diff --git a/main.cpp b/main.cpp index dfb4038..74ebc4c 100644 --- a/main.cpp +++ b/main.cpp @@ -12,7 +12,9 @@ microscopeStruct* SCOPE; #include "fieldslice.h" #include "fileout.h" -#include "options.h" +//#include "options.h" +#include "arguments.h" +#include "rts/tools/arguments.h" #include "montecarlo.h" #include "rts/math/point.h" #include "rts/math/spherical_bessel.h" @@ -29,28 +31,42 @@ microscopeStruct* SCOPE; #include "qtMainDialog.h" bool gui = false; +#ifdef _WIN32 +bool ansi = false; +#else +bool ansi = true; +#endif + fileoutStruct gFileOut; bool verbose = false; using namespace std; -int cbessjyva(double v,complex z,double &vm,complex*cjv, +int cbessjyva(double v,complex z,double &vm,complex*cjv, complex*cyv,complex*cjvp,complex*cyvp); int main(int argc, char *argv[]) { + //arguments test + rts::arglist args; + SetArguments(args); - //benchtest planewave class - rts::vector k(1, 0, 0); - rts::vector E(2, 2, 0); - planewave P(k, E); - - std::cout<init(); diff --git a/nearfield.h b/nearfield.h index 56f9120..b682c0c 100644 --- a/nearfield.h +++ b/nearfield.h @@ -8,7 +8,7 @@ #include "sphere.h" #include -#define EPSILON_FLOAT 0.000001 +#define EPSILON_FLOAT 0.000001f //This structure stores values relevant to creating the near field struct nearfieldStruct @@ -29,7 +29,7 @@ struct nearfieldStruct bsVector k; //cartesian coordinates, normalized bsPoint focus; - //slice position and orientation in world space + //slice position and orientation in world space rts::quad pos; //slices for the focused field diff --git a/nfScalarUf.cu b/nfScalarUf.cu index 261c3d8..4a2aa27 100644 --- a/nfScalarUf.cu +++ b/nfScalarUf.cu @@ -3,8 +3,8 @@ #include "rts/math/legendre.h" #include #include "rts/cuda/error.h" -#include "rts/cuda/timer.h" - +#include "rts/cuda/timer.h" + //Incident field for a single plane wave __global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR) { @@ -33,15 +33,15 @@ __global__ void gpuScalarUfp(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, p //get the rtsPoint in world space and then the r vector bsPoint p = ABCD(u, v); bsVector r = p - f; - //ptype d = r.len(); - - ptype k_dot_r = kmag * k.dot(r); - bsComplex d(0, k_dot_r); - + //ptype d = r.len(); + + ptype k_dot_r = kmag * k.dot(r); + bsComplex d(0, k_dot_r); + Uf[i] = exp(d) * A; } - + //Incident field for a focused point source __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, ptype A, bsRect ABCD, int uR, int vR, ptype cosAlpha, ptype cosBeta, int nl, ptype j_conv = 1.4) { @@ -70,11 +70,11 @@ __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, pt //get the rtsPoint in world space and then the r vector bsPoint p = ABCD(u, v); bsVector r = p - f; - ptype d = r.len(); - if(d < EPSILON_FLOAT) - { - Uf[i] = A * 2 * PI * (cosAlpha - cosBeta); - return; + ptype d = r.len(); + if(d < EPSILON_FLOAT) + { + Uf[i] = A * 2 * PI * (cosAlpha - cosBeta); + return; } //get info for the light direction and frequency @@ -94,10 +94,10 @@ __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, pt ptype P[2]; //get the angle between k and r (light direction and position vector) ptype cosTheta; - cosTheta = k.dot(r); - - //deal with the degenerate case where r == 0 - //if(isnan(cosTheta)) + cosTheta = k.dot(r); + + //deal with the degenerate case where r == 0 + //if(isnan(cosTheta)) // cosTheta = 0; rts::init_legendre(cosTheta, P[0], P[1]); @@ -162,12 +162,12 @@ __global__ void gpuScalarUf(bsComplex* Uf, bsVector k, ptype kmag, bsPoint f, pt void nearfieldStruct::scalarUf() { - - gpuStartTimer(); + + gpuStartTimer(); //create one thread for each pixel of the field slice dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); - dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); + dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); //if we are computing a plane wave, call the gpuScalarUfp function if(planeWave) @@ -176,15 +176,15 @@ void nearfieldStruct::scalarUf() } //otherwise compute the condenser info and create a focused field else - { - //pre-compute the cosine of the obscuration and objective angles - //cout<<"Condenser angle in: "<>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1], cosAlpha, cosBeta, m); - } - - t_Uf = gpuStopTimer(); -} + } + + t_Uf = gpuStopTimer(); +} diff --git a/nfScalarUfLut.cu b/nfScalarUfLut.cu index 71abd73..31c866d 100644 --- a/nfScalarUfLut.cu +++ b/nfScalarUfLut.cu @@ -1,7 +1,7 @@ #include "nearfield.h" #include "rts/math/legendre.h" -#include "rts/cuda/error.h" +#include "rts/cuda/error.h" #include "rts/cuda/timer.h" texture texJ; @@ -25,100 +25,100 @@ __global__ void gpuScalarUfLut(bsComplex* Uf, bsRect ABCD, int uR, int vR, bsPoi */ - //get the current coordinate in the plane slice - int iu = blockIdx.x * blockDim.x + threadIdx.x; - int iv = blockIdx.y * blockDim.y + threadIdx.y; - - //make sure that the thread indices are in-bounds - if(iu >= uR || iv >= vR) return; - - //compute the index (easier access to the scalar field array) - int i = iv*uR + iu; - - //compute the parameters for u and v - ptype u = (ptype)iu / (uR); - ptype v = (ptype)iv / (vR); - - - - //get the rtsPoint in world space and then the r vector - bsPoint p = ABCD(u, v); - bsVector r = p - f; + //get the current coordinate in the plane slice + int iu = blockIdx.x * blockDim.x + threadIdx.x; + int iv = blockIdx.y * blockDim.y + threadIdx.y; + + //make sure that the thread indices are in-bounds + if(iu >= uR || iv >= vR) return; + + //compute the index (easier access to the scalar field array) + int i = iv*uR + iu; + + //compute the parameters for u and v + ptype u = (ptype)iu / (uR); + ptype v = (ptype)iv / (vR); + + + + //get the rtsPoint in world space and then the r vector + bsPoint p = ABCD(u, v); + bsVector r = p - f; ptype d = r.len(); if(d == 0) { Uf[i] = A * 2 * PI * (cosAlpha - cosBeta); return; - } - - //get info for the light direction and frequency - r = r.norm(); - - //compute the imaginary factor i^l - bsComplex im = bsComplex(0, 1); - bsComplex il = bsComplex(1, 0); - - //Legendre functions are computed dynamically to save memory - //initialize the Legendre functions - - ptype P[2]; - //get the angle between k and r (light direction and position vector) - ptype cosTheta; + } + + //get info for the light direction and frequency + r = r.norm(); + + //compute the imaginary factor i^l + bsComplex im = bsComplex(0, 1); + bsComplex il = bsComplex(1, 0); + + //Legendre functions are computed dynamically to save memory + //initialize the Legendre functions + + ptype P[2]; + //get the angle between k and r (light direction and position vector) + ptype cosTheta; cosTheta = k.dot(r); - rts::init_legendre(cosTheta, P[0], P[1]); - - //initialize legendre functions for the cassegrain angles - ptype Palpha[3]; - rts::init_legendre(cosAlpha, Palpha[0], Palpha[1]); - Palpha[2] = 1; - - ptype Pbeta[3]; - rts::init_legendre(cosBeta, Pbeta[0], Pbeta[1]); - Pbeta[2] = 1; - - //for each order l - bsComplex sumUf(0.0, 0.0); - ptype jl = 0.0; + rts::init_legendre(cosTheta, P[0], P[1]); + + //initialize legendre functions for the cassegrain angles + ptype Palpha[3]; + rts::init_legendre(cosAlpha, Palpha[0], Palpha[1]); + Palpha[2] = 1; + + ptype Pbeta[3]; + rts::init_legendre(cosBeta, Pbeta[0], Pbeta[1]); + Pbeta[2] = 1; + + //for each order l + bsComplex sumUf(0, 0); + ptype jl = 0; ptype Pl; - ptype di = ( (d - dmin)/(dmax - dmin) ) * (dR - 1); - for(int l = 0; l<=nl; l++) - { - jl = tex2D(texJ, l + 0.5, di + 0.5); - if(l==0) - Pl = P[0]; - else if(l==1) - { - Pl = P[1]; - - //adjust the cassegrain Legendre function - Palpha[2] = Palpha[0]; - rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); - Pbeta[2] = Pbeta[0]; - rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); - } - else - { - rts::shift_legendre(l, cosTheta, P[0], P[1]); - - Pl = P[1]; - - //adjust the cassegrain outer Legendre function - Palpha[2] = Palpha[0]; - rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); - Pbeta[2] = Pbeta[0]; - rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); - } - + ptype di = ( (d - dmin)/(dmax - dmin) ) * (dR - 1); + for(int l = 0; l<=nl; l++) + { + jl = tex2D(texJ, l + 0.5f, di + 0.5f); + if(l==0) + Pl = P[0]; + else if(l==1) + { + Pl = P[1]; + + //adjust the cassegrain Legendre function + Palpha[2] = Palpha[0]; + rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); + Pbeta[2] = Pbeta[0]; + rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); + } + else + { + rts::shift_legendre(l, cosTheta, P[0], P[1]); + + Pl = P[1]; + + //adjust the cassegrain outer Legendre function + Palpha[2] = Palpha[0]; + rts::shift_legendre(l+1, cosAlpha, Palpha[0], Palpha[1]); + Pbeta[2] = Pbeta[0]; + rts::shift_legendre(l+1, cosBeta, Pbeta[0], Pbeta[1]); + } + sumUf += il * jl * Pl * (Palpha[1] - Palpha[2] - Pbeta[1] + Pbeta[2]); //sumUf += jl; - - il *= im; - } - + + il *= im; + } + Uf[i] = sumUf * 2 * PI * A; - //Uf[i] = u; + //Uf[i] = u; //return; } @@ -159,23 +159,23 @@ void nearfieldStruct::scalarUfLut() HANDLE_ERROR( cudaMemcpy2DToArray(arrayJ, 0, 0, j, (m+1)*sizeof(float), (m+1)*sizeof(float), dR, cudaMemcpyHostToDevice)); //----------------Compute the focused field - //create one thread for each pixel of the field slice - dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); + //create one thread for each pixel of the field slice + dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); dim3 dimGrid((Uf.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (Uf.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); - - //if we are computing a plane wave, call the gpuScalarUfp function - if(planeWave) - { - gpuScalarUfp<<>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1]); - } - //otherwise compute the condenser info and create a focused field - else + + //if we are computing a plane wave, call the gpuScalarUfp function + if(planeWave) + { + gpuScalarUfp<<>>(Uf.x_hat, k, 2 * PI / lambda, focus, A, pos, Uf.R[0], Uf.R[1]); + } + //otherwise compute the condenser info and create a focused field + else { //pre-compute the cosine of the obscuration and objective angles ptype cosAlpha = cos(asin(condenser[0])); - ptype cosBeta = cos(asin(condenser[1])); - //compute the scalar Uf field (this will be in the x_hat channel of Uf) - gpuScalarUfLut<<>>(Uf.x_hat, pos, Uf.R[0], Uf.R[1], focus, k, A, cosAlpha, cosBeta, m, d_min, d_max, dR); + ptype cosBeta = cos(asin(condenser[1])); + //compute the scalar Uf field (this will be in the x_hat channel of Uf) + gpuScalarUfLut<<>>(Uf.x_hat, pos, Uf.R[0], Uf.R[1], focus, k, A, cosAlpha, cosBeta, m, d_min, d_max, dR); } diff --git a/nfScalarUpLut.cu b/nfScalarUpLut.cu index 9a2c0a1..9fbbdad 100644 --- a/nfScalarUpLut.cu +++ b/nfScalarUpLut.cu @@ -3,30 +3,30 @@ #include "rts/math/legendre.h" #include #include "rts/cuda/error.h" -#include "rts/cuda/timer.h" - -texture texUsp; -texture texUip; - +#include "rts/cuda/timer.h" + +texture texUsp; +texture texUip; + __global__ void gpuScalarUpLut(bsComplex* Us, bsVector* k, int nk, ptype kmag, ptype a, ptype dmin, ptype dmax, bsPoint f, bsPoint ps, ptype A, bsRect ABCD, int uR, int vR, int dR, int aR, int thetaR) -{ - /*This function uses Monte-Carlo integration to sample a texture-based LUT describing the scattered field - produced by a plane wave through a sphere. The MC sampling is used to approximate a focused field. - - Us = final scattered field - k = list of incoming plane waves (Monte-Carlo samples) - nk = number of incoming MC samples - kmag= magnitude of the incoming field 2pi/lambda - dmin= minimum distance of the Usp texture - dmax= maximum distance of the Usp texture - f = position of the focus - ps = position of the sphere - A = total amplitude of the incident field arriving at the focal spot - ABCD= rectangle representing the field slice - uR = resolution of the field slice in the u direction - vR = resolution of the field slice in the v direction - dR = resolution of the Usp texture in the d direction - thetaR= resolution of the Usp texture in the theta direction +{ + /*This function uses Monte-Carlo integration to sample a texture-based LUT describing the scattered field + produced by a plane wave through a sphere. The MC sampling is used to approximate a focused field. + + Us = final scattered field + k = list of incoming plane waves (Monte-Carlo samples) + nk = number of incoming MC samples + kmag= magnitude of the incoming field 2pi/lambda + dmin= minimum distance of the Usp texture + dmax= maximum distance of the Usp texture + f = position of the focus + ps = position of the sphere + A = total amplitude of the incident field arriving at the focal spot + ABCD= rectangle representing the field slice + uR = resolution of the field slice in the u direction + vR = resolution of the field slice in the v direction + dR = resolution of the Usp texture in the d direction + thetaR= resolution of the Usp texture in the theta direction */ //get the current coordinate in the plane slice @@ -46,47 +46,47 @@ __global__ void gpuScalarUpLut(bsComplex* Us, bsVector* k, int nk, ptype kmag, p //get the rtsPoint in world space and then the r vector bsPoint p = ABCD(u, v); bsVector r = p - ps; - ptype d = r.len(); - float di = ( (d - max(a, dmin))/(dmax - max(a, dmin)) ) * (dR - 1); - float ai = ( (d - dmin)/(a - dmin)) * (aR - 1); - - bsComplex sumUs(0, 0); - //for each plane wave in the wave list - for(int iw = 0; iw < nk; iw++) - { - //normalize the direction vectors and find their inner product - r = r.norm(); - ptype cos_theta = k[iw].dot(r); - if(cos_theta < -1) - cos_theta = -1; - if(cos_theta > 1) - cos_theta = 1; - float thetai = ( acos(cos_theta) / PI ) * (thetaR - 1); - - //compute the phase factor for spheres that are not at the origin - bsVector c = ps - f; - bsComplex phase = exp(bsComplex(0, kmag * k[iw].dot(c))); - - //compute the internal field if we are inside a sphere + ptype d = r.len(); + float di = ( (d - max(a, dmin))/(dmax - max(a, dmin)) ) * (dR - 1); + float ai = ( (d - dmin)/(a - dmin)) * (aR - 1); + + bsComplex sumUs(0, 0); + //for each plane wave in the wave list + for(int iw = 0; iw < nk; iw++) + { + //normalize the direction vectors and find their inner product + r = r.norm(); + ptype cos_theta = k[iw].dot(r); + if(cos_theta < -1) + cos_theta = -1; + if(cos_theta > 1) + cos_theta = 1; + float thetai = ( acos(cos_theta) / PI ) * (thetaR - 1); + + //compute the phase factor for spheres that are not at the origin + bsVector c = ps - f; + bsComplex phase = exp(bsComplex(0, kmag * k[iw].dot(c))); + + //compute the internal field if we are inside a sphere if(d < a) { - float2 Uip = tex2D(texUip, ai + 0.5, thetai + 0.5); - sumUs += (1.0/nk) * A * phase * bsComplex(Uip.x, Uip.y); - } - //otherwise compute the scattered field - else - { - float2 Usp = tex2D(texUsp, di + 0.5, thetai + 0.5); - sumUs += (1.0/nk) * A * phase * bsComplex(Usp.x, Usp.y); - } - - } - - Us[i] += sumUs; -} - -void nearfieldStruct::scalarUpLut() -{ + float2 Uip = tex2D(texUip, ai + 0.5f, thetai + 0.5f); + sumUs += (1.0f/nk) * A * phase * bsComplex(Uip.x, Uip.y); + } + //otherwise compute the scattered field + else + { + float2 Usp = tex2D(texUsp, di + 0.5f, thetai + 0.5f); + sumUs += (1.0f/nk) * A * phase * bsComplex(Usp.x, Usp.y); + } + + } + + Us[i] += sumUs; +} + +void nearfieldStruct::scalarUpLut() +{ //get the number of spheres int nSpheres = sVector.size(); @@ -103,90 +103,90 @@ void nearfieldStruct::scalarUpLut() //create one thread for each pixel of the field slice dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK); dim3 dimGrid((U.R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (U.R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK); - - //copy Monte-Carlo samples to the GPU and determine the incident amplitude (plane-wave specific stuff) - bsVector* gpuk; - int nWaves; - ptype subA; - if(planeWave) - { - nWaves = 1; - HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) ) ); - HANDLE_ERROR(cudaMemcpy( gpuk, &k, sizeof(bsVector), cudaMemcpyHostToDevice)); - subA = A; - } - else - { - nWaves = inWaves.size(); - HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) * nWaves ) ); - HANDLE_ERROR(cudaMemcpy( gpuk, &inWaves[0], sizeof(bsVector) * nWaves, cudaMemcpyHostToDevice)); - //compute the amplitude that makes it through the condenser - subA = 2 * PI * A * ( (1 - cos(asin(condenser[1]))) - (1 - cos(asin(condenser[0]))) ); - } + + //copy Monte-Carlo samples to the GPU and determine the incident amplitude (plane-wave specific stuff) + bsVector* gpuk; + int nWaves; + ptype subA; + if(planeWave) + { + nWaves = 1; + HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) ) ); + HANDLE_ERROR(cudaMemcpy( gpuk, &k, sizeof(bsVector), cudaMemcpyHostToDevice)); + subA = A; + } + else + { + nWaves = inWaves.size(); + HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) * nWaves ) ); + HANDLE_ERROR(cudaMemcpy( gpuk, &inWaves[0], sizeof(bsVector) * nWaves, cudaMemcpyHostToDevice)); + //compute the amplitude that makes it through the condenser + subA = 2 * PI * A * ( (1 - cos(asin(condenser[1]))) - (1 - cos(asin(condenser[0]))) ); + } //for each sphere for(int s = 0; s>>(U.x_hat, - gpuk, + gpuk, nWaves, - 2 * PI / lambda, - sVector[s].a, - sVector[s].d_min, + 2 * PI / lambda, + sVector[s].a, + sVector[s].d_min, sVector[s].d_max, focus, - sVector[s].p, - subA, + sVector[s].p, + subA, pos, U.R[0], - U.R[1], - dR, - aR, - thetaR); - - cudaFreeArray(arrayUsp); - cudaFreeArray(arrayUip); - - } - + U.R[1], + dR, + aR, + thetaR); + + cudaFreeArray(arrayUsp); + cudaFreeArray(arrayUip); + + } + //store the time to compute the scattered field - t_Us = gpuStopTimer(); - - //free monte-carlo samples - cudaFree(gpuk); - -} + t_Us = gpuStopTimer(); + + //free monte-carlo samples + cudaFree(gpuk); + +} diff --git a/nfScalarUs.cu b/nfScalarUs.cu index 813fb0b..2aeb227 100644 --- a/nfScalarUs.cu +++ b/nfScalarUs.cu @@ -7,47 +7,47 @@ __device__ bsComplex calc_Us(ptype kd, ptype cos_theta, int Nl, bsComplex* B) { - //initialize the spherical Bessel functions - ptype j[2]; - rts::init_sbesselj(kd, j); - ptype y[2]; - rts::init_sbessely(kd, y); - - //initialize the Legendre polynomial - ptype P[2]; - rts::init_legendre(cos_theta, P[0], P[1]); - - //initialize the spherical Hankel function - bsComplex h((ptype)0, (ptype)0); - - //initialize the result - bsComplex Us((ptype)0, (ptype)0); - - //for each order up to Nl - for(int l=0; l<=Nl; l++) - { - if(l == 0) - { - h.r = j[0]; - h.i = y[0]; - Us += B[0] * h * P[0]; - } - else - { - //shift the bessel functions and legendre polynomials - if(l > 1) - { - rts::shift_legendre(l, cos_theta, P[0], P[1]); - rts::shift_sbesselj(l, kd, j); - rts::shift_sbessely(l, kd, y); - } - - h.r = j[1]; - h.i = y[1]; - Us += B[l] * h * P[1]; - - - } + //initialize the spherical Bessel functions + ptype j[2]; + rts::init_sbesselj(kd, j); + ptype y[2]; + rts::init_sbessely(kd, y); + + //initialize the Legendre polynomial + ptype P[2]; + rts::init_legendre(cos_theta, P[0], P[1]); + + //initialize the spherical Hankel function + bsComplex h((ptype)0, (ptype)0); + + //initialize the result + bsComplex Us((ptype)0, (ptype)0); + + //for each order up to Nl + for(int l=0; l<=Nl; l++) + { + if(l == 0) + { + h.r = j[0]; + h.i = y[0]; + Us += B[0] * h * P[0]; + } + else + { + //shift the bessel functions and legendre polynomials + if(l > 1) + { + rts::shift_legendre(l, cos_theta, P[0], P[1]); + rts::shift_sbesselj(l, kd, j); + rts::shift_sbessely(l, kd, y); + } + + h.r = j[1]; + h.i = y[1]; + Us += B[l] * h * P[1]; + + + } } return Us; } @@ -59,41 +59,41 @@ __device__ bsComplex calc_Ui(bsComplex knd, ptype cos_theta, int Nl, bsComplex* bsComplex Ui((ptype)0, (ptype)0); //deal with rtsPoints near zero - if(real(knd) < EPSILON_FLOAT) - { - //for(int l=0; l(knd, j); - - //initialize the Legendre polynomial - ptype P[2]; - rts::init_legendre(cos_theta, P[0], P[1]); - - //for each order up to Nl - for(int l=0; l<=Nl; l++) - { - if(l == 0) - { - Ui += A[0] * j[0] * P[0]; - } - else - { - //shift the bessel functions and legendre polynomials - if(l > 1) - { - rts::shift_legendre(l, cos_theta, P[0], P[1]); - rts::shift_sbesselj(l, knd, j); - } - - Ui += A[l] * j[1] * P[1]; - - - } + //initialize the spherical Bessel functions + bsComplex j[2]; + rts::init_sbesselj(knd, j); + + //initialize the Legendre polynomial + ptype P[2]; + rts::init_legendre(cos_theta, P[0], P[1]); + + //for each order up to Nl + for(int l=0; l<=Nl; l++) + { + if(l == 0) + { + Ui += A[0] * j[0] * P[0]; + } + else + { + //shift the bessel functions and legendre polynomials + if(l > 1) + { + rts::shift_legendre(l, cos_theta, P[0], P[1]); + rts::shift_sbesselj(l, knd, j); + } + + Ui += A[l] * j[1] * P[1]; + + + } } return Ui; } @@ -118,39 +118,39 @@ __global__ void gpuScalarUsp(bsComplex* Us, bsVector* k, int nk, ptype kmag, bsP //get the rtsPoint in world space and then the r vector bsPoint p = ABCD(u, v); bsVector r = p - ps; - ptype d = r.len(); - - bsComplex sumUs(0, 0); - //for each plane wave in the wave list - for(int iw = 0; iw < nk; iw++) - { - //normalize the direction vectors and find their inner product - r = r.norm(); - ptype cos_theta = k[iw].dot(r); - - //compute the phase factor for spheres that are not at the origin - bsVector c = ps - f; - bsComplex phase = exp(bsComplex(0, kmag * k[iw].dot(c))); - - //compute the internal field if we are inside a sphere + ptype d = r.len(); + + bsComplex sumUs(0, 0); + //for each plane wave in the wave list + for(int iw = 0; iw < nk; iw++) + { + //normalize the direction vectors and find their inner product + r = r.norm(); + ptype cos_theta = k[iw].dot(r); + + //compute the phase factor for spheres that are not at the origin + bsVector c = ps - f; + bsComplex phase = exp(bsComplex(0, kmag * k[iw].dot(c))); + + //compute the internal field if we are inside a sphere if(d <= a) { bsComplex knd = kmag * d * n; - sumUs += (1.0/nk) * A * phase * calc_Ui(knd, cos_theta, Nl, Alpha); - } - //otherwise compute the scattered field - else - { - //compute the argument for the spherical Hankel function - ptype kd = kmag * d; - sumUs += (1.0/nk) * A * phase * calc_Us(kd, cos_theta, Nl, Beta); - } - - } - - Us[i] += sumUs; - - + sumUs += (1.0f/nk) * A * phase * calc_Ui(knd, cos_theta, Nl, Alpha); + } + //otherwise compute the scattered field + else + { + //compute the argument for the spherical Hankel function + ptype kd = kmag * d; + sumUs += (1.0f/nk) * A * phase * calc_Us(kd, cos_theta, Nl, Beta); + } + + } + + Us[i] += sumUs; + + } void nearfieldStruct::scalarUs() @@ -190,17 +190,17 @@ void nearfieldStruct::scalarUs() HANDLE_ERROR(cudaMemcpy(gpuB, &sVector[s].B[0], (Nl+1) * sizeof(bsComplex), cudaMemcpyHostToDevice)); HANDLE_ERROR(cudaMemcpy(gpuA, &sVector[s].A[0], (Nl+1) * sizeof(bsComplex), cudaMemcpyHostToDevice)); - //if we are computing a plane wave, call the gpuScalarUfp function - sphere S = sVector[s]; - bsVector* gpuk; + //if we are computing a plane wave, call the gpuScalarUfp function + sphere S = sVector[s]; + bsVector* gpuk; if(planeWave) - { - //if this is a single plane wave, assume it goes along direction k (copy the k vector to the GPU) - HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) )); + { + //if this is a single plane wave, assume it goes along direction k (copy the k vector to the GPU) + HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) )); HANDLE_ERROR(cudaMemcpy( gpuk, &k, sizeof(bsVector), cudaMemcpyHostToDevice)); gpuScalarUsp<<>>(U.x_hat, - gpuk, + gpuk, 1, 2 * PI / lambda, focus, @@ -213,20 +213,20 @@ void nearfieldStruct::scalarUs() A, pos, U.R[0], - U.R[1]); + U.R[1]); HANDLE_ERROR(cudaFree(gpuk)); - } - //otherwise copy all of the monte-carlo samples to the GPU and compute - else - { - HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) * inWaves.size() )); - HANDLE_ERROR(cudaMemcpy( gpuk, &inWaves[0], sizeof(bsVector) * inWaves.size(), cudaMemcpyHostToDevice)); - - //compute the amplitude that makes it through the condenser - ptype subA = 2 * PI * A * ( (1 - cos(asin(condenser[1]))) - (1 - cos(asin(condenser[0]))) ); - + } + //otherwise copy all of the monte-carlo samples to the GPU and compute + else + { + HANDLE_ERROR(cudaMalloc( (void**)&gpuk, sizeof(bsVector) * inWaves.size() )); + HANDLE_ERROR(cudaMemcpy( gpuk, &inWaves[0], sizeof(bsVector) * inWaves.size(), cudaMemcpyHostToDevice)); + + //compute the amplitude that makes it through the condenser + ptype subA = 2 * PI * A * ( (1 - cos(asin(condenser[1]))) - (1 - cos(asin(condenser[0]))) ); + gpuScalarUsp<<>>(U.x_hat, - gpuk, + gpuk, inWaves.size(), 2 * PI / lambda, focus, @@ -239,17 +239,17 @@ void nearfieldStruct::scalarUs() subA, pos, U.R[0], - U.R[1]); - HANDLE_ERROR(cudaFree(gpuk)); - - - } - - //free memory for scattering coefficients - HANDLE_ERROR(cudaFree(gpuA)); + U.R[1]); + HANDLE_ERROR(cudaFree(gpuk)); + + + } + + //free memory for scattering coefficients + HANDLE_ERROR(cudaFree(gpuA)); HANDLE_ERROR(cudaFree(gpuB)); - } - + } + //store the time to compute the scattered field t_Us = gpuStopTimer(); diff --git a/options.h b/options.h deleted file mode 100644 index 3aa53ea..0000000 --- a/options.h +++ /dev/null @@ -1,609 +0,0 @@ -//AnyOption for command-line processing -//#include "anyoption.h" - -#include "rts/optics/material.h" - -#include "nearfield.h" -#include "microscope.h" -#include "rts/visualization/colormap.h" -#include "fileout.h" -//extern nearfieldStruct* NF; -extern microscopeStruct* SCOPE; -extern fileoutStruct gFileOut; - -//default values -#include "defaults.h" - -#include -#include -#include -#include -using namespace std; - -#include -namespace po = boost::program_options; - -extern bool verbose; -extern bool gui; - - - -static void lNearfield(po::variables_map vm) -{ - //test to see if we are running a vector field simulation - bool vectorField = false; - if(vm.count("vector")) - vectorField = true; - SCOPE->scalarSim = !vectorField; - - //test to see if we are simulating a plane wave - bool planeWave = DEFAULT_PLANEWAVE; - if(vm.count("plane-wave")) - planeWave = !planeWave; - SCOPE->nf.planeWave = planeWave; - - //get the incident field amplitude - SCOPE->nf.A = vm["amplitude"].as(); - - //get the condenser parameters - SCOPE->nf.condenser[0] = DEFAULT_CONDENSER_MIN; - SCOPE->nf.condenser[1] = DEFAULT_CONDENSER_MAX; - - if(vm.count("condenser")) - { - vector cparams = vm["condenser"].as< vector >(); - - if(cparams.size() == 1) - SCOPE->nf.condenser[1] = cparams[0]; - else - { - SCOPE->nf.condenser[0] = cparams[0]; - SCOPE->nf.condenser[1] = cparams[1]; - } - } - - - //get the focal rtsPoint position - SCOPE->nf.focus[0] = DEFAULT_FOCUS_X; - SCOPE->nf.focus[1] = DEFAULT_FOCUS_Y; - SCOPE->nf.focus[2] = DEFAULT_FOCUS_Z; - if(vm.count("focus")) - { - vector fpos = vm["focus"].as< vector >(); - if(fpos.size() != 3) - { - cout<<"BIMSIM Error - the incident focal point is incorrectly specified; it must have three components."<nf.focus[0] = fpos[0]; - SCOPE->nf.focus[1] = fpos[1]; - SCOPE->nf.focus[2] = fpos[2]; - } - - //get the incident light direction (k-vector) - bsVector spherical(1, 0, 0); - - //if a k-vector is specified - if(vm.count("k")) - { - vector kvec = vm["k"].as< vector >(); - if(kvec.size() != 2) - { - cout<<"BIMSIM Error - k-vector is not specified correctly: it must contain two elements"<nf.k = spherical.sph2cart(); - - - //incident field order - SCOPE->nf.m = vm["field-order"].as(); - - //number of Monte-Carlo samples - SCOPE->nf.nWaves = vm["samples"].as(); - - //random number seed for Monte-Carlo samples - if(vm.count("seed")) - srand(vm["seed"].as()); - - - -} - - -static void loadOutputParams(po::variables_map vm) -{ - //append simulation results to previous binary files - gFileOut.append = DEFAULT_APPEND; - if(vm.count("append")) - gFileOut.append = true; - - //image parameters - //component of the field to be saved - std::string fieldStr; - fieldStr = vm["output-type"].as(); - - if(fieldStr == "magnitude") - gFileOut.field = fileoutStruct::fieldMag; - else if(fieldStr == "intensity") - gFileOut.field = fileoutStruct::fieldIntensity; - else if(fieldStr == "polarization") - gFileOut.field = fileoutStruct::fieldPolar; - else if(fieldStr == "imaginary") - gFileOut.field = fileoutStruct::fieldImag; - else if(fieldStr == "real") - gFileOut.field = fileoutStruct::fieldReal; - else if(fieldStr == "angular-spectrum") - gFileOut.field = fileoutStruct::fieldAngularSpectrum; - - - //image file names - gFileOut.intFile = vm["intensity"].as(); - gFileOut.absFile = vm["absorbance"].as(); - gFileOut.transFile = vm["transmittance"].as(); - gFileOut.nearFile = vm["near-field"].as(); - gFileOut.farFile = vm["far-field"].as(); - - //colormap - std::string cmapStr; - cmapStr = vm["colormap"].as(); - if(cmapStr == "brewer") - gFileOut.colormap = rts::cmBrewer; - else if(cmapStr == "gray") - gFileOut.colormap = rts::cmGrayscale; - else - cout<<"color-map value not recognized (using default): "<nf.lut_us = false; - SCOPE->nf.lut_uf = false; - } - else if(vm.count("recursive-us")) - { - SCOPE->nf.lut_us = false; - } - else if(vm.count("lut-uf")) - { - SCOPE->nf.lut_uf = true; - } - - //gui - if(vm.count("gui")) - gui = true; -} - -void lWavelength(po::variables_map vm) -{ - //load the wavelength - if(vm.count("nu")) - { - //wavelength is given in wavenumber - transform and flag - SCOPE->nf.lambda = 10000/vm["nu"].as(); - gFileOut.wavenumber = true; - } - //otherwise we are using lambda = wavelength - else - { - SCOPE->nf.lambda = vm["lambda"].as(); - gFileOut.wavenumber = false; - } -} - -static void lSpheres(string sphereList) -{ - /*This function loads a list of sphere given in the string sphereList - The format is: - x y z a m - where - x, y, z = sphere position (required) - a = sphere radius (required) - m = material ID (optional) */ - - std::stringstream ss(sphereList); - - while(!ss.eof()) - { - //create a new sphere - sphere newS; - - //get the sphere data - ss>>newS.p[0]; - ss>>newS.p[1]; - ss>>newS.p[2]; - ss>>newS.a; - - if(ss.peek() != '\n') - ss>>newS.iMaterial; - - //add the new sphere to the sphere vector - SCOPE->nf.sVector.push_back(newS); - - //ignore the rest of the line - ss.ignore(std::numeric_limits::max(), '\n'); - - //check out the next element (this should set the EOF error flag) - ss.peek(); - } -} - -void lSpheres(po::variables_map vm) -{ - //if a sphere is specified at the command line - if(vm.count("spheres")) - { - //convert the sphere to a string - vector sdesc = vm["spheres"].as< vector >(); - - //compute the number of spheres specified - unsigned int nS; - if(sdesc.size() <= 5) - nS = 1; - else - { - //if the number of parameters is divisible by 4, compute the number of spheres - if(sdesc.size() % 5 == 0) - nS = sdesc.size() / 5; - else - { - cout<<"BIMSIM Error: Invalid number of sphere parameters."< filenames = vm["sphere-file"].as< vector >(); - //load each file - for(unsigned int iS=0; iS(ifs)), std::istreambuf_iterator()); - - //load the list of spheres from a string - lSpheres(instr); - } - } - - //make sure the appropriate materials are loaded - unsigned int nS = SCOPE->nf.sVector.size(); - - //for each sphere - for(unsigned int s = 0; snf.sVector[s].iMaterial + 1 > SCOPE->nf.mVector.size()) - { - //otherwise output an error - cout<<"BIMSIM Error - A material is not loaded for sphere "<nf.sVector[s].iMaterial + 1<nf.mVector.size()< matVec = vm["materials"].as< vector >(); - if(matVec.size() == 1) - { - rts::material newM(SCOPE->nf.lambda, matVec[0], 0); - SCOPE->nf.mVector.push_back(newM); - } - else if(matVec.size() %2 != 0) - { - cout<<"BIMSim Error: materials must be specified in n, k pairs"< newM(SCOPE->nf.lambda, matVec[i], matVec[i+1]); - SCOPE->nf.mVector.push_back(newM); - } - } - } - - //if file names are specified, load the materials - if(vm.count("material-file")) - { - vector filenames = vm["material-file"].as< vector >(); - for(unsigned int i=0; i(ifs)), std::istreambuf_iterator()); - - //load the list of spheres from a string - rts::material newM(filenames[i].c_str()); - //newM.fromStr(instr, ""); - SCOPE->nf.mVector.push_back(newM); - } - } - -} - -static void lOptics(po::variables_map vm) -{ - SCOPE->objective[0] = DEFAULT_OBJECTIVE_MIN; - SCOPE->objective[1] = DEFAULT_OBJECTIVE_MAX; - if(vm.count("objective")) - { - vector oparams = vm["objective"].as< vector >(); - - if(oparams.size() == 1) - SCOPE->objective[1] = oparams[0]; - else - { - SCOPE->objective[0] = oparams[0]; - SCOPE->objective[1] = oparams[1]; - } - } -} - -static void lImagePlane(po::variables_map vm) -{ - bsPoint pMin(DEFAULT_PLANE_MIN_X, DEFAULT_PLANE_MIN_Y, DEFAULT_PLANE_MIN_Z); - bsPoint pMax(DEFAULT_PLANE_MAX_X, DEFAULT_PLANE_MAX_Y, DEFAULT_PLANE_MAX_Z); - bsVector normal(DEFAULT_PLANE_NORM_X, DEFAULT_PLANE_NORM_Y, DEFAULT_PLANE_NORM_Z); - - //set the default values for the slice position and orientation - if(vm.count("plane-lower-left") && vm.count("plane-upper-right") && vm.count("plane-normal")) - { - vector ll = vm["plane-lower-left"].as< vector >(); - if(ll.size() != 3) - { - cout<<"BIMSIM Error - The lower-left corner of the image plane is incorrectly specified."< ur = vm["plane-lower-left"].as< vector >(); - if(ur.size() != 3) - { - cout<<"BIMSIM Error - The upper-right corner of the image plane is incorrectly specified."< norm = vm["plane-lower-left"].as< vector >(); - if(norm.size() != 3) - { - cout<<"BIMSIM Error - The normal of the image plane is incorrectly specified."< xy = vm["xy"].as< vector >(); - if(xy.size() >= 1) - s = xy[0]; - if(xy.size() >= 2) - pos = xy[1]; - - //calculate the plane corners and normal based on the size and position - pMin = bsPoint(-s/2, -s/2, pos); - pMax = bsPoint(s/2, s/2, pos); - normal = bsVector(0, 0, 1); - } - else if(vm.count("xz")) - { - //default plane size in microns - ptype size = DEFAULT_PLANE_SIZE; - ptype pos = DEFAULT_PLANE_POSITION; - - vector xz = vm["xz"].as< vector >(); - if(xz.size() >= 1) - size = xz[0]; - if(xz.size() >= 2) - pos = xz[1]; - - //calculate the plane corners and normal based on the size and position - pMin = bsPoint(-size/2, pos, -size/2); - pMax = bsPoint(size/2, pos, size/2); - normal = bsVector(0, -1, 0); - } - else if(vm.count("yz")) - { - //default plane size in microns - ptype size = DEFAULT_PLANE_SIZE; - ptype pos = DEFAULT_PLANE_POSITION; - - vector yz = vm["yz"].as< vector >(); - if(yz.size() >= 1) - size = yz[0]; - if(yz.size() >= 2) - pos = yz[1]; - - //calculate the plane corners and normal based on the size and position - pMin = bsPoint(pos, -size/2, -size/2); - pMax = bsPoint(pos, size/2, size/2); - normal = bsVector(1, 0, 0); - } - SCOPE->setPos(pMin, pMax, normal); - - //resolution - SCOPE->setRes(vm["resolution"].as(), - vm["resolution"].as(), - vm["padding"].as(), - vm["supersample"].as()); - - - - - - SCOPE->setNearfield(); -} - -static void OutputOptions() -{ - cout<toStr(); - - cout<<"# of source points: "<focalPoints.size()< test; -static void SetOptions(po::options_description &desc) -{ - desc.add_options() - ("help", "prints this help") - ("gui", "run using the Qt GUI") - ("verbose", "verbose output\n\nOutput Parameters\n--------------------------") - - ("vector", "run a vector field simulation") - ("intensity", po::value()->default_value(DEFAULT_INTENSITY_FILE), "output measured intensity (filename)") - ("absorbance", po::value()->default_value(DEFAULT_ABSORBANCE_FILE), "output measured absorbance (filename)") - ("transmittance", po::value()->default_value(DEFAULT_TRANSMITTANCE_FILE), "output measured transmittance (filename)") - ("far-field", po::value()->default_value(DEFAULT_FAR_FILE), "output far-field at detector (filename)") - ("near-field", po::value()->default_value(DEFAULT_NEAR_FILE), "output field at focal plane (filename)") - ("extended-source", po::value()->default_value(DEFAULT_EXTENDED_SOURCE), "image of source at focus (filename)") - ("output-type", po::value()->default_value(DEFAULT_FIELD_TYPE), "output field value:\n magnitude, polarization, real, imaginary, angular-spectrum") - ("colormap", po::value()->default_value(DEFAULT_COLORMAP), "colormap: gray, brewer") - ("append", "append result to an existing file\n (binary files only)\n\nSphere Parameters\n--------------------------") - - ("spheres", po::value< vector >()->multitoken(), "sphere position: x y z a m") - ("sphere-file", po::value< vector >()->multitoken(), "sphere file:\n [x y z radius material]") - ("materials", po::value< vector >()->multitoken(), "refractive indices as n, k pairs:\n ex. -m n0 k0 n1 k1 n2 k2") - ("material-file", po::value< vector >()->multitoken(), "material file:\n [lambda n k]\n\nOptics\n--------------------------") - - ("lambda", po::value()->default_value(DEFAULT_LAMBDA), "incident wavelength") - ("nu", po::value(), "incident frequency (in cm^-1)\n(if specified, lambda is ignored)") - ("k", po::value< vector >()->multitoken(), "k-vector direction: -k theta phi\n theta = [0 2*pi], phi = [0 pi]") - ("amplitude", po::value()->default_value(DEFAULT_AMPLITUDE), "incident field amplitude") - ("condenser", po::value< vector >()->multitoken(), "condenser numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout") - ("objective", po::value< vector >()->multitoken(), "objective numerical aperature\nA pair of values can be used to specify an inner obscuration: -c NAin NAout") - ("focus", po::value< vector >()->multitoken(), "focal position for the incident point source\n (default = --focus 0 0 0)") - ("plane-wave", "simulates an incident plane wave\n\n\nImaging Parameters\n--------------------------") - - ("resolution", po::value()->default_value(DEFAULT_SLICE_RES), "resolution of the detector") - ("plane-lower-left", po::value< vector >()->multitoken(), "lower-left position of the image plane") - ("plane-upper-right", po::value< vector >()->multitoken(), "upper-right position of the image plane") - ("plane-normal", po::value< vector >()->multitoken(), "normal for the image plane") - ("xy", po::value< vector >()->multitoken(), "specify an x-y image plane\n (standard microscope)") - ("xz", po::value< vector >()->multitoken(), "specify a x-z image plane\n (cross-section of the focal volume)") - ("yz", po::value< vector >()->multitoken(), "specify a y-z image plane\n (cross-section of the focal volume)\n\nSampling Parameters\n--------------------------") - - ("samples", po::value()->default_value(DEFAULT_SAMPLES), "Monte-Carlo samples used to compute Us") - ("padding", po::value()->default_value(DEFAULT_PADDING), "FFT padding for the objective bandpass") - ("supersample", po::value()->default_value(DEFAULT_SUPERSAMPLE), "super-sampling rate for the detector field") - ("field-order", po::value()->default_value(DEFAULT_FIELD_ORDER), "order of the incident field") - ("seed", po::value(), "seed for the Monte-Carlo random number generator") - ("recursive", "evaluate all Bessel functions recursively\n") - ("recursive-us", "evaluate scattered-field Bessel functions recursively\n") - ("lut-uf", "evaluate the focused-field using a look-up table\n") - ; -} - -static void LoadParameters(int argc, char *argv[]) -{ - //create an option description - po::options_description desc("BimSim arguments"); - - //fill it with options - SetOptions(desc); - - po::variables_map vm; - po::store(po::parse_command_line(argc, argv, desc, po::command_line_style::unix_style ^ po::command_line_style::allow_short), vm); - po::notify(vm); - - - //load flags (help, verbose output) - lFlags(vm, desc); - - //load the wavelength - lWavelength(vm); - - //load materials - lMaterials(vm); - - //load the sphere data - lSpheres(vm); - - //load the optics - lOptics(vm); - - //load the position and orientation of the image plane - lImagePlane(vm); - - - lNearfield(vm); - - loadOutputParams(vm); - - //if an extended source will be used - if(vm["extended-source"].as() != "") - { - //load the point sources - string filename = vm["extended-source"].as(); - SCOPE->LoadExtendedSource(filename); - - } - - - - - -} diff --git a/planewave.h b/planewave.h index 57437f0..5e1fb43 100644 --- a/planewave.h +++ b/planewave.h @@ -38,7 +38,7 @@ class planewave } // refract will bend the wave vector k to correspond to the normalized vector v - refract(rts::vector v) + void refract(rts::vector v) { //make sure that v is normalized v = v.norm(); @@ -46,25 +46,25 @@ class planewave } - std::string toStr() - { - std::stringstream ss; - + std::string toStr() + { + std::stringstream ss; + ss<<"k = "< @@ -69,7 +70,7 @@ struct sphere void calcNl(ptype lambda) { - Nl = ceil( (2 * PI * a) / lambda + 4 * pow( (2 * PI * a) / lambda, 1.0/3.0) + 2); + Nl = ceil( (2 * PI * a) / lambda + 4 * pow( (2 * PI * a) / lambda, (ptype)(1.0/3.0)) + 2); } //compute the scattering coefficients -- libgit2 0.21.4