From 0ef519a4a62c88829271fa984c25fd465b707890 Mon Sep 17 00:00:00 2001 From: David Mayerich Date: Thu, 7 Aug 2014 14:16:43 -0500 Subject: [PATCH] optimized material class, still have to add KK --- math/complex.h | 27 +++++++++++++++++++++++++++ math/function.h | 57 +++++++++++++++++++++++++++++++++++++++++++++------------ optics/material.h | 663 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- optics/mirst-1d.cuh | 19 +++++++++++-------- ui/arguments.h | 17 ++++++++++------- 5 files changed, 199 insertions(+), 584 deletions(-) diff --git a/math/complex.h b/math/complex.h index 5430956..2473ff8 100644 --- a/math/complex.h +++ b/math/complex.h @@ -294,6 +294,19 @@ struct complex return false; } + CUDA_CALLABLE bool operator<(complex rhs){ + return abs() < rhs.abs(); + } + CUDA_CALLABLE bool operator<=(complex rhs){ + return abs() <= rhs.abs(); + } + CUDA_CALLABLE bool operator>(complex rhs){ + return abs() > rhs.abs(); + } + CUDA_CALLABLE bool operator >=(complex rhs){ + return abs() >= rhs.abs(); + } + //CASTING operators template < typename otherT > operator complex() @@ -469,6 +482,20 @@ std::ostream& operator<<(std::ostream& os, rts::complex x) return os; } +template +std::istream& operator>>(std::istream& is, rts::complex& x) +{ + A r, i; + r = i = 0; //initialize the real and imaginary parts to zero + is>>r; //parse + is>>i; + + x.real(r); //assign the parsed values to x + x.imag(i); + + return is; //return the stream +} + //#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 //template using rtsComplex = rts::complex; //#endif diff --git a/math/function.h b/math/function.h index 11655be..97a9c69 100644 --- a/math/function.h +++ b/math/function.h @@ -2,6 +2,7 @@ #define RTS_FUNCTION_H #include +#include namespace rts{ @@ -9,8 +10,12 @@ namespace rts{ template class function { +protected: + std::vector X; std::vector Y; + Ty range[2]; + Ty bounding[2]; //comparison function for searching lambda static bool findCeiling(Tx ax, Tx bx){ @@ -21,11 +26,6 @@ class function void process_string(std::string s){ std::stringstream ss(s); - //std::string test; - //std::getline(ss, test); - //std::cout<>x; //read the x value lstream>>y; //read the y value - //std::cout< range[1]) range[1] = y; + } + public: + function(){ + range[0] = range[1] = NAN; + bounding[0] = bounding[1] = 0; + } //linear interpolation Ty linear(Tx x) const { - if(X.size() == 0) return (Ty)0; //return zero if the function is empty + if(X.size() == 0) return bounding[0]; //return zero if the function is empty + //return bounding values if x is outside the sample domain + if(x < X[0]) return bounding[0]; + if(x > X.back()) return bounding[1]; unsigned int N = X.size(); //number of sample points @@ -92,6 +107,8 @@ public: { unsigned int N = X.size(); //number of sample points + update_range(y); //update the range of the function + if(N == 0 || X[N-1] < x){ X.push_back(x); Y.push_back(y); @@ -172,14 +189,30 @@ public: } - + //output a string description of the function (used for console output and debugging) std::string str(){ - stringstream ss; + unsigned int N = X.size(); + + std::stringstream ss; + ss<<"# of samples: "< #include @@ -9,596 +9,145 @@ #include #include #include "rts/math/complex.h" +#include "rts/math/constants.h" #include "rts/math/function.h" -#define PI 3.14159f - namespace rts{ - enum field_type {field_microns, field_wavenumber, field_n, field_k, field_A, field_ignore}; - - //conversion functions - - //convert wavenumber to lambda - template - static T _wn(T inverse_cm) - { - return (T)10000.0/inverse_cm; - } +//Material class - default representation for the material property is the refractive index (RI) +template +class material : public function< T, complex >{ - template - static T _2wn(T lambda) - { - return (T)10000.0/lambda; - } +public: + enum wave_property{microns, inverse_cm}; + enum material_property{ri, absorbance}; - //convert absorbance to k - template - static T _A(T absorbance, T lambda) - { - return (absorbance * lambda) / (4 * PI); - } - template - static T _2A(T k, T lambda) - { - return (4 * PI * k)/lambda; - } +private: - //define the dispersion as a single wavelength/refractive index pair - template - struct refIndex - { - //wavelength (in microns) - T lambda; - complex n; - }; - - template - struct entryType - { - //list of value types per entry - std::vector valueList; - - entryType(std::string format) - { - //location of the end of a parameter - size_t e; - - //string storing a token - std::string token; - - do - { - //find the end of the first parameter - e = format.find_first_of(','); - - //get the substring up to the comma - token = format.substr(0, e); - - //turn the token into a val_type - if(token == "microns") - valueList.push_back(field_microns); - else if(token == "wavenumber") - valueList.push_back(field_wavenumber); - else if(token == "n") - valueList.push_back(field_n); - else if(token == "k") - valueList.push_back(field_k); - else if(token == "A") - valueList.push_back(field_A); - else - valueList.push_back(field_ignore); - - //remove the first token from the format string - format = format.substr(e+1, format.length()-1); - }while(e != std::string::npos); - - - - } - - void addValue(field_type value) - { - valueList.push_back(value); - } - - refIndex inputEntry(std::string line, T scaleA = 1.0) - { - T val; - std::stringstream ss(line); - - //create a new refractive index - refIndex newRI; - - - //read the entry from an input string - for(unsigned int i=0; i '9') - { - //cout<<"bad char"<>val; - //cout< material) - { - //std::string result; - std::stringstream ss; - - //for each field in the entry - for(int i=0; i 0) - ss<<"\t"; - //store the value in the appropriate location - switch(valueList[i]) - { - case field_microns: - ss< >::X; + using function< T, complex >::Y; + using function< T, complex >::insert; + using function< T, complex >::bounding; - } - return ss.str(); - } + std::string name; //name for the material (defaults to file name) + void process_header(std::string str, wave_property& wp, material_property& mp){ - }; - - - //a material is a list of refractive index values - template - class material - { - std::string name; - //dispersion (refractive index as a function of wavelength) - //std::vector< refIndex > dispersion; - function< T, complex > dispersion; - - //average refractive index (approximately 1.4) - T n0; - - /*void add(refIndex ri) - { - //refIndex converted = convert(ri, measurement); - dispersion.push_back(ri); - }*/ - - //comparison function for sorting - static bool compare(refIndex a, refIndex b) - { - return (a.lambda < b.lambda); - } - - //comparison function for searching lambda - /*static bool findCeiling(refIndex a, refIndex b) - { - return (a.lambda > b.lambda); - }*/ - - public: - void add(T lambda, complex n) - { - dispersion.insert(lambda, n); - } - - std::string getName() - { - return name; - } - void setName(std::string n) - { - name = n; - } - T getN0() - { - return n0; - } - void setN0(T n) - { - n0 = n; - } + std::stringstream ss(str); //create a stream from the data string + std::string line; + std::getline(ss, line); //get the first line as a string + while(line[0] == '#'){ //continue looping while the line is a comment - void setM(function< T, complex > m) - { - dispersion = m; - } - unsigned int nSamples() - { - return dispersion.size(); - } - material computeN(T _n0, unsigned int n_samples = 0, T pf = 2) - { - /* This function computes the real part of the refractive index - from the imaginary part. The Hilbert transform is required. I - use an FFT in order to simplify this, so either the FFTW or CUFFT - packages are required. CUFFT is used if this file is passed through - a CUDA compiler. Otherwise, FFTW is used if available. - */ - - n0 = _n0; - - int N; - if(n_samples) - N = n_samples; - else - N = dispersion.size(); - - -#ifdef FFTW_AVAILABLE - //allocate memory for the FFT - complex* Chi2 = (complex*)fftw_malloc(sizeof(complex) * N * pf); - complex* Chi2FFT = (complex*)fftw_malloc(sizeof(complex) * N * pf); - complex* Chi1 = (complex*)fftw_malloc(sizeof(complex) * N * pf); - - //create an FFT plan for the forward and inverse transforms - fftw_plan planForward, planInverse; - planForward = fftw_plan_dft_1d(N*pf, (fftw_complex*)Chi2, (fftw_complex*)Chi2FFT, FFTW_FORWARD, FFTW_ESTIMATE); - planInverse = fftw_plan_dft_1d(N*pf, (fftw_complex*)Chi2FFT, (fftw_complex*)Chi1, FFTW_BACKWARD, FFTW_ESTIMATE); - - float k, alpha; - T chi_temp; - - //the spectrum will be re-sampled in uniform values of wavenumber - T nuMin = _2wn(dispersion.back().lambda); - T nuMax = _2wn(dispersion.front().lambda); - T dnu = (nuMax - nuMin)/(N-1); - T lambda, tlambda; - for(int i=0; i nMin = dispersion.back(); - //complex nMax = dispersion.front(); - T a; - for(int i=N; i>prop; - Chi2[i] = 0.0;//Chi2[N-1]; + if(prop == "X"){ + std::string wp_name; + lstream>>wp_name; + if(wp_name == "microns") wp = microns; + else if(wp_name == "inverse_cm") wp = inverse_cm; } - - //perform the FFT - fftw_execute(planForward); - - //perform the Hilbert transform in the Fourier domain - complex j(0, 1); - for(int i=0; i N/2, multiply by -i - else - Chi2FFT[i] *= -j; + else if(prop == "Y"){ + std::string mp_name; + lstream>>mp_name; + if(mp_name == "ri") mp = ri; + else if(mp_name == "absorbance") mp = absorbance; } - //execute the inverse Fourier transform (completing the Hilbert transform) - fftw_execute(planInverse); - - //divide the Chi1 values by N - for(int i=0; i newM; - newM.dispersion.clear(); - refIndex ri; - for(int i=0; i >::process_string(str); + } - ofstream outN("n.txt"); - for(int i=0; i(); -#endif - } + void init(){ + bounding[0] = bounding[1] = rts::complex(1, 0); + } - material(T lambda = 1, T n = 1, T k = 0) - { - dispersion.insert(lambda, complex(n, k)); - /*//create a default refractive index - refIndex def; - def.lambda = lambda; - def.n.real(n); - def.n.imag(k); - add(def); - */ - //set n0 - n0 = 0; - } - - material(std::string filename, std::string format = "", T scaleA = 1.0) - { - fromFile(filename, format); - n0 = 0; - } - - void fromFile(std::string filename, std::string format = "", T scaleA = 1.0) - { - name = filename; - //clear any previous values - dispersion = rts::function< T, complex >(); - - //load the file into a string - std::ifstream ifs(filename.c_str()); - - std::string line; - - if(!ifs.is_open()) - { - std::cout<<"Material Error -- file not found: "<(ifs)), std::istreambuf_iterator()); - fromStr(instr, format, scaleA); - } - - void fromStr(std::string str, std::string format = "", T scaleA = 1.0) - { - //create a string stream to process the input data - std::stringstream ss(str); - - //this string will be read line-by-line (where each line is an entry) - std::string line; - - //if the format is not provided, see if it is in the file, otherwise use a default - if(format == "") - { - //set a default format of "lambda,n,k" - format = "microns,n,k"; - - //see if the first line is a comment - char c = ss.peek(); - if(c == '#') - { - //get the first line - getline(ss, line); - //look for a bracket, denoting the start of a format string - int istart = line.find('['); - if(istart != std::string::npos) - { - //look for a bracket terminating the format string - int iend = line.find(']'); - if(iend != std::string::npos) - { - //read the string between the brackets - format = line.substr(istart+1, iend - istart - 1); - } - } - } - - } - - entryType entry(format); - - //std::cout<<"Loading material with format: "<::compare); - } - - //convert the material to a string - std::string toStr(std::string format = "microns,n,k", bool reverse_order = false) - { - std::stringstream ss; - entryType entry(format); - - if(reverse_order) - { - for(int l=dispersion.size() - 1; l>=0; l--) - { - if(l < dispersion.size() - 1) ss< 0) ss<& operator[](unsigned int i) - { - return dispersion[i]; +public: - } + material(std::string filename, wave_property wp, material_property mp){ + name = filename; + load(filename, wp, mp); + } - complex getN(T l) - { - //return complex(1.0, 0.0); - complex ri = dispersion.linear(l) + n0; - return ri; - } + material(std::string filename){ + name = filename; + load(filename); + } - function > getF() - { - return dispersion + complex(n0, 0.0); - } + material(){ + init(); + } - //returns the scattering efficiency at wavelength l - complex eta(T l) - { - //get the refractive index - complex ri = getN(l); + complex getN(T lambda){ + return function< T, complex >::linear(lambda); + } - //convert the refractive index to scattering efficiency - return ri*ri - 1.0; + void load(std::string filename, wave_property wp, material_property mp){ - } - //interpolate the given lambda value and return the index of refraction - complex operator()(T l) - { - return getN(l); - } + //load the file as a function + function< T, complex >::load(filename); + } + void load(std::string filename){ + + wave_property wp = inverse_cm; + material_property mp = ri; + //turn the file into a string + std::ifstream t(filename.c_str()); //open the file as a stream + + if(!t){ + std::cout<<"ERROR: Couldn't open the material file '"<(t)), + std::istreambuf_iterator()); + + //process the header information + process_header(str, wp, mp); + + //convert units + if(wp == inverse_cm) + from_inverse_cm(); + //set the bounding values + bounding[0] = Y[0]; + bounding[1] = Y.back(); + } + std::string str(){ + std::stringstream ss; + ss< >::str(); + return ss.str(); + } + std::string get_name(){ + return name; + } - }; -} //end namespace rts + void set_name(std::string str){ + name = str; + } -template -std::ostream& operator<<(std::ostream& os, rts::material m) -{ - os<* dest, complex* ri, //T l = w * ri[inu].real(); //complex e(0.0, -2 * PI * fz * (zf[inu] + zR/2 - l/2.0)); - complex e(0, -2 * PI * fz * (zf[inu] + opl/2)); + complex e(0, -2 * rtsPI * fz * (zf[inu] + opl/2)); complex eta = ri[inu] * ri[inu] - 1; @@ -59,7 +60,7 @@ __global__ void gpu_mirst1d_layer_fft(complex* dest, complex* ri, if(ifz == 0) dest[i] += opl * exp(e) * eta * src[inu]; else - dest[i] += opl * exp(e) * eta * src[inu] * sin(PI * fz * opl) / (PI * fz * opl); + dest[i] += opl * exp(e) * eta * src[inu] * sin(rtsPI * fz * opl) / (rtsPI * fz * opl); } template @@ -139,7 +140,7 @@ __global__ void gpu_mirst1d_apply_filter(complex* sampleFFT, T* lambda, //compute the final filter T mirst = 0; if(fz != 0) - mirst = 2 * PI * r * s * Q * (1/fz); + mirst = 2 * rtsPI * r * s * Q * (1/fz); sampleFFT[i] *= mirst; @@ -246,7 +247,7 @@ private: //store the refractive index and source profile in a CPU array for(int inu=0; inu), cudaMemcpyHostToDevice )); //save the source profile to the GPU @@ -320,6 +321,7 @@ public: pad = padding; NA[0] = 0; NA[1] = 0.8; + S = 0; } //add a layer, thickness = microns @@ -397,16 +399,17 @@ public: std::string str(){ stringstream ss; - + ss<<"1D MIRST Simulation========================="<