codebase / stimlib

Commit 0ef519a4a62c88829271fa984c25fd465b707890

Authored by David Mayerich
1 parent 7c246292

optimized material class, still have to add KK

Showing 5 changed files with 199 additions and 584 deletions   Show diff stats
Inline Side-by-side
@@ -294,6 +294,19 @@ struct complex @@ -294,6 +294,19 @@ struct complex
294 return false; 294 return false;
295 } 295 }
296 296
  297 + CUDA_CALLABLE bool operator<(complex<T> rhs){
  298 + return abs() < rhs.abs();
  299 + }
  300 + CUDA_CALLABLE bool operator<=(complex<T> rhs){
  301 + return abs() <= rhs.abs();
  302 + }
  303 + CUDA_CALLABLE bool operator>(complex<T> rhs){
  304 + return abs() > rhs.abs();
  305 + }
  306 + CUDA_CALLABLE bool operator >=(complex<T> rhs){
  307 + return abs() >= rhs.abs();
  308 + }
  309 +
297 //CASTING operators 310 //CASTING operators
298 template < typename otherT > 311 template < typename otherT >
299 operator complex<otherT>() 312 operator complex<otherT>()
@@ -469,6 +482,20 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, rts::complex&lt;A&gt; x) @@ -469,6 +482,20 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, rts::complex&lt;A&gt; x)
469 return os; 482 return os;
470 } 483 }
471 484
  485 +template<class A>
  486 +std::istream& operator>>(std::istream& is, rts::complex<A>& x)
  487 +{
  488 + A r, i;
  489 + r = i = 0; //initialize the real and imaginary parts to zero
  490 + is>>r; //parse
  491 + is>>i;
  492 +
  493 + x.real(r); //assign the parsed values to x
  494 + x.imag(i);
  495 +
  496 + return is; //return the stream
  497 +}
  498 +
472 //#if __GNUC__ > 3 && __GNUC_MINOR__ > 7 499 //#if __GNUC__ > 3 && __GNUC_MINOR__ > 7
473 //template<class T> using rtsComplex = rts::complex<T>; 500 //template<class T> using rtsComplex = rts::complex<T>;
474 //#endif 501 //#endif
@@ -2,6 +2,7 @@ @@ -2,6 +2,7 @@
2 #define RTS_FUNCTION_H 2 #define RTS_FUNCTION_H
3 3
4 #include <string> 4 #include <string>
  5 +#include <cmath>
5 6
6 namespace rts{ 7 namespace rts{
7 8
@@ -9,8 +10,12 @@ namespace rts{ @@ -9,8 +10,12 @@ namespace rts{
9 template <class Tx, class Ty> 10 template <class Tx, class Ty>
10 class function 11 class function
11 { 12 {
  13 +protected:
  14 +
12 std::vector<Tx> X; 15 std::vector<Tx> X;
13 std::vector<Ty> Y; 16 std::vector<Ty> Y;
  17 + Ty range[2];
  18 + Ty bounding[2];
14 19
15 //comparison function for searching lambda 20 //comparison function for searching lambda
16 static bool findCeiling(Tx ax, Tx bx){ 21 static bool findCeiling(Tx ax, Tx bx){
@@ -21,11 +26,6 @@ class function @@ -21,11 +26,6 @@ class function
21 void process_string(std::string s){ 26 void process_string(std::string s){
22 std::stringstream ss(s); 27 std::stringstream ss(s);
23 28
24 - //std::string test;  
25 - //std::getline(ss, test);  
26 - //std::cout<<test;  
27 - //exit(1);  
28 -  
29 Tx x; 29 Tx x;
30 Ty y; 30 Ty y;
31 std::string line; 31 std::string line;
@@ -40,21 +40,36 @@ class function @@ -40,21 +40,36 @@ class function
40 lstream>>x; //read the x value 40 lstream>>x; //read the x value
41 lstream>>y; //read the y value 41 lstream>>y; //read the y value
42 42
43 - //std::cout<<x<<", "<<y<<std::endl;  
44 -  
45 if(ss.eof()) break; 43 if(ss.eof()) break;
46 insert(x, y); //insert the read value into the function 44 insert(x, y); //insert the read value into the function
47 -  
48 } 45 }
49 } 46 }
50 47
  48 + void update_range(Ty y){
  49 + //if the function is empty, just set the range to the given value
  50 + if(X.size() == 0){
  51 + range[0] = range[1] = y;
  52 + }
  53 +
  54 + //otherwise compare the values and set them appropriately
  55 + if(y < range[0]) range[0] = y;
  56 + if(y > range[1]) range[1] = y;
  57 + }
  58 +
51 59
52 public: 60 public:
  61 + function(){
  62 + range[0] = range[1] = NAN;
  63 + bounding[0] = bounding[1] = 0;
  64 + }
53 65
54 //linear interpolation 66 //linear interpolation
55 Ty linear(Tx x) const 67 Ty linear(Tx x) const
56 { 68 {
57 - if(X.size() == 0) return (Ty)0; //return zero if the function is empty 69 + if(X.size() == 0) return bounding[0]; //return zero if the function is empty
  70 + //return bounding values if x is outside the sample domain
  71 + if(x < X[0]) return bounding[0];
  72 + if(x > X.back()) return bounding[1];
58 73
59 unsigned int N = X.size(); //number of sample points 74 unsigned int N = X.size(); //number of sample points
60 75
@@ -92,6 +107,8 @@ public: @@ -92,6 +107,8 @@ public:
92 { 107 {
93 unsigned int N = X.size(); //number of sample points 108 unsigned int N = X.size(); //number of sample points
94 109
  110 + update_range(y); //update the range of the function
  111 +
95 if(N == 0 || X[N-1] < x){ 112 if(N == 0 || X[N-1] < x){
96 X.push_back(x); 113 X.push_back(x);
97 Y.push_back(y); 114 Y.push_back(y);
@@ -172,14 +189,30 @@ public: @@ -172,14 +189,30 @@ public:
172 } 189 }
173 190
174 191
175 - 192 + //output a string description of the function (used for console output and debugging)
176 std::string str(){ 193 std::string str(){
177 - stringstream ss; 194 + unsigned int N = X.size();
  195 +
  196 + std::stringstream ss;
  197 + ss<<"# of samples: "<<N<<std::endl;
  198 +
  199 + if(N == 0) return ss.str();
  200 +
  201 + ss<<"domain: ["<<X[0]<<", "<<X[X.size()-1]<<"]"<<std::endl;
  202 + ss<<"range: ["<<range[0]<<", "<<range[1]<<"]"<<std::endl;
  203 +
  204 +
  205 + return ss.str();
  206 +
  207 + }
  208 + std::string str_vals(){
  209 + std::stringstream ss;
178 210
179 unsigned int N = X.size(); //number of sample points 211 unsigned int N = X.size(); //number of sample points
180 //output each function value 212 //output each function value
181 for(unsigned int i = 0; i<N; i++){ 213 for(unsigned int i = 0; i<N; i++){
182 - ss<<X[i]<<", "<<Y[i]<<std::endl; 214 + if(i != 0) ss<<std::endl;
  215 + ss<<X[i]<<", "<<Y[i];
183 } 216 }
184 217
185 return ss.str(); 218 return ss.str();
1 -#ifndef MATERIALSTRUCT_H  
2 -#define MATERIALSTRUCT_H 1 +#ifndef RTS_MATERIAL_H
  2 +#define RTS_MATERIAL_H
3 3
4 #include <vector> 4 #include <vector>
5 #include <ostream> 5 #include <ostream>
@@ -9,596 +9,145 @@ @@ -9,596 +9,145 @@
9 #include <algorithm> 9 #include <algorithm>
10 #include <sstream> 10 #include <sstream>
11 #include "rts/math/complex.h" 11 #include "rts/math/complex.h"
  12 +#include "rts/math/constants.h"
12 #include "rts/math/function.h" 13 #include "rts/math/function.h"
13 14
14 -#define PI 3.14159f  
15 -  
16 namespace rts{ 15 namespace rts{
17 16
18 - enum field_type {field_microns, field_wavenumber, field_n, field_k, field_A, field_ignore};  
19 -  
20 - //conversion functions  
21 -  
22 - //convert wavenumber to lambda  
23 - template <class T>  
24 - static T _wn(T inverse_cm)  
25 - {  
26 - return (T)10000.0/inverse_cm;  
27 - } 17 +//Material class - default representation for the material property is the refractive index (RI)
  18 +template<typename T>
  19 +class material : public function< T, complex<T> >{
28 20
29 - template <class T>  
30 - static T _2wn(T lambda)  
31 - {  
32 - return (T)10000.0/lambda;  
33 - } 21 +public:
  22 + enum wave_property{microns, inverse_cm};
  23 + enum material_property{ri, absorbance};
34 24
35 - //convert absorbance to k  
36 - template <class T>  
37 - static T _A(T absorbance, T lambda)  
38 - {  
39 - return (absorbance * lambda) / (4 * PI);  
40 - }  
41 - template <class T>  
42 - static T _2A(T k, T lambda)  
43 - {  
44 - return (4 * PI * k)/lambda;  
45 - } 25 +private:
46 26
47 - //define the dispersion as a single wavelength/refractive index pair  
48 - template <class T>  
49 - struct refIndex  
50 - {  
51 - //wavelength (in microns)  
52 - T lambda;  
53 - complex<T> n;  
54 - };  
55 -  
56 - template <class T>  
57 - struct entryType  
58 - {  
59 - //list of value types per entry  
60 - std::vector<field_type> valueList;  
61 -  
62 - entryType(std::string format)  
63 - {  
64 - //location of the end of a parameter  
65 - size_t e;  
66 -  
67 - //string storing a token  
68 - std::string token;  
69 -  
70 - do  
71 - {  
72 - //find the end of the first parameter  
73 - e = format.find_first_of(',');  
74 -  
75 - //get the substring up to the comma  
76 - token = format.substr(0, e);  
77 -  
78 - //turn the token into a val_type  
79 - if(token == "microns")  
80 - valueList.push_back(field_microns);  
81 - else if(token == "wavenumber")  
82 - valueList.push_back(field_wavenumber);  
83 - else if(token == "n")  
84 - valueList.push_back(field_n);  
85 - else if(token == "k")  
86 - valueList.push_back(field_k);  
87 - else if(token == "A")  
88 - valueList.push_back(field_A);  
89 - else  
90 - valueList.push_back(field_ignore);  
91 -  
92 - //remove the first token from the format string  
93 - format = format.substr(e+1, format.length()-1);  
94 - }while(e != std::string::npos);  
95 -  
96 -  
97 -  
98 - }  
99 -  
100 - void addValue(field_type value)  
101 - {  
102 - valueList.push_back(value);  
103 - }  
104 -  
105 - refIndex<T> inputEntry(std::string line, T scaleA = 1.0)  
106 - {  
107 - T val;  
108 - std::stringstream ss(line);  
109 -  
110 - //create a new refractive index  
111 - refIndex<T> newRI;  
112 -  
113 -  
114 - //read the entry from an input string  
115 - for(unsigned int i=0; i<valueList.size(); i++)  
116 - {  
117 -  
118 -  
119 - while(ss.peek() < '0' || ss.peek() > '9')  
120 - {  
121 - //cout<<"bad char"<<endl;  
122 - ss.ignore();  
123 - }  
124 -  
125 - //retrieve the value  
126 - ss>>val;  
127 - //cout<<val<<endl;  
128 -  
129 - //store the value in the appropriate location  
130 - switch(valueList[i])  
131 - {  
132 - case field_microns:  
133 - newRI.lambda = val;  
134 - break;  
135 - case field_wavenumber:  
136 - newRI.lambda = _wn(val);  
137 - break;  
138 - case field_n:  
139 - newRI.n.real(val);  
140 - break;  
141 - case field_k:  
142 - newRI.n.imag(val);  
143 - break;  
144 - case field_A:  
145 - newRI.n.imag(_A(val * scaleA, newRI.lambda));  
146 - break;  
147 - default:  
148 - break;  
149 - }  
150 - }  
151 -  
152 - //return the refractive index associated with the entry  
153 - return newRI;  
154 -  
155 - }  
156 -  
157 - std::string outputEntry(refIndex<T> material)  
158 - {  
159 - //std::string result;  
160 - std::stringstream ss;  
161 -  
162 - //for each field in the entry  
163 - for(int i=0; i<valueList.size(); i++)  
164 - {  
165 - if(i > 0)  
166 - ss<<"\t";  
167 - //store the value in the appropriate location  
168 - switch(valueList[i])  
169 - {  
170 - case field_microns:  
171 - ss<<material.lambda;  
172 - break;  
173 - case field_wavenumber:  
174 - ss<<_2wn(material.lambda);  
175 - break;  
176 - case field_n:  
177 - ss<<material.n.real();  
178 - break;  
179 - case field_k:  
180 - ss<<material.n.imag();  
181 - break;  
182 - case field_A:  
183 - ss<<_2A(material.n.imag(), material.lambda);  
184 - break;  
185 - } 27 + using function< T, complex<T> >::X;
  28 + using function< T, complex<T> >::Y;
  29 + using function< T, complex<T> >::insert;
  30 + using function< T, complex<T> >::bounding;
186 31
187 - }  
188 - return ss.str();  
189 - } 32 + std::string name; //name for the material (defaults to file name)
190 33
  34 + void process_header(std::string str, wave_property& wp, material_property& mp){
191 35
192 - };  
193 -  
194 -  
195 - //a material is a list of refractive index values  
196 - template <class T>  
197 - class material  
198 - {  
199 - std::string name;  
200 - //dispersion (refractive index as a function of wavelength)  
201 - //std::vector< refIndex<T> > dispersion;  
202 - function< T, complex<T> > dispersion;  
203 -  
204 - //average refractive index (approximately 1.4)  
205 - T n0;  
206 -  
207 - /*void add(refIndex<T> ri)  
208 - {  
209 - //refIndex<T> converted = convert(ri, measurement);  
210 - dispersion.push_back(ri);  
211 - }*/  
212 -  
213 - //comparison function for sorting  
214 - static bool compare(refIndex<T> a, refIndex<T> b)  
215 - {  
216 - return (a.lambda < b.lambda);  
217 - }  
218 -  
219 - //comparison function for searching lambda  
220 - /*static bool findCeiling(refIndex<T> a, refIndex<T> b)  
221 - {  
222 - return (a.lambda > b.lambda);  
223 - }*/  
224 -  
225 - public:  
226 - void add(T lambda, complex<T> n)  
227 - {  
228 - dispersion.insert(lambda, n);  
229 - }  
230 -  
231 - std::string getName()  
232 - {  
233 - return name;  
234 - }  
235 - void setName(std::string n)  
236 - {  
237 - name = n;  
238 - }  
239 - T getN0()  
240 - {  
241 - return n0;  
242 - }  
243 - void setN0(T n)  
244 - {  
245 - n0 = n;  
246 - } 36 + std::stringstream ss(str); //create a stream from the data string
  37 + std::string line;
  38 + std::getline(ss, line); //get the first line as a string
  39 + while(line[0] == '#'){ //continue looping while the line is a comment
247 40
248 - void setM(function< T, complex<T> > m)  
249 - {  
250 - dispersion = m;  
251 - }  
252 - unsigned int nSamples()  
253 - {  
254 - return dispersion.size();  
255 - }  
256 - material<T> computeN(T _n0, unsigned int n_samples = 0, T pf = 2)  
257 - {  
258 - /* This function computes the real part of the refractive index  
259 - from the imaginary part. The Hilbert transform is required. I  
260 - use an FFT in order to simplify this, so either the FFTW or CUFFT  
261 - packages are required. CUFFT is used if this file is passed through  
262 - a CUDA compiler. Otherwise, FFTW is used if available.  
263 - */  
264 -  
265 - n0 = _n0;  
266 -  
267 - int N;  
268 - if(n_samples)  
269 - N = n_samples;  
270 - else  
271 - N = dispersion.size();  
272 -  
273 -  
274 -#ifdef FFTW_AVAILABLE  
275 - //allocate memory for the FFT  
276 - complex<T>* Chi2 = (complex<T>*)fftw_malloc(sizeof(complex<T>) * N * pf);  
277 - complex<T>* Chi2FFT = (complex<T>*)fftw_malloc(sizeof(complex<T>) * N * pf);  
278 - complex<T>* Chi1 = (complex<T>*)fftw_malloc(sizeof(complex<T>) * N * pf);  
279 -  
280 - //create an FFT plan for the forward and inverse transforms  
281 - fftw_plan planForward, planInverse;  
282 - planForward = fftw_plan_dft_1d(N*pf, (fftw_complex*)Chi2, (fftw_complex*)Chi2FFT, FFTW_FORWARD, FFTW_ESTIMATE);  
283 - planInverse = fftw_plan_dft_1d(N*pf, (fftw_complex*)Chi2FFT, (fftw_complex*)Chi1, FFTW_BACKWARD, FFTW_ESTIMATE);  
284 -  
285 - float k, alpha;  
286 - T chi_temp;  
287 -  
288 - //the spectrum will be re-sampled in uniform values of wavenumber  
289 - T nuMin = _2wn(dispersion.back().lambda);  
290 - T nuMax = _2wn(dispersion.front().lambda);  
291 - T dnu = (nuMax - nuMin)/(N-1);  
292 - T lambda, tlambda;  
293 - for(int i=0; i<N; i++)  
294 - {  
295 - //go from back-to-front (wavenumber is the inverse of wavelength)  
296 - lambda = _wn(nuMax - i * dnu);  
297 -  
298 - //compute the frequency  
299 - k = 2 * PI / (lambda);  
300 -  
301 - //get the absorbance  
302 - alpha = getN(lambda).imag() * (2 * k);  
303 -  
304 - //compute chi2  
305 - Chi2[i] = -alpha * (n0 / k);  
306 - } 41 + std::stringstream lstream(line); //create a stream from the line
  42 + lstream.ignore(); //ignore the first character ('#')
307 43
308 - //use linear interpolation between the start and end points to pad the spectrum  
309 - //complex<T> nMin = dispersion.back();  
310 - //complex<T> nMax = dispersion.front();  
311 - T a;  
312 - for(int i=N; i<N*pf; i++)  
313 - {  
314 - //a = (T)(i-N)/(T)(N*pf - N);  
315 - //Chi2[i] = a * Chi2[0] + ((T)1 - a) * Chi2[N-1]; 44 + std::string prop; //get the property name
  45 + lstream>>prop;
316 46
317 - Chi2[i] = 0.0;//Chi2[N-1]; 47 + if(prop == "X"){
  48 + std::string wp_name;
  49 + lstream>>wp_name;
  50 + if(wp_name == "microns") wp = microns;
  51 + else if(wp_name == "inverse_cm") wp = inverse_cm;
318 } 52 }
319 -  
320 - //perform the FFT  
321 - fftw_execute(planForward);  
322 -  
323 - //perform the Hilbert transform in the Fourier domain  
324 - complex<T> j(0, 1);  
325 - for(int i=0; i<N*pf; i++)  
326 - {  
327 - //if w = 0, set the DC component to zero  
328 - if(i == 0)  
329 - Chi2FFT[i] *= (T)0.0;  
330 - //if w <0, multiply by i  
331 - else if(i < N*pf/2.0)  
332 - Chi2FFT[i] *= j;  
333 - //if i > N/2, multiply by -i  
334 - else  
335 - Chi2FFT[i] *= -j; 53 + else if(prop == "Y"){
  54 + std::string mp_name;
  55 + lstream>>mp_name;
  56 + if(mp_name == "ri") mp = ri;
  57 + else if(mp_name == "absorbance") mp = absorbance;
336 } 58 }
337 59
338 - //execute the inverse Fourier transform (completing the Hilbert transform)  
339 - fftw_execute(planInverse);  
340 -  
341 - //divide the Chi1 values by N  
342 - for(int i=0; i<N*pf; i++)  
343 - Chi1[i] /= (T)(N*pf);  
344 -  
345 - //create a new material  
346 - material<T> newM;  
347 - newM.dispersion.clear();  
348 - refIndex<T> ri;  
349 - for(int i=0; i<N; i++)  
350 - {  
351 - ri.lambda = _wn(nuMax - i * dnu);  
352 - ri.n.real(Chi1[i].real() / (2 * n0) + n0);  
353 - ri.n.imag(getN(ri.lambda).imag());  
354 -  
355 - newM.dispersion.push_back(ri);  
356 - }  
357 -  
358 -  
359 - //dispersion[i].n.real(Chi1[i].real() / (2 * n0) + n0);  
360 -  
361 -  
362 - /*//output the Absorbance value  
363 - ofstream outOrig("origN.txt");  
364 - for(int i=0; i<N; i++)  
365 - outOrig<<dispersion[i].lambda<<" "<<dispersion[i].n.real()<<endl;  
366 -  
367 - //output the Chi2 value  
368 - ofstream outChi2("chi2.txt");  
369 - for(int i=0; i<N; i++)  
370 - outChi2<<dispersion[i].lambda<<" "<<Chi2[i].real()<<endl;  
371 -  
372 - //output the Fourier transform  
373 - ofstream outFFT("chi2_FFT.txt");  
374 - for(int i=0; i<N; i++)  
375 - {  
376 - float mag = std::sqrt( std::pow(Chi2FFT[i].real(), 2.0) + std::pow(Chi2FFT[i].imag(), 2.0));  
377 - outFFT<<dispersion[i].lambda<<" "<<mag<<endl;  
378 - } 60 + std::getline(ss, line); //get the next line
  61 + }
379 62
380 - //output the computed Chi1 value  
381 - ofstream outChi1("chi1.txt");  
382 - for(int i=0; i<N; i++)  
383 - {  
384 - outChi1<<dispersion[i].lambda<<" "<<Chi1[i].real()<<" "<<Chi1[i].imag()<<endl;  
385 - } 63 + function< T, rts::complex<T> >::process_string(str);
  64 + }
386 65
387 - ofstream outN("n.txt");  
388 - for(int i=0; i<N; i++)  
389 - outN<<dispersion[i].lambda<<" "<<Chi1[i].real() / (2 * n0) + n0<<endl;*/ 66 + void from_inverse_cm(){
  67 + //convert inverse centimeters to wavelength (in microns)
  68 + for(unsigned int i=0; i<X.size(); i++)
  69 + X[i] = 10000 / X[i];
390 70
  71 + //reverse the function array
  72 + std::reverse(X.begin(), X.end());
  73 + std::reverse(Y.begin(), Y.end());
391 74
392 - //de-allocate memory  
393 - fftw_destroy_plan(planForward);  
394 - fftw_destroy_plan(planInverse);  
395 - fftw_free(Chi2);  
396 - fftw_free(Chi2FFT);  
397 - fftw_free(Chi1); 75 + }
398 76
399 - return newM;  
400 -#else  
401 - std::cout<<"MATERIAL Error: Must have FFTW in order to compute Kramers-Kronig."<<std::endl;  
402 - return material<T>();  
403 -#endif  
404 - } 77 + void init(){
  78 + bounding[0] = bounding[1] = rts::complex<T>(1, 0);
  79 + }
405 80
406 - material(T lambda = 1, T n = 1, T k = 0)  
407 - {  
408 - dispersion.insert(lambda, complex<T>(n, k));  
409 - /*//create a default refractive index  
410 - refIndex<T> def;  
411 - def.lambda = lambda;  
412 - def.n.real(n);  
413 - def.n.imag(k);  
414 - add(def);  
415 - */  
416 - //set n0  
417 - n0 = 0;  
418 - }  
419 -  
420 - material(std::string filename, std::string format = "", T scaleA = 1.0)  
421 - {  
422 - fromFile(filename, format);  
423 - n0 = 0;  
424 - }  
425 -  
426 - void fromFile(std::string filename, std::string format = "", T scaleA = 1.0)  
427 - {  
428 - name = filename;  
429 - //clear any previous values  
430 - dispersion = rts::function< T, complex<T> >();  
431 -  
432 - //load the file into a string  
433 - std::ifstream ifs(filename.c_str());  
434 -  
435 - std::string line;  
436 -  
437 - if(!ifs.is_open())  
438 - {  
439 - std::cout<<"Material Error -- file not found: "<<filename<<std::endl;  
440 - exit(1);  
441 - }  
442 -  
443 - //process the file as a string  
444 - std::string instr((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());  
445 - fromStr(instr, format, scaleA);  
446 - }  
447 -  
448 - void fromStr(std::string str, std::string format = "", T scaleA = 1.0)  
449 - {  
450 - //create a string stream to process the input data  
451 - std::stringstream ss(str);  
452 -  
453 - //this string will be read line-by-line (where each line is an entry)  
454 - std::string line;  
455 -  
456 - //if the format is not provided, see if it is in the file, otherwise use a default  
457 - if(format == "")  
458 - {  
459 - //set a default format of "lambda,n,k"  
460 - format = "microns,n,k";  
461 -  
462 - //see if the first line is a comment  
463 - char c = ss.peek();  
464 - if(c == '#')  
465 - {  
466 - //get the first line  
467 - getline(ss, line);  
468 - //look for a bracket, denoting the start of a format string  
469 - int istart = line.find('[');  
470 - if(istart != std::string::npos)  
471 - {  
472 - //look for a bracket terminating the format string  
473 - int iend = line.find(']');  
474 - if(iend != std::string::npos)  
475 - {  
476 - //read the string between the brackets  
477 - format = line.substr(istart+1, iend - istart - 1);  
478 - }  
479 - }  
480 - }  
481 -  
482 - }  
483 -  
484 - entryType<T> entry(format);  
485 -  
486 - //std::cout<<"Loading material with format: "<<format<<std::endl;  
487 -  
488 - while(!ss.eof())  
489 - {  
490 - //read a line from the string  
491 - getline(ss, line);  
492 -  
493 - //if the line is not a comment, process it  
494 - if(line[0] != '#')  
495 - {  
496 - //load the entry and add it to the dispersion list  
497 - add(entry.inputEntry(line, scaleA).lambda, entry.inputEntry(line, scaleA).n);  
498 - }  
499 - //generally have to peek to trigger the eof flag  
500 - ss.peek();  
501 - }  
502 -  
503 - //sort the vector by lambda  
504 - //sort(dispersion.begin(), dispersion.end(), &material<T>::compare);  
505 - }  
506 -  
507 - //convert the material to a string  
508 - std::string toStr(std::string format = "microns,n,k", bool reverse_order = false)  
509 - {  
510 - std::stringstream ss;  
511 - entryType<T> entry(format);  
512 -  
513 - if(reverse_order)  
514 - {  
515 - for(int l=dispersion.size() - 1; l>=0; l--)  
516 - {  
517 - if(l < dispersion.size() - 1) ss<<std::endl;  
518 - ss<<entry.outputEntry(dispersion[l]);  
519 - }  
520 81
521 - }  
522 - else  
523 - {  
524 - for(unsigned int l=0; l<dispersion.size(); l++)  
525 - {  
526 - if(l > 0) ss<<std::endl;  
527 - ss<<entry.outputEntry(dispersion[l]);  
528 - }  
529 - }  
530 -  
531 - return ss.str();  
532 - }  
533 -  
534 - void save(std::string filename, std::string format = "microns,n,k", bool reverse_order = false)  
535 - {  
536 - std::ofstream outfile(filename.c_str());  
537 - outfile<<"#material file saved as [" + format + "]"<<std::endl;  
538 - outfile<<toStr(format, reverse_order)<<std::endl;  
539 -  
540 - }  
541 -  
542 - //convert between wavelength and wavenumber  
543 - /*void nu2lambda(T s = (T)1)  
544 - {  
545 - for(int i=0; i<dispersion.size(); i++)  
546 - dispersion[i].lambda = s/dispersion[i].lambda;  
547 - }  
548 -  
549 - void lambda2nu(T s = (T)1)  
550 - {  
551 - for(int i=0; i<dispersion.size(); i++)  
552 - dispersion[i].lambda = s/dispersion[i].lambda;  
553 - }*/  
554 -  
555 -  
556 - refIndex<T>& operator[](unsigned int i)  
557 - {  
558 - return dispersion[i]; 82 +public:
559 83
560 - } 84 + material(std::string filename, wave_property wp, material_property mp){
  85 + name = filename;
  86 + load(filename, wp, mp);
  87 + }
561 88
562 - complex<T> getN(T l)  
563 - {  
564 - //return complex<T>(1.0, 0.0);  
565 - complex<T> ri = dispersion.linear(l) + n0;  
566 - return ri;  
567 - } 89 + material(std::string filename){
  90 + name = filename;
  91 + load(filename);
  92 + }
568 93
569 - function<T, complex<T> > getF()  
570 - {  
571 - return dispersion + complex<T>(n0, 0.0);  
572 - } 94 + material(){
  95 + init();
  96 + }
573 97
574 - //returns the scattering efficiency at wavelength l  
575 - complex<T> eta(T l)  
576 - {  
577 - //get the refractive index  
578 - complex<T> ri = getN(l); 98 + complex<T> getN(T lambda){
  99 + return function< T, complex<T> >::linear(lambda);
  100 + }
579 101
580 - //convert the refractive index to scattering efficiency  
581 - return ri*ri - 1.0; 102 + void load(std::string filename, wave_property wp, material_property mp){
582 103
583 - }  
584 - //interpolate the given lambda value and return the index of refraction  
585 - complex<T> operator()(T l)  
586 - {  
587 - return getN(l);  
588 - } 104 + //load the file as a function
  105 + function< T, complex<T> >::load(filename);
  106 + }
589 107
  108 + void load(std::string filename){
  109 +
  110 + wave_property wp = inverse_cm;
  111 + material_property mp = ri;
  112 + //turn the file into a string
  113 + std::ifstream t(filename.c_str()); //open the file as a stream
  114 +
  115 + if(!t){
  116 + std::cout<<"ERROR: Couldn't open the material file '"<<filename<<"'"<<std::endl;
  117 + exit(1);
  118 + }
  119 + std::string str((std::istreambuf_iterator<char>(t)),
  120 + std::istreambuf_iterator<char>());
  121 +
  122 + //process the header information
  123 + process_header(str, wp, mp);
  124 +
  125 + //convert units
  126 + if(wp == inverse_cm)
  127 + from_inverse_cm();
  128 + //set the bounding values
  129 + bounding[0] = Y[0];
  130 + bounding[1] = Y.back();
  131 + }
  132 + std::string str(){
  133 + std::stringstream ss;
  134 + ss<<name<<std::endl;
  135 + ss<<function< T, complex<T> >::str();
  136 + return ss.str();
  137 + }
  138 + std::string get_name(){
  139 + return name;
  140 + }
590 141
591 - };  
592 -} //end namespace rts 142 + void set_name(std::string str){
  143 + name = str;
  144 + }
593 145
594 -template <typename T>  
595 -std::ostream& operator<<(std::ostream& os, rts::material<T> m)  
596 -{  
597 - os<<m.toStr(); 146 +};
598 147
599 - return os;  
600 } 148 }
601 149
602 150
603 151
  152 +
604 #endif 153 #endif
optics/mirst-1d.cuh
  Show/Hide comments   View file @0ef519a
1 #include "../optics/material.h" 1 #include "../optics/material.h"
2 #include "../math/complexfield.cuh" 2 #include "../math/complexfield.cuh"
  3 +#include "../math/constants.h"
3 4
4 #include "cufft.h" 5 #include "cufft.h"
5 6
@@ -51,7 +52,7 @@ __global__ void gpu_mirst1d_layer_fft(complex&lt;T&gt;* dest, complex&lt;T&gt;* ri, @@ -51,7 +52,7 @@ __global__ void gpu_mirst1d_layer_fft(complex&lt;T&gt;* dest, complex&lt;T&gt;* ri,
51 52
52 //T l = w * ri[inu].real(); 53 //T l = w * ri[inu].real();
53 //complex<T> e(0.0, -2 * PI * fz * (zf[inu] + zR/2 - l/2.0)); 54 //complex<T> e(0.0, -2 * PI * fz * (zf[inu] + zR/2 - l/2.0));
54 - complex<T> e(0, -2 * PI * fz * (zf[inu] + opl/2)); 55 + complex<T> e(0, -2 * rtsPI * fz * (zf[inu] + opl/2));
55 56
56 complex<T> eta = ri[inu] * ri[inu] - 1; 57 complex<T> eta = ri[inu] * ri[inu] - 1;
57 58
@@ -59,7 +60,7 @@ __global__ void gpu_mirst1d_layer_fft(complex&lt;T&gt;* dest, complex&lt;T&gt;* ri, @@ -59,7 +60,7 @@ __global__ void gpu_mirst1d_layer_fft(complex&lt;T&gt;* dest, complex&lt;T&gt;* ri,
59 if(ifz == 0) 60 if(ifz == 0)
60 dest[i] += opl * exp(e) * eta * src[inu]; 61 dest[i] += opl * exp(e) * eta * src[inu];
61 else 62 else
62 - dest[i] += opl * exp(e) * eta * src[inu] * sin(PI * fz * opl) / (PI * fz * opl); 63 + dest[i] += opl * exp(e) * eta * src[inu] * sin(rtsPI * fz * opl) / (rtsPI * fz * opl);
63 } 64 }
64 65
65 template<typename T> 66 template<typename T>
@@ -139,7 +140,7 @@ __global__ void gpu_mirst1d_apply_filter(complex&lt;T&gt;* sampleFFT, T* lambda, @@ -139,7 +140,7 @@ __global__ void gpu_mirst1d_apply_filter(complex&lt;T&gt;* sampleFFT, T* lambda,
139 //compute the final filter 140 //compute the final filter
140 T mirst = 0; 141 T mirst = 0;
141 if(fz != 0) 142 if(fz != 0)
142 - mirst = 2 * PI * r * s * Q * (1/fz); 143 + mirst = 2 * rtsPI * r * s * Q * (1/fz);
143 144
144 sampleFFT[i] *= mirst; 145 sampleFFT[i] *= mirst;
145 146
@@ -246,7 +247,7 @@ private: @@ -246,7 +247,7 @@ private:
246 //store the refractive index and source profile in a CPU array 247 //store the refractive index and source profile in a CPU array
247 for(int inu=0; inu<S; inu++){ 248 for(int inu=0; inu<S; inu++){
248 //save the refractive index to the GPU 249 //save the refractive index to the GPU
249 - ri = matlist[n].getN(lambdas[inu]); 250 + ri = matlist[n].getN(lambdas[inu]);
250 HANDLE_ERROR(cudaMemcpy( gpuRi + inu, &ri, sizeof(complex<T>), cudaMemcpyHostToDevice )); 251 HANDLE_ERROR(cudaMemcpy( gpuRi + inu, &ri, sizeof(complex<T>), cudaMemcpyHostToDevice ));
251 252
252 //save the source profile to the GPU 253 //save the source profile to the GPU
@@ -320,6 +321,7 @@ public: @@ -320,6 +321,7 @@ public:
320 pad = padding; 321 pad = padding;
321 NA[0] = 0; 322 NA[0] = 0;
322 NA[1] = 0.8; 323 NA[1] = 0.8;
  324 + S = 0;
323 } 325 }
324 326
325 //add a layer, thickness = microns 327 //add a layer, thickness = microns
@@ -397,16 +399,17 @@ public: @@ -397,16 +399,17 @@ public:
397 std::string str(){ 399 std::string str(){
398 400
399 stringstream ss; 401 stringstream ss;
400 - 402 + ss<<"1D MIRST Simulation========================="<<std::endl;
401 ss<<"z-axis resolution: "<<Z<<std::endl; 403 ss<<"z-axis resolution: "<<Z<<std::endl;
  404 + ss<<"simulation domain: ["<<lambdas[0]<<", "<<lambdas.back()<<"]"<<std::endl;
402 ss<<"number of wavelengths: "<<lambdas.size()<<std::endl; 405 ss<<"number of wavelengths: "<<lambdas.size()<<std::endl;
403 ss<<"padding: "<<pad<<std::endl; 406 ss<<"padding: "<<pad<<std::endl;
404 ss<<"sample thickness: "<<thickness()<<" um"<<std::endl; 407 ss<<"sample thickness: "<<thickness()<<" um"<<std::endl;
405 ss<<"dz: "<<dz()<<" um"<<std::endl; 408 ss<<"dz: "<<dz()<<" um"<<std::endl;
406 ss<<std::endl; 409 ss<<std::endl;
407 - ss<<"layers: "<<layers.size()<<" -----------------"<<std::endl; 410 + ss<<layers.size()<<" layers-------------"<<std::endl;
408 for(unsigned int l=0; l<layers.size(); l++) 411 for(unsigned int l=0; l<layers.size(); l++)
409 - ss<<"layer "<<l<<": "<<layers[l]<<" um"<<" ("<<matlist[l].getName()<<")"<<std::endl; 412 + ss<<"layer "<<l<<": "<<layers[l]<<" um"<<"---------"<<std::endl<<matlist[l].str()<<std::endl;
410 413
411 return ss.str(); 414 return ss.str();
412 415
@@ -417,4 +420,4 @@ public: @@ -417,4 +420,4 @@ public:
417 420
418 }; 421 };
419 422
420 -}  
421 \ No newline at end of file 423 \ No newline at end of file
  424 +}
@@ -196,15 +196,15 @@ namespace rts{ @@ -196,15 +196,15 @@ namespace rts{
196 if(d == 0) 196 if(d == 0)
197 ss<<desc[0]; 197 ss<<desc[0];
198 else 198 else
199 - ss<<std::endl<<setfill(' ')<<setw(width)<<" "<<desc[d]; 199 + ss<<std::endl<<std::setfill(' ')<<std::setw(width)<<" "<<desc[d];
200 200
201 } 201 }
202 202
203 //output the range in the right column 203 //output the range in the right column
204 if(range != "" && ansi) 204 if(range != "" && ansi)
205 - ss<<std::endl<<setfill(' ')<<setw(width)<<" "<<" "<<std::string("setw(width)<<" "<<" "<<std::string("\033[1;36m") + range + "\033[0m";33[1;36m") + range + "setw(width)<<" "<<" "<<std::string("\033[1;36m") + range + "\033[0m";33[0m"; 205 + ss<<std::endl<<std::setfill(' ')<<std::setw(width)<<" "<<" "<<std::string("setw(width)<<" "<<" "<<std::string("\033[1;36m") + range + "\033[0m";33[1;36m") + range + "setw(width)<<" "<<" "<<std::string("\033[1;36m") + range + "\033[0m";33[0m";
206 else if(range != "") 206 else if(range != "")
207 - ss<<std::endl<<setfill(' ')<<setw(width)<<" "<<" "<<range; 207 + ss<<std::endl<<std::setfill(' ')<<std::setw(width)<<" "<<" "<<range;
208 208
209 return ss.str(); 209 return ss.str();
210 } 210 }
@@ -253,7 +253,10 @@ namespace rts{ @@ -253,7 +253,10 @@ namespace rts{
253 253
254 public: 254 public:
255 255
256 - arglist(){ col_width = 0; } 256 + arglist(){
  257 + col_width = 0;
  258 + ansi = true;
  259 + }
257 260
258 void set_ansi(bool b) 261 void set_ansi(bool b)
259 { 262 {
@@ -295,9 +298,9 @@ namespace rts{ @@ -295,9 +298,9 @@ namespace rts{
295 if(si != -1 && a == sections[si].index) 298 if(si != -1 && a == sections[si].index)
296 { 299 {
297 if(ansi) 300 if(ansi)
298 - ss<<std::endl<<std::left<<setfill('=')<<setw(col_width)<<std::string("setw(col_width)<<std::string("\033[1;31m") + sections[si].name<<"\033[0m"<<std::endl;33[1;31m") + sections[si].name<<"setw(col_width)<<std::string("\033[1;31m") + sections[si].name<<"\033[0m"<<std::endl;33[0m"<<std::endl; 301 + ss<<std::endl<<std::left<<std::setfill('=')<<std::setw(col_width)<<std::string("setw(col_width)<<std::string("\033[1;31m") + sections[si].name<<"\033[0m"<<std::endl;33[1;31m") + sections[si].name<<"setw(col_width)<<std::string("\033[1;31m") + sections[si].name<<"\033[0m"<<std::endl;33[0m"<<std::endl;
299 else 302 else
300 - ss<<std::endl<<std::left<<setfill('=')<<setw(col_width)<<sections[si].name<<std::endl; 303 + ss<<std::endl<<std::left<<std::setfill('=')<<std::setw(col_width)<<sections[si].name<<std::endl;
301 si++; 304 si++;
302 if(si == sections.size()) si = -1; 305 if(si == sections.size()) si = -1;
303 } 306 }
@@ -326,7 +329,7 @@ namespace rts{ @@ -326,7 +329,7 @@ namespace rts{
326 { 329 {
327 args[i].set(_value); 330 args[i].set(_value);
328 //adjust the column width if necessary 331 //adjust the column width if necessary
329 - col_width = max(col_width, args[i].col_width()); 332 + col_width = std::max(col_width, args[i].col_width());
330 } 333 }
331 else 334 else
332 std::cout<<"ERROR - Argument not recognized: "<<_name<<std::endl; 335 std::cout<<"ERROR - Argument not recognized: "<<_name<<std::endl;