Commit b6179de6d7efffaef298858c1a56f9079f21a89a

Authored by dmayerich
1 parent 967e46a1

added scripts for spectral simulation

@@ -75,10 +75,10 @@ @@ -75,10 +75,10 @@
75 75
76 #define DEFAULT_INTENSITY_FILE "testappend" 76 #define DEFAULT_INTENSITY_FILE "testappend"
77 #define DEFAULT_TRANSMITTANCE_FILE "" 77 #define DEFAULT_TRANSMITTANCE_FILE ""
78 -#define DEFAULT_ABSORBANCE_FILE "" 78 +#define DEFAULT_ABSORBANCE_FILE "out_a"
79 #define DEFAULT_NEAR_FILE "out_n.bmp" 79 #define DEFAULT_NEAR_FILE "out_n.bmp"
80 #define DEFAULT_FAR_FILE "out_f.bmp" 80 #define DEFAULT_FAR_FILE "out_f.bmp"
81 -#define DEFAULT_EXTENDED_SOURCE "" 81 +#define DEFAULT_EXTENDED_SOURCE "einstein_small.jpg"
82 #define DEFAULT_FIELD_TYPE "magnitude" 82 #define DEFAULT_FIELD_TYPE "magnitude"
83 #define DEFAULT_FORMAT fileoutStruct::formatImage 83 #define DEFAULT_FORMAT fileoutStruct::formatImage
84 #define DEFAULT_COLORMAP "brewer" 84 #define DEFAULT_COLORMAP "brewer"
@@ -113,7 +113,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope) @@ -113,7 +113,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope)
113 scalarslice I = scope->getIntensity(); 113 scalarslice I = scope->getIntensity();
114 114
115 if(is_binary(intFile)) 115 if(is_binary(intFile))
116 - I.toEnvi(intFile, scope->nf.lambda, append); 116 + {
  117 + if(wavenumber)
  118 + I.toEnvi(intFile, 10000.0f/scope->nf.lambda, append);
  119 + else
  120 + I.toEnvi(intFile, scope->nf.lambda, append);
  121 + }
117 else 122 else
118 I.toImage(intFile); 123 I.toImage(intFile);
119 } 124 }
@@ -123,7 +128,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope) @@ -123,7 +128,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope)
123 scalarslice I = scope->getAbsorbance(); 128 scalarslice I = scope->getAbsorbance();
124 129
125 if(is_binary(absFile)) 130 if(is_binary(absFile))
126 - I.toEnvi(absFile, scope->nf.lambda, append); 131 + {
  132 + if(wavenumber)
  133 + I.toEnvi(absFile, 10000.0f/scope->nf.lambda, append);
  134 + else
  135 + I.toEnvi(absFile, scope->nf.lambda, append);
  136 + }
127 else 137 else
128 I.toImage(absFile); 138 I.toImage(absFile);
129 } 139 }
@@ -133,7 +143,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope) @@ -133,7 +143,12 @@ void fileoutStruct::saveDetector(microscopeStruct* scope)
133 scalarslice I = scope->getTransmittance(); 143 scalarslice I = scope->getTransmittance();
134 144
135 if(is_binary(transFile)) 145 if(is_binary(transFile))
136 - I.toEnvi(transFile, scope->nf.lambda, append); 146 + {
  147 + if(wavenumber)
  148 + I.toEnvi(transFile, 10000.0f/scope->nf.lambda, append);
  149 + else
  150 + I.toEnvi(transFile, scope->nf.lambda, append);
  151 + }
137 else 152 else
138 I.toImage(transFile); 153 I.toImage(transFile);
139 } 154 }
@@ -28,13 +28,15 @@ struct fileoutStruct{ @@ -28,13 +28,15 @@ struct fileoutStruct{
28 28
29 field_type field; 29 field_type field;
30 30
  31 + //flag for output in wavenumber
  32 + bool wavenumber;
  33 +
31 //image_source source; 34 //image_source source;
32 35
33 //color map info 36 //color map info
34 rts::colormap::colormapType colormap; 37 rts::colormap::colormapType colormap;
35 ptype colorMax; 38 ptype colorMax;
36 39
37 -  
38 void Save(microscopeStruct* scope); 40 void Save(microscopeStruct* scope);
39 void Simulate(microscopeStruct* scope); 41 void Simulate(microscopeStruct* scope);
40 42
@@ -42,9 +44,10 @@ struct fileoutStruct{ @@ -42,9 +44,10 @@ struct fileoutStruct{
42 bool is_binary(std::string filename); 44 bool is_binary(std::string filename);
43 void saveNearField(nearfieldStruct* nf); 45 void saveNearField(nearfieldStruct* nf);
44 void saveFarField(microscopeStruct* scope); 46 void saveFarField(microscopeStruct* scope);
45 - void saveDetector(microscopeStruct* scope);  
46 - 47 + void saveDetector(microscopeStruct* scope);
47 48
48 }; 49 };
49 50
  51 +
  52 +
50 #endif 53 #endif
@@ -138,14 +138,14 @@ static void loadMaterials(po::variables_map vm) @@ -138,14 +138,14 @@ static void loadMaterials(po::variables_map vm)
138 138
139 for(int i=0; i<matVec.size(); i+=2) 139 for(int i=0; i<matVec.size(); i+=2)
140 { 140 {
141 - rts::material<ptype> newM(vm["lambda"].as<ptype>(), matVec[i], matVec[i+1]); 141 + rts::material<ptype> newM(SCOPE->nf.lambda, matVec[i], matVec[i+1]);
142 SCOPE->nf.mVector.push_back(newM); 142 SCOPE->nf.mVector.push_back(newM);
143 } 143 }
144 } 144 }
145 else 145 else
146 { 146 {
147 //add the command line material as the default (material 0) 147 //add the command line material as the default (material 0)
148 - rts::material<ptype> newM(vm["lambda"].as<ptype>(), vm["n"].as<ptype>(), vm["k"].as<ptype>()); 148 + rts::material<ptype> newM(SCOPE->nf.lambda, vm["n"].as<ptype>(), vm["k"].as<ptype>());
149 SCOPE->nf.mVector.push_back(newM); 149 SCOPE->nf.mVector.push_back(newM);
150 } 150 }
151 151
@@ -178,7 +178,7 @@ static void loadNearfieldParams(po::variables_map vm) @@ -178,7 +178,7 @@ static void loadNearfieldParams(po::variables_map vm)
178 SCOPE->nf.planeWave = planeWave; 178 SCOPE->nf.planeWave = planeWave;
179 179
180 //get the wavelength 180 //get the wavelength
181 - SCOPE->nf.lambda = vm["lambda"].as<ptype>(); 181 + //SCOPE->nf.lambda = vm["lambda"].as<ptype>();
182 182
183 //get the incident field amplitude 183 //get the incident field amplitude
184 SCOPE->nf.A = vm["amplitude"].as<ptype>(); 184 SCOPE->nf.A = vm["amplitude"].as<ptype>();
@@ -328,6 +328,7 @@ static void SetOptions(po::options_description &amp;desc) @@ -328,6 +328,7 @@ static void SetOptions(po::options_description &amp;desc)
328 ("material-file,M", po::value< vector<string> >()->multitoken(), "material file:\n [lambda n k]") 328 ("material-file,M", po::value< vector<string> >()->multitoken(), "material file:\n [lambda n k]")
329 ("materials", po::value< vector<ptype> >()->multitoken(), "materials specified using n, k pairs:\n ex. --materials n1 k1 n2 k2\n (if used --n and --k are ignored)") 329 ("materials", po::value< vector<ptype> >()->multitoken(), "materials specified using n, k pairs:\n ex. --materials n1 k1 n2 k2\n (if used --n and --k are ignored)")
330 ("lambda,l", po::value<ptype>()->default_value(DEFAULT_LAMBDA), "incident wavelength") 330 ("lambda,l", po::value<ptype>()->default_value(DEFAULT_LAMBDA), "incident wavelength")
  331 + ("nu", po::value<ptype>(), "incident frequency (in cm^-1)\n(if specified, lambda is ignored)")
331 ("theta,t", po::value<ptype>()->default_value(DEFAULT_K_THETA), "light direction (polar coords)") 332 ("theta,t", po::value<ptype>()->default_value(DEFAULT_K_THETA), "light direction (polar coords)")
332 ("phi,p", po::value<ptype>()->default_value(DEFAULT_K_PHI)) 333 ("phi,p", po::value<ptype>()->default_value(DEFAULT_K_PHI))
333 ("fx", po::value<ptype>()->default_value(DEFAULT_FOCUS_X), "incident focal point") 334 ("fx", po::value<ptype>()->default_value(DEFAULT_FOCUS_X), "incident focal point")
@@ -374,6 +375,20 @@ static void LoadParameters(int argc, char *argv[]) @@ -374,6 +375,20 @@ static void LoadParameters(int argc, char *argv[])
374 exit(1); 375 exit(1);
375 } 376 }
376 377
  378 + //load the wavelength
  379 + if(vm.count("nu"))
  380 + {
  381 + //wavelength is given in wavenumber - transform and flag
  382 + SCOPE->nf.lambda = 10000/vm["nu"].as<ptype>();
  383 + gFileOut.wavenumber = true;
  384 + }
  385 + //otherwise we are using lambda = wavelength
  386 + else
  387 + {
  388 + SCOPE->nf.lambda = vm["lambda"].as<ptype>();
  389 + gFileOut.wavenumber = false;
  390 + }
  391 +
377 //load spheres 392 //load spheres
378 loadSpheres(vm); 393 loadSpheres(vm);
379 394
1 -#!/usr/bin/python3  
2 -import subprocess  
3 -  
4 -command = "bimsim"  
5 -  
6 -#images  
7 -intImage = "out_i.bmp"  
8 -absImage = "out_a"  
9 -#detector specs  
10 -dsize = 208  
11 -dsample = 6.5  
12 -res = int(dsize / dsample)  
13 -sampling = 4  
14 -#incident field  
15 -#nu = 800  
16 -source = "central.bmp"  
17 -order = 100  
18 -#sphere  
19 -x = 0  
20 -y = 0  
21 -z = 0  
22 -a = 5  
23 -#spectral samples  
24 -nuStart = 800  
25 -nuStop = 4000  
26 -nuStep = 10  
27 -iters = int((nuStop - nuStart) / nuStep)  
28 -  
29 -  
30 -#set the position of the image plane  
31 -command += " -u " + str(-dsize/2)  
32 -command += " -v " + str(-dsize/2)  
33 -command += " -w " + str(a)  
34 -command += " -U " + str(dsize/2)  
35 -command += " -V " + str(dsize/2)  
36 -command += " -W " + str(a)  
37 -command += " --plane-norm-x " + str(0)  
38 -command += " --plane-norm-y " + str(0)  
39 -command += " --plane-norm-z " + str(1)  
40 -command += " -R " + str(res)  
41 -command += " --supersample " + str(sampling)  
42 -command += " -X " + source  
43 -command += " --field-order " + str(order)  
44 -command += " -I " + intImage  
45 -command += " -A " + absImage  
46 -command += " --append"  
47 -command += " -x " + str(x)  
48 -command += " -y " + str(y)  
49 -command += " -z " + str(z)  
50 -command += " -r " + str(a)  
51 -  
52 -for inu in range(0, iters):  
53 - print("Iteration # " + str(inu + 1) + "/" + str(iters))  
54 - nu = nuStart + inu * nuStep  
55 - lam = 10000.0/nu  
56 -  
57 - runcommand = command + " -l " + str(lam)  
58 - print(runcommand)  
59 - subprocess.call(runcommand, shell=True)  
60 -  
61 -#print("Hello world!")  
62 \ No newline at end of file 1 \ No newline at end of file
  2 +#!/usr/bin/python3
  3 +
  4 +
  5 +import subprocess, sys
  6 +
  7 +nuStart = 800
  8 +nuStop = 4000
  9 +nuStep = 10
  10 +
  11 +if len(sys.argv) > 1:
  12 + nuStart = sys.argv[1]
  13 +if len(sys.argv) > 2:
  14 + nuStop = sys.argv[2]
  15 +if len(sys.argv) > 3:
  16 + nuStep = sys.argv[3]
  17 +
  18 +command = "bimsim"
  19 +
  20 +#images
  21 +intImage = "out_i.bmp"
  22 +absImage = "out_a"
  23 +#detector specs
  24 +dsize = 208
  25 +dsample = 1
  26 +res = int(dsize / dsample)
  27 +sampling = 1
  28 +padding = 1
  29 +#incident field
  30 +source = ""
  31 +order = 100
  32 +mc = 400
  33 +#sphere
  34 +x = 0
  35 +y = 0
  36 +z = 0
  37 +a = 5
  38 +n = 1.4
  39 +k = 0.0
  40 +#spectral samples
  41 +iters = int((nuStop - nuStart) / nuStep) + 1
  42 +#optics
  43 +NAin = 0.2
  44 +NAout = 0.5
  45 +
  46 +
  47 +#set the position of the image plane
  48 +command += " -u " + str(-dsize/2)
  49 +command += " -v " + str(-dsize/2)
  50 +command += " -w " + str(a)
  51 +command += " -U " + str(dsize/2)
  52 +command += " -V " + str(dsize/2)
  53 +command += " -W " + str(a)
  54 +command += " --plane-norm-x " + str(0)
  55 +command += " --plane-norm-y " + str(0)
  56 +command += " --plane-norm-z " + str(1)
  57 +command += " -R " + str(res)
  58 +command += " --supersample " + str(sampling)
  59 +
  60 +if source != "":
  61 + command += " -X " + source
  62 +command += " --field-order " + str(order)
  63 +command += " -d " + str(padding)
  64 +command += " -s " + str(mc)
  65 +command += " -I " + intImage
  66 +command += " -A " + absImage
  67 +command += " --append"
  68 +#sphere
  69 +command += " -x " + str(x)
  70 +command += " -y " + str(y)
  71 +command += " -z " + str(z)
  72 +command += " -r " + str(a)
  73 +command += " -n " + str(n)
  74 +command += " -k " + str(k)
  75 +#set the optics
  76 +command += " -c " + str(NAin)
  77 +command += " -C " + str(NAout)
  78 +command += " -o " + str(NAin)
  79 +command += " -O " + str(NAout)
  80 +
  81 +for inu in range(0, iters):
  82 + print("Iteration # " + str(inu + 1) + "/" + str(iters))
  83 + nu = nuStart + inu * nuStep
  84 +
  85 + runcommand = command + " --nu " + str(nu)
  86 + print(runcommand)
  87 + subprocess.call(runcommand, shell=True)
  88 +
  89 +#print("Hello world!")