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

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

Added a new template "TotalTemplate?.pxt" that includes the package loader as part of the macros menu. from there, everything can be loaded (but nothing unloaded).

Many changes to SASCALC and MultiScatter? to make them correct, and play nice together. An XOP version is functional, but has not yet been added, and won't be until the MC simulation is checked for correctness. Without incoherent scattering, it looks good.

File size: 22.6 KB
RevLine 
[429]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.
[430]9// A lot of the setup, wave creation, and post-calculations are done in SASCALC->ReCalculateInten()
[429]10//
11//
12// - Most importantly, this needs to be checked for correctness of the MC simulation
[430]13// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation
[429]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"?
[430]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.
[429]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.
[430]27// X if no MC desired, still use the selected model
28// X better display of MC results on panel
[429]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??
[430]32// - make sure that the ratio of scattering coherent/incoherent is properly adjusted for the sample composition
33//   or the volume fraction of solvent.
[429]34//
[430]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// - 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// - 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//
[429]51
52
53
54//////////
55//    PROGRAM Monte_SANS
56//    PROGRAM simulates multiple SANS.
57//       revised 2/12/99  JGB
58//            added calculation of random deviate, and 2D 10/2008 SRK
59
60//    N1 = NUMBER OF INCIDENT NEUTRONS.
61//    N2 = NUMBER INTERACTED IN THE SAMPLE.
62//    N3 = NUMBER ABSORBED.
63//    THETA = SCATTERING ANGLE.
64
65//        IMON = 'Enter number of neutrons to use in simulation.'
66//        NUM_BINS = 'Enter number of THETA BINS TO use. (<5000).'
67//        R1 = 'Enter beam radius. (cm)'
68//        R2 = 'Enter sample radius. (cm)'
69//        thick = 'Enter sample thickness. (cm)'
70//        wavelength = 'Enter neutron wavelength. (A)'
71//        R0 = 'Enter sphere radius. (A)'
72//
73
[430]74Function Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results)
75        WAVE inputWave,ran_dev,nt,j1,j2,nn,MC_linear_data,results
[429]76
[430]77        Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
[429]78        Variable NUM_BINS,N_INDEX
79        Variable RHO,SIGSAS,SIGABS_0
80        Variable ii,jj,IND,idum,INDEX,IR,NQ
81        Variable qmax,theta_max,R_DAB,R0,BOUND,I0,q0,zpow
82        Variable N1,N2,N3,DTH,zz,tt,SIG_SINGLE,xx,yy,PHI,UU,SIG
83        Variable THETA,Ran,ll,D_OMEGA,RR,Tabs,Ttot,I1_sumI
84        //in John's implementation, he dimensioned the indices of the arrays to begin
85        // at 0, making things much easier for me...
86        //DIMENSION  NT(0:5000),J1(0:5000),J2(0:5000),NN(0:100)
[430]87
[429]88        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat
89        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single
90        Variable DONE,FIND_THETA,err            //used as logicals
91
[430]92        Variable Vx,Vy,Vz,Theta_z,qq
[429]93        Variable Sig_scat,Sig_abs,Ratio,Sig_total
94        Variable isOn=0,testQ,testPhi,xPixel,yPixel
95       
[430]96        imon = inputWave[0]
97        r1 = inputWave[1]
98        r2 = inputWave[2]
99        xCtr = inputWave[3]
100        yCtr = inputWave[4]
101        sdd = inputWave[5]
102        pixSize = inputWave[6]
103        thick = inputWave[7]
104        wavelength = inputWave[8]
105        sig_incoh = inputWave[9]
106        sig_sas = inputWave[10]
107       
[429]108//      SetRandomSeed 0.1               //to get a reproduceable sequence
109
110//scattering power and maximum qvalue to bin
111//      zpow = .1               //scattering power, calculated below
[430]112        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
[429]113        sigabs_0 = 0.0          // ignore absorption cross section/wavelength [1/(cm A)]
114        N_INDEX = 50            // maximum number of scattering events per neutron
115        num_bins = 200          //number of 1-D bins (not really used)
116       
117// my additions - calculate the randome deviate function as needed
118// and calculate the scattering power from the model function
119//
120        Variable left = leftx(ran_dev)
121        Variable delta = deltax(ran_dev)
122       
123//c       total SAS cross-section
124//
125// input a test value for the incoherent scattering from water
126//
127//      sig_incoh = 5.6                 //[1/cm] as calculated for H2O, now a parameter
128//
129//      SIG_SAS = zpow/thick
130        zpow = sig_sas*thick                    //since I now calculate the sig_sas from the model
131        SIG_ABS = SIGABS_0 * WAVElength
132        SIG_TOTAL =SIG_ABS + SIG_SAS + sig_incoh
133//      Print "The TOTAL XSECTION. (CM-1) is ",sig_total
134//      Print "The TOTAL SAS XSECTION. (CM-1) is ",sig_sas
135        results[0] = sig_total
136        results[1] = sig_sas
137//      RATIO = SIG_ABS / SIG_TOTAL
138        RATIO = sig_incoh / SIG_TOTAL
[430]139//!!!! the ratio is not yet properly weighted for the volume fractions of each component!!!
[429]140       
141//c       assuming theta = sin(theta)...OK
142        theta_max = wavelength*qmax/(2*pi)
143//C     SET Theta-STEP SIZE.
144        DTH = Theta_max/NUM_BINS
145//      Print "theta bin size = dth = ",dth
146
147//C     INITIALIZE COUNTERS.
148        N1 = 0
149        N2 = 0
150        N3 = 0
151
152//C     INITIALIZE ARRAYS.
153        j1 = 0
154        j2 = 0
155        nt = 0
156        nn=0
157       
158//C     MONITOR LOOP - looping over the number of incedent neutrons
159//note that zz, is the z-position in the sample - NOT the scattering power
160
161// NOW, start the loop, throwing neutrons at the sample.
162        do
163                Vx = 0.0                        // Initialize direction vector.
164                Vy = 0.0
165                Vz = 1.0
166               
167                Theta = 0.0             //      Initialize scattering angle.
168                Phi = 0.0                       //      Intialize azimuthal angle.
169                N1 += 1                 //      Increment total number neutrons counter.
170                DONE = 0                        //      True when neutron is absorbed or when  scattered out of the sample.
171                INDEX = 0                       //      Set counter for number of scattering events.
172                zz = 0.0                        //      Set entering dimension of sample.
173               
174                do                                      //      Makes sure position is within circle.
175                        ran = abs(enoise(1))            //[0,1]
176                        xx = 2.0*R1*(Ran-0.5)           //X beam position of neutron entering sample.
177                        ran = abs(enoise(1))            //[0,1]
178                        yy = 2.0*R1*(Ran-0.5)           //Y beam position ...
179                        RR = SQRT(xx*xx+yy*yy)          //Radial position of neutron in incident beam.
180                while(rr>r1)
181
182                do    //Scattering Loop, will exit when "done" == 1
183                                // keep scattering multiple times until the neutron exits the sample
184                        ran = abs(enoise(1))            //[0,1]  RANDOM NUMBER FOR DETERMINING PATH LENGTH
185                        ll = PATH_len(ran,Sig_total)
186                        //Determine new scattering direction vector.
[430]187                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function
188
[429]189                        //X,Y,Z-POSITION OF SCATTERING EVENT.
190                        xx += ll*vx
191                        yy += ll*vy
192                        zz += ll*vz
193                        RR = sqrt(xx*xx+yy*yy)          //radial position of scattering event.
194
195                        //Check whether interaction occurred within sample volume.
[430]196                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
[429]197                                //NEUTRON INTERACTED.
198                                INDEX += 1                      //Increment counter of scattering events.
199                                IF(INDEX == 1)
200                                        N2 += 1                 //Increment # of scat. neutrons
201                                Endif
202                                ran = abs(enoise(1))            //[0,1]
203                                //Split neutron interactions into scattering and absorption events
204                                IF(ran > ratio )                //C             NEUTRON SCATTERED coherently
205                                        FIND_THETA = 0                  //false
206                                        DO
[430]207                                                //ran = abs(enoise(1))          //[0,1]
[429]208                                                //theta = Scat_angle(Ran,R_DAB,wavelength)      // CHOOSE DAB ANGLE -- this is 2Theta
209                                                //Q0 = 2*PI*THETA/WAVElength                                    // John chose theta, calculated Q
210
211                                                // pick a q-value from the deviate function
212                                                // pnt2x truncates the point to an integer before returning the x
213                                                // so get it from the wave scaling instead
214                                                Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta
215                                                theta = Q0/2/Pi*wavelength              //SAS approximation
216                                               
[430]217                                                //Print "q0, theta = ",q0,theta
218                                               
219                                                FIND_THETA = 1          //always accept
220
[429]221                                        while(!find_theta)
222                                        ran = abs(enoise(1))            //[0,1]
223                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
224                                ELSE
225                                        //NEUTRON scattered incoherently
226                   // N3 += 1
227                  // DONE = 1
228                  // phi and theta are random over the entire sphere of scattering
229                       
230                        ran = abs(enoise(1))            //[0,1]
231                                        theta = pi*ran
232                       
233                        ran = abs(enoise(1))            //[0,1]
234                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
235                                ENDIF           //(ran > ratio)
236                        ELSE
237                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
238                                DONE = 1                //done = true, will exit from loop
239                                //Increment #scattering events array
240                                If (index <= N_Index)
241                                        NN[INDEX] += 1
242                                Endif
243                                //IF (VZ > 1.0)         // FIX INVALID ARGUMENT
244                                        //VZ = 1.0 - 1.2e-7
245                                //ENDIF
246                                Theta_z = acos(Vz)              // Angle WITH respect to z axis.
247                                testQ = 2*pi*sin(theta_z)/wavelength
248                               
249                                // pick a random phi angle, and see if it lands on the detector
250                                // since the scattering is isotropic, I can safely pick a new, random value
251                                // this would not be true if simulating anisotropic scattering.
[430]252                                testPhi = abs(enoise(1))*2*Pi
[429]253                                // is it on the detector?       
254                                FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
[430]255                               
256//                              if(xPixel != xCtr && yPixel != yCtr)
257//                                      Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel
258//                              endif
259                               
[429]260                                if(xPixel != -1 && yPixel != -1)
261                                        isOn += 1
262                                        //if(index==1)  // only the single scattering events
263                                                MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering
264                                        //endif
265                                endif
266
267                                If(theta_z < theta_max)
268                                        //Choose index for scattering angle array.
269                                        //IND = NINT(THETA_z/DTH + 0.4999999)
270                                        ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint()
271                                        NT[ind] += 1                    //Increment bin for angle.
272                                        //Increment angle array for single scattering events.
273                                        IF(INDEX == 1)
274                                                j1[ind] += 1
275                                        Endif
276                                        //Increment angle array for double scattering events.
277                                        IF (INDEX == 2)
278                                                j2[ind] += 1
279                                        Endif
280                                EndIf
281                               
282                        ENDIF
283                while (!done)
284        while(n1 < imon)
285
286//      Print "Monte Carlo Done"
287        trans_th = exp(-sig_total*thick)
288//      TRANS_exp = (N1-N2) / N1                        // Transmission
289        // dsigma/domega assuming isotropic scattering, with no absorption.
290//      Print "trans_exp = ",trans_exp
291//      Print "total # of neutrons reaching 2D detector",isOn
292//      Print "fraction of incident neutrons reaching detector ",isOn/iMon
293        results[2] = isOn
294        results[3] = isOn/iMon
295       
296                               
297// OUTPUT of the 1D data, not necessary now since I want the 2D
298//      Make/O/N=(num_bins) qvals,int_ms,sig_ms,int_sing,int_doub
299//      ii=1
300//      Print "binning"
301//      do
302//              //CALCULATE MEAN THETA IN BIN.
303//              THETA_z = (ii-0.5)*DTH                  // Mean scattering angle of bin.
304//              //Solid angle of Ith bin.
305//              D_OMEGA = 2*PI*ABS( COS(ii*DTH) - COS((ii-1)*DTH) )
306//              //SOLID ANGLE CORRECTION: YIELDING CROSS-SECTION.
307//              dS_dW = NT[ii]/(N1*THICK*Trans_th*D_OMEGA)
308//              SIG = NT[ii]^0.5/(N1*THICK*Trans_th*D_OMEGA)
309//              ds_dw_single = J1[ii]/(N1*THICK*Trans_th*D_OMEGA)
310//              ds_dw_double = J2[ii]/(N1*THICK*Trans_th*D_OMEGA)
311//              //Deviation from isotropic model.
312//              qq = 4*pi*sin(0.5*theta_z)/wavelength
313//              qvals[ii-1] = qq
314//              int_ms[ii-1] = dS_dW
315//              sig_ms[ii-1] = sig
316//              int_sing[ii-1] = ds_dw_single
317//              int_doub[ii-1] = ds_dw_double
318//              //140     WRITE (7,145) qq,dS_dW,SIG,ds_dw_single,ds_dw_double
319//              ii+=1
320//      while(ii<=num_bins)
321//      Print "done binning"
322
323
324//        Write(7,*)
325//        Write(7,*) '#Times Sc.   #Events    '
326//      DO 150 I = 1,N_INDEX
327//150    WRITE (7,146) I,NN(I)
328//146   Format (I5,T10,F8.0)
329
330///       Write(7,171) N1
331//        Write(7,172) N2
332//        Write(7,173) N3
333//        Write(7,174) N2-N3
334
335//171     Format('Total number of neutrons:         N1= ',E10.5)
336///172     Format('Number of neutrons that interact: N2= ',E10.5)
337//173     Format('Number of absorption events:      N3= ',E10.5)
338//174     Format('# of neutrons that scatter out:(N2-N3)= ',E10.5)
339
340//      Print "Total number of neutrons = ",N1
341//      Print "Total number of neutrons that interact = ",N2
342//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
343        results[4] = N2
344        results[5] = sum(j1,-inf,inf)/N2
345       
346        Tabs = (N1-N3)/N1
347        Ttot = (N1-N2)/N1
348        I1_sumI = NN[0]/(N2-N3)
349//      Print "Tabs = ",Tabs
350//      Print "Transmitted neutrons = ",Ttot
351        results[6] = Ttot
352//      Print "I1 / all I1 = ", I1_sumI
353//      Print "DONE!"
354
355End
356////////        end of main function for calculating multiple scattering
357
358
359// returns the random deviate as a wave
360// and the total SAS cross-section [1/cm] sig_sas
361Function        CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
362        FUNCREF SANSModelAAO_MCproto func
363        WAVE coef
364        Variable lam
365        String outWave
366        Variable &SASxs
367
368        Variable nPts_ran=10000,qu
369        qu = 4*pi/lam           
370       
371        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"
372        WAVE Gq = root:Packages:NIST:SAS:gQ
373        WAVE xw = root:Packages:NIST:SAS:xw
[430]374        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
[429]375        xw=x                                                                                            //for the AAO
376        func(coef,Gq,xw)                                                                        //call as AAO
377
378//      Gq = x*Gq                                                                                                       // SAS approximation
379        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact
380        //
381        Integrate/METH=1 Gq/D=Gq_INT
382       
383        SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1]
384       
385        Gq_INT /= Gq_INT[nPts_ran-1]
386       
387        Duplicate/O Gq_INT $outWave
388
389        return(0)
390End
391
392Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
393        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
394
395        Variable theta,dy,dx,qx,qy
396        //decompose to qx,qy
397        qx = testQ*cos(testPhi)
398        qy = testQ*sin(testPhi)
399
400        //convert qx,qy to pixel locations relative to # of pixels x, y from center
401        theta = 2*asin(qy*lam/4/pi)
402        dy = sdd*tan(theta)
403        yPixel = round(yCtr + dy/pixSize)
404       
405        theta = 2*asin(qx*lam/4/pi)
406        dx = sdd*tan(theta)
407        xPixel = round(xCtr + dx/pixSize)
408
409        //if on detector, return xPix and yPix values, otherwise -1
410        if(yPixel > 127 || yPixel < 0)
411                yPixel = -1
412        endif
413        if(xPixel > 127 || xPixel < 0)
414                xPixel = -1
415        endif
416       
417        return(0)
418End
419
420Function MC_CheckFunctionAndCoef(funcStr,coefStr)
421        String funcStr,coefStr
422       
423        SVAR/Z listStr=root:Packages:NIST:coefKWStr
424        if(SVAR_Exists(listStr) == 1)
425                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
426                if(cmpstr("",properCoefStr)==0)
427                        return(0)               //false, no match found, so properCoefStr is returned null
428                endif
429                if(cmpstr(coefStr,properCoefStr)==0)
430                        return(1)               //true, the coef is the correct match
431                endif
432        endif
433        return(0)                       //false, wrong coef
434End
435
436Function/S MC_getFunctionCoef(funcStr)
437        String funcStr
438
439        SVAR/Z listStr=root:Packages:NIST:coefKWStr
440        String coefStr=""
441        if(SVAR_Exists(listStr) == 1)
442                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
443        endif
444        return(coefStr)
445End
446
447Function SANSModelAAO_MCproto(w,yw,xw)
448        Wave w,yw,xw
449       
450        Print "in SANSModelAAO_MCproto function"
451        return(1)
452end
453
454Function/S MC_FunctionPopupList()
455        String list,tmp
456        list = FunctionList("*",";","KIND:10")          //get every user defined curve fit function
457
458        //now start to remove everything the user doesn't need to see...
459               
460        tmp = FunctionList("*_proto",";","KIND:10")             //prototypes
461        list = RemoveFromList(tmp, list  ,";")
462        //prototypes that show up if GF is loaded
463        list = RemoveFromList("GFFitFuncTemplate", list)
464        list = RemoveFromList("GFFitAllAtOnceTemplate", list)
465        list = RemoveFromList("NewGlblFitFunc", list)
466        list = RemoveFromList("NewGlblFitFuncAllAtOnce", list)
467        list = RemoveFromList("GlobalFitFunc", list)
468        list = RemoveFromList("GlobalFitAllAtOnce", list)
469        list = RemoveFromList("GFFitAAOStructTemplate", list)
470        list = RemoveFromList("NewGF_SetXWaveInList", list)
471        list = RemoveFromList("NewGlblFitFuncAAOStruct", list)
472       
473        // more to remove as a result of 2D/Gizmo
474        list = RemoveFromList("A_WMRunLessThanDelta", list)
475        list = RemoveFromList("WMFindNaNValue", list)
476        list = RemoveFromList("WM_Make3DBarChartParametricWave", list)
477        list = RemoveFromList("UpdateQxQy2Mat", list)
478        list = RemoveFromList("MakeBSMask", list)
479       
480        // MOTOFIT/GenFit bits
481        tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;"
482        list = RemoveFromList(tmp, list  ,";")
483
484        // SANS Reduction bits
485        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;"
486        list = RemoveFromList(tmp, list  ,";")
[430]487        list = RemoveFromList("Monte_SANS", list)
[429]488
489        tmp = FunctionList("f*",";","NPARAMS:2")                //point calculations
490        list = RemoveFromList(tmp, list  ,";")
491       
492        tmp = FunctionList("fSmear*",";","NPARAMS:3")           //smeared dependency calculations
493        list = RemoveFromList(tmp, list  ,";")
494       
495//      tmp = FunctionList("*X",";","KIND:4")           //XOPs, but these shouldn't show up if KIND:10 is used initially
496//      Print "X* = ",tmp
497//      print " "
498//      list = RemoveFromList(tmp, list  ,";")
499       
500        //non-fit functions that I can't seem to filter out
501        list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";")
502
503        if(strlen(list)==0)
504                list = "No functions plotted"
505        endif
506       
507        list = SortList(list)
508        return(list)
[430]509End             
[429]510
511
512//Function Scat_Angle(Ran,R_DAB,wavelength)
513//      Variable Ran,r_dab,wavelength
514//
515//      Variable qq,arg,theta
516//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
517//      arg = qq*wavelength/(4*pi)
518//      If (arg < 1.0)
519//              theta = 2.*asin(arg)
520//      else
521//              theta = pi
522//      endif
523//      Return (theta)
524//End
525
526//calculates new direction (xyz) from an old direction
527//theta and phi don't change
[430]528Function NewDirection(vx,vy,vz,theta,phi)
529        Variable &vx,&vy,&vz
[429]530        Variable theta,phi
531       
532        Variable err=0,vx0,vy0,vz0
[430]533        Variable nx,ny,mag_xy,tx,ty,tz
[429]534       
535        //store old direction vector
536        vx0 = vx
537        vy0 = vy
538        vz0 = vz
539       
540        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
541        if(mag_xy < 1e-12)
542                //old vector lies along beam direction
543                nx = 0
544                ny = 1
545                tx = 1
546                ty = 0
547                tz = 0
548        else
549                Nx = -Vy0 / Mag_XY
550                Ny = Vx0 / Mag_XY
551                Tx = -Vz0*Vx0 / Mag_XY
552                Ty = -Vz0*Vy0 / Mag_XY
553                Tz = Mag_XY
554        endif
555       
556        //new scattered direction vector
557        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
558        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
559        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
560       
561        Return(err)
562End
563
564Function path_len(aval,sig_tot)
565        Variable aval,sig_tot
566       
567        Variable retval
568       
[430]569        retval = -1*ln(1-aval)/sig_tot
[429]570       
571        return(retval)
572End
573
574Window MC_SASCALC() : Panel
575        PauseUpdate; Silent 1           // building window...
[430]576        NewPanel /W=(787,44,1088,563)  /N=MC_SASCALC as "SANS Simulator"
[429]577        CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation"
578        CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo
579        SetVariable MC_setvar0,pos={11,38},size={144,15},bodyWidth=80,title="# of neutrons"
580        SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
581        SetVariable MC_setvar0_1,pos={11,121},size={131,15},bodyWidth=60,title="Thickness (cm)"
582        SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
583        SetVariable MC_setvar0_2,pos={11,93},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
584        SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
585        SetVariable MC_setvar0_3,pos={11,149},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
586        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
587        PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function"
588        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
[430]589        Button MC_button0,pos={17,455},size={50,20},proc=MC_DoItButtonProc,title="Do It"
590        Button MC_button1,pos={15,484},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
[429]591       
592        SetDataFolder root:Packages:NIST:SAS:
593        Edit/W=(13,174,284,435)/HOST=# results_desc,results
594        ModifyTable width(Point)=0,width(results_desc)=150
595        SetDataFolder root:
596        RenameWindow #,T_results
597        SetActiveSubwindow ##
598EndMacro
599
600Function MC_ModelPopMenuProc(pa) : PopupMenuControl
601        STRUCT WMPopupAction &pa
602
603        switch( pa.eventCode )
604                case 2: // mouse up
605                        Variable popNum = pa.popNum
606                        String popStr = pa.popStr
607                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
608                        gStr = popStr
609                       
610                        break
611        endswitch
612
613        return 0
614End
[430]615
616Function MC_DoItButtonProc(ba) : ButtonControl
617        STRUCT WMButtonAction &ba
618
619        switch( ba.eventCode )
620                case 2: // mouse up
621                        // click code here
622                        ReCalculateInten(1)
623                        break
624        endswitch
625
626        return 0
627End
628
629
630Function MC_Display2DButtonProc(ba) : ButtonControl
631        STRUCT WMButtonAction &ba
632
633        switch( ba.eventCode )
634                case 2: // mouse up
635                        // click code here
636                        Execute "ChangeDisplay(\"SAS\")"
637                        break
638        endswitch
639
640        return 0
641End
642
643/////UNUSED, testing routines that have note been updated to work with SASCALC
644//
645//Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr)
646//      Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1
647//      String funcStr
648//      Prompt funcStr, "Pick the model function", popup,       MC_FunctionPopupList()
649//     
650//      String coefStr = MC_getFunctionCoef(funcStr)
651//      Variable pixSize = 0.5          // can't have 11 parameters in macro!
652//     
653//      if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
654//              Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
655//      endif
656//     
657//      Make/O/D/N=10 root:Packages:NIST:SAS:results
658//      Make/O/T/N=10 root:Packages:NIST:SAS:results_desc
659//     
660//      RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results)
661//     
662//End
663//
664//Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr)
665//      Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh
666//      String funcStr,coefStr
667//      WAVE results
668//     
669//      FUNCREF SANSModelAAO_MCproto func=$funcStr
670//     
671//      Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results)
672//End
673//
674////// END UNUSED BLOCK
Note: See TracBrowser for help on using the repository browser.