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

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

Added new model with Yun Liu's TwoYukawa? structure factor. A one-yukawa version is also added, but only in the XOP. Note that the SQ calculation is AAO and NOT threadsafe, unlike all other fitting functions. The SQ calculation has not yet been added to any form factors.

Corrected Fractal_PolyCore to return P*S rather than just P.

Only allow a 2D simulation to be saved inn 2D VAX format if it is RAW counts, not absolute rescaled data. Otherwise the DP values are mangled on conversion to VAX integers. In the future, the data may be "unscaled" if there is an outcry.

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