source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf @ 757

Last change on this file since 757 was 757, checked in by srkline, 12 years ago

Fixed bug re-introduced by using the global ASCII code for Angstrom. Now it will use a default "A" if the global does not exist, rather than throwin up an error.

For 2D MC simulation, now if the data was scaled to ABS, the 2D data is "un-scaled" back to raw counts before saving as VAX format. A few counts may be lost here and there during the conversion, but less than 0.01% in my limited testing.

File size: 73.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4
5//
6// Monte Carlo simulator for SASCALC
7// October 2008 SRK
8//
9// This code simulates the scattering for a selected model based on the instrument configuration
10// This code is based directly on John Barker's code, which includes multiple scattering effects.
11// A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten()
12//
13//
14// at the end of the procedure fils is a *very* simple example of scripting for unattended simulations
15// - not for the casual user at all.
16
17
18
19// the RNG issue is really not worth the effort. multiple copies with different RNG is as good as I need. Plus,
20// whatever XOP crashing was happining during threading is really unlikely to be from the RNG
21//
22// -- June 2010 - calls from different threads to the same RNG really seems to cause a crash. Probably as soon
23//                              as the different threads try to call at the same time. Found this out by accident doing the
24//                              wavelength spread. Each thread called ran3 at that point, and the crash came quickly. Went
25//                              away immediately when I kept the ran calls consistent and isolated within threads.
26//
27// *** look into erand48() as the (pseudo) random number generator (it's a standard c-lib function, at least on unix)
28//     and is apparantly thread safe. drand48() returns values [0.0,1.0)
29//http://qnxcs.unomaha.edu/help/product/neutrino/lib_ref/e/erand48.html
30//
31
32
33// - Why am I off by a factor of 2.7 - 3.7 (MC too high) relative to real data?
34//   I need to include efficiency (70%?) - do I knock these off before the simulation or do I
35//    really simulate that some fraction of neutrons on the detector don't actually get counted?
36//   Is the flux estimate up-to-date? !! Flux estimates at NG3 are out-of-date....
37// - my simulated transmission is larger than what is measured, even after correcting for the quartz cell.
38//   Why? Do I need to include absorption? Just inherent problems with incoherent cross sections?
39
40// - Most importantly, this needs to be checked for correctness of the MC simulation
41// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation
42// X why does my integrated tau not match up with John's analytical calculations? where are the assumptions?
43// X get rid of all small angle assumptions - to make sure that the calculation is correct at all angles
44
45//
46// X at the larger angles, is the "flat" detector being properly accounted for - in terms of
47//   the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector
48//   so that what I see is already "corrected"?
49// X the MC will, of course benefit greatly from being XOPized. Maybe think about parallel implementation
50//   by allowing the data arrays to accumulate. First pass at the XOP is done. Not pretty, not the speediest (5.8x)
51//   but it is functional. Add spinCursor() for long calculations. See WaveAccess XOP example.
52// X the background parameter for the model MUST be zero, or the integration for scattering
53//    power will be incorrect. (now the LAST point in a copy of the coef wave is set to zero, only for the rad_dev calculation
54// X fully use the SASCALC input, most importantly, flux on sample.
55// X if no MC desired, still use the selected model
56// X better display of MC results on panel
57// X settings for "count for X seconds" or "how long to 1E6 cts on detector" (but 1E6 is typically too many counts...)
58// X warn of projected simulation time
59// - add quartz window scattering to the simulation somehow
60// -?- do smeared models make any sense?? Yes, John agrees that they do, and may be used in a more realistic simulation
61//   -?- but the random deviate can't be properly calculated...
62// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition
63//   or the volume fraction of solvent.
64//
65// X add to the results the fraction of coherently scattered neutrons that are singly scattered, different than
66//   the overall fraction of singly scattered, and maybe more important to know.
67//
68// X change the fraction reaching the detector to exclude those that don't interact. These transmitted neutrons
69//   aren't counted. Is the # that interact a better number?
70//
71// - do we want to NOT offset the data by a multiplicative factor as it is "frozen" , so that the
72//   effects on the absolute scale can be seen?
73//
74// X why is "pure" incoherent scattering giving me a q^-1 slope, even with the detector all the way back?
75// -NO- can I speed up by assuming everything interacts? This would compromise the ability to calculate multiple scattering
76// X ask John how to verify what is going on
77// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up.
78//   a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable
79//   unless something else can be done. Using a log-spacing of points doesn't seem to help, and it introduces a lot of
80//   other problems. Not the way to go.
81// - if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it
82//   should always default to...
83//
84
85// --- TO ADD ---
86// X- wavelength distribution = another RNG to select the wavelength
87// ---- done Jun 2010, approximating the wavelength distribution as a Gaussian, based on the triangular
88//                      FWHM. Wavelength distribution added to XOP too, and now very accurately matches the shape of the 1D
89//                      simulation.
90//
91//
92// X- quartz windows (an empirical model?? or measure some real data - power Law + background)
93// X- blocked beam (measure this too, and have some empirical model for this too - Broad Peak)
94// --- Done (mostly). quartz cell and blocked beam have been added empirically, giving the count rate and predicted
95//     scattering. Count time for the simulated scattering is the same as the sample. The simulated EC
96//               data can be plotted, but only by hand right now. EC and blocked beam are combined.
97//
98// -- divergence / size of the incoming beam. Currently everything is parallel, and anything that is transmitted
99//              simply ends up in (xCtr,yCtr), and the "real" profile of the beam is not captured.
100
101
102
103
104
105// setting the flag to 1 == 2D simulation
106// any other value for flag == 1D simulation
107//
108// must remember to close/reopen the simulation control panel to get the correct panel
109//
110Function Set_2DMonteCarlo_Flag(value)
111        Variable value
112       
113        NVAR flag=root:Packages:NIST:SAS:gDoMonteCarlo
114        flag=value
115        return(0)
116end
117
118// threaded call to the main function, adds up the individual runs, and returns what is to be displayed
119// results is calculated and sent back for display
120Function Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
121        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
122
123        //initialize ran1 in the XOP by passing a negative integer
124        // does nothing in the Igor code
125        Duplicate/O results retWave
126
127        Variable NNeutron=inputWave[0]
128        Variable i,nthreads= ThreadProcessorCount
129
130// make sure that the XOP exists if we are going to thread     
131#if exists("Monte_SANSX4")
132        //OK
133        if(nthreads>4)          //only support 4 processors until I can figure out how to properly thread the XOP and to loop it
134                                                        //AND - just use 4 threads rather than the 8 (4 + 4 hyperthread?) my quad-core reports.
135                nthreads=4
136        endif
137#else
138        nthreads = 1
139#endif
140       
141//      nthreads = 1
142        NVAR mt=root:myGlobals:gThreadGroupID
143        mt = ThreadGroupCreate(nthreads)
144        NVAR gInitTime = root:Packages:NIST:SAS:gRanDateTime            //time that SASCALC was started
145//      Print "thread group ID = ",mt
146       
147        inputWave[0] = NNeutron/nthreads                //split up the number of neutrons
148       
149        for(i=0;i<nthreads;i+=1)
150                Duplicate/O nt $("nt"+num2istr(i))              //new instance for each thread
151                Duplicate/O j1 $("j1"+num2istr(i))
152                Duplicate/O j2 $("j2"+num2istr(i))
153                Duplicate/O nn $("nn"+num2istr(i))
154                Duplicate/O linear_data $("linear_data"+num2istr(i))
155                Duplicate/O retWave $("retWave"+num2istr(i))
156                Duplicate/O inputWave $("inputWave"+num2istr(i))
157                Duplicate/O ran_dev $("ran_dev"+num2istr(i))
158                               
159                // ?? I need explicit wave references?
160                // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach.
161                // more likely there is something bad going on in the XOP code.
162                if(i==0)
163                        WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0
164                        retWave0 = 0            //clear the return wave
165                        retWave0[0] = -1*trunc(datetime-gInitTime)              //to initialize ran3
166                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0)
167                        Print "started thread 1"
168                endif
169                if(i==1)
170                        WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1
171                        retWave1 = 0                    //clear the return wave
172                        retWave1[0] = -1*trunc(datetime-gInitTime-2)            //to initialize ran1
173                        ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1)
174                        Print "started thread 2"
175                endif
176                if(i==2)
177                        WAVE inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2
178                        retWave2[0] = -1*trunc(datetime-gInitTime-3)            //to initialize ran3a
179                        ThreadStart mt,i,Monte_SANS_W3(inputWave2,ran_dev2,nt2,j12,j22,nn2,linear_data2,retWave2)
180                        Print "started thread 3"
181                endif
182                if(i==3)
183                        WAVE inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3
184                        retWave3[0] = -1*trunc(datetime-gInitTime-4)            //to initialize ran1a
185                        ThreadStart mt,i,Monte_SANS_W4(inputWave3,ran_dev3,nt3,j13,j23,nn3,linear_data3,retWave3)
186                        Print "started thread 4"
187                endif
188        endfor
189
190// wait until done
191        do
192                variable tgs= ThreadGroupWait(mt,100)
193        while( tgs != 0 )
194        variable dummy= ThreadGroupRelease(mt)
195        mt=0
196        Print "done with all threads"
197
198        // calculate all of the bits for the results
199        if(nthreads == 1)
200                nt = nt0                // add up each instance
201                j1 = j10
202                j2 = j20
203                nn = nn0
204                linear_data = linear_data0
205                retWave = retWave0
206        endif
207        if(nthreads == 2)
208                nt = nt0+nt1            // add up each instance
209                j1 = j10+j11
210                j2 = j20+j21
211                nn = nn0+nn1
212                linear_data = linear_data0+linear_data1
213                retWave = retWave0+retWave1
214        endif
215        if(nthreads == 3)
216                nt = nt0+nt1+nt2                // add up each instance
217                j1 = j10+j11+j12
218                j2 = j20+j21+j22
219                nn = nn0+nn1+nn2
220                linear_data = linear_data0+linear_data1+linear_data2
221                retWave = retWave0+retWave1+retWave2
222        endif
223        if(nthreads == 4)
224                nt = nt0+nt1+nt2+nt3            // add up each instance
225                j1 = j10+j11+j12+j13
226                j2 = j20+j21+j22+j23
227                nn = nn0+nn1+nn2+nn3
228                linear_data = linear_data0+linear_data1+linear_data2+linear_data3
229                retWave = retWave0+retWave1+retWave2+retWave3
230        endif
231       
232        // fill up the results wave
233        Variable xc,yc
234        xc=inputWave[3]
235        yc=inputWave[4]
236        results[0] = inputWave[9]+inputWave[10]         //total XS
237        results[1] = inputWave[10]                                              //SAS XS
238        results[2] = retWave[1]                                                 //number that interact n2
239        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
240        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
241        results[5] = retWave[4]/retWave[7]                                              //single coherent fraction
242        results[6] = retWave[5]/retWave[7]                              //multiple coherent fraction
243        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
244        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
245       
246        return(0)
247End
248
249// worker function for threads, does nothing except switch between XOP and Igor versions
250//
251// uses ran3
252ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
253        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
254       
255#if exists("Monte_SANSX")
256        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
257#else
258        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
259#endif
260
261        return (0)
262End
263
264// worker function for threads, does nothing except switch between XOP and Igor versions
265//
266// uses ran1
267ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
268        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
269       
270#if exists("Monte_SANSX2")
271        Monte_SANSX2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
272#else
273//      Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
274#endif
275
276        return (0)
277End
278
279// uses ran3a
280ThreadSafe Function Monte_SANS_W3(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
281        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
282       
283#if exists("Monte_SANSX3")
284        Monte_SANSX3(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
285#else
286//      Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
287#endif
288
289        return (0)
290End
291
292// uses ran1a
293ThreadSafe Function Monte_SANS_W4(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
294        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
295       
296#if exists("Monte_SANSX4")
297        Monte_SANSX4(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
298#else
299//      Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
300#endif
301
302        return (0)
303End
304
305
306
307// NON-threaded call to the main function returns what is to be displayed
308// results is calculated and sent back for display
309Function Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
310        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results
311
312        //initialize ran1 in the XOP by passing a negative integer
313        // does nothing in the Igor code, enoise is already initialized
314        Duplicate/O results retWave
315        WAVE retWave
316        retWave[0] = -1*abs(trunc(100000*enoise(1)))
317       
318#if exists("Monte_SANSX")
319        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
320#else
321        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,retWave)
322#endif
323
324        // fill up the results wave
325        Variable xc,yc
326        xc=inputWave[3]
327        yc=inputWave[4]
328        results[0] = inputWave[9]+inputWave[10]         //total XS
329        results[1] = inputWave[10]                                              //SAS XS
330        results[2] = retWave[1]                                                 //number that interact n2
331        results[3] = retWave[2] - linear_data[xc][yc]                           //# reaching detector minus Q(0)
332        results[4] = retWave[3]/retWave[1]                              //avg# times scattered
333        results[5] = retWave[4]/retWave[7]                                              //single coherent fraction
334        results[6] = retWave[5]/retWave[7]                              //double coherent fraction
335        results[7] = retWave[6]/retWave[1]                              //multiple scatter fraction
336        results[8] = (retWave[0]-retWave[1])/retWave[0]                 //transmitted fraction
337       
338        return(0)
339End
340
341
342
343//////////
344//    PROGRAM Monte_SANS
345//    PROGRAM simulates multiple SANS.
346//       revised 2/12/99  JGB
347//            added calculation of random deviate, and 2D 10/2008 SRK
348
349//    N1 = NUMBER OF INCIDENT NEUTRONS.
350//    N2 = NUMBER INTERACTED IN THE SAMPLE.
351//    N3 = NUMBER ABSORBED.
352//    THETA = SCATTERING ANGLE.
353
354//        IMON = 'Enter number of neutrons to use in simulation.'
355//        NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).'
356//        R1 = 'Enter beam radius. (cm)'
357//        R2 = 'Enter sample radius. (cm)'
358//        thick = 'Enter sample thickness. (cm)'
359//        wavelength = 'Enter neutron wavelength. (A)'
360//        R0 = 'Enter sphere radius. (A)'
361//
362
363ThreadSafe Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)
364        WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results
365
366        Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
367        Variable NUM_BINS,N_INDEX
368        Variable RHO,SIGSAS,SIGABS_0
369        Variable ii,jj,IND,idum,INDEX,IR,NQ
370        Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow
371        Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG
372        Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI
373        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat
374        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single
375        Variable DONE,FIND_THETA,err            //used as logicals
376
377        Variable Vx,Vy,Vz,Theta_z,qq
378        Variable Sig_scat,Sig_abs,Ratio,Sig_total
379        Variable isOn=0,testQ,testPhi,xPixel,yPixel
380        Variable NSingleIncoherent,NSingleCoherent,NScatterEvents,incoherentEvent,coherentEvent
381        Variable NDoubleCoherent,NMultipleScatter,countIt,detEfficiency
382        Variable NMultipleCoherent,NCoherentEvents
383        Variable deltaLam,v1,v2,currWavelength,rsq,fac          //for simulating wavelength distribution
384       
385        // don't set to other than one here. Detector efficiency is handled outside, only passing the number of
386        // countable neutrons to any of the simulation functions (n=imon*eff)
387        detEfficiency = 1.0             //70% counting efficiency = 0.7
388       
389        imon = inputWave[0]
390        r1 = inputWave[1]
391        r2 = inputWave[2]
392        xCtr = inputWave[3]
393        yCtr = inputWave[4]
394        sdd = inputWave[5]
395        pixSize = inputWave[6]
396        thick = inputWave[7]
397        wavelength = inputWave[8]
398        sig_incoh = inputWave[9]
399        sig_sas = inputWave[10]
400        deltaLam = inputWave[11]
401       
402//      SetRandomSeed 0.1               //to get a reproduceable sequence
403
404//scattering power and maximum qvalue to bin
405//      zpow = .1               //scattering power, calculated below
406        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
407        sigabs_0 = 0.0          // ignore absorption cross section/wavelength [1/(cm A)]
408        N_INDEX = 50            // maximum number of scattering events per neutron
409        num_bins = 200          //number of 1-D bins (not really used)
410       
411// my additions - calculate the random deviate function as needed
412// and calculate the scattering power from the model function (passed in as a wave)
413//
414        Variable left = leftx(ran_dev)
415        Variable delta = deltax(ran_dev)
416       
417//c       total SAS cross-section
418//      SIG_SAS = zpow/thick
419        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model
420        SIG_ABS = SIGABS_0 * WAVElength
421        sig_abs = 0.0           //cm-1
422        SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh
423//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
424//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
425//      results[0] = sig_total          //assign these after everything's done
426//      results[1] = sig_sas
427//      variable ratio1,ratio2
428//      ratio1 = sig_abs/sig_total
429//      ratio2 = sig_incoh/sig_total
430//      // 0->ratio1 = abs
431//      // ratio1 -> ratio2 = incoh
432//      // > ratio2 = coherent
433        RATIO = sig_incoh / SIG_TOTAL
434       
435//c       assuming theta = sin(theta)...OK
436        theta_max = wavelength*qmax/(2*pi)
437//C     SET Theta-STEP SIZE.
438        DTH = Theta_max/NUM_BINS
439//      Print "theta bin size = dth = ",dth
440
441//C     INITIALIZE COUNTERS.
442        N1 = 0
443        N2 = 0
444        N3 = 0
445        NSingleIncoherent = 0
446        NSingleCoherent = 0
447        NDoubleCoherent = 0
448        NMultipleScatter = 0
449        NScatterEvents = 0
450        NMultipleCoherent = 0
451        NCoherentEvents = 0
452
453//C     INITIALIZE ARRAYS.
454        j1 = 0
455        j2 = 0
456        nt = 0
457        nn=0
458       
459//C     MONITOR LOOP - looping over the number of incedent neutrons
460//note that zz, is the z-position in the sample - NOT the scattering power
461
462// NOW, start the loop, throwing neutrons at the sample.
463        do
464                Vx = 0.0                        // Initialize direction vector.
465                Vy = 0.0
466                Vz = 1.0
467               
468                Theta = 0.0             //      Initialize scattering angle.
469                Phi = 0.0                       //      Intialize azimuthal angle.
470                N1 += 1                 //      Increment total number neutrons counter.
471                DONE = 0                        //      True when neutron is scattered out of the sample.
472                INDEX = 0                       //      Set counter for number of scattering events.
473                zz = 0.0                        //      Set entering dimension of sample.
474                incoherentEvent = 0
475                coherentEvent = 0
476               
477                do                                      //      Makes sure position is within circle.
478                        ran = abs(enoise(1))            //[0,1]
479                        xx = 2.0*R1*(Ran-0.5)           //X beam position of neutron entering sample.
480                        ran = abs(enoise(1))            //[0,1]
481                        yy = 2.0*R1*(Ran-0.5)           //Y beam position ...
482                        RR = SQRT(xx*xx+yy*yy)          //Radial position of neutron in incident beam.
483                while(rr>r1)
484
485                //pick the wavelength out of the wavelength spread, approximate as a gaussian
486                // from NR - pg 288. Needs random # from [0,1]. del is deltaLam/lam (as FWHM) and the
487                // 2.35 converts to a gaussian std dev.
488                do
489                        v1=2.0*abs(enoise(1))-1.0
490                        v2=2.0*abs(enoise(1))-1.0
491                        rsq=v1*v1+v2*v2
492                while (rsq >= 1.0 || rsq == 0.0)
493                fac=sqrt(-2.0*log(rsq)/rsq)
494               
495//              gset=v1*fac             //technically, I'm throwing away one of the two values
496               
497                currWavelength = (v2*fac)*deltaLam*wavelength/2.35 + wavelength
498               
499               
500                do    //Scattering Loop, will exit when "done" == 1
501                                // keep scattering multiple times until the neutron exits the sample
502                        ran = abs(enoise(1))            //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
503                        ll = PATH_len(ran,Sig_total)
504                        //Determine new scattering direction vector.
505                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function
506
507                        //X,Y,Z-POSITION OF SCATTERING EVENT.
508                        xx += ll*vx
509                        yy += ll*vy
510                        zz += ll*vz
511                        RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event.
512
513                        //Check whether interaction occurred within sample volume.
514                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
515                                //NEUTRON INTERACTED.
516                                ran = abs(enoise(1))            //[0,1]
517                               
518//                              if(ran<ratio1)
519//                                      //absorption event
520//                                      n3 +=1
521//                                      done=1
522//                              else
523
524                                INDEX += 1                      //Increment counter of scattering events.
525                                IF(INDEX == 1)
526                                        N2 += 1                 //Increment # of scat. neutrons
527                                Endif
528                                //Split neutron interactions into scattering and absorption events
529//                              IF(ran > (ratio1+ratio2) )              //C             NEUTRON SCATTERED coherently
530                                IF(ran > ratio)         //C             NEUTRON SCATTERED coherently
531                                        coherentEvent += 1
532                                        FIND_THETA = 0                  //false
533                                        DO
534                                                //ran = abs(enoise(1))          //[0,1]
535                                                //theta = Scat_angle(Ran,R_DAB,wavelength)      // CHOOSE DAB ANGLE -- this is 2Theta
536                                                //Q0 = 2*PI*THETA/WAVElength                                    // John chose theta, calculated Q
537
538                                                // pick a q-value from the deviate function
539                                                // pnt2x truncates the point to an integer before returning the x
540                                                // so get it from the wave scaling instead
541                                                Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta
542                                                theta = Q0/2/Pi*currWavelength          //SAS approximation. 1% error at theta=30 deg (theta/2=15deg)
543                                               
544                                                //Print "q0, theta = ",q0,theta
545                                               
546                                                FIND_THETA = 1          //always accept
547
548                                        while(!find_theta)
549                                        ran = abs(enoise(1))            //[0,1]
550                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
551                                ELSE
552                                        //NEUTRON scattered incoherently
553                   // N3 += 1
554                  // DONE = 1
555                  // phi and theta are random over the entire sphere of scattering
556                  // !can't just choose random theta and phi, won't be random over sphere solid angle
557                        incoherentEvent += 1
558                       
559                        ran = abs(enoise(1))            //[0,1]
560                                        theta = acos(2*ran-1)           
561                       
562                        ran = abs(enoise(1))            //[0,1]
563                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
564                                ENDIF           //(ran > ratio)
565//                              endif           // event was absorption
566                        ELSE
567                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
568                                DONE = 1                //done = true, will exit from loop
569                               
570//                              countIt = 1
571//                              if(abs(enoise(1)) > detEfficiency)              //efficiency of 70% wired @top
572//                                      countIt = 0                                     //detector does not register
573//                              endif
574                               
575                                //Increment #scattering events array
576                                If (index <= N_Index)
577                                        NN[INDEX] += 1
578                                Endif
579                               
580                                if(index != 0)          //the neutron interacted at least once, figure out where it ends up
581
582                                        Theta_z = acos(Vz)              // Angle WITH respect to z axis.
583                                        testQ = 2*pi*sin(theta_z)/currWavelength
584                                       
585                                        // pick a random phi angle, and see if it lands on the detector
586                                        // since the scattering is isotropic, I can safely pick a new, random value
587                                        // this would not be true if simulating anisotropic scattering.
588                                        //testPhi = abs(enoise(1))*2*Pi
589                                        testPhi = MC_FindPhi(Vx,Vy)             //use the exiting phi value as defined by Vx and Vy
590                                       
591                                        // is it on the detector?       
592                                        FindPixel(testQ,testPhi,currWavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
593                                       
594                                        if(xPixel != -1 && yPixel != -1)
595                                                //if(index==1)  // only the single scattering events
596                                                        MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering
597                                                //endif
598                                                        isOn += 1               // neutron that lands on detector
599                                        endif
600       
601                                        If(theta_z < theta_max)
602                                                //Choose index for scattering angle array.
603                                                //IND = NINT(THETA_z/DTH + 0.4999999)
604                                                ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint()
605                                                NT[ind] += 1                    //Increment bin for angle.
606                                                //Increment angle array for single scattering events.
607                                                IF(INDEX == 1)
608                                                        j1[ind] += 1
609                                                Endif
610                                                //Increment angle array for double scattering events.
611                                                IF (INDEX == 2)
612                                                        j2[ind] += 1
613                                                Endif
614                                        EndIf
615                                       
616                                        // increment all of the counters now since done==1 here and I'm sure to exit and get another neutron
617                                        NScatterEvents += index         //total number of scattering events
618                                        if(index == 1 && incoherentEvent == 1)
619                                                NSingleIncoherent += 1
620                                        endif
621                                        if(index == 1 && coherentEvent == 1)
622                                                NSingleCoherent += 1
623                                        endif
624                                        if(index == 2 && coherentEvent == 1 && incoherentEvent == 0)
625                                                NDoubleCoherent += 1
626                                        endif
627                                        if(index > 1)
628                                                NMultipleScatter += 1
629                                        endif
630                                        if(coherentEvent >= 1 && incoherentEvent == 0)
631                                                NCoherentEvents += 1
632                                        endif
633                                        if(coherentEvent > 1 && incoherentEvent == 0)
634                                                NMultipleCoherent += 1
635                                        endif
636                                       
637                                       
638                                        //Print "n1,index (x,y) = ",n1,index, xpixel,ypixel
639                                       
640                                else    // if neutron escaped without interacting
641                               
642                                        // then it must be a transmitted neutron
643                                        // don't need to calculate, just increment the proper counters
644                                       
645                                        MC_linear_data[xCtr+xx/pixsize][yCtr+yy/pixsize] += 1
646                                        isOn += 1
647                                        nt[0] += 1
648                                       
649                                endif           //if interacted
650                        ENDIF
651                while (!done)
652        while(n1 < imon)
653
654//      Print "Monte Carlo Done"
655        results[0] = n1
656        results[1] = n2
657        results[2] = isOn
658        results[3] = NScatterEvents             //sum of # of times that neutrons scattered (coh+incoh)
659        results[4] = NSingleCoherent            //# of events that are single, coherent
660        results[5] = NMultipleCoherent  //# of scattered neutrons that are coherently scattered more than once
661        results[6] = NMultipleScatter           //# of scattered neutrons that are scattered more than once (coh and/or incoh)
662        results[7] = NCoherentEvents            //# of scattered neutrons that are scattered coherently one or more times
663       
664//      Print "# absorbed = ",n3
665
666//      trans_th = exp(-sig_total*thick)
667//      TRANS_exp = (N1-N2) / N1                        // Transmission
668        // dsigma/domega assuming isotropic scattering, with no absorption.
669//      Print "trans_exp = ",trans_exp
670//      Print "total # of neutrons reaching 2D detector",isOn
671//      Print "fraction of incident neutrons reaching detector ",isOn/iMon
672       
673//      Print "Total number of neutrons = ",N1
674//      Print "Total number of neutrons that interact = ",N2
675//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
676//      results[2] = N2                                         //number that scatter
677//      results[3] = isOn - MC_linear_data[xCtr][yCtr]                  //# scattered reaching detector minus zero angle
678
679       
680//      Tabs = (N1-N3)/N1
681//      Ttot = (N1-N2)/N1
682//      I1_sumI = NN[0]/(N2-N3)
683//      Print "Tabs = ",Tabs
684//      Print "Transmitted neutrons = ",Ttot
685//      results[8] = Ttot
686//      Print "I1 / all I1 = ", I1_sumI
687
688End
689////////        end of main function for calculating multiple scattering
690
691
692// returns the random deviate as a wave
693// and the total SAS cross-section [1/cm] sig_sas
694Function CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
695        FUNCREF SANSModelAAO_MCproto func
696        WAVE coef
697        Variable lam
698        String outWave
699        Variable &SASxs
700
701        Variable nPts_ran=10000,qu
702        qu = 4*pi/lam           
703       
704// hard-wired into the Simulation directory rather than the SAS folder.
705// plotting resolution-smeared models won't work any other way
706        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated"
707        WAVE Gq = root:Simulation:gQ
708        WAVE xw = root:Simulation:xw
709        SetScale/I x (0+1e-4),qu*(1-1e-10),"", Gq,xw                    //don't start at zero or run up all the way to qu to avoid numerical errors
710
711/// if all of the coefficients are well-behaved, then the last point is the background
712// and I can set it to zero here (only for the calculation)
713        Duplicate/O coef,tmp_coef
714        Variable num=numpnts(coef)
715        tmp_coef[num-1] = 0
716       
717        xw=x                                                                                            //for the AAO
718        func(tmp_coef,Gq,xw)                                                                    //call as AAO
719
720//      Gq = x*Gq                                                                                                       // SAS approximation
721        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact
722        //
723        //
724        Integrate/METH=1 Gq/D=Gq_INT
725       
726//      SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1]                 //if the approximation is used
727        SASxs = lam*Gq_INT[nPts_ran-1]
728       
729        Gq_INT /= Gq_INT[nPts_ran-1]
730       
731        Duplicate/O Gq_INT $outWave
732
733        return(0)
734End
735
736// returns the random deviate as a wave
737// and the total SAS cross-section [1/cm] sig_sas
738//
739// uses a log spacing of x for better coverage
740// downside is that it doesn't use built-in integrate, which is automatically cumulative
741//
742// --- Currently does not work - since the usage of the random deviate in the MC routine is based on the
743// wave properties of ran_dev, that is it must have the proper scaling and be equally spaced.
744//
745// -- not really sure that this will solve any of the problems with some functions (notably those with power-laws)
746// giving unreasonably large SAS cross sections. (>>10)
747//
748Function CalculateRandomDeviate_log(func,coef,lam,outWave,SASxs)
749        FUNCREF SANSModelAAO_MCproto func
750        WAVE coef
751        Variable lam
752        String outWave
753        Variable &SASxs
754
755        Variable nPts_ran=1000,qu,qmin,ii
756        qmin=1e-5
757        qu = 4*pi/lam           
758
759// hard-wired into the Simulation directory rather than the SAS folder.
760// plotting resolution-smeared models won't work any other way
761        Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw             // if these waves are 1000 pts, the results are "pixelated"
762        WAVE Gq = root:Simulation:gQ
763        WAVE xw = root:Simulation:xw
764//      SetScale/I x (0+1e-4),qu*(1-1e-10),"", Gq,xw                    //don't start at zero or run up all the way to qu to avoid numerical errors
765        xw =  alog(log(qmin) + x*((log(qu)-log(qmin))/nPts_ran))
766
767/// if all of the coefficients are well-behaved, then the last point is the background
768// and I can set it to zero here (only for the calculation)
769        Duplicate/O coef,tmp_coef
770        Variable num=numpnts(coef)
771        tmp_coef[num-1] = 0
772       
773        func(tmp_coef,Gq,xw)                                                                    //call as AAO
774        Gq = Gq*sin(2*asin(xw/qu))/sqrt(1-(xw/qu))                      // exact
775
776       
777        Duplicate/O Gq Gq_INT
778        Gq_INT = 0
779        for(ii=0;ii<nPts_ran;ii+=1)
780                Gq_INT[ii] = AreaXY(xw,Gq,qmin,xw[ii])
781        endfor
782       
783        SASxs = lam*Gq_INT[nPts_ran-1]
784       
785        Gq_INT /= Gq_INT[nPts_ran-1]
786       
787        Duplicate/O Gq_INT $outWave
788
789        return(0)
790End
791
792ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
793        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
794
795        Variable theta,dy,dx,qx,qy
796        //decompose to qx,qy
797        qx = testQ*cos(testPhi)
798        qy = testQ*sin(testPhi)
799
800        //convert qx,qy to pixel locations relative to # of pixels x, y from center
801        theta = 2*asin(qy*lam/4/pi)
802        dy = sdd*tan(theta)
803        yPixel = round(yCtr + dy/pixSize)
804       
805        theta = 2*asin(qx*lam/4/pi)
806        dx = sdd*tan(theta)
807        xPixel = round(xCtr + dx/pixSize)
808
809        NVAR pixelsX = root:myGlobals:gNPixelsX
810        NVAR pixelsY = root:myGlobals:gNPixelsY
811       
812        //if on detector, return xPix and yPix values, otherwise -1
813        if(yPixel > pixelsY || yPixel < 0)
814                yPixel = -1
815        endif
816        if(xPixel > pixelsX || xPixel < 0)
817                xPixel = -1
818        endif
819       
820        return(0)
821End
822
823Function MC_CheckFunctionAndCoef(funcStr,coefStr)
824        String funcStr,coefStr
825       
826        SVAR/Z listStr=root:Packages:NIST:coefKWStr
827        if(SVAR_Exists(listStr) == 1)
828                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
829                if(cmpstr("",properCoefStr)==0)
830                        return(0)               //false, no match found, so properCoefStr is returned null
831                endif
832                if(cmpstr(coefStr,properCoefStr)==0)
833                        return(1)               //true, the coef is the correct match
834                endif
835        endif
836        return(0)                       //false, wrong coef
837End
838
839Function/S MC_getFunctionCoef(funcStr)
840        String funcStr
841
842        SVAR/Z listStr=root:Packages:NIST:coefKWStr
843        String coefStr=""
844        if(SVAR_Exists(listStr) == 1)
845                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
846        endif
847        return(coefStr)
848End
849
850Function SANSModelAAO_MCproto(w,yw,xw)
851        Wave w,yw,xw
852       
853        Print "in SANSModelAAO_MCproto function"
854        return(1)
855end
856
857Function/S MC_FunctionPopupList()
858        String list,tmp
859        list = User_FunctionPopupList()
860       
861        //simplify the display, forcing smeared calculations behind the scenes
862        tmp = FunctionList("Smear*",";","NPARAMS:1")            //smeared dependency calculations
863        list = RemoveFromList(tmp, list,";")
864
865
866        if(strlen(list)==0)
867                list = "No functions plotted"
868        endif
869       
870        list = SortList(list)
871       
872        return(list)
873End             
874
875
876//Function Scat_Angle(Ran,R_DAB,wavelength)
877//      Variable Ran,r_dab,wavelength
878//
879//      Variable qq,arg,theta
880//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
881//      arg = qq*wavelength/(4*pi)
882//      If (arg < 1.0)
883//              theta = 2.*asin(arg)
884//      else
885//              theta = pi
886//      endif
887//      Return (theta)
888//End
889
890//calculates new direction (xyz) from an old direction
891//theta and phi don't change
892ThreadSafe Function NewDirection(vx,vy,vz,theta,phi)
893        Variable &vx,&vy,&vz
894        Variable theta,phi
895       
896        Variable err=0,vx0,vy0,vz0
897        Variable nx,ny,mag_xy,tx,ty,tz
898       
899        //store old direction vector
900        vx0 = vx
901        vy0 = vy
902        vz0 = vz
903       
904        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
905        if(mag_xy < 1e-12)
906                //old vector lies along beam direction
907                nx = 0
908                ny = 1
909                tx = 1
910                ty = 0
911                tz = 0
912        else
913                Nx = -Vy0 / Mag_XY
914                Ny = Vx0 / Mag_XY
915                Tx = -Vz0*Vx0 / Mag_XY
916                Ty = -Vz0*Vy0 / Mag_XY
917                Tz = Mag_XY
918        endif
919       
920        //new scattered direction vector
921        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
922        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
923        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
924       
925        Return(err)
926End
927
928ThreadSafe Function path_len(aval,sig_tot)
929        Variable aval,sig_tot
930       
931        Variable retval
932       
933        retval = -1*ln(1-aval)/sig_tot
934       
935        return(retval)
936End
937
938// globals are initialized in SASCALC.ipf
939// coordinates if I make this part of the panel - but this breaks other things...
940//
941//Proc MC_SASCALC()
942////    PauseUpdate; Silent 1           // building window...
943//
944////    NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator"
945//      SetVariable MC_setvar0,pos={491,73},size={144,15},bodyWidth=80,title="# of neutrons"
946//      SetVariable MC_setvar0,format="%5.4g"
947//      SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
948//      SetVariable MC_setvar0_1,pos={491,119},size={131,15},bodyWidth=60,title="Thickness (cm)"
949//      SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
950//      SetVariable MC_setvar0_2,pos={491,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
951//      SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
952//      SetVariable MC_setvar0_3,pos={491,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
953//      SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
954//      PopupMenu MC_popup0,pos={476,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
955//      PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
956//      Button MC_button0,pos={480,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
957//      Button MC_button1,pos={644,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
958//      SetVariable setvar0_3,pos={568,484},size={50,20},disable=1
959//      GroupBox group0,pos={478,42},size={267,130},title="Monte Carlo"
960//      SetVariable cntVar,pos={653,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)"
961//      SetVariable cntVar,format="%d"
962//      SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime
963//     
964//      String fldrSav0= GetDataFolder(1)
965//      SetDataFolder root:Packages:NIST:SAS:
966//      Edit/W=(476,217,746,450)/HOST=#  results_desc,results
967//      ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
968//      SetDataFolder fldrSav0
969//      RenameWindow #,T_results
970//      SetActiveSubwindow ##
971//EndMacro
972
973// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right
974// for various reasons...
975Proc MC_SASCALC()
976
977        // when opening the panel, set the raw counts check to 1
978        root:Packages:NIST:SAS:gRawCounts = 1
979       
980        PauseUpdate; Silent 1           // building window...
981        NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator"
982        DoWindow/C MC_SASCALC
983       
984        SetVariable MC_setvar0,pos={26,73},size={144,15},title="# of neutrons"
985        SetVariable MC_setvar0,format="%5.4g"
986        SetVariable MC_setvar0,limits={0,inf,100},value= root:Packages:NIST:SAS:gImon
987        SetVariable MC_setvar0_1,pos={26,119},size={140,15},title="Thickness (cm)"
988        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
989        SetVariable MC_setvar0_2,pos={26,96},size={165,15},title="Incoherent XS (1/cm)"
990        SetVariable MC_setvar0_2,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
991        SetVariable MC_setvar0_3,pos={26,142},size={155,15},title="Sample Radius (cm)"
992        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
993        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function"
994        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
995        Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation"
996        Button MC_button0,fColor=(3,52428,1)
997        Button MC_button1,pos={17,208},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
998        Button MC_button3,pos={210,94},size={25,20},proc=showIncohXSHelp,title="?"
999        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1
1000        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo"
1001        SetVariable cntVar,pos={185,73},size={90,15},proc=CountTimeSetVarProc,title="time(s)"
1002        SetVariable cntVar,format="%d"
1003        SetVariable cntVar,limits={1,3600,1},value= root:Packages:NIST:SAS:gCntTime
1004        Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX"
1005        CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts
1006        CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1007        CheckBox check0_2,pos={216,199+19},size={60,14},title="Beam Stop in",variable= root:Packages:NIST:SAS:gBeamStopIn
1008        CheckBox check0_3,pos={216,199+2*19},size={60,14},title="use XOP",variable= root:Packages:NIST:SAS:gUse_MC_XOP
1009       
1010        String fldrSav0= GetDataFolder(1)
1011        SetDataFolder root:Packages:NIST:SAS:
1012        Edit/W=(344,23,606,248)/HOST=#  results_desc,results
1013        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150
1014        SetDataFolder fldrSav0
1015        RenameWindow #,T_results
1016        SetActiveSubwindow ##
1017EndMacro
1018
1019
1020//
1021Proc showIncohXSHelp(ctrlName): ButtonControl
1022        String ctrlName
1023        DisplayHelpTopic/K=1/Z "Approximate Incoherent Cross Section"
1024        if(V_flag !=0)
1025                DoAlert 0,"The SANS Simulation Help file could not be found"
1026        endif
1027end
1028
1029
1030Function CountTimeSetVarProc(sva) : SetVariableControl
1031        STRUCT WMSetVariableAction &sva
1032
1033        switch( sva.eventCode )
1034                case 1: // mouse up
1035                case 2: // Enter key
1036                case 3: // Live update
1037                        Variable dval = sva.dval
1038
1039                        // get the neutron flux, multiply, and reset the global for # neutrons
1040                        NVAR imon=root:Packages:NIST:SAS:gImon
1041                        imon = dval*beamIntensity()
1042                       
1043                        break
1044        endswitch
1045
1046        return 0
1047End
1048
1049
1050Function MC_ModelPopMenuProc(pa) : PopupMenuControl
1051        STRUCT WMPopupAction &pa
1052
1053        switch( pa.eventCode )
1054                case 2: // mouse up
1055                        Variable popNum = pa.popNum
1056                        String popStr = pa.popStr
1057                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1058                        gStr = popStr
1059                       
1060                        break
1061        endswitch
1062
1063        return 0
1064End
1065
1066Function MC_DoItButtonProc(ba) : ButtonControl
1067        STRUCT WMButtonAction &ba
1068
1069        switch( ba.eventCode )
1070                case 2: // mouse up
1071                        // click code here
1072                        NVAR doMC = root:Packages:NIST:SAS:gDoMonteCarlo
1073                        doMC = 1
1074                        ReCalculateInten(1)
1075                        doMC = 0                //so the next time won't be MC
1076                        break
1077        endswitch
1078
1079        return 0
1080End
1081
1082
1083Function MC_Display2DButtonProc(ba) : ButtonControl
1084        STRUCT WMButtonAction &ba
1085
1086        switch( ba.eventCode )
1087                case 2: // mouse up
1088                        // click code here
1089                        Execute "ChangeDisplay(\"SAS\")"
1090                        break
1091        endswitch
1092
1093        return 0
1094End
1095
1096// after a 2d data image is averaged in the usual way, take the waves and generate a "fake" folder of the 1d
1097// data, to appear as if it was loaded from a real data file.
1098//
1099// ---- use FakeUSANSDataFolder() if you want to fake a 1D USANS data set ----
1100//
1101Function        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,dataFolder)
1102        WAVE qval,aveint,sigave,sigmaQ,qbar,fSubs
1103        String dataFolder
1104
1105        String baseStr=dataFolder
1106        if(DataFolderExists("root:"+baseStr))
1107                SetDataFolder $("root:"+baseStr)
1108        else
1109                NewDataFolder/S $("root:"+baseStr)
1110        endif
1111
1112        ////overwrite the existing data, if it exists
1113        Duplicate/O qval, $(baseStr+"_q")
1114        Duplicate/O aveint, $(baseStr+"_i")
1115        Duplicate/O sigave, $(baseStr+"_s")
1116
1117
1118        // make a resolution matrix for SANS data
1119        Variable np=numpnts(qval)
1120        Make/D/O/N=(np,4) $(baseStr+"_res")
1121        Wave res=$(baseStr+"_res")
1122       
1123        res[][0] = sigmaQ[p]            //sigQ
1124        res[][1] = qBar[p]              //qBar
1125        res[][2] = fSubS[p]             //fShad
1126        res[][3] = qval[p]              //Qvalues
1127       
1128        // keep a copy of everything in SAS too... the smearing wrapper function looks for
1129        // data in folders based on waves it is passed - an I lose control of that
1130        Duplicate/O res, $("root:Packages:NIST:SAS:"+baseStr+"_res")
1131        Duplicate/O qval,  $("root:Packages:NIST:SAS:"+baseStr+"_q")
1132        Duplicate/O aveint,  $("root:Packages:NIST:SAS:"+baseStr+"_i")
1133        Duplicate/O sigave,  $("root:Packages:NIST:SAS:"+baseStr+"_s")
1134       
1135        //clean up             
1136        SetDataFolder root:
1137       
1138End
1139
1140// writes out a VAX binary data file
1141// automatically generates a name
1142// will prompt for the sample label
1143//
1144// currently hard-wired for SAS data folder
1145//
1146// A later call to Write_VAXRaw_Data() will check for the simulation data, and
1147// convert ABS simulated data to raw counts if necessary
1148//
1149Function SaveAsVAXButtonProc(ctrlName,[runIndex,simLabel])
1150        String ctrlName
1151        Variable runIndex
1152        String simLabel
1153
1154       
1155        //first, check to be sure that the data is RAW counts before trying to save the VAX format
1156        // the easy answer is to abort, but I possibly could unscale the data and get back to integer
1157        // counts.
1158        // A later call to Write_VAXRaw_Data() will check for the simulation data, and
1159        // convert ABS simulated data to raw counts if necessary
1160//      NVAR isRAW = root:Packages:NIST:SAS:gRawCounts
1161//      if(!isRAW)
1162//              Abort "The simulation must be in RAW counts for it to be saved as RAW VAX format"
1163//      endif
1164       
1165        // if default parameters were passed in, use them
1166        // if not, set them to "bad" values so that the user will be prompted later     
1167        NVAR autoSaveIndex = root:Packages:NIST:SAS:gAutoSaveIndex
1168        SVAR autoSaveLabel = root:Packages:NIST:SAS:gAutoSaveLabel
1169       
1170        // Determine if the optional parameters were supplied
1171        if( ParamIsDefault(runIndex))           //==1 if parameter was NOT specified
1172                print "runIndex not specified"
1173                autoSaveIndex=0                                 // 0 == bad value, test for this later
1174        else
1175                autoSaveIndex=runIndex
1176        endif
1177       
1178        if( ParamIsDefault(simLabel))           //==1 if parameter was NOT specified
1179                print "simLabel not specified"
1180                autoSaveLabel=""                                        // "" == bad value, test for this later
1181        else
1182                autoSaveLabel=simLabel
1183        endif
1184       
1185        String fullpath="",destStr=""
1186        Variable refnum
1187       
1188        fullpath = Write_RawData_File("SAS","",0)
1189       
1190        // write out the results into a text file
1191        destStr = "root:Packages:NIST:SAS:"
1192        SetDataFolder $destStr
1193
1194        WAVE results=results
1195        WAVE/T results_desc=results_desc
1196       
1197        //check each wave
1198        If(!(WaveExists(results)))
1199                Abort "results DNExist WriteVAXData()"
1200        Endif
1201        If(!(WaveExists(results_desc)))
1202                Abort "results_desc DNExist WriteVAXData()"
1203        Endif
1204       
1205        NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1206        String str = ""
1207        sprintf str,"%30s\t\t%g seconds\r","MonteCarlo Simulation time = ",actSimTime
1208               
1209        Open refNum as fullpath+".txt"
1210                wfprintf refNum, "%30s\t\t%g\r",results_desc,results
1211                FBinWrite refNum,str
1212                FStatus refNum
1213                FSetPos refNum,V_logEOF
1214        Close refNum
1215       
1216        ///////////////////////////////
1217       
1218        // could also automatically do the average here, but probably not worth the the effort...
1219       
1220        SetDataFolder root:
1221       
1222        return(0)
1223End
1224
1225// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1226// and qmin and qmax
1227//
1228//
1229// still some question of the corners and number of pixels per q-bin
1230Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1231        wave ran_dev
1232        Variable Qmin,Qmax
1233       
1234        Variable r1,r2,frac
1235        r1=x2pnt(ran_dev,Qmin)
1236        r2=x2pnt(ran_dev,Qmax)
1237       
1238        // no normalization needed - the full q-range is defined as [0,1]
1239        frac = ran_dev[r2] - ran_dev[r1]
1240       
1241        return frac
1242End
1243
1244
1245/// called in SASCALC:ReCalculateInten()
1246Function Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1247        String funcStr
1248        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1249
1250        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
1251        WAVE rw=root:Packages:NIST:SAS:realsRead
1252       
1253// Try to nicely exit from a threading error, if possible
1254        Variable err=0
1255        if(!exists("root:myGlobals:gThreadGroupID"))
1256                Variable/G root:myGlobals:gThreadGroupID=0
1257        endif
1258        NVAR mt=root:myGlobals:gThreadGroupID
1259
1260        if(mt!=0)       //there was an error with the stopping of the threads, possibly user abort
1261                err = ThreadGroupRelease(mt)
1262                Print "threading err = ",err
1263                if(err == 0)
1264                        // all *should* be OK
1265                else
1266                        return(0)
1267                endif
1268        endif
1269
1270        NVAR imon = root:Packages:NIST:SAS:gImon
1271        NVAR thick = root:Packages:NIST:SAS:gThick
1272        NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1273        NVAR r2 = root:Packages:NIST:SAS:gR2
1274
1275        // do the simulation here, or not
1276        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength,deltaLam
1277        String coefStr,abortStr,str
1278
1279        r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
1280        xCtr = rw[16]
1281        yCtr = rw[17]
1282        sdd = rw[18]*100                //conver header of [m] to [cm]
1283        pixSize = rw[10]/10             // convert pix size in mm to cm
1284        wavelength = rw[26]
1285        deltaLam = rw[27]
1286        coefStr = MC_getFunctionCoef(funcStr)
1287       
1288        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1289                doMonteCarlo = 0                //we're getting out now, reset the flag so we don't continually end up here
1290                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1291        endif
1292       
1293        Variable sig_sas
1294//              FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)         //a wrapper for the structure version
1295        FUNCREF SANSModelAAO_MCproto func=$(funcStr)            //unsmeared
1296        WAVE results = root:Packages:NIST:SAS:results
1297        WAVE linear_data = root:Packages:NIST:SAS:linear_data
1298        WAVE data = root:Packages:NIST:SAS:data
1299
1300        results = 0
1301        linear_data = 0
1302       
1303        CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1304        if(sig_sas > 100)
1305                sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1306                Abort abortStr
1307        endif
1308       
1309        WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1310       
1311        Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1312        Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1313        Make/O/D/N=15 root:Packages:NIST:SAS:inputWave=0
1314       
1315        WAVE nt = root:Packages:NIST:SAS:nt
1316        WAVE j1 = root:Packages:NIST:SAS:j1
1317        WAVE j2 = root:Packages:NIST:SAS:j2
1318        WAVE nn = root:Packages:NIST:SAS:nn
1319        WAVE inputWave = root:Packages:NIST:SAS:inputWave
1320
1321        inputWave[0] = imon
1322        inputWave[1] = r1
1323        inputWave[2] = r2
1324        inputWave[3] = xCtr
1325        inputWave[4] = yCtr
1326        inputWave[5] = sdd
1327        inputWave[6] = pixSize
1328        inputWave[7] = thick
1329        inputWave[8] = wavelength
1330        inputWave[9] = sig_incoh
1331        inputWave[10] = sig_sas
1332        inputWave[11] = deltaLam
1333//      inputWave[] 12-14 are currently unused
1334
1335        linear_data = 0         //initialize
1336
1337        Variable t0,trans
1338       
1339        // get a time estimate, and give the user a chance to exit if they're unsure.
1340        t0 = stopMStimer(-2)
1341        inputWave[0] = 1000
1342        NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP                //if zero, will use non-threaded Igor code
1343       
1344        if(useXOP)
1345                //use a single thread, otherwise time is dominated by overhead
1346                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1347        else
1348                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1349        endif
1350       
1351        t0 = (stopMSTimer(-2) - t0)*1e-6
1352        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation)
1353        t0 *= 2         //empirical correction
1354        Print "Estimated Simulation time (s) = ",t0
1355       
1356// to correct for detector efficiency, send only the fraction of neutrons that are actually counted     
1357        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1358        NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1359        NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
1360
1361        inputWave[0] = imon     * detectorEff                   //reset number of input neutrons before full simulation
1362       
1363        if(t0>SimTimeWarn)
1364                sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1365                DoAlert 1,str
1366                if(V_flag == 2)
1367                        doMonteCarlo = 0
1368                        reCalculateInten(1)             //come back in and do the smeared calculation
1369                        return(0)
1370                endif
1371        endif
1372       
1373        linear_data = 0         //initialize
1374// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1375// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1376// use a different rng. This works. 1.75x speedup.     
1377        t0 = stopMStimer(-2)
1378
1379        if(useXOP)
1380                Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1381        else
1382                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1383        endif
1384       
1385        t0 = (stopMSTimer(-2) - t0)*1e-6
1386        Printf  "MC sim time = %g seconds\r",t0
1387        actSimTime = t0
1388       
1389        trans = results[8]                      //(n1-n2)/n1
1390        if(trans == 0)
1391                trans = 1
1392        endif
1393
1394//      Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1395       
1396//              linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike
1397//              Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1398
1399        // or simulate a beamstop
1400        NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn              //if zero, beam stop is "out", as in a transmission measurement
1401       
1402        Variable rad=beamstopDiam()/2           //beamstop radius in cm
1403        if(MC_BS_in)
1404                rad /= 0.5                              //convert cm to pixels
1405                rad += 0.                                       // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1406                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1407                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1408                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1
1409               
1410                linear_data *= tmp_mask
1411        endif
1412       
1413        results[9] = sum(linear_data,-inf,inf)
1414        //              Print "counts on detector not behind beamstop = ",results[9]
1415       
1416        // convert to absolute scale
1417        Variable kappa          //,beaminten = beamIntensity()
1418//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten)
1419        kappa = thick*(pixSize/sdd)^2*trans*iMon
1420       
1421        //use kappa to get back to counts => linear_data = round(linear_data*kappa)
1422        Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1423       
1424        NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1425        if(!rawCts)                     //go ahead and do the abs scaling to the linear_data
1426                linear_data = linear_data / kappa
1427                linear_data /= detectorEff
1428        endif
1429       
1430        // add a signature to the data file to tag as a simulation
1431        linear_data[0][0] = 1
1432        linear_data[2][0] = 1
1433        linear_data[0][2] = 1
1434        linear_data[1][1] = 1
1435        linear_data[2][2] = 1
1436        linear_data[1][0] = 0
1437        linear_data[2][1] = 0
1438        linear_data[0][1] = 0
1439        linear_data[1][2] = 0
1440
1441        linear_data[0][3] = 0
1442        linear_data[1][3] = 0
1443        linear_data[2][3] = 0
1444        linear_data[3][3] = 0
1445        linear_data[3][2] = 0
1446        linear_data[3][1] = 0
1447        linear_data[3][0] = 0
1448
1449                       
1450        data = linear_data
1451        // re-average the 2D data
1452        S_CircularAverageTo1D("SAS")
1453       
1454        // put the new result into the simulation folder
1455        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")     
1456                               
1457        // simulate the empty cell scattering, only in 1D
1458        Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs)
1459        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1460        Print "Sample Simulation (2D) CR = ",results[9]/ctTime
1461        if(WinType("SANS_Data") ==1)
1462                Execute "ChangeDisplay(\"SAS\")"                //equivalent to pressing "Show 2D"
1463        endif
1464
1465        return(0)
1466end
1467
1468//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1469ThreadSafe Function MC_FindPhi(vx,vy)
1470        variable vx,vy
1471       
1472        variable phi
1473       
1474        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1475       
1476        // special cases
1477        if(vx==0 && vy > 0)
1478                return(pi/2)
1479        endif
1480        if(vx==0 && vy < 0)
1481                return(3*pi/2)
1482        endif
1483        if(vx >= 0 && vy == 0)
1484                return(0)
1485        endif
1486        if(vx < 0 && vy == 0)
1487                return(pi)
1488        endif
1489       
1490       
1491        if(vx > 0 && vy > 0)
1492                return(phi)
1493        endif
1494        if(vx < 0 && vy > 0)
1495                return(phi + pi)
1496        endif
1497        if(vx < 0 && vy < 0)
1498                return(phi + pi)
1499        endif
1500        if( vx > 0 && vy < 0)
1501                return(phi + 2*pi)
1502        endif
1503       
1504        return(phi)
1505end
1506
1507
1508
1509
1510
1511Proc Sim_1D_Panel()
1512        PauseUpdate; Silent 1           // building window...
1513        NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1514        DoWindow/C Sim_1D_Panel
1515       
1516        SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1517        SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1518        SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1519        SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1520        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1521        SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1522
1523        SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1524        SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1525        SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1526
1527        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1528        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1529        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1530        Button MC_button0,fColor=(3,52428,1)
1531        Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data"
1532        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1533        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1534        CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1535        CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1536       
1537// a box for the results
1538        GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1539        ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1540        ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1541        ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1542        ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1543        ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1544        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1545        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1546        ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1547        ValDisplay valdisp0_4,pos={326,145},size={220,13},title="Multiple Coherent Scattering"
1548        ValDisplay valdisp0_4,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_MultScattFraction
1549        // set the flags here -- do the simulation, but not 2D
1550       
1551        root:Packages:NIST:SAS:doSimulation     = 1     // == 1 if 1D simulated data, 0 if other from the checkbox
1552        root:Packages:NIST:SAS:gDoMonteCarlo     = 0  // == 1 if 2D MonteCarlo set by hidden flag
1553
1554       
1555EndMacro
1556
1557Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1558        STRUCT WMSetVariableAction &sva
1559
1560        switch( sva.eventCode )
1561                case 1: // mouse up
1562                case 2: // Enter key
1563                case 3: // Live update
1564                        Variable dval = sva.dval
1565
1566                        ReCalculateInten(1)
1567                       
1568                        break
1569        endswitch
1570
1571        return 0
1572End
1573
1574Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1575        STRUCT WMSetVariableAction &sva
1576
1577        switch( sva.eventCode )
1578                case 1: // mouse up
1579                case 2: // Enter key
1580                case 3: // Live update
1581                        Variable dval = sva.dval
1582
1583                        ReCalculateInten(1)
1584                       
1585                        break
1586        endswitch
1587
1588        return 0
1589End
1590
1591Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1592        STRUCT WMSetVariableAction &sva
1593
1594        switch( sva.eventCode )
1595                case 1: // mouse up
1596                case 2: // Enter key
1597                case 3: // Live update
1598                        Variable dval = sva.dval
1599
1600                        ReCalculateInten(1)
1601                       
1602                        break
1603        endswitch
1604
1605        return 0
1606End
1607
1608
1609Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1610        STRUCT WMPopupAction &pa
1611
1612        switch( pa.eventCode )
1613                case 2: // mouse up
1614                        Variable popNum = pa.popNum
1615                        String popStr = pa.popStr
1616                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1617                        gStr = popStr
1618                       
1619                        break
1620        endswitch
1621
1622        return 0
1623End
1624
1625
1626Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1627        STRUCT WMButtonAction &ba
1628
1629        switch( ba.eventCode )
1630                case 2: // mouse up
1631               
1632                        ReCalculateInten(1)
1633                       
1634                        break
1635        endswitch
1636
1637        return 0
1638End
1639
1640
1641//
1642//
1643//
1644Function Save_1DSimData(ctrlName) : ButtonControl
1645        String ctrlName
1646
1647        String type="SAS",fullpath=""
1648        Variable dialog=1               //=1 will present dialog for name
1649       
1650        String destStr=""
1651        destStr = "root:Packages:NIST:"+type
1652       
1653        Variable refNum
1654        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
1655        String fname,ave="C",hdrStr1="",hdrStr2=""
1656        Variable step=1
1657       
1658        If(1)
1659                //setup a "fake protocol" wave, sice I have no idea of the current state of the data
1660                Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol
1661                Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol"
1662                SVAR funcStr = root:Packages:NIST:SAS:gFuncStr
1663                String junk="****SIMULATED DATA****"
1664
1665                //stick in the fake protocol...
1666                NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1667                NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts                       //summed counts (simulated)
1668                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate
1669                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt
1670                NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1671       
1672                SIMProtocol[0] = "*** SIMULATED DATA from "+funcStr+" ***"
1673                SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime)
1674                SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts)
1675                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR)
1676                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat)
1677                SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat)
1678                SIMProtocol[6] = ""
1679                SIMProtocol[7] = "*** SIMULATED DATA ***"
1680                //set the global
1681                String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol"
1682        Endif
1683       
1684       
1685        NVAR useXMLOutput = root:Packages:NIST:gXML_Write
1686       
1687        if (useXMLOutput == 1)
1688                WriteXMLWaves_W_Protocol(type,"",1)
1689        else
1690                WriteWaves_W_Protocol(type,"",1)                //"" is an empty path, 1 will force a dialog
1691        endif
1692       
1693//     
1694//      //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
1695//      WAVE intw=$(destStr + ":integersRead")
1696//      WAVE rw=$(destStr + ":realsRead")
1697//      WAVE/T textw=$(destStr + ":textRead")
1698//      WAVE qvals =$(destStr + ":qval")
1699//      WAVE inten=$(destStr + ":aveint")
1700//      WAVE sig=$(destStr + ":sigave")
1701//      WAVE qbar = $(destStr + ":QBar")
1702//      WAVE sigmaq = $(destStr + ":SigmaQ")
1703//      WAVE fsubs = $(destStr + ":fSubS")
1704//
1705//      SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
1706//      Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
1707//     
1708//      //check each wave
1709//      If(!(WaveExists(intw)))
1710//              Abort "intw DNExist Save_1DSimData()"
1711//      Endif
1712//      If(!(WaveExists(rw)))
1713//              Abort "rw DNExist Save_1DSimData()"
1714//      Endif
1715//      If(!(WaveExists(textw)))
1716//              Abort "textw DNExist Save_1DSimData()"
1717//      Endif
1718//      If(!(WaveExists(qvals)))
1719//              Abort "qvals DNExist Save_1DSimData()"
1720//      Endif
1721//      If(!(WaveExists(inten)))
1722//              Abort "inten DNExist Save_1DSimData()"
1723//      Endif
1724//      If(!(WaveExists(sig)))
1725//              Abort "sig DNExist Save_1DSimData()"
1726//      Endif
1727//      If(!(WaveExists(qbar)))
1728//              Abort "qbar DNExist Save_1DSimData()"
1729//      Endif
1730//      If(!(WaveExists(sigmaq)))
1731//              Abort "sigmaq DNExist Save_1DSimData()"
1732//      Endif
1733//      If(!(WaveExists(fsubs)))
1734//              Abort "fsubs DNExist Save_1DSimData()"
1735//      Endif
1736//      If(!(WaveExists(proto)))
1737//              Abort "current protocol wave DNExist Save_1DSimData()"
1738//      Endif
1739//
1740//      //strings can be too long to print-- must trim to 255 chars
1741//      Variable ii,num=8
1742//      Make/O/T/N=(num) tempShortProto
1743//      for(ii=0;ii<num;ii+=1)
1744//              tempShortProto[ii] = (proto[ii])[0,240]
1745//      endfor
1746//     
1747//      if(dialog)
1748//              PathInfo/S catPathName
1749//              fullPath = DoSaveFileDialog("Save data as")
1750//              If(cmpstr(fullPath,"")==0)
1751//                      //user cancel, don't write out a file
1752//                      Close/A
1753//                      Abort "no data file was written"
1754//              Endif
1755//              //Print "dialog fullpath = ",fullpath
1756//      Endif
1757//     
1758//      NVAR monCt = root:Packages:NIST:SAS:gImon
1759//      NVAR thick = root:Packages:NIST:SAS:gThick
1760//      NVAR trans = root:Packages:NIST:SAS:gSamTrans                   //for 1D, default value
1761//     
1762//
1763//     
1764//      hdrStr1 = num2str(monCt)+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
1765//      hdrStr1 += "     "+num2str(trans)+"     "+num2str(thick) + ave +"   "+num2str(step) + "\r\n"
1766//
1767//      hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
1768//      hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+"ORNL  " + "\r\n"
1769//     
1770//      //actually open the file here
1771//      Open refNum as fullpath
1772//     
1773//      //write out the standard header information
1774//      fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +"  "+ time())
1775//      fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA"
1776//      fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
1777//      fprintf refnum,hdrStr1
1778//      fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
1779//      fprintf refnum,hdrStr2
1780////    fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
1781//
1782//      //insert protocol information here
1783//      //-1 list of sample files
1784//      //0 - bkg
1785//      //1 - emp
1786//      //2 - div
1787//      //3 - mask
1788//      //4 - abs params c2-c5
1789//      //5 - average params
1790//      fprintf refnum, "SAM: %s\r\n",tempShortProto[0]
1791//      fprintf refnum, "BGD: %s\r\n",tempShortProto[1]
1792//      fprintf refnum, "EMP: %s\r\n",tempShortProto[2]
1793//      fprintf refnum, "DIV: %s\r\n",tempShortProto[3]
1794//      fprintf refnum, "MASK: %s\r\n",tempShortProto[4]
1795//      fprintf refnum, "ABS: %s\r\n",tempShortProto[5]
1796//      fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6]
1797//     
1798//      //write out the data columns
1799//      fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"
1800//      wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
1801//     
1802//      Close refnum
1803//     
1804//      SetDataFolder root:             //(redundant)
1805//     
1806//      //write confirmation of write operation to history area
1807//      Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
1808//      KillWaves/Z tempShortProto
1809//
1810//      //clear the stuff that was created for case of saving files
1811//      If(1)
1812//              Killwaves/Z root:myGlobals:Protocols:SIMProtocol
1813//              String/G root:myGlobals:Protocols:gProtoStr = ""
1814//      Endif
1815//     
1816       
1817        return(0)
1818       
1819End
1820
1821
1822/// called in SASCALC:ReCalculateInten()
1823Function Simulate_1D(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1824        String funcStr
1825        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1826
1827        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1828        String coefStr,abortStr,str     
1829
1830        FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
1831        FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
1832        coefStr = MC_getFunctionCoef(funcStr)
1833       
1834        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1835                Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
1836        endif
1837       
1838        Wave inten=$"root:Simulation:Simulation_i"              // this will exist and send the smeared calculation to the corect DF
1839       
1840        // the resolution-smeared intensity is calculated, including the incoherent background
1841        func($coefStr,inten,qval)
1842
1843        NVAR imon = root:Packages:NIST:SAS:gImon
1844        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1845        NVAR thick = root:Packages:NIST:SAS:gThick
1846        NVAR trans = root:Packages:NIST:SAS:gSamTrans
1847        NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
1848        NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
1849        NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
1850        NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
1851        NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1852        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1853       
1854        WAVE rw=root:Packages:NIST:SAS:realsRead
1855        WAVE nCells=root:Packages:NIST:SAS:nCells                               
1856                                       
1857        pixSize = rw[10]/10             // convert pix size in mm to cm
1858        sdd = rw[18]*100                //convert header of [m] to [cm]
1859        wavelength = rw[26]             // in 1/A
1860       
1861        imon = beamIntensity()*ctTime
1862       
1863        // calculate the scattering cross section simply to be able to estimate the transmission
1864        Variable sig_sas=0
1865       
1866        // remember that the random deviate is the coherent portion ONLY - the incoherent background is
1867        // subtracted before the calculation.
1868        CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
1869       
1870//                              if(sig_sas > 100)
1871//                                      sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1872//                              endif
1873
1874        // calculate the multiple scattering fraction for display (10/2009)
1875        Variable ii,nMax=10,tau
1876        mScat=0
1877        tau = thick*sig_sas
1878        // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
1879        // neutrons out of those that were scattered
1880        for(ii=2;ii<nMax;ii+=1)
1881                mScat += tau^(ii)/factorial(ii)
1882//              print tau^(ii)/factorial(ii)
1883        endfor
1884        estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
1885        mscat *= (estTrans)/(1-estTrans)
1886
1887        if(mScat > 0.1)         //  /Z to supress error if this was a 1D calc with the 2D panel open
1888                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=(65535,32768,32768)
1889        else
1890                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=0
1891        endif
1892
1893
1894        Print "Sig_sas = ",sig_sas
1895       
1896        Duplicate/O qval prob_i,countsInAnnulus
1897       
1898        // not needed - nCells takes care of this when the error is correctly calculated
1899//                              Duplicate/O qval circle_fraction,rval,nCells_expected
1900//                              rval = sdd*tan(2*asin(qval*wavelength/4/pi))            //radial distance in cm
1901//                              nCells_expected = 2*pi*rval/pixSize                                     //does this need to be an integer?
1902//                              circle_fraction = nCells / nCells_expected
1903       
1904                               
1905//                              prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten                       //probability of a neutron in q-bin(i) that has nCells
1906        prob_i = trans*thick*(pixSize/sdd)^2*inten                      //probability of a neutron in q-bin(i)
1907       
1908        Variable P_on = sum(prob_i,-inf,inf)
1909        Print "P_on = ",P_on
1910       
1911//                              fracScat = P_on
1912        fracScat = 1-estTrans
1913       
1914//                              aveint = (Imon)*prob_i / circle_fraction / nCells_expected
1915// added correction for detector efficiency, since SASCALC is flux on sample
1916        aveint = (Imon)*prob_i*detectorEff
1917
1918        countsInAnnulus = aveint*nCells
1919        SimDetCts = sum(countsInAnnulus,-inf,inf)
1920        estDetCR = SimDetCts/ctTime
1921       
1922       
1923        NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1924        NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1925       
1926        // this is where the number of cells comes in - the calculation of the error bars
1927        // sigma[i] = SUM(sigma[ij]^2) / nCells^2
1928        // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
1929        // then...
1930        sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99
1931       
1932        // add in random error in aveint based on the sigave
1933        if(addNoise)
1934                aveint += gnoise(sigave)
1935        endif
1936
1937        // signature in the standard deviation, do this after the noise is added
1938        // start at 10 to be out of the beamstop (makes for nicer plotting)
1939        // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
1940        sigave[10,50;10] = 10*sigave[p]
1941
1942        // convert to absolute scale, remembering to un-correct for the detector efficiency
1943        if(doABS)
1944                Variable kappa = thick*(pixSize/sdd)^2*trans*iMon
1945                aveint /= kappa
1946                sigave /= kappa
1947                aveint /= detectorEff
1948                sigave /= detectorEff
1949        endif
1950                               
1951       
1952        Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs)
1953        Print "Sample Simulation (1D) CR = ",estDetCR
1954       
1955        return(0)
1956End
1957
1958/// for testing only
1959Function testProbability(sas,thick)
1960        Variable sas,thick
1961       
1962        Variable tau,trans,p2p,p3p,p4p
1963       
1964        tau = sas*thick
1965        trans = exp(-tau)
1966       
1967        Print "tau = ",tau
1968        Print "trans = ",trans
1969       
1970        p2p = tau^2/factorial(2)*trans/(1-trans)
1971        p3p = tau^3/factorial(3)*trans/(1-trans)
1972        p4p = tau^4/factorial(4)*trans/(1-trans)
1973       
1974        Print "double scattering = ",p2p
1975        Print "triple scattering = ",p3p
1976        Print "quadruple scattering = ",p4p
1977       
1978        Variable ii,nMax=10,mScat=0
1979        for(ii=2;ii<nMax;ii+=1)
1980                mScat += tau^(ii)/factorial(ii)
1981//              print tau^(ii)/factorial(ii)
1982        endfor
1983        mscat *= (Trans)/(1-Trans)
1984
1985        Print "Total fraction of multiple scattering = ",mScat
1986
1987        return(mScat)
1988End
1989
1990
1991
1992//
1993// -- empirical simulation of the scattering from an empty quartz cell + background (combined)
1994// - there is little difference vs. the empty cell alone.
1995//
1996// - data was fit to the TwoLevel model, which fits rather nicely
1997//
1998Function Simulate_1D_EmptyCell(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1999        String funcStr
2000        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
2001
2002        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
2003        String coefStr,abortStr,str     
2004
2005        FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
2006        FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
2007       
2008        Make/O/D root:Packages:NIST:SAS:coef_Empty = {1,1.84594,714.625,5e-08,2.63775,0.0223493,3.94009,0.0153754,1.72127,0}
2009        WAVE coefW = root:Packages:NIST:SAS:coef_Empty
2010       
2011        Wave samInten=$"root:Simulation:Simulation_i"           // this will exist and send the smeared calculation to the corect DF
2012        Duplicate/O samInten, root:Simulation:Simulation_EC_i
2013        Wave inten_EC=$"root:Simulation:Simulation_EC_i"
2014
2015        // the resolution-smeared intensity of the empty cell
2016        func(coefW,inten_EC,qval)
2017
2018        NVAR imon = root:Packages:NIST:SAS:gImon
2019        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
2020//      NVAR thick = root:Packages:NIST:SAS:gThick
2021        NVAR trans = root:Packages:NIST:SAS:gSamTrans
2022//      NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
2023//      NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
2024//      NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
2025//      NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
2026//      NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
2027        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
2028
2029//      use local variables here for the Empty cell - maybe use globals later, if I really want to save them
2030// - here, just print them out for now
2031        Variable SimDetCts,estDetCR,fracScat,estTrans,mScat,thick
2032       
2033// for two 1/16" quartz windows, thick = 0.32 cm
2034        thick = 0.32
2035       
2036        WAVE rw=root:Packages:NIST:SAS:realsRead
2037        WAVE nCells=root:Packages:NIST:SAS:nCells                               
2038                                       
2039        pixSize = rw[10]/10             // convert pix size in mm to cm
2040        sdd = rw[18]*100                //convert header of [m] to [cm]
2041        wavelength = rw[26]             // in 1/A
2042       
2043        imon = beamIntensity()*ctTime
2044       
2045        // calculate the scattering cross section simply to be able to estimate the transmission
2046        Variable sig_sas=0
2047       
2048        // remember that the random deviate is the coherent portion ONLY - the incoherent background is
2049        // subtracted before the calculation.
2050        CalculateRandomDeviate(funcUnsmeared,coefW,wavelength,"root:Packages:NIST:SAS:ran_dev_EC",sig_sas)
2051
2052        // calculate the multiple scattering fraction for display (10/2009)
2053        Variable ii,nMax=10,tau
2054        mScat=0
2055        tau = thick*sig_sas
2056        // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
2057        // neutrons out of those that were scattered
2058        for(ii=2;ii<nMax;ii+=1)
2059                mScat += tau^(ii)/factorial(ii)
2060//              print tau^(ii)/factorial(ii)
2061        endfor
2062        estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
2063        mscat *= (estTrans)/(1-estTrans)
2064
2065       
2066        Duplicate/O qval prob_i_EC,countsInAnnulus_EC
2067       
2068        prob_i_EC = trans*thick*(pixSize/sdd)^2*inten_EC                        //probability of a neutron in q-bin(i)
2069       
2070        Variable P_on = sum(prob_i_EC,-inf,inf)
2071//      Print "P_on = ",P_on
2072       
2073        fracScat = 1-estTrans
2074       
2075// added correction for detector efficiency, since SASCALC is flux on sample
2076        Duplicate/O aveint root:Packages:NIST:SAS:aveint_EC,root:Packages:NIST:SAS:sigave_EC
2077        WAVE aveint_EC = root:Packages:NIST:SAS:aveint_EC
2078        WAVE sigave_EC = root:Packages:NIST:SAS:sigave_EC
2079        aveint_EC = (Imon)*prob_i_EC*detectorEff
2080
2081        countsInAnnulus_EC = aveint_EC*nCells
2082        SimDetCts = sum(countsInAnnulus_EC,-inf,inf)
2083        estDetCR = SimDetCts/ctTime
2084       
2085//      Print "Empty Cell Sig_sas = ",sig_sas
2086        Print "Empty Cell Count Rate : ",estDetCR
2087       
2088        NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
2089        NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
2090       
2091        // this is where the number of cells comes in - the calculation of the error bars
2092        // sigma[i] = SUM(sigma[ij]^2) / nCells^2
2093        // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
2094        // then...
2095        sigave_EC = sqrt(aveint_EC/nCells)              // corrected based on John's memo, from 8/9/99
2096       
2097        // add in random error in aveint based on the sigave
2098        if(addNoise)
2099                aveint_EC += gnoise(sigave_EC)
2100        endif
2101
2102        // signature in the standard deviation, do this after the noise is added
2103        // start at 10 to be out of the beamstop (makes for nicer plotting)
2104        // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
2105        sigave_EC[10,50;10] = 10*sigave_EC[p]
2106
2107        // convert to absolute scale, remembering to un-correct for the detector efficiency
2108        if(doABS)
2109                Variable kappa = thick*(pixSize/sdd)^2*trans*iMon
2110                aveint_EC /= kappa
2111                sigave_EC /= kappa
2112                aveint_EC /= detectorEff
2113                sigave_EC /= detectorEff
2114        endif
2115                               
2116                               
2117        return(0)
2118End
2119
2120
2121// instead of including the Beaucage model in everything, keep a local copy here
2122
2123//AAO version, uses XOP if available
2124// simply calls the original single point calculation with
2125// a wave assignment (this will behave nicely if given point ranges)
2126Function TwoLevel_EC(cw,yw,xw)
2127        Wave cw,yw,xw
2128       
2129#if exists("TwoLevelX")
2130        yw = TwoLevelX(cw,xw)
2131#else
2132        yw = fTwoLevel_EC(cw,xw)
2133#endif
2134        return(0)
2135End
2136
2137Function fTwoLevel_EC(w,x)
2138        Wave w
2139        Variable x
2140       
2141        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
2142        Variable erf1,erf2,prec=1e-15,scale
2143       
2144        //Rsub = Rs
2145        scale = w[0]
2146        G1 = w[1]       //equivalent to I(0)
2147        Rg1 = w[2]
2148        B1 = w[3]
2149        Pow1 = w[4]
2150        G2 = w[5]
2151        Rg2 = w[6]
2152        B2 = w[7]
2153        Pow2 = w[8]
2154        bkg = w[9]
2155       
2156        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
2157        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
2158        //Print erf1
2159       
2160        ans = G1*exp(-x*x*Rg1*Rg1/3)
2161        ans += B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
2162        ans += G2*exp(-x*x*Rg2*Rg2/3)
2163        ans += B2*(erf2^3/x)^Pow2
2164       
2165        if(x == 0)
2166                ans = G1 + G2
2167        endif
2168       
2169        ans *= scale
2170        ans += bkg
2171       
2172        Return(ans)
2173End
2174
2175
2176Function SmearedTwoLevel_EC(s)
2177        Struct ResSmearAAOStruct &s
2178
2179//      the name of your unsmeared model (AAO) is the first argument
2180        Smear_Model_20(TwoLevel_EC,s.coefW,s.xW,s.yW,s.resW)
2181
2182        return(0)
2183End
2184
2185//wrapper to calculate the smeared model as an AAO-Struct
2186// fills the struct and calls the ususal function with the STRUCT parameter
2187//
2188// used only for the dependency, not for fitting
2189//
2190Function fSmearedTwoLevel_EC(coefW,yW,xW)
2191        Wave coefW,yW,xW
2192       
2193        String str = getWavesDataFolder(yW,0)
2194        String DF="root:"+str+":"
2195       
2196        WAVE resW = $(DF+str+"_res")
2197       
2198        STRUCT ResSmearAAOStruct fs
2199        WAVE fs.coefW = coefW   
2200        WAVE fs.yW = yW
2201        WAVE fs.xW = xW
2202        WAVE fs.resW = resW
2203       
2204        Variable err
2205        err = SmearedTwoLevel_EC(fs)
2206       
2207        return (0)
2208End
2209
2210
2211
2212
2213//// this is a very simple example of how to script the MC simulation to run unattended
2214//
2215//  you need to supply for each "run":  the run index (you increment manually)
2216//                                                                                              the sample label (as a string)
2217//
2218// changing the various configuration paramters will have to be done on a case-by-case basis
2219// looking into SASCALC to see what is really changed,
2220// or the configuration parameters of the MC_SASCALC panel
2221//
2222//
2223//Function Script_2DMC()
2224//
2225//
2226//      NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
2227//      SimTimeWarn = 36000                     //sets the threshold for the warning dialog to 10 hours
2228//      STRUCT WMButtonAction ba
2229//      ba.eventCode = 2                        //fake mouse click on button
2230//     
2231//      NVAR detDist = root:Packages:NIST:SAS:gDetDist
2232//
2233//      detDist = 200           //set directly in cm
2234//      MC_DoItButtonProc(ba)
2235//      SaveAsVAXButtonProc("",runIndex=105,simLabel="this is run 105, SDD = 200")
2236//     
2237//      detDist = 300           //set directly in cm
2238//      MC_DoItButtonProc(ba)
2239//      SaveAsVAXButtonProc("",runIndex=106,simLabel="this is run 106, SDD = 300")
2240//
2241//      detDist = 400           //set directly in cm
2242//      MC_DoItButtonProc(ba)
2243//      SaveAsVAXButtonProc("",runIndex=107,simLabel="this is run 107, SDD = 400")
2244//     
2245//
2246// SimTimeWarn = 10             //back to 10 seconds for manual operation
2247//      return(0)
2248//end
Note: See TracBrowser for help on using the repository browser.