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

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

the 1d output of the simulations now is aware of the useXML flag. It also writes out the name of the function that was used for the simulation.

File size: 73.3 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//
1146Function SaveAsVAXButtonProc(ctrlName,[runIndex,simLabel])
1147        String ctrlName
1148        Variable runIndex
1149        String simLabel
1150
1151       
1152        // if default parameters were passed in, use them
1153        // if not, set them to "bad" values so that the user will be prompted later     
1154        NVAR autoSaveIndex = root:Packages:NIST:SAS:gAutoSaveIndex
1155        SVAR autoSaveLabel = root:Packages:NIST:SAS:gAutoSaveLabel
1156       
1157        // Determine if the optional parameters were supplied
1158        if( ParamIsDefault(runIndex))           //==1 if parameter was NOT specified
1159                print "runIndex not specified"
1160                autoSaveIndex=0                                 // 0 == bad value, test for this later
1161        else
1162                autoSaveIndex=runIndex
1163        endif
1164       
1165        if( ParamIsDefault(simLabel))           //==1 if parameter was NOT specified
1166                print "simLabel not specified"
1167                autoSaveLabel=""                                        // "" == bad value, test for this later
1168        else
1169                autoSaveLabel=simLabel
1170        endif
1171       
1172        String fullpath="",destStr=""
1173        Variable refnum
1174       
1175        fullpath = Write_RawData_File("SAS","",0)
1176       
1177        // write out the results into a text file
1178        destStr = "root:Packages:NIST:SAS:"
1179        SetDataFolder $destStr
1180
1181        WAVE results=results
1182        WAVE/T results_desc=results_desc
1183       
1184        //check each wave
1185        If(!(WaveExists(results)))
1186                Abort "results DNExist WriteVAXData()"
1187        Endif
1188        If(!(WaveExists(results_desc)))
1189                Abort "results_desc DNExist WriteVAXData()"
1190        Endif
1191       
1192        NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1193        String str = ""
1194        sprintf str,"%30s\t\t%g seconds\r","MonteCarlo Simulation time = ",actSimTime
1195               
1196        Open refNum as fullpath+".txt"
1197                wfprintf refNum, "%30s\t\t%g\r",results_desc,results
1198                FBinWrite refNum,str
1199                FStatus refNum
1200                FSetPos refNum,V_logEOF
1201        Close refNum
1202       
1203        ///////////////////////////////
1204       
1205        // could also automatically do the average here, but probably not worth the the effort...
1206       
1207        SetDataFolder root:
1208       
1209        return(0)
1210End
1211
1212// calculates the fraction of the scattering that reaches the detector, given the random deviate function
1213// and qmin and qmax
1214//
1215//
1216// still some question of the corners and number of pixels per q-bin
1217Function FractionReachingDetector(ran_dev,Qmin,Qmax)
1218        wave ran_dev
1219        Variable Qmin,Qmax
1220       
1221        Variable r1,r2,frac
1222        r1=x2pnt(ran_dev,Qmin)
1223        r2=x2pnt(ran_dev,Qmax)
1224       
1225        // no normalization needed - the full q-range is defined as [0,1]
1226        frac = ran_dev[r2] - ran_dev[r1]
1227       
1228        return frac
1229End
1230
1231
1232/// called in SASCALC:ReCalculateInten()
1233Function Simulate_2D_MC(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1234        String funcStr
1235        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1236
1237        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if 2D MonteCarlo set by hidden flag
1238        WAVE rw=root:Packages:NIST:SAS:realsRead
1239       
1240// Try to nicely exit from a threading error, if possible
1241        Variable err=0
1242        if(!exists("root:myGlobals:gThreadGroupID"))
1243                Variable/G root:myGlobals:gThreadGroupID=0
1244        endif
1245        NVAR mt=root:myGlobals:gThreadGroupID
1246
1247        if(mt!=0)       //there was an error with the stopping of the threads, possibly user abort
1248                err = ThreadGroupRelease(mt)
1249                Print "threading err = ",err
1250                if(err == 0)
1251                        // all *should* be OK
1252                else
1253                        return(0)
1254                endif
1255        endif
1256
1257        NVAR imon = root:Packages:NIST:SAS:gImon
1258        NVAR thick = root:Packages:NIST:SAS:gThick
1259        NVAR sig_incoh = root:Packages:NIST:SAS:gSig_incoh
1260        NVAR r2 = root:Packages:NIST:SAS:gR2
1261
1262        // do the simulation here, or not
1263        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength,deltaLam
1264        String coefStr,abortStr,str
1265
1266        r1 = rw[24]/2/10                // sample diameter convert diam in [mm] to radius in cm
1267        xCtr = rw[16]
1268        yCtr = rw[17]
1269        sdd = rw[18]*100                //conver header of [m] to [cm]
1270        pixSize = rw[10]/10             // convert pix size in mm to cm
1271        wavelength = rw[26]
1272        deltaLam = rw[27]
1273        coefStr = MC_getFunctionCoef(funcStr)
1274       
1275        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1276                doMonteCarlo = 0                //we're getting out now, reset the flag so we don't continually end up here
1277                Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
1278        endif
1279       
1280        Variable sig_sas
1281//              FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)         //a wrapper for the structure version
1282        FUNCREF SANSModelAAO_MCproto func=$(funcStr)            //unsmeared
1283        WAVE results = root:Packages:NIST:SAS:results
1284        WAVE linear_data = root:Packages:NIST:SAS:linear_data
1285        WAVE data = root:Packages:NIST:SAS:data
1286
1287        results = 0
1288        linear_data = 0
1289       
1290        CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS)
1291        if(sig_sas > 100)
1292                sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1293                Abort abortStr
1294        endif
1295       
1296        WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev"
1297       
1298        Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0
1299        Make/O/D/N=100 root:Packages:NIST:SAS:nn=0
1300        Make/O/D/N=15 root:Packages:NIST:SAS:inputWave=0
1301       
1302        WAVE nt = root:Packages:NIST:SAS:nt
1303        WAVE j1 = root:Packages:NIST:SAS:j1
1304        WAVE j2 = root:Packages:NIST:SAS:j2
1305        WAVE nn = root:Packages:NIST:SAS:nn
1306        WAVE inputWave = root:Packages:NIST:SAS:inputWave
1307
1308        inputWave[0] = imon
1309        inputWave[1] = r1
1310        inputWave[2] = r2
1311        inputWave[3] = xCtr
1312        inputWave[4] = yCtr
1313        inputWave[5] = sdd
1314        inputWave[6] = pixSize
1315        inputWave[7] = thick
1316        inputWave[8] = wavelength
1317        inputWave[9] = sig_incoh
1318        inputWave[10] = sig_sas
1319        inputWave[11] = deltaLam
1320//      inputWave[] 12-14 are currently unused
1321
1322        linear_data = 0         //initialize
1323
1324        Variable t0,trans
1325       
1326        // get a time estimate, and give the user a chance to exit if they're unsure.
1327        t0 = stopMStimer(-2)
1328        inputWave[0] = 1000
1329        NVAR useXOP = root:Packages:NIST:SAS:gUse_MC_XOP                //if zero, will use non-threaded Igor code
1330       
1331        if(useXOP)
1332                //use a single thread, otherwise time is dominated by overhead
1333                Monte_SANS_NotThreaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1334        else
1335                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1336        endif
1337       
1338        t0 = (stopMSTimer(-2) - t0)*1e-6
1339        t0 *= imon/1000/ThreadProcessorCount                    //projected time, in seconds (using threads for the calculation)
1340        t0 *= 2         //empirical correction
1341        Print "Estimated Simulation time (s) = ",t0
1342       
1343// to correct for detector efficiency, send only the fraction of neutrons that are actually counted     
1344        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1345        NVAR actSimTime = root:Packages:NIST:SAS:g_actSimTime
1346        NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
1347
1348        inputWave[0] = imon     * detectorEff                   //reset number of input neutrons before full simulation
1349       
1350        if(t0>SimTimeWarn)
1351                sprintf str,"The simulation will take approximately %d seconds.\r- Proceed?",t0
1352                DoAlert 1,str
1353                if(V_flag == 2)
1354                        doMonteCarlo = 0
1355                        reCalculateInten(1)             //come back in and do the smeared calculation
1356                        return(0)
1357                endif
1358        endif
1359       
1360        linear_data = 0         //initialize
1361// threading crashes!! - there must be some operation in the XOP that is not threadSafe. What, I don't know...
1362// I think it's the ran() calls, being "non-reentrant". So the XOP now defines two separate functions, that each
1363// use a different rng. This works. 1.75x speedup.     
1364        t0 = stopMStimer(-2)
1365
1366        if(useXOP)
1367                Monte_SANS_Threaded(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1368        else
1369                Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results)
1370        endif
1371       
1372        t0 = (stopMSTimer(-2) - t0)*1e-6
1373        Printf  "MC sim time = %g seconds\r",t0
1374        actSimTime = t0
1375       
1376        trans = results[8]                      //(n1-n2)/n1
1377        if(trans == 0)
1378                trans = 1
1379        endif
1380
1381//      Print "counts on detector, including transmitted = ",sum(linear_data,-inf,inf)
1382       
1383//              linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike
1384//              Print "counts on detector not transmitted = ",sum(linear_data,-inf,inf)
1385
1386        // or simulate a beamstop
1387        NVAR MC_BS_in = root:Packages:NIST:SAS:gBeamStopIn              //if zero, beam stop is "out", as in a transmission measurement
1388       
1389        Variable rad=beamstopDiam()/2           //beamstop radius in cm
1390        if(MC_BS_in)
1391                rad /= 0.5                              //convert cm to pixels
1392                rad += 0.                                       // (no - it cuts off the low Q artificially) add an extra pixel to each side to account for edge
1393                Duplicate/O linear_data,root:Packages:NIST:SAS:tmp_mask//,root:Packages:NIST:SAS:MC_linear_data
1394                WAVE tmp_mask = root:Packages:NIST:SAS:tmp_mask
1395                tmp_mask = (sqrt((p-xCtr)^2+(q-yCtr)^2) < rad) ? 0 : 1          //behind beamstop = 0, away = 1
1396               
1397                linear_data *= tmp_mask
1398        endif
1399       
1400        results[9] = sum(linear_data,-inf,inf)
1401        //              Print "counts on detector not behind beamstop = ",results[9]
1402       
1403        // convert to absolute scale
1404        Variable kappa          //,beaminten = beamIntensity()
1405//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*trans*(iMon/beaminten)
1406        kappa = thick*(pixSize/sdd)^2*trans*iMon
1407       
1408        //use kappa to get back to counts => linear_data = round(linear_data*kappa)
1409        Note/K linear_data ,"KAPPA="+num2str(kappa)+";"
1410       
1411        NVAR rawCts = root:Packages:NIST:SAS:gRawCounts
1412        if(!rawCts)                     //go ahead and do the abs scaling to the linear_data
1413                linear_data = linear_data / kappa
1414                linear_data /= detectorEff
1415        endif
1416       
1417        // add a signature to the data file to tag as a simulation
1418        linear_data[0][0] = 1
1419        linear_data[2][0] = 1
1420        linear_data[0][2] = 1
1421        linear_data[1][1] = 1
1422        linear_data[2][2] = 1
1423        linear_data[1][0] = 0
1424        linear_data[2][1] = 0
1425        linear_data[0][1] = 0
1426        linear_data[1][2] = 0
1427
1428        linear_data[0][3] = 0
1429        linear_data[1][3] = 0
1430        linear_data[2][3] = 0
1431        linear_data[3][3] = 0
1432        linear_data[3][2] = 0
1433        linear_data[3][1] = 0
1434        linear_data[3][0] = 0
1435
1436                       
1437        data = linear_data
1438        // re-average the 2D data
1439        S_CircularAverageTo1D("SAS")
1440       
1441        // put the new result into the simulation folder
1442        Fake1DDataFolder(qval,aveint,sigave,sigmaQ,qbar,fSubs,"Simulation")     
1443                               
1444        // simulate the empty cell scattering, only in 1D
1445        Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs)
1446        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1447        Print "Sample Simulation (2D) CR = ",results[9]/ctTime
1448        if(WinType("SANS_Data") ==1)
1449                Execute "ChangeDisplay(\"SAS\")"                //equivalent to pressing "Show 2D"
1450        endif
1451
1452        return(0)
1453end
1454
1455//phi is defined from +x axis, proceeding CCW around [0,2Pi]
1456ThreadSafe Function MC_FindPhi(vx,vy)
1457        variable vx,vy
1458       
1459        variable phi
1460       
1461        phi = atan(vy/vx)               //returns a value from -pi/2 to pi/2
1462       
1463        // special cases
1464        if(vx==0 && vy > 0)
1465                return(pi/2)
1466        endif
1467        if(vx==0 && vy < 0)
1468                return(3*pi/2)
1469        endif
1470        if(vx >= 0 && vy == 0)
1471                return(0)
1472        endif
1473        if(vx < 0 && vy == 0)
1474                return(pi)
1475        endif
1476       
1477       
1478        if(vx > 0 && vy > 0)
1479                return(phi)
1480        endif
1481        if(vx < 0 && vy > 0)
1482                return(phi + pi)
1483        endif
1484        if(vx < 0 && vy < 0)
1485                return(phi + pi)
1486        endif
1487        if( vx > 0 && vy < 0)
1488                return(phi + 2*pi)
1489        endif
1490       
1491        return(phi)
1492end
1493
1494
1495
1496
1497
1498Proc Sim_1D_Panel()
1499        PauseUpdate; Silent 1           // building window...
1500        NewPanel /W=(92,556,713,818)/K=1 as "1D SANS Simulator"
1501        DoWindow/C Sim_1D_Panel
1502       
1503        SetVariable cntVar,pos={26,68},size={160,15},title="Counting time(s)",format="%d"
1504        SetVariable cntVar,limits={1,36000,10},value= root:Packages:NIST:SAS:gCntTime
1505        SetVariable cntVar, proc=Sim_1D_CountTimeSetVarProc
1506        SetVariable MC_setvar0_1,pos={26,91},size={160,15},title="Thickness (cm)"
1507        SetVariable MC_setvar0_1,limits={0,inf,0.1},value= root:Packages:NIST:SAS:gThick
1508        SetVariable MC_setvar0_1, proc=Sim_1D_SamThickSetVarProc
1509
1510        SetVariable MC_setvar0_3,pos={26,114},size={160,15},title="Sample Transmission"
1511        SetVariable MC_setvar0_3,limits={0,1,0.01},value= root:Packages:NIST:SAS:gSamTrans
1512        SetVariable MC_setvar0_3, proc=Sim_1D_SamTransSetVarProc
1513
1514        PopupMenu MC_popup0,pos={13,13},size={165,20},proc=Sim_1D_ModelPopMenuProc,title="Model Function"
1515        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
1516        Button MC_button0,pos={17,181},size={130,20},proc=Sim_1D_DoItButtonProc,title="Do 1D Simulation"
1517        Button MC_button0,fColor=(3,52428,1)
1518        Button MC_button1,pos={17,211},size={150,20},proc=Save_1DSimData,title="Save Simulated Data"
1519        GroupBox group0,pos={15,42},size={280,130},title="Sample Setup"
1520        CheckBox check0_1,pos={216,179},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset
1521        CheckBox check0_2,pos={216,199},size={60,14},title="Abs scale?",variable= root:Packages:NIST:SAS:g_1D_DoABS
1522        CheckBox check0_3,pos={216,219},size={60,14},title="Noise?",variable= root:Packages:NIST:SAS:g_1D_AddNoise
1523       
1524// a box for the results
1525        GroupBox group1,pos={314,23},size={277,163},title="Simulation Results"
1526        ValDisplay valdisp0,pos={326,48},size={220,13},title="Total detector counts"
1527        ValDisplay valdisp0,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DTotCts
1528        ValDisplay valdisp0_1,pos={326,72},size={220,13},title="Estimated count rate (1/s)"
1529        ValDisplay valdisp0_1,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstDetCR
1530        ValDisplay valdisp0_2,pos={326,96},size={220,13},title="Fraction of beam scattered"
1531        ValDisplay valdisp0_2,limits={0,0,0},barmisc={0,1000},value= root:Packages:NIST:SAS:g_1DFracScatt
1532        ValDisplay valdisp0_3,pos={326,121},size={220,13},title="Estimated transmission"
1533        ValDisplay valdisp0_3,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_1DEstTrans
1534        ValDisplay valdisp0_4,pos={326,145},size={220,13},title="Multiple Coherent Scattering"
1535        ValDisplay valdisp0_4,limits={0,0,0},barmisc={0,1000},value=root:Packages:NIST:SAS:g_MultScattFraction
1536        // set the flags here -- do the simulation, but not 2D
1537       
1538        root:Packages:NIST:SAS:doSimulation     = 1     // == 1 if 1D simulated data, 0 if other from the checkbox
1539        root:Packages:NIST:SAS:gDoMonteCarlo     = 0  // == 1 if 2D MonteCarlo set by hidden flag
1540
1541       
1542EndMacro
1543
1544Function Sim_1D_CountTimeSetVarProc(sva) : SetVariableControl
1545        STRUCT WMSetVariableAction &sva
1546
1547        switch( sva.eventCode )
1548                case 1: // mouse up
1549                case 2: // Enter key
1550                case 3: // Live update
1551                        Variable dval = sva.dval
1552
1553                        ReCalculateInten(1)
1554                       
1555                        break
1556        endswitch
1557
1558        return 0
1559End
1560
1561Function Sim_1D_SamThickSetVarProc(sva) : SetVariableControl
1562        STRUCT WMSetVariableAction &sva
1563
1564        switch( sva.eventCode )
1565                case 1: // mouse up
1566                case 2: // Enter key
1567                case 3: // Live update
1568                        Variable dval = sva.dval
1569
1570                        ReCalculateInten(1)
1571                       
1572                        break
1573        endswitch
1574
1575        return 0
1576End
1577
1578Function Sim_1D_SamTransSetVarProc(sva) : SetVariableControl
1579        STRUCT WMSetVariableAction &sva
1580
1581        switch( sva.eventCode )
1582                case 1: // mouse up
1583                case 2: // Enter key
1584                case 3: // Live update
1585                        Variable dval = sva.dval
1586
1587                        ReCalculateInten(1)
1588                       
1589                        break
1590        endswitch
1591
1592        return 0
1593End
1594
1595
1596Function Sim_1D_ModelPopMenuProc(pa) : PopupMenuControl
1597        STRUCT WMPopupAction &pa
1598
1599        switch( pa.eventCode )
1600                case 2: // mouse up
1601                        Variable popNum = pa.popNum
1602                        String popStr = pa.popStr
1603                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
1604                        gStr = popStr
1605                       
1606                        break
1607        endswitch
1608
1609        return 0
1610End
1611
1612
1613Function Sim_1D_DoItButtonProc(ba) : ButtonControl
1614        STRUCT WMButtonAction &ba
1615
1616        switch( ba.eventCode )
1617                case 2: // mouse up
1618               
1619                        ReCalculateInten(1)
1620                       
1621                        break
1622        endswitch
1623
1624        return 0
1625End
1626
1627
1628//
1629//
1630//
1631Function Save_1DSimData(ctrlName) : ButtonControl
1632        String ctrlName
1633
1634        String type="SAS",fullpath=""
1635        Variable dialog=1               //=1 will present dialog for name
1636       
1637        String destStr=""
1638        destStr = "root:Packages:NIST:"+type
1639       
1640        Variable refNum
1641        String formatStr = "%15.4g %15.4g %15.4g %15.4g %15.4g %15.4g\r\n"
1642        String fname,ave="C",hdrStr1="",hdrStr2=""
1643        Variable step=1
1644       
1645        If(1)
1646                //setup a "fake protocol" wave, sice I have no idea of the current state of the data
1647                Make/O/T/N=8 root:myGlobals:Protocols:SIMProtocol
1648                Wave/T SIMProtocol = $"root:myGlobals:Protocols:SIMProtocol"
1649                SVAR funcStr = root:Packages:NIST:SAS:gFuncStr
1650                String junk="****SIMULATED DATA****"
1651
1652                //stick in the fake protocol...
1653                NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1654                NVAR totalCts = root:Packages:NIST:SAS:g_1DTotCts                       //summed counts (simulated)
1655                NVAR detCR = root:Packages:NIST:SAS:g_1DEstDetCR                // estimated detector count rate
1656                NVAR fractScat = root:Packages:NIST:SAS:g_1DFracScatt
1657                NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1658       
1659                SIMProtocol[0] = "*** SIMULATED DATA from "+funcStr+" ***"
1660                SIMProtocol[1] = "\tCounting time (s) = "+num2str(ctTime)
1661                SIMProtocol[2] = "\tTotal detector counts = "+num2str(totalCts)
1662                SIMProtocol[3] = "\tDetector countrate (1/s) = "+num2str(detCR)
1663                SIMProtocol[4] = "\tFraction of beam scattered coherently = "+num2str(fractScat)
1664                SIMProtocol[5] = "\tFraction of multiple coherent scattering = "+num2str(mScat)
1665                SIMProtocol[6] = ""
1666                SIMProtocol[7] = "*** SIMULATED DATA ***"
1667                //set the global
1668                String/G root:myGlobals:Protocols:gProtoStr = "SIMProtocol"
1669        Endif
1670       
1671       
1672        NVAR useXMLOutput = root:Packages:NIST:gXML_Write
1673       
1674        if (useXMLOutput == 1)
1675                WriteXMLWaves_W_Protocol(type,"",1)
1676        else
1677                WriteWaves_W_Protocol(type,"",1)                //"" is an empty path, 1 will force a dialog
1678        endif
1679       
1680//     
1681//      //*****these waves MUST EXIST, or IGOR Pro will crash, with a type 2 error****
1682//      WAVE intw=$(destStr + ":integersRead")
1683//      WAVE rw=$(destStr + ":realsRead")
1684//      WAVE/T textw=$(destStr + ":textRead")
1685//      WAVE qvals =$(destStr + ":qval")
1686//      WAVE inten=$(destStr + ":aveint")
1687//      WAVE sig=$(destStr + ":sigave")
1688//      WAVE qbar = $(destStr + ":QBar")
1689//      WAVE sigmaq = $(destStr + ":SigmaQ")
1690//      WAVE fsubs = $(destStr + ":fSubS")
1691//
1692//      SVAR gProtoStr = root:myGlobals:Protocols:gProtoStr
1693//      Wave/T proto=$("root:myGlobals:Protocols:"+gProtoStr)
1694//     
1695//      //check each wave
1696//      If(!(WaveExists(intw)))
1697//              Abort "intw DNExist Save_1DSimData()"
1698//      Endif
1699//      If(!(WaveExists(rw)))
1700//              Abort "rw DNExist Save_1DSimData()"
1701//      Endif
1702//      If(!(WaveExists(textw)))
1703//              Abort "textw DNExist Save_1DSimData()"
1704//      Endif
1705//      If(!(WaveExists(qvals)))
1706//              Abort "qvals DNExist Save_1DSimData()"
1707//      Endif
1708//      If(!(WaveExists(inten)))
1709//              Abort "inten DNExist Save_1DSimData()"
1710//      Endif
1711//      If(!(WaveExists(sig)))
1712//              Abort "sig DNExist Save_1DSimData()"
1713//      Endif
1714//      If(!(WaveExists(qbar)))
1715//              Abort "qbar DNExist Save_1DSimData()"
1716//      Endif
1717//      If(!(WaveExists(sigmaq)))
1718//              Abort "sigmaq DNExist Save_1DSimData()"
1719//      Endif
1720//      If(!(WaveExists(fsubs)))
1721//              Abort "fsubs DNExist Save_1DSimData()"
1722//      Endif
1723//      If(!(WaveExists(proto)))
1724//              Abort "current protocol wave DNExist Save_1DSimData()"
1725//      Endif
1726//
1727//      //strings can be too long to print-- must trim to 255 chars
1728//      Variable ii,num=8
1729//      Make/O/T/N=(num) tempShortProto
1730//      for(ii=0;ii<num;ii+=1)
1731//              tempShortProto[ii] = (proto[ii])[0,240]
1732//      endfor
1733//     
1734//      if(dialog)
1735//              PathInfo/S catPathName
1736//              fullPath = DoSaveFileDialog("Save data as")
1737//              If(cmpstr(fullPath,"")==0)
1738//                      //user cancel, don't write out a file
1739//                      Close/A
1740//                      Abort "no data file was written"
1741//              Endif
1742//              //Print "dialog fullpath = ",fullpath
1743//      Endif
1744//     
1745//      NVAR monCt = root:Packages:NIST:SAS:gImon
1746//      NVAR thick = root:Packages:NIST:SAS:gThick
1747//      NVAR trans = root:Packages:NIST:SAS:gSamTrans                   //for 1D, default value
1748//     
1749//
1750//     
1751//      hdrStr1 = num2str(monCt)+"  "+num2str(rw[26])+"       "+num2str(rw[19])+"     "+num2str(rw[18])
1752//      hdrStr1 += "     "+num2str(trans)+"     "+num2str(thick) + ave +"   "+num2str(step) + "\r\n"
1753//
1754//      hdrStr2 = num2str(rw[16])+"  "+num2str(rw[17])+"  "+num2str(rw[23])+"    "+num2str(rw[24])+"    "
1755//      hdrStr2 += num2str(rw[25])+"    "+num2str(rw[27])+"    "+num2str(rw[21])+"    "+"ORNL  " + "\r\n"
1756//     
1757//      //actually open the file here
1758//      Open refNum as fullpath
1759//     
1760//      //write out the standard header information
1761//      fprintf refnum,"FILE: %s\t\t CREATED: %s\r\n","SIMULATED DATA",(date() +"  "+ time())
1762//      fprintf refnum,"LABEL: %s\r\n","SIMULATED DATA"
1763//      fprintf refnum,"MON CNT   LAMBDA   DET ANG   DET DIST   TRANS   THICK   AVE   STEP\r\n"
1764//      fprintf refnum,hdrStr1
1765//      fprintf refnum,"BCENT(X,Y)   A1(mm)   A2(mm)   A1A2DIST(m)   DL/L   BSTOP(mm)   DET_TYP \r\n"
1766//      fprintf refnum,hdrStr2
1767////    fprintf refnum,headerFormat,rw[0],rw[26],rw[19],rw[18],rw[4],rw[5],ave,step
1768//
1769//      //insert protocol information here
1770//      //-1 list of sample files
1771//      //0 - bkg
1772//      //1 - emp
1773//      //2 - div
1774//      //3 - mask
1775//      //4 - abs params c2-c5
1776//      //5 - average params
1777//      fprintf refnum, "SAM: %s\r\n",tempShortProto[0]
1778//      fprintf refnum, "BGD: %s\r\n",tempShortProto[1]
1779//      fprintf refnum, "EMP: %s\r\n",tempShortProto[2]
1780//      fprintf refnum, "DIV: %s\r\n",tempShortProto[3]
1781//      fprintf refnum, "MASK: %s\r\n",tempShortProto[4]
1782//      fprintf refnum, "ABS: %s\r\n",tempShortProto[5]
1783//      fprintf refnum, "Average Choices: %s\r\n",tempShortProto[6]
1784//     
1785//      //write out the data columns
1786//      fprintf refnum,"The 6 columns are | Q (1/A) | I(Q) (1/cm) | std. dev. I(Q) (1/cm) | sigmaQ | meanQ | ShadowFactor|\r\n"
1787//      wfprintf refnum, formatStr, qvals,inten,sig,sigmaq,qbar,fsubs
1788//     
1789//      Close refnum
1790//     
1791//      SetDataFolder root:             //(redundant)
1792//     
1793//      //write confirmation of write operation to history area
1794//      Print "Averaged File written: ", GetFileNameFromPathNoSemi(fullPath)
1795//      KillWaves/Z tempShortProto
1796//
1797//      //clear the stuff that was created for case of saving files
1798//      If(1)
1799//              Killwaves/Z root:myGlobals:Protocols:SIMProtocol
1800//              String/G root:myGlobals:Protocols:gProtoStr = ""
1801//      Endif
1802//     
1803       
1804        return(0)
1805       
1806End
1807
1808
1809/// called in SASCALC:ReCalculateInten()
1810Function Simulate_1D(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1811        String funcStr
1812        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1813
1814        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1815        String coefStr,abortStr,str     
1816
1817        FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
1818        FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
1819        coefStr = MC_getFunctionCoef(funcStr)
1820       
1821        if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
1822                Abort "Function and coefficients do not match. You must plot the unsmeared function before simulation."
1823        endif
1824       
1825        Wave inten=$"root:Simulation:Simulation_i"              // this will exist and send the smeared calculation to the corect DF
1826       
1827        // the resolution-smeared intensity is calculated, including the incoherent background
1828        func($coefStr,inten,qval)
1829
1830        NVAR imon = root:Packages:NIST:SAS:gImon
1831        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
1832        NVAR thick = root:Packages:NIST:SAS:gThick
1833        NVAR trans = root:Packages:NIST:SAS:gSamTrans
1834        NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
1835        NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
1836        NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
1837        NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
1838        NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
1839        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
1840       
1841        WAVE rw=root:Packages:NIST:SAS:realsRead
1842        WAVE nCells=root:Packages:NIST:SAS:nCells                               
1843                                       
1844        pixSize = rw[10]/10             // convert pix size in mm to cm
1845        sdd = rw[18]*100                //convert header of [m] to [cm]
1846        wavelength = rw[26]             // in 1/A
1847       
1848        imon = beamIntensity()*ctTime
1849       
1850        // calculate the scattering cross section simply to be able to estimate the transmission
1851        Variable sig_sas=0
1852       
1853        // remember that the random deviate is the coherent portion ONLY - the incoherent background is
1854        // subtracted before the calculation.
1855        CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
1856       
1857//                              if(sig_sas > 100)
1858//                                      sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas
1859//                              endif
1860
1861        // calculate the multiple scattering fraction for display (10/2009)
1862        Variable ii,nMax=10,tau
1863        mScat=0
1864        tau = thick*sig_sas
1865        // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
1866        // neutrons out of those that were scattered
1867        for(ii=2;ii<nMax;ii+=1)
1868                mScat += tau^(ii)/factorial(ii)
1869//              print tau^(ii)/factorial(ii)
1870        endfor
1871        estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
1872        mscat *= (estTrans)/(1-estTrans)
1873
1874        if(mScat > 0.1)         //  /Z to supress error if this was a 1D calc with the 2D panel open
1875                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=(65535,32768,32768)
1876        else
1877                ValDisplay/Z valdisp0_4 win=Sim_1D_Panel,labelBack=0
1878        endif
1879
1880
1881        Print "Sig_sas = ",sig_sas
1882       
1883        Duplicate/O qval prob_i,countsInAnnulus
1884       
1885        // not needed - nCells takes care of this when the error is correctly calculated
1886//                              Duplicate/O qval circle_fraction,rval,nCells_expected
1887//                              rval = sdd*tan(2*asin(qval*wavelength/4/pi))            //radial distance in cm
1888//                              nCells_expected = 2*pi*rval/pixSize                                     //does this need to be an integer?
1889//                              circle_fraction = nCells / nCells_expected
1890       
1891                               
1892//                              prob_i = trans*thick*nCells*(pixSize/sdd)^2*inten                       //probability of a neutron in q-bin(i) that has nCells
1893        prob_i = trans*thick*(pixSize/sdd)^2*inten                      //probability of a neutron in q-bin(i)
1894       
1895        Variable P_on = sum(prob_i,-inf,inf)
1896        Print "P_on = ",P_on
1897       
1898//                              fracScat = P_on
1899        fracScat = 1-estTrans
1900       
1901//                              aveint = (Imon)*prob_i / circle_fraction / nCells_expected
1902// added correction for detector efficiency, since SASCALC is flux on sample
1903        aveint = (Imon)*prob_i*detectorEff
1904
1905        countsInAnnulus = aveint*nCells
1906        SimDetCts = sum(countsInAnnulus,-inf,inf)
1907        estDetCR = SimDetCts/ctTime
1908       
1909       
1910        NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
1911        NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
1912       
1913        // this is where the number of cells comes in - the calculation of the error bars
1914        // sigma[i] = SUM(sigma[ij]^2) / nCells^2
1915        // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
1916        // then...
1917        sigave = sqrt(aveint/nCells)            // corrected based on John's memo, from 8/9/99
1918       
1919        // add in random error in aveint based on the sigave
1920        if(addNoise)
1921                aveint += gnoise(sigave)
1922        endif
1923
1924        // signature in the standard deviation, do this after the noise is added
1925        // start at 10 to be out of the beamstop (makes for nicer plotting)
1926        // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
1927        sigave[10,50;10] = 10*sigave[p]
1928
1929        // convert to absolute scale, remembering to un-correct for the detector efficiency
1930        if(doABS)
1931                Variable kappa = thick*(pixSize/sdd)^2*trans*iMon
1932                aveint /= kappa
1933                sigave /= kappa
1934                aveint /= detectorEff
1935                sigave /= detectorEff
1936        endif
1937                               
1938       
1939        Simulate_1D_EmptyCell("TwoLevel_EC",aveint,qval,sigave,sigmaq,qbar,fsubs)
1940        Print "Sample Simulation (1D) CR = ",estDetCR
1941       
1942        return(0)
1943End
1944
1945/// for testing only
1946Function testProbability(sas,thick)
1947        Variable sas,thick
1948       
1949        Variable tau,trans,p2p,p3p,p4p
1950       
1951        tau = sas*thick
1952        trans = exp(-tau)
1953       
1954        Print "tau = ",tau
1955        Print "trans = ",trans
1956       
1957        p2p = tau^2/factorial(2)*trans/(1-trans)
1958        p3p = tau^3/factorial(3)*trans/(1-trans)
1959        p4p = tau^4/factorial(4)*trans/(1-trans)
1960       
1961        Print "double scattering = ",p2p
1962        Print "triple scattering = ",p3p
1963        Print "quadruple scattering = ",p4p
1964       
1965        Variable ii,nMax=10,mScat=0
1966        for(ii=2;ii<nMax;ii+=1)
1967                mScat += tau^(ii)/factorial(ii)
1968//              print tau^(ii)/factorial(ii)
1969        endfor
1970        mscat *= (Trans)/(1-Trans)
1971
1972        Print "Total fraction of multiple scattering = ",mScat
1973
1974        return(mScat)
1975End
1976
1977
1978
1979//
1980// -- empirical simulation of the scattering from an empty quartz cell + background (combined)
1981// - there is little difference vs. the empty cell alone.
1982//
1983// - data was fit to the TwoLevel model, which fits rather nicely
1984//
1985Function Simulate_1D_EmptyCell(funcStr,aveint,qval,sigave,sigmaq,qbar,fsubs)
1986        String funcStr
1987        WAVE aveint,qval,sigave,sigmaq,qbar,fsubs
1988
1989        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength
1990        String coefStr,abortStr,str     
1991
1992        FUNCREF SANSModelAAO_MCproto func=$("fSmeared"+funcStr)                 //a wrapper for the structure version
1993        FUNCREF SANSModelAAO_MCproto funcUnsmeared=$(funcStr)           //unsmeared
1994       
1995        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}
1996        WAVE coefW = root:Packages:NIST:SAS:coef_Empty
1997       
1998        Wave samInten=$"root:Simulation:Simulation_i"           // this will exist and send the smeared calculation to the corect DF
1999        Duplicate samInten, root:Simulation:Simulation_EC_i
2000        Wave inten_EC=$"root:Simulation:Simulation_EC_i"
2001
2002        // the resolution-smeared intensity of the empty cell
2003        func(coefW,inten_EC,qval)
2004
2005        NVAR imon = root:Packages:NIST:SAS:gImon
2006        NVAR ctTime = root:Packages:NIST:SAS:gCntTime
2007//      NVAR thick = root:Packages:NIST:SAS:gThick
2008        NVAR trans = root:Packages:NIST:SAS:gSamTrans
2009//      NVAR SimDetCts = root:Packages:NIST:SAS:g_1DTotCts                      //summed counts (simulated)
2010//      NVAR estDetCR = root:Packages:NIST:SAS:g_1DEstDetCR                     // estimated detector count rate
2011//      NVAR fracScat = root:Packages:NIST:SAS:g_1DFracScatt            // fraction of beam captured on detector
2012//      NVAR estTrans = root:Packages:NIST:SAS:g_1DEstTrans             // estimated transmission of sample
2013//      NVAR mScat = root:Packages:NIST:SAS:g_MultScattFraction
2014        NVAR detectorEff = root:Packages:NIST:SAS:g_detectorEff
2015
2016//      use local variables here for the Empty cell - maybe use globals later, if I really want to save them
2017// - here, just print them out for now
2018        Variable SimDetCts,estDetCR,fracScat,estTrans,mScat,thick
2019       
2020// for two 1/16" quartz windows, thick = 0.32 cm
2021        thick = 0.32
2022       
2023        WAVE rw=root:Packages:NIST:SAS:realsRead
2024        WAVE nCells=root:Packages:NIST:SAS:nCells                               
2025                                       
2026        pixSize = rw[10]/10             // convert pix size in mm to cm
2027        sdd = rw[18]*100                //convert header of [m] to [cm]
2028        wavelength = rw[26]             // in 1/A
2029       
2030        imon = beamIntensity()*ctTime
2031       
2032        // calculate the scattering cross section simply to be able to estimate the transmission
2033        Variable sig_sas=0
2034       
2035        // remember that the random deviate is the coherent portion ONLY - the incoherent background is
2036        // subtracted before the calculation.
2037        CalculateRandomDeviate(funcUnsmeared,coefW,wavelength,"root:Packages:NIST:SAS:ran_dev_EC",sig_sas)
2038
2039        // calculate the multiple scattering fraction for display (10/2009)
2040        Variable ii,nMax=10,tau
2041        mScat=0
2042        tau = thick*sig_sas
2043        // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
2044        // neutrons out of those that were scattered
2045        for(ii=2;ii<nMax;ii+=1)
2046                mScat += tau^(ii)/factorial(ii)
2047//              print tau^(ii)/factorial(ii)
2048        endfor
2049        estTrans = exp(-1*thick*sig_sas)                //thickness and sigma both in units of cm
2050        mscat *= (estTrans)/(1-estTrans)
2051
2052       
2053        Duplicate/O qval prob_i_EC,countsInAnnulus_EC
2054       
2055        prob_i_EC = trans*thick*(pixSize/sdd)^2*inten_EC                        //probability of a neutron in q-bin(i)
2056       
2057        Variable P_on = sum(prob_i_EC,-inf,inf)
2058//      Print "P_on = ",P_on
2059       
2060        fracScat = 1-estTrans
2061       
2062// added correction for detector efficiency, since SASCALC is flux on sample
2063        Duplicate/O aveint root:Packages:NIST:SAS:aveint_EC,root:Packages:NIST:SAS:sigave_EC
2064        WAVE aveint_EC = root:Packages:NIST:SAS:aveint_EC
2065        WAVE sigave_EC = root:Packages:NIST:SAS:sigave_EC
2066        aveint_EC = (Imon)*prob_i_EC*detectorEff
2067
2068        countsInAnnulus_EC = aveint_EC*nCells
2069        SimDetCts = sum(countsInAnnulus_EC,-inf,inf)
2070        estDetCR = SimDetCts/ctTime
2071       
2072//      Print "Empty Cell Sig_sas = ",sig_sas
2073        Print "Empty Cell Count Rate : ",estDetCR
2074       
2075        NVAR doABS = root:Packages:NIST:SAS:g_1D_DoABS
2076        NVAR addNoise = root:Packages:NIST:SAS:g_1D_AddNoise
2077       
2078        // this is where the number of cells comes in - the calculation of the error bars
2079        // sigma[i] = SUM(sigma[ij]^2) / nCells^2
2080        // and since in the simulation, SUM(sigma[ij]^2) = nCells*sigma[ij]^2 = nCells*Inten
2081        // then...
2082        sigave_EC = sqrt(aveint_EC/nCells)              // corrected based on John's memo, from 8/9/99
2083       
2084        // add in random error in aveint based on the sigave
2085        if(addNoise)
2086                aveint_EC += gnoise(sigave_EC)
2087        endif
2088
2089        // signature in the standard deviation, do this after the noise is added
2090        // start at 10 to be out of the beamstop (makes for nicer plotting)
2091        // end at 50 to leave the natural statistics at the end of the set (may have a total of 80+ points if no offset)
2092        sigave_EC[10,50;10] = 10*sigave_EC[p]
2093
2094        // convert to absolute scale, remembering to un-correct for the detector efficiency
2095        if(doABS)
2096                Variable kappa = thick*(pixSize/sdd)^2*trans*iMon
2097                aveint_EC /= kappa
2098                sigave_EC /= kappa
2099                aveint_EC /= detectorEff
2100                sigave_EC /= detectorEff
2101        endif
2102                               
2103                               
2104        return(0)
2105End
2106
2107
2108// instead of including the Beaucage model in everything, keep a local copy here
2109
2110//AAO version, uses XOP if available
2111// simply calls the original single point calculation with
2112// a wave assignment (this will behave nicely if given point ranges)
2113Function TwoLevel_EC(cw,yw,xw)
2114        Wave cw,yw,xw
2115       
2116#if exists("TwoLevelX")
2117        yw = TwoLevelX(cw,xw)
2118#else
2119        yw = fTwoLevel_EC(cw,xw)
2120#endif
2121        return(0)
2122End
2123
2124Function fTwoLevel_EC(w,x)
2125        Wave w
2126        Variable x
2127       
2128        Variable ans,G1,Rg1,B1,G2,Rg2,B2,Pow1,Pow2,bkg
2129        Variable erf1,erf2,prec=1e-15,scale
2130       
2131        //Rsub = Rs
2132        scale = w[0]
2133        G1 = w[1]       //equivalent to I(0)
2134        Rg1 = w[2]
2135        B1 = w[3]
2136        Pow1 = w[4]
2137        G2 = w[5]
2138        Rg2 = w[6]
2139        B2 = w[7]
2140        Pow2 = w[8]
2141        bkg = w[9]
2142       
2143        erf1 = erf( (x*Rg1/sqrt(6)) ,prec)
2144        erf2 = erf( (x*Rg2/sqrt(6)) ,prec)
2145        //Print erf1
2146       
2147        ans = G1*exp(-x*x*Rg1*Rg1/3)
2148        ans += B1*exp(-x*x*Rg2*Rg2/3)*(erf1^3/x)^Pow1
2149        ans += G2*exp(-x*x*Rg2*Rg2/3)
2150        ans += B2*(erf2^3/x)^Pow2
2151       
2152        if(x == 0)
2153                ans = G1 + G2
2154        endif
2155       
2156        ans *= scale
2157        ans += bkg
2158       
2159        Return(ans)
2160End
2161
2162
2163Function SmearedTwoLevel_EC(s)
2164        Struct ResSmearAAOStruct &s
2165
2166//      the name of your unsmeared model (AAO) is the first argument
2167        Smear_Model_20(TwoLevel_EC,s.coefW,s.xW,s.yW,s.resW)
2168
2169        return(0)
2170End
2171
2172//wrapper to calculate the smeared model as an AAO-Struct
2173// fills the struct and calls the ususal function with the STRUCT parameter
2174//
2175// used only for the dependency, not for fitting
2176//
2177Function fSmearedTwoLevel_EC(coefW,yW,xW)
2178        Wave coefW,yW,xW
2179       
2180        String str = getWavesDataFolder(yW,0)
2181        String DF="root:"+str+":"
2182       
2183        WAVE resW = $(DF+str+"_res")
2184       
2185        STRUCT ResSmearAAOStruct fs
2186        WAVE fs.coefW = coefW   
2187        WAVE fs.yW = yW
2188        WAVE fs.xW = xW
2189        WAVE fs.resW = resW
2190       
2191        Variable err
2192        err = SmearedTwoLevel_EC(fs)
2193       
2194        return (0)
2195End
2196
2197
2198
2199
2200//// this is a very simple example of how to script the MC simulation to run unattended
2201//
2202//  you need to supply for each "run":  the run index (you increment manually)
2203//                                                                                              the sample label (as a string)
2204//
2205// changing the various configuration paramters will have to be done on a case-by-case basis
2206// looking into SASCALC to see what is really changed,
2207// or the configuration parameters of the MC_SASCALC panel
2208//
2209//
2210//Function Script_2DMC()
2211//
2212//
2213//      NVAR SimTimeWarn = root:Packages:NIST:SAS:g_SimTimeWarn
2214//      SimTimeWarn = 36000                     //sets the threshold for the warning dialog to 10 hours
2215//      STRUCT WMButtonAction ba
2216//      ba.eventCode = 2                        //fake mouse click on button
2217//     
2218//      NVAR detDist = root:Packages:NIST:SAS:gDetDist
2219//
2220//      detDist = 200           //set directly in cm
2221//      MC_DoItButtonProc(ba)
2222//      SaveAsVAXButtonProc("",runIndex=105,simLabel="this is run 105, SDD = 200")
2223//     
2224//      detDist = 300           //set directly in cm
2225//      MC_DoItButtonProc(ba)
2226//      SaveAsVAXButtonProc("",runIndex=106,simLabel="this is run 106, SDD = 300")
2227//
2228//      detDist = 400           //set directly in cm
2229//      MC_DoItButtonProc(ba)
2230//      SaveAsVAXButtonProc("",runIndex=107,simLabel="this is run 107, SDD = 400")
2231//     
2232//
2233// SimTimeWarn = 10             //back to 10 seconds for manual operation
2234//      return(0)
2235//end
Note: See TracBrowser for help on using the repository browser.