1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | |
---|
3 | // |
---|
4 | // Monte Carlo simulator for SASCALC |
---|
5 | // October 2008 SRK |
---|
6 | // |
---|
7 | // This code simulates the scattering for a selected model based on the instrument configuration |
---|
8 | // This code is based directly on John Barker's code, which includes multiple scattering effects. |
---|
9 | // A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten() |
---|
10 | // |
---|
11 | // |
---|
12 | // - Most importantly, this needs to be checked for correctness of the MC simulation |
---|
13 | // X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation |
---|
14 | // - why does my integrated tau not match up with John's analytical calculations? where are the assumptions? |
---|
15 | // - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles |
---|
16 | // - what is magical about Qu? Is this an assumpution? |
---|
17 | |
---|
18 | // - at the larger angles, is the "flat" detector being properly accounted for - in terms of |
---|
19 | // the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector |
---|
20 | // so that what I see is already "corrected"? |
---|
21 | // X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation |
---|
22 | // by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x) |
---|
23 | // but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example. |
---|
24 | // - the background parameter for the model MUST be zero, or the integration for scattering |
---|
25 | // power will be incorrect. |
---|
26 | // - fully use the SASCALC input, most importantly, flux on sample. |
---|
27 | // X if no MC desired, still use the selected model |
---|
28 | // X better display of MC results on panel |
---|
29 | // - settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) |
---|
30 | // - add quartz window scattering to the simulation somehow |
---|
31 | // - do smeared models make any sense?? |
---|
32 | // - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition |
---|
33 | // or the volume fraction of solvent. |
---|
34 | // |
---|
35 | // - add to the results the fraction of coherently scattered neutrons that are singly scattered, different than |
---|
36 | // the overall fraction of singly scattered, and maybe more important to know. |
---|
37 | // |
---|
38 | // - change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons |
---|
39 | // aren't counted. Is the # that interact a better number? |
---|
40 | // |
---|
41 | // - do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the |
---|
42 | // effects on the absolute scale can be seen? |
---|
43 | // |
---|
44 | // - why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back? |
---|
45 | // - can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering |
---|
46 | // - ask John how to verify what is going on |
---|
47 | // - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up. |
---|
48 | // a warning has been added - but the models are better fixed with the limiting value. |
---|
49 | // |
---|
50 | // |
---|
51 | |
---|
52 | |
---|
53 | |
---|
54 | ////////// |
---|
55 | // PROGRAM Monte_SANS |
---|
56 | // PROGRAM simulates multiple SANS. |
---|
57 | // revised 2/12/99 JGB |
---|
58 | // added calculation of random deviate, and 2D 10/2008 SRK |
---|
59 | |
---|
60 | // N1 = NUMBER OF INCIDENT NEUTRONS. |
---|
61 | // N2 = NUMBER INTERACTED IN THE SAMPLE. |
---|
62 | // N3 = NUMBER ABSORBED. |
---|
63 | // THETA = SCATTERING ANGLE. |
---|
64 | |
---|
65 | // IMON = 'Enter number of neutrons to use in simulation.' |
---|
66 | // NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).' |
---|
67 | // R1 = 'Enter beam radius. (cm)' |
---|
68 | // R2 = 'Enter sample radius. (cm)' |
---|
69 | // thick = 'Enter sample thickness. (cm)' |
---|
70 | // wavelength = 'Enter neutron wavelength. (A)' |
---|
71 | // R0 = 'Enter sphere radius. (A)' |
---|
72 | // |
---|
73 | |
---|
74 | Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results) |
---|
75 | WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results |
---|
76 | |
---|
77 | Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas |
---|
78 | Variable NUM_BINS,N_INDEX |
---|
79 | Variable RHO,SIGSAS,SIGABS_0 |
---|
80 | Variable ii,jj,IND,idum,INDEX,IR,NQ |
---|
81 | Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow |
---|
82 | Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG |
---|
83 | Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI |
---|
84 | //in John's implementation, he dimensioned the indices of the arrays to begin |
---|
85 | // at 0, making things much easier for me... |
---|
86 | //DIMENSION NT(0:5000),J1(0:5000),J2(0:5000),NN(0:100) |
---|
87 | |
---|
88 | Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat |
---|
89 | Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single |
---|
90 | Variable DONE,FIND_THETA,err //used as logicals |
---|
91 | |
---|
92 | Variable Vx,Vy,Vz,Theta_z,qq |
---|
93 | Variable Sig_scat,Sig_abs,Ratio,Sig_total |
---|
94 | Variable isOn=0,testQ,testPhi,xPixel,yPixel |
---|
95 | |
---|
96 | imon = inputWave[0] |
---|
97 | r1 = inputWave[1] |
---|
98 | r2 = inputWave[2] |
---|
99 | xCtr = inputWave[3] |
---|
100 | yCtr = inputWave[4] |
---|
101 | sdd = inputWave[5] |
---|
102 | pixSize = inputWave[6] |
---|
103 | thick = inputWave[7] |
---|
104 | wavelength = inputWave[8] |
---|
105 | sig_incoh = inputWave[9] |
---|
106 | sig_sas = inputWave[10] |
---|
107 | |
---|
108 | // SetRandomSeed 0.1 //to get a reproduceable sequence |
---|
109 | |
---|
110 | //scattering power and maximum qvalue to bin |
---|
111 | // zpow = .1 //scattering power, calculated below |
---|
112 | qmax = 4*pi/wavelength //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value) |
---|
113 | sigabs_0 = 0.0 // ignore absorption cross section/wavelength [1/(cm A)] |
---|
114 | N_INDEX = 50 // maximum number of scattering events per neutron |
---|
115 | num_bins = 200 //number of 1-D bins (not really used) |
---|
116 | |
---|
117 | // my additions - calculate the randome deviate function as needed |
---|
118 | // and calculate the scattering power from the model function |
---|
119 | // |
---|
120 | Variable left = leftx(ran_dev) |
---|
121 | Variable delta = deltax(ran_dev) |
---|
122 | |
---|
123 | //c total SAS cross-section |
---|
124 | // |
---|
125 | // input a test value for the incoherent scattering from water |
---|
126 | // |
---|
127 | // sig_incoh = 5.6 //[1/cm] as calculated for H2O, now a parameter |
---|
128 | // |
---|
129 | // SIG_SAS = zpow/thick |
---|
130 | zpow = sig_sas*thick //since I now calculate the sig_sas from the model |
---|
131 | SIG_ABS = SIGABS_0 * WAVElength |
---|
132 | SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh |
---|
133 | // Print "The TOTAL XSECTION. (CM-1) is ",sig_total |
---|
134 | // Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas |
---|
135 | results[0] = sig_total |
---|
136 | results[1] = sig_sas |
---|
137 | // RATIO = SIG_ABS / SIG_TOTAL |
---|
138 | RATIO = sig_incoh / SIG_TOTAL |
---|
139 | //!!!! the ratio is not yet properly weighted for the volume fractions of each component!!! |
---|
140 | |
---|
141 | //c assuming theta = sin(theta)...OK |
---|
142 | theta_max = wavelength*qmax/(2*pi) |
---|
143 | //C SET Theta-STEP SIZE. |
---|
144 | DTH = Theta_max/NUM_BINS |
---|
145 | // Print "theta bin size = dth = ",dth |
---|
146 | |
---|
147 | //C INITIALIZE COUNTERS. |
---|
148 | N1 = 0 |
---|
149 | N2 = 0 |
---|
150 | N3 = 0 |
---|
151 | |
---|
152 | //C INITIALIZE ARRAYS. |
---|
153 | j1 = 0 |
---|
154 | j2 = 0 |
---|
155 | nt = 0 |
---|
156 | nn=0 |
---|
157 | |
---|
158 | //C MONITOR LOOP - looping over the number of incedent neutrons |
---|
159 | //note that zz, is the z-position in the sample - NOT the scattering power |
---|
160 | |
---|
161 | // NOW, start the loop, throwing neutrons at the sample. |
---|
162 | do |
---|
163 | Vx = 0.0 // Initialize direction vector. |
---|
164 | Vy = 0.0 |
---|
165 | Vz = 1.0 |
---|
166 | |
---|
167 | Theta = 0.0 // Initialize scattering angle. |
---|
168 | Phi = 0.0 // Intialize azimuthal angle. |
---|
169 | N1 += 1 // Increment total number neutrons counter. |
---|
170 | DONE = 0 // True when neutron is absorbed or when scattered out of the sample. |
---|
171 | INDEX = 0 // Set counter for number of scattering events. |
---|
172 | zz = 0.0 // Set entering dimension of sample. |
---|
173 | |
---|
174 | do // Makes sure position is within circle. |
---|
175 | ran = abs(enoise(1)) //[0,1] |
---|
176 | xx = 2.0*R1*(Ran-0.5) //X beam position of neutron entering sample. |
---|
177 | ran = abs(enoise(1)) //[0,1] |
---|
178 | yy = 2.0*R1*(Ran-0.5) //Y beam position ... |
---|
179 | RR = SQRT(xx*xx+yy*yy) //Radial position of neutron in incident beam. |
---|
180 | while(rr>r1) |
---|
181 | |
---|
182 | do //Scattering Loop, will exit when "done" == 1 |
---|
183 | // keep scattering multiple times until the neutron exits the sample |
---|
184 | ran = abs(enoise(1)) //[0,1] RANDOM NUMBER FOR DETERMINING PATH LENGTH |
---|
185 | ll = PATH_len(ran,Sig_total) |
---|
186 | //Determine new scattering direction vector. |
---|
187 | err = NewDirection(vx,vy,vz,Theta,Phi) //vx,vy,vz is updated, theta, phi unchanged by function |
---|
188 | |
---|
189 | //X,Y,Z-POSITION OF SCATTERING EVENT. |
---|
190 | xx += ll*vx |
---|
191 | yy += ll*vy |
---|
192 | zz += ll*vz |
---|
193 | RR = sqrt(xx*xx+yy*yy) //radial position of scattering event. |
---|
194 | |
---|
195 | //Check whether interaction occurred within sample volume. |
---|
196 | IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) |
---|
197 | //NEUTRON INTERACTED. |
---|
198 | INDEX += 1 //Increment counter of scattering events. |
---|
199 | IF(INDEX == 1) |
---|
200 | N2 += 1 //Increment # of scat. neutrons |
---|
201 | Endif |
---|
202 | ran = abs(enoise(1)) //[0,1] |
---|
203 | //Split neutron interactions into scattering and absorption events |
---|
204 | IF(ran > ratio ) //C NEUTRON SCATTERED coherently |
---|
205 | FIND_THETA = 0 //false |
---|
206 | DO |
---|
207 | //ran = abs(enoise(1)) //[0,1] |
---|
208 | //theta = Scat_angle(Ran,R_DAB,wavelength) // CHOOSE DAB ANGLE -- this is 2Theta |
---|
209 | //Q0 = 2*PI*THETA/WAVElength // John chose theta, calculated Q |
---|
210 | |
---|
211 | // pick a q-value from the deviate function |
---|
212 | // pnt2x truncates the point to an integer before returning the x |
---|
213 | // so get it from the wave scaling instead |
---|
214 | Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta |
---|
215 | theta = Q0/2/Pi*wavelength //SAS approximation |
---|
216 | |
---|
217 | //Print "q0, theta = ",q0,theta |
---|
218 | |
---|
219 | FIND_THETA = 1 //always accept |
---|
220 | |
---|
221 | while(!find_theta) |
---|
222 | ran = abs(enoise(1)) //[0,1] |
---|
223 | PHI = 2.0*PI*Ran //Chooses azimuthal scattering angle. |
---|
224 | ELSE |
---|
225 | //NEUTRON scattered incoherently |
---|
226 | // N3 += 1 |
---|
227 | // DONE = 1 |
---|
228 | // phi and theta are random over the entire sphere of scattering |
---|
229 | // !can't just choose random theta and phi, won't be random over sphere solid angle |
---|
230 | |
---|
231 | ran = abs(enoise(1)) //[0,1] |
---|
232 | theta = acos(2*ran-1) |
---|
233 | |
---|
234 | ran = abs(enoise(1)) //[0,1] |
---|
235 | PHI = 2.0*PI*Ran //Chooses azimuthal scattering angle. |
---|
236 | ENDIF //(ran > ratio) |
---|
237 | ELSE |
---|
238 | //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere |
---|
239 | DONE = 1 //done = true, will exit from loop |
---|
240 | //Increment #scattering events array |
---|
241 | If (index <= N_Index) |
---|
242 | NN[INDEX] += 1 |
---|
243 | Endif |
---|
244 | //IF (VZ > 1.0) // FIX INVALID ARGUMENT |
---|
245 | //VZ = 1.0 - 1.2e-7 |
---|
246 | //ENDIF |
---|
247 | Theta_z = acos(Vz) // Angle WITH respect to z axis. |
---|
248 | testQ = 2*pi*sin(theta_z)/wavelength |
---|
249 | |
---|
250 | // pick a random phi angle, and see if it lands on the detector |
---|
251 | // since the scattering is isotropic, I can safely pick a new, random value |
---|
252 | // this would not be true if simulating anisotropic scattering. |
---|
253 | testPhi = abs(enoise(1))*2*Pi |
---|
254 | // is it on the detector? |
---|
255 | FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) |
---|
256 | |
---|
257 | // if(xPixel != xCtr && yPixel != yCtr) |
---|
258 | // Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel |
---|
259 | // endif |
---|
260 | |
---|
261 | if(xPixel != -1 && yPixel != -1) |
---|
262 | isOn += 1 |
---|
263 | //if(index==1) // only the single scattering events |
---|
264 | MC_linear_data[xPixel][yPixel] += 1 //this is the total scattering, including multiple scattering |
---|
265 | //endif |
---|
266 | endif |
---|
267 | |
---|
268 | If(theta_z < theta_max) |
---|
269 | //Choose index for scattering angle array. |
---|
270 | //IND = NINT(THETA_z/DTH + 0.4999999) |
---|
271 | ind = round(THETA_z/DTH + 0.4999999) //round is eqivalent to nint() |
---|
272 | NT[ind] += 1 //Increment bin for angle. |
---|
273 | //Increment angle array for single scattering events. |
---|
274 | IF(INDEX == 1) |
---|
275 | j1[ind] += 1 |
---|
276 | Endif |
---|
277 | //Increment angle array for double scattering events. |
---|
278 | IF (INDEX == 2) |
---|
279 | j2[ind] += 1 |
---|
280 | Endif |
---|
281 | EndIf |
---|
282 | |
---|
283 | ENDIF |
---|
284 | while (!done) |
---|
285 | while(n1 < imon) |
---|
286 | |
---|
287 | // Print "Monte Carlo Done" |
---|
288 | trans_th = exp(-sig_total*thick) |
---|
289 | // TRANS_exp = (N1-N2) / N1 // Transmission |
---|
290 | // dsigma/domega assuming isotropic scattering, with no absorption. |
---|
291 | // Print "trans_exp = ",trans_exp |
---|
292 | // Print "total # of neutrons reaching 2D detector",isOn |
---|
293 | // Print "fraction of incident neutrons reaching detector ",isOn/iMon |
---|
294 | results[2] = isOn |
---|
295 | results[3] = isOn/iMon |
---|
296 | |
---|
297 | |
---|
298 | // OUTPUT of the 1D data, not necessary now since I want the 2D |
---|
299 | // Make/O/N=(num_bins) qvals,int_ms,sig_ms,int_sing,int_doub |
---|
300 | // ii=1 |
---|
301 | // Print "binning" |
---|
302 | // do |
---|
303 | // //CALCULATE MEAN THETA IN BIN. |
---|
304 | // THETA_z = (ii-0.5)*DTH // Mean scattering angle of bin. |
---|
305 | // //Solid angle of Ith bin. |
---|
306 | // D_OMEGA = 2*PI*ABS( COS(ii*DTH) - COS((ii-1)*DTH) ) |
---|
307 | // //SOLID ANGLE CORRECTION: YIELDING CROSS-SECTION. |
---|
308 | // dS_dW = NT[ii]/(N1*THICK*Trans_th*D_OMEGA) |
---|
309 | // SIG = NT[ii]^0.5/(N1*THICK*Trans_th*D_OMEGA) |
---|
310 | // ds_dw_single = J1[ii]/(N1*THICK*Trans_th*D_OMEGA) |
---|
311 | // ds_dw_double = J2[ii]/(N1*THICK*Trans_th*D_OMEGA) |
---|
312 | // //Deviation from isotropic model. |
---|
313 | // qq = 4*pi*sin(0.5*theta_z)/wavelength |
---|
314 | // qvals[ii-1] = qq |
---|
315 | // int_ms[ii-1] = dS_dW |
---|
316 | // sig_ms[ii-1] = sig |
---|
317 | // int_sing[ii-1] = ds_dw_single |
---|
318 | // int_doub[ii-1] = ds_dw_double |
---|
319 | // //140 WRITE (7,145) qq,dS_dW,SIG,ds_dw_single,ds_dw_double |
---|
320 | // ii+=1 |
---|
321 | // while(ii<=num_bins) |
---|
322 | // Print "done binning" |
---|
323 | |
---|
324 | |
---|
325 | // Write(7,*) |
---|
326 | // Write(7,*) '#Times Sc. #Events ' |
---|
327 | // DO 150 I = 1,N_INDEX |
---|
328 | //150 WRITE (7,146) I,NN(I) |
---|
329 | //146 Format (I5,T10,F8.0) |
---|
330 | |
---|
331 | /// Write(7,171) N1 |
---|
332 | // Write(7,172) N2 |
---|
333 | // Write(7,173) N3 |
---|
334 | // Write(7,174) N2-N3 |
---|
335 | |
---|
336 | //171 Format('Total number of neutrons: N1= ',E10.5) |
---|
337 | ///172 Format('Number of neutrons that interact: N2= ',E10.5) |
---|
338 | //173 Format('Number of absorption events: N3= ',E10.5) |
---|
339 | //174 Format('# of neutrons that scatter out:(N2-N3)= ',E10.5) |
---|
340 | |
---|
341 | // Print "Total number of neutrons = ",N1 |
---|
342 | // Print "Total number of neutrons that interact = ",N2 |
---|
343 | // Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2 |
---|
344 | results[4] = N2 |
---|
345 | results[5] = sum(j1,-inf,inf)/N2 |
---|
346 | |
---|
347 | Tabs = (N1-N3)/N1 |
---|
348 | Ttot = (N1-N2)/N1 |
---|
349 | I1_sumI = NN[0]/(N2-N3) |
---|
350 | // Print "Tabs = ",Tabs |
---|
351 | // Print "Transmitted neutrons = ",Ttot |
---|
352 | results[6] = Ttot |
---|
353 | // Print "I1 / all I1 = ", I1_sumI |
---|
354 | // Print "DONE!" |
---|
355 | |
---|
356 | End |
---|
357 | //////// end of main function for calculating multiple scattering |
---|
358 | |
---|
359 | |
---|
360 | // returns the random deviate as a wave |
---|
361 | // and the total SAS cross-section [1/cm] sig_sas |
---|
362 | Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs) |
---|
363 | FUNCREF SANSModelAAO_MCproto func |
---|
364 | WAVE coef |
---|
365 | Variable lam |
---|
366 | String outWave |
---|
367 | Variable &SASxs |
---|
368 | |
---|
369 | Variable nPts_ran=10000,qu |
---|
370 | qu = 4*pi/lam |
---|
371 | |
---|
372 | Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw // if these waves are 1000 pts, the results are "pixelated" |
---|
373 | WAVE Gq = root:Packages:NIST:SAS:gQ |
---|
374 | WAVE xw = root:Packages:NIST:SAS:xw |
---|
375 | SetScale/I x (0+1e-6),qu*(1-1e-10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors |
---|
376 | xw=x //for the AAO |
---|
377 | func(coef,Gq,xw) //call as AAO |
---|
378 | |
---|
379 | // Gq = x*Gq // SAS approximation |
---|
380 | Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu)) // exact |
---|
381 | // |
---|
382 | Integrate/METH=1 Gq/D=Gq_INT |
---|
383 | |
---|
384 | SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1] |
---|
385 | |
---|
386 | Gq_INT /= Gq_INT[nPts_ran-1] |
---|
387 | |
---|
388 | Duplicate/O Gq_INT $outWave |
---|
389 | |
---|
390 | return(0) |
---|
391 | End |
---|
392 | |
---|
393 | Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) |
---|
394 | Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel |
---|
395 | |
---|
396 | Variable theta,dy,dx,qx,qy |
---|
397 | //decompose to qx,qy |
---|
398 | qx = testQ*cos(testPhi) |
---|
399 | qy = testQ*sin(testPhi) |
---|
400 | |
---|
401 | //convert qx,qy to pixel locations relative to # of pixels x, y from center |
---|
402 | theta = 2*asin(qy*lam/4/pi) |
---|
403 | dy = sdd*tan(theta) |
---|
404 | yPixel = round(yCtr + dy/pixSize) |
---|
405 | |
---|
406 | theta = 2*asin(qx*lam/4/pi) |
---|
407 | dx = sdd*tan(theta) |
---|
408 | xPixel = round(xCtr + dx/pixSize) |
---|
409 | |
---|
410 | //if on detector, return xPix and yPix values, otherwise -1 |
---|
411 | if(yPixel > 127 || yPixel < 0) |
---|
412 | yPixel = -1 |
---|
413 | endif |
---|
414 | if(xPixel > 127 || xPixel < 0) |
---|
415 | xPixel = -1 |
---|
416 | endif |
---|
417 | |
---|
418 | return(0) |
---|
419 | End |
---|
420 | |
---|
421 | Function MC_CheckFunctionAndCoef(funcStr,coefStr) |
---|
422 | String funcStr,coefStr |
---|
423 | |
---|
424 | SVAR/Z listStr=root:Packages:NIST:coefKWStr |
---|
425 | if(SVAR_Exists(listStr) == 1) |
---|
426 | String properCoefStr = StringByKey(funcStr, listStr ,"=",";",0) |
---|
427 | if(cmpstr("",properCoefStr)==0) |
---|
428 | return(0) //false, no match found, so properCoefStr is returned null |
---|
429 | endif |
---|
430 | if(cmpstr(coefStr,properCoefStr)==0) |
---|
431 | return(1) //true, the coef is the correct match |
---|
432 | endif |
---|
433 | endif |
---|
434 | return(0) //false, wrong coef |
---|
435 | End |
---|
436 | |
---|
437 | Function/S MC_getFunctionCoef(funcStr) |
---|
438 | String funcStr |
---|
439 | |
---|
440 | SVAR/Z listStr=root:Packages:NIST:coefKWStr |
---|
441 | String coefStr="" |
---|
442 | if(SVAR_Exists(listStr) == 1) |
---|
443 | coefStr = StringByKey(funcStr, listStr ,"=",";",0) |
---|
444 | endif |
---|
445 | return(coefStr) |
---|
446 | End |
---|
447 | |
---|
448 | Function SANSModelAAO_MCproto(w,yw,xw) |
---|
449 | Wave w,yw,xw |
---|
450 | |
---|
451 | Print "in SANSModelAAO_MCproto function" |
---|
452 | return(1) |
---|
453 | end |
---|
454 | |
---|
455 | Function/S MC_FunctionPopupList() |
---|
456 | String list,tmp |
---|
457 | list = FunctionList("*",";","KIND:10") //get every user defined curve fit function |
---|
458 | |
---|
459 | //now start to remove everything the user doesn't need to see... |
---|
460 | |
---|
461 | tmp = FunctionList("*_proto",";","KIND:10") //prototypes |
---|
462 | list = RemoveFromList(tmp, list ,";") |
---|
463 | //prototypes that show up if GF is loaded |
---|
464 | list = RemoveFromList("GFFitFuncTemplate", list) |
---|
465 | list = RemoveFromList("GFFitAllAtOnceTemplate", list) |
---|
466 | list = RemoveFromList("NewGlblFitFunc", list) |
---|
467 | list = RemoveFromList("NewGlblFitFuncAllAtOnce", list) |
---|
468 | list = RemoveFromList("GlobalFitFunc", list) |
---|
469 | list = RemoveFromList("GlobalFitAllAtOnce", list) |
---|
470 | list = RemoveFromList("GFFitAAOStructTemplate", list) |
---|
471 | list = RemoveFromList("NewGF_SetXWaveInList", list) |
---|
472 | list = RemoveFromList("NewGlblFitFuncAAOStruct", list) |
---|
473 | |
---|
474 | // more to remove as a result of 2D/Gizmo |
---|
475 | list = RemoveFromList("A_WMRunLessThanDelta", list) |
---|
476 | list = RemoveFromList("WMFindNaNValue", list) |
---|
477 | list = RemoveFromList("WM_Make3DBarChartParametricWave", list) |
---|
478 | list = RemoveFromList("UpdateQxQy2Mat", list) |
---|
479 | list = RemoveFromList("MakeBSMask", list) |
---|
480 | |
---|
481 | // MOTOFIT/GenFit bits |
---|
482 | tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;" |
---|
483 | list = RemoveFromList(tmp, list ,";") |
---|
484 | |
---|
485 | // SANS Reduction bits |
---|
486 | tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;" |
---|
487 | list = RemoveFromList(tmp, list ,";") |
---|
488 | list = RemoveFromList("Monte_SANS", list) |
---|
489 | |
---|
490 | tmp = FunctionList("f*",";","NPARAMS:2") //point calculations |
---|
491 | list = RemoveFromList(tmp, list ,";") |
---|
492 | |
---|
493 | tmp = FunctionList("fSmear*",";","NPARAMS:3") //smeared dependency calculations |
---|
494 | list = RemoveFromList(tmp, list ,";") |
---|
495 | |
---|
496 | // tmp = FunctionList("*X",";","KIND:4") //XOPs, but these shouldn't show up if KIND:10 is used initially |
---|
497 | // Print "X* = ",tmp |
---|
498 | // print " " |
---|
499 | // list = RemoveFromList(tmp, list ,";") |
---|
500 | |
---|
501 | //non-fit functions that I can't seem to filter out |
---|
502 | list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";") |
---|
503 | |
---|
504 | if(strlen(list)==0) |
---|
505 | list = "No functions plotted" |
---|
506 | endif |
---|
507 | |
---|
508 | list = SortList(list) |
---|
509 | return(list) |
---|
510 | End |
---|
511 | |
---|
512 | |
---|
513 | //Function Scat_Angle(Ran,R_DAB,wavelength) |
---|
514 | // Variable Ran,r_dab,wavelength |
---|
515 | // |
---|
516 | // Variable qq,arg,theta |
---|
517 | // qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) ) |
---|
518 | // arg = qq*wavelength/(4*pi) |
---|
519 | // If (arg < 1.0) |
---|
520 | // theta = 2.*asin(arg) |
---|
521 | // else |
---|
522 | // theta = pi |
---|
523 | // endif |
---|
524 | // Return (theta) |
---|
525 | //End |
---|
526 | |
---|
527 | //calculates new direction (xyz) from an old direction |
---|
528 | //theta and phi don't change |
---|
529 | Function NewDirection(vx,vy,vz,theta,phi) |
---|
530 | Variable &vx,&vy,&vz |
---|
531 | Variable theta,phi |
---|
532 | |
---|
533 | Variable err=0,vx0,vy0,vz0 |
---|
534 | Variable nx,ny,mag_xy,tx,ty,tz |
---|
535 | |
---|
536 | //store old direction vector |
---|
537 | vx0 = vx |
---|
538 | vy0 = vy |
---|
539 | vz0 = vz |
---|
540 | |
---|
541 | mag_xy = sqrt(vx0*vx0 + vy0*vy0) |
---|
542 | if(mag_xy < 1e-12) |
---|
543 | //old vector lies along beam direction |
---|
544 | nx = 0 |
---|
545 | ny = 1 |
---|
546 | tx = 1 |
---|
547 | ty = 0 |
---|
548 | tz = 0 |
---|
549 | else |
---|
550 | Nx = -Vy0 / Mag_XY |
---|
551 | Ny = Vx0 / Mag_XY |
---|
552 | Tx = -Vz0*Vx0 / Mag_XY |
---|
553 | Ty = -Vz0*Vy0 / Mag_XY |
---|
554 | Tz = Mag_XY |
---|
555 | endif |
---|
556 | |
---|
557 | //new scattered direction vector |
---|
558 | Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0 |
---|
559 | Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0 |
---|
560 | Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0 |
---|
561 | |
---|
562 | Return(err) |
---|
563 | End |
---|
564 | |
---|
565 | Function path_len(aval,sig_tot) |
---|
566 | Variable aval,sig_tot |
---|
567 | |
---|
568 | Variable retval |
---|
569 | |
---|
570 | retval = -1*ln(1-aval)/sig_tot |
---|
571 | |
---|
572 | return(retval) |
---|
573 | End |
---|
574 | |
---|
575 | Window MC_SASCALC() : Panel |
---|
576 | PauseUpdate; Silent 1 // building window... |
---|
577 | NewPanel /W=(787,44,1088,563) /N=MC_SASCALC as "SANS Simulator" |
---|
578 | CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation" |
---|
579 | CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo |
---|
580 | SetVariable MC_setvar0,pos={11,38},size={144,15},bodyWidth=80,title="# of neutrons" |
---|
581 | SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon |
---|
582 | SetVariable MC_setvar0_1,pos={11,121},size={131,15},bodyWidth=60,title="Thickness (cm)" |
---|
583 | SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick |
---|
584 | SetVariable MC_setvar0_2,pos={11,93},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" |
---|
585 | SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh |
---|
586 | SetVariable MC_setvar0_3,pos={11,149},size={150,15},bodyWidth=60,title="Sample Radius (cm)" |
---|
587 | SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2 |
---|
588 | PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function" |
---|
589 | PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" |
---|
590 | Button MC_button0,pos={17,455},size={50,20},proc=MC_DoItButtonProc,title="Do It" |
---|
591 | Button MC_button1,pos={15,484},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" |
---|
592 | |
---|
593 | SetDataFolder root:Packages:NIST:SAS: |
---|
594 | Edit/W=(13,174,284,435)/HOST=# results_desc,results |
---|
595 | ModifyTable width(Point)=0,width(results_desc)=150 |
---|
596 | SetDataFolder root: |
---|
597 | RenameWindow #,T_results |
---|
598 | SetActiveSubwindow ## |
---|
599 | EndMacro |
---|
600 | |
---|
601 | Function MC_ModelPopMenuProc(pa) : PopupMenuControl |
---|
602 | STRUCT WMPopupAction &pa |
---|
603 | |
---|
604 | switch( pa.eventCode ) |
---|
605 | case 2: // mouse up |
---|
606 | Variable popNum = pa.popNum |
---|
607 | String popStr = pa.popStr |
---|
608 | SVAR gStr = root:Packages:NIST:SAS:gFuncStr |
---|
609 | gStr = popStr |
---|
610 | |
---|
611 | break |
---|
612 | endswitch |
---|
613 | |
---|
614 | return 0 |
---|
615 | End |
---|
616 | |
---|
617 | Function MC_DoItButtonProc(ba) : ButtonControl |
---|
618 | STRUCT WMButtonAction &ba |
---|
619 | |
---|
620 | switch( ba.eventCode ) |
---|
621 | case 2: // mouse up |
---|
622 | // click code here |
---|
623 | ReCalculateInten(1) |
---|
624 | break |
---|
625 | endswitch |
---|
626 | |
---|
627 | return 0 |
---|
628 | End |
---|
629 | |
---|
630 | |
---|
631 | Function MC_Display2DButtonProc(ba) : ButtonControl |
---|
632 | STRUCT WMButtonAction &ba |
---|
633 | |
---|
634 | switch( ba.eventCode ) |
---|
635 | case 2: // mouse up |
---|
636 | // click code here |
---|
637 | Execute "ChangeDisplay(\"SAS\")" |
---|
638 | break |
---|
639 | endswitch |
---|
640 | |
---|
641 | return 0 |
---|
642 | End |
---|
643 | |
---|
644 | /////UNUSED, testing routines that have note been updated to work with SASCALC |
---|
645 | // |
---|
646 | //Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr) |
---|
647 | // Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1 |
---|
648 | // String funcStr |
---|
649 | // Prompt funcStr, "Pick the model function", popup, MC_FunctionPopupList() |
---|
650 | // |
---|
651 | // String coefStr = MC_getFunctionCoef(funcStr) |
---|
652 | // Variable pixSize = 0.5 // can't have 11 parameters in macro! |
---|
653 | // |
---|
654 | // if(!MC_CheckFunctionAndCoef(funcStr,coefStr)) |
---|
655 | // Abort "The coefficients and function type do not match. Please correct the selections in the popup menus." |
---|
656 | // endif |
---|
657 | // |
---|
658 | // Make/O/D/N=10 root:Packages:NIST:SAS:results |
---|
659 | // Make/O/T/N=10 root:Packages:NIST:SAS:results_desc |
---|
660 | // |
---|
661 | // RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results) |
---|
662 | // |
---|
663 | //End |
---|
664 | // |
---|
665 | //Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr) |
---|
666 | // Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh |
---|
667 | // String funcStr,coefStr |
---|
668 | // WAVE results |
---|
669 | // |
---|
670 | // FUNCREF SANSModelAAO_MCproto func=$funcStr |
---|
671 | // |
---|
672 | // Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results) |
---|
673 | //End |
---|
674 | // |
---|
675 | ////// END UNUSED BLOCK |
---|