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

Last change on this file was 1011, checked in by srkline, 6 years ago

typo in #if statement

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