#ifndef MATERIALSTRUCT_H #define MATERIALSTRUCT_H #include #include #include #include #include #include #include #include "rts/math/complex.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; } template static T _2wn(T lambda) { return (T)10000.0/lambda; } //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; } //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< 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; } 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 j(0, 1); for(int i=0; i N/2, multiply by -i else Chi2FFT[i] *= -j; } //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(); #endif } 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]; } complex getN(T l) { //return complex(1.0, 0.0); complex ri = dispersion.linear(l) + n0; return ri; } function > getF() { return dispersion + complex(n0, 0.0); } //returns the scattering efficiency at wavelength l complex eta(T l) { //get the refractive index complex ri = getN(l); //convert the refractive index to scattering efficiency return ri*ri - 1.0; } //interpolate the given lambda value and return the index of refraction complex operator()(T l) { return getN(l); } }; } //end namespace rts template std::ostream& operator<<(std::ostream& os, rts::material m) { os<