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

Last change on this file since 434 was 434, checked in by srkline, 14 years ago

More tweaking of the MonteCarlo?, trying in vain to thread the XOP. As is, works fine with the XOP and single threading.

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