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

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

A variety of fixes to make some procedures compatible with the new Packages:NIST subfolder structure

Also, some additional changes to include the proper procedures when loading overall packages

File size: 22.7 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??
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// - 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//
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
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
76
77        Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas
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)
87
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
92        Variable Vx,Vy,Vz,Theta_z,qq
93        Variable Sig_scat,Sig_abs,Ratio,Sig_total
94        Variable isOn=0,testQ,testPhi,xPixel,yPixel
95       
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       
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
112        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value)
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
139//!!!! the ratio is not yet properly weighted for the volume fractions of each component!!!
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.
187                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function
188
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.
196                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2))
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
207                                                //ran = abs(enoise(1))          //[0,1]
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                                               
217                                                //Print "q0, theta = ",q0,theta
218                                               
219                                                FIND_THETA = 1          //always accept
220
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                  // !can't just choose random theta and phi, won't be random over sphere solid angle
230                       
231                        ran = abs(enoise(1))            //[0,1]
232                                        theta = acos(2*ran-1)           
233                       
234                        ran = abs(enoise(1))            //[0,1]
235                                        PHI = 2.0*PI*Ran                        //Chooses azimuthal scattering angle.
236                                ENDIF           //(ran > ratio)
237                        ELSE
238                                //NEUTRON ESCAPES FROM SAMPLE -- bin it somewhere
239                                DONE = 1                //done = true, will exit from loop
240                                //Increment #scattering events array
241                                If (index <= N_Index)
242                                        NN[INDEX] += 1
243                                Endif
244                                //IF (VZ > 1.0)         // FIX INVALID ARGUMENT
245                                        //VZ = 1.0 - 1.2e-7
246                                //ENDIF
247                                Theta_z = acos(Vz)              // Angle WITH respect to z axis.
248                                testQ = 2*pi*sin(theta_z)/wavelength
249                               
250                                // pick a random phi angle, and see if it lands on the detector
251                                // since the scattering is isotropic, I can safely pick a new, random value
252                                // this would not be true if simulating anisotropic scattering.
253                                testPhi = abs(enoise(1))*2*Pi
254                                // is it on the detector?       
255                                FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
256                               
257//                              if(xPixel != xCtr && yPixel != yCtr)
258//                                      Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel
259//                              endif
260                               
261                                if(xPixel != -1 && yPixel != -1)
262                                        isOn += 1
263                                        //if(index==1)  // only the single scattering events
264                                                MC_linear_data[xPixel][yPixel] += 1             //this is the total scattering, including multiple scattering
265                                        //endif
266                                endif
267
268                                If(theta_z < theta_max)
269                                        //Choose index for scattering angle array.
270                                        //IND = NINT(THETA_z/DTH + 0.4999999)
271                                        ind = round(THETA_z/DTH + 0.4999999)            //round is eqivalent to nint()
272                                        NT[ind] += 1                    //Increment bin for angle.
273                                        //Increment angle array for single scattering events.
274                                        IF(INDEX == 1)
275                                                j1[ind] += 1
276                                        Endif
277                                        //Increment angle array for double scattering events.
278                                        IF (INDEX == 2)
279                                                j2[ind] += 1
280                                        Endif
281                                EndIf
282                               
283                        ENDIF
284                while (!done)
285        while(n1 < imon)
286
287//      Print "Monte Carlo Done"
288        trans_th = exp(-sig_total*thick)
289//      TRANS_exp = (N1-N2) / N1                        // Transmission
290        // dsigma/domega assuming isotropic scattering, with no absorption.
291//      Print "trans_exp = ",trans_exp
292//      Print "total # of neutrons reaching 2D detector",isOn
293//      Print "fraction of incident neutrons reaching detector ",isOn/iMon
294        results[2] = isOn
295        results[3] = isOn/iMon
296       
297                               
298// OUTPUT of the 1D data, not necessary now since I want the 2D
299//      Make/O/N=(num_bins) qvals,int_ms,sig_ms,int_sing,int_doub
300//      ii=1
301//      Print "binning"
302//      do
303//              //CALCULATE MEAN THETA IN BIN.
304//              THETA_z = (ii-0.5)*DTH                  // Mean scattering angle of bin.
305//              //Solid angle of Ith bin.
306//              D_OMEGA = 2*PI*ABS( COS(ii*DTH) - COS((ii-1)*DTH) )
307//              //SOLID ANGLE CORRECTION: YIELDING CROSS-SECTION.
308//              dS_dW = NT[ii]/(N1*THICK*Trans_th*D_OMEGA)
309//              SIG = NT[ii]^0.5/(N1*THICK*Trans_th*D_OMEGA)
310//              ds_dw_single = J1[ii]/(N1*THICK*Trans_th*D_OMEGA)
311//              ds_dw_double = J2[ii]/(N1*THICK*Trans_th*D_OMEGA)
312//              //Deviation from isotropic model.
313//              qq = 4*pi*sin(0.5*theta_z)/wavelength
314//              qvals[ii-1] = qq
315//              int_ms[ii-1] = dS_dW
316//              sig_ms[ii-1] = sig
317//              int_sing[ii-1] = ds_dw_single
318//              int_doub[ii-1] = ds_dw_double
319//              //140     WRITE (7,145) qq,dS_dW,SIG,ds_dw_single,ds_dw_double
320//              ii+=1
321//      while(ii<=num_bins)
322//      Print "done binning"
323
324
325//        Write(7,*)
326//        Write(7,*) '#Times Sc.   #Events    '
327//      DO 150 I = 1,N_INDEX
328//150    WRITE (7,146) I,NN(I)
329//146   Format (I5,T10,F8.0)
330
331///       Write(7,171) N1
332//        Write(7,172) N2
333//        Write(7,173) N3
334//        Write(7,174) N2-N3
335
336//171     Format('Total number of neutrons:         N1= ',E10.5)
337///172     Format('Number of neutrons that interact: N2= ',E10.5)
338//173     Format('Number of absorption events:      N3= ',E10.5)
339//174     Format('# of neutrons that scatter out:(N2-N3)= ',E10.5)
340
341//      Print "Total number of neutrons = ",N1
342//      Print "Total number of neutrons that interact = ",N2
343//      Print "Fraction of singly scattered neutrons = ",sum(j1,-inf,inf)/N2
344        results[4] = N2
345        results[5] = sum(j1,-inf,inf)/N2
346       
347        Tabs = (N1-N3)/N1
348        Ttot = (N1-N2)/N1
349        I1_sumI = NN[0]/(N2-N3)
350//      Print "Tabs = ",Tabs
351//      Print "Transmitted neutrons = ",Ttot
352        results[6] = Ttot
353//      Print "I1 / all I1 = ", I1_sumI
354//      Print "DONE!"
355
356End
357////////        end of main function for calculating multiple scattering
358
359
360// returns the random deviate as a wave
361// and the total SAS cross-section [1/cm] sig_sas
362Function        CalculateRandomDeviate(func,coef,lam,outWave,SASxs)
363        FUNCREF SANSModelAAO_MCproto func
364        WAVE coef
365        Variable lam
366        String outWave
367        Variable &SASxs
368
369        Variable nPts_ran=10000,qu
370        qu = 4*pi/lam           
371       
372        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"
373        WAVE Gq = root:Packages:NIST:SAS:gQ
374        WAVE xw = root:Packages:NIST:SAS:xw
375        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
376        xw=x                                                                                            //for the AAO
377        func(coef,Gq,xw)                                                                        //call as AAO
378
379//      Gq = x*Gq                                                                                                       // SAS approximation
380        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact
381        //
382        Integrate/METH=1 Gq/D=Gq_INT
383       
384        SASxs = lam*lam/2/pi*Gq_INT[nPts_ran-1]
385       
386        Gq_INT /= Gq_INT[nPts_ran-1]
387       
388        Duplicate/O Gq_INT $outWave
389
390        return(0)
391End
392
393Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel)
394        Variable testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,&xPixel,&yPixel
395
396        Variable theta,dy,dx,qx,qy
397        //decompose to qx,qy
398        qx = testQ*cos(testPhi)
399        qy = testQ*sin(testPhi)
400
401        //convert qx,qy to pixel locations relative to # of pixels x, y from center
402        theta = 2*asin(qy*lam/4/pi)
403        dy = sdd*tan(theta)
404        yPixel = round(yCtr + dy/pixSize)
405       
406        theta = 2*asin(qx*lam/4/pi)
407        dx = sdd*tan(theta)
408        xPixel = round(xCtr + dx/pixSize)
409
410        //if on detector, return xPix and yPix values, otherwise -1
411        if(yPixel > 127 || yPixel < 0)
412                yPixel = -1
413        endif
414        if(xPixel > 127 || xPixel < 0)
415                xPixel = -1
416        endif
417       
418        return(0)
419End
420
421Function MC_CheckFunctionAndCoef(funcStr,coefStr)
422        String funcStr,coefStr
423       
424        SVAR/Z listStr=root:Packages:NIST:coefKWStr
425        if(SVAR_Exists(listStr) == 1)
426                String properCoefStr = StringByKey(funcStr, listStr  ,"=",";",0)
427                if(cmpstr("",properCoefStr)==0)
428                        return(0)               //false, no match found, so properCoefStr is returned null
429                endif
430                if(cmpstr(coefStr,properCoefStr)==0)
431                        return(1)               //true, the coef is the correct match
432                endif
433        endif
434        return(0)                       //false, wrong coef
435End
436
437Function/S MC_getFunctionCoef(funcStr)
438        String funcStr
439
440        SVAR/Z listStr=root:Packages:NIST:coefKWStr
441        String coefStr=""
442        if(SVAR_Exists(listStr) == 1)
443                coefStr = StringByKey(funcStr, listStr  ,"=",";",0)
444        endif
445        return(coefStr)
446End
447
448Function SANSModelAAO_MCproto(w,yw,xw)
449        Wave w,yw,xw
450       
451        Print "in SANSModelAAO_MCproto function"
452        return(1)
453end
454
455Function/S MC_FunctionPopupList()
456        String list,tmp
457        list = FunctionList("*",";","KIND:10")          //get every user defined curve fit function
458
459        //now start to remove everything the user doesn't need to see...
460               
461        tmp = FunctionList("*_proto",";","KIND:10")             //prototypes
462        list = RemoveFromList(tmp, list  ,";")
463        //prototypes that show up if GF is loaded
464        list = RemoveFromList("GFFitFuncTemplate", list)
465        list = RemoveFromList("GFFitAllAtOnceTemplate", list)
466        list = RemoveFromList("NewGlblFitFunc", list)
467        list = RemoveFromList("NewGlblFitFuncAllAtOnce", list)
468        list = RemoveFromList("GlobalFitFunc", list)
469        list = RemoveFromList("GlobalFitAllAtOnce", list)
470        list = RemoveFromList("GFFitAAOStructTemplate", list)
471        list = RemoveFromList("NewGF_SetXWaveInList", list)
472        list = RemoveFromList("NewGlblFitFuncAAOStruct", list)
473       
474        // more to remove as a result of 2D/Gizmo
475        list = RemoveFromList("A_WMRunLessThanDelta", list)
476        list = RemoveFromList("WMFindNaNValue", list)
477        list = RemoveFromList("WM_Make3DBarChartParametricWave", list)
478        list = RemoveFromList("UpdateQxQy2Mat", list)
479        list = RemoveFromList("MakeBSMask", list)
480       
481        // MOTOFIT/GenFit bits
482        tmp = "GEN_allatoncefitfunc;GEN_fitfunc;GetCheckBoxesState;MOTO_GFFitAllAtOnceTemplate;MOTO_GFFitFuncTemplate;MOTO_NewGF_SetXWaveInList;MOTO_NewGlblFitFunc;MOTO_NewGlblFitFuncAllAtOnce;"
483        list = RemoveFromList(tmp, list  ,";")
484
485        // SANS Reduction bits
486        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;"
487        list = RemoveFromList(tmp, list  ,";")
488        list = RemoveFromList("Monte_SANS", list)
489
490        tmp = FunctionList("f*",";","NPARAMS:2")                //point calculations
491        list = RemoveFromList(tmp, list  ,";")
492       
493        tmp = FunctionList("fSmear*",";","NPARAMS:3")           //smeared dependency calculations
494        list = RemoveFromList(tmp, list  ,";")
495       
496//      tmp = FunctionList("*X",";","KIND:4")           //XOPs, but these shouldn't show up if KIND:10 is used initially
497//      Print "X* = ",tmp
498//      print " "
499//      list = RemoveFromList(tmp, list  ,";")
500       
501        //non-fit functions that I can't seem to filter out
502        list = RemoveFromList("BinaryHS_PSF11;BinaryHS_PSF12;BinaryHS_PSF22;EllipCyl_Integrand;PP_Inner;PP_Outer;Phi_EC;TaE_Inner;TaE_Outer;",list,";")
503
504        if(strlen(list)==0)
505                list = "No functions plotted"
506        endif
507       
508        list = SortList(list)
509        return(list)
510End             
511
512
513//Function Scat_Angle(Ran,R_DAB,wavelength)
514//      Variable Ran,r_dab,wavelength
515//
516//      Variable qq,arg,theta
517//      qq = 1. / ( R_DAB*sqrt(1.0/Ran - 1.0) )
518//      arg = qq*wavelength/(4*pi)
519//      If (arg < 1.0)
520//              theta = 2.*asin(arg)
521//      else
522//              theta = pi
523//      endif
524//      Return (theta)
525//End
526
527//calculates new direction (xyz) from an old direction
528//theta and phi don't change
529Function NewDirection(vx,vy,vz,theta,phi)
530        Variable &vx,&vy,&vz
531        Variable theta,phi
532       
533        Variable err=0,vx0,vy0,vz0
534        Variable nx,ny,mag_xy,tx,ty,tz
535       
536        //store old direction vector
537        vx0 = vx
538        vy0 = vy
539        vz0 = vz
540       
541        mag_xy = sqrt(vx0*vx0 + vy0*vy0)
542        if(mag_xy < 1e-12)
543                //old vector lies along beam direction
544                nx = 0
545                ny = 1
546                tx = 1
547                ty = 0
548                tz = 0
549        else
550                Nx = -Vy0 / Mag_XY
551                Ny = Vx0 / Mag_XY
552                Tx = -Vz0*Vx0 / Mag_XY
553                Ty = -Vz0*Vy0 / Mag_XY
554                Tz = Mag_XY
555        endif
556       
557        //new scattered direction vector
558        Vx = cos(phi)*sin(theta)*Tx + sin(phi)*sin(theta)*Nx + cos(theta)*Vx0
559        Vy = cos(phi)*sin(theta)*Ty + sin(phi)*sin(theta)*Ny + cos(theta)*Vy0
560        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0
561       
562        Return(err)
563End
564
565Function path_len(aval,sig_tot)
566        Variable aval,sig_tot
567       
568        Variable retval
569       
570        retval = -1*ln(1-aval)/sig_tot
571       
572        return(retval)
573End
574
575Window MC_SASCALC() : Panel
576        PauseUpdate; Silent 1           // building window...
577        NewPanel /W=(787,44,1088,563)  /N=MC_SASCALC as "SANS Simulator"
578        CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation"
579        CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo
580        SetVariable MC_setvar0,pos={11,38},size={144,15},bodyWidth=80,title="# of neutrons"
581        SetVariable MC_setvar0,limits={-inf,inf,100},value= root:Packages:NIST:SAS:gImon
582        SetVariable MC_setvar0_1,pos={11,121},size={131,15},bodyWidth=60,title="Thickness (cm)"
583        SetVariable MC_setvar0_1,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gThick
584        SetVariable MC_setvar0_2,pos={11,93},size={149,15},bodyWidth=60,title="Incoherent XS (cm)"
585        SetVariable MC_setvar0_2,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh
586        SetVariable MC_setvar0_3,pos={11,149},size={150,15},bodyWidth=60,title="Sample Radius (cm)"
587        SetVariable MC_setvar0_3,limits={-inf,inf,0.1},value= root:Packages:NIST:SAS:gR2
588        PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function"
589        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()"
590        Button MC_button0,pos={17,455},size={50,20},proc=MC_DoItButtonProc,title="Do It"
591        Button MC_button1,pos={15,484},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D"
592       
593        SetDataFolder root:Packages:NIST:SAS:
594        Edit/W=(13,174,284,435)/HOST=# results_desc,results
595        ModifyTable width(Point)=0,width(results_desc)=150
596        SetDataFolder root:
597        RenameWindow #,T_results
598        SetActiveSubwindow ##
599EndMacro
600
601Function MC_ModelPopMenuProc(pa) : PopupMenuControl
602        STRUCT WMPopupAction &pa
603
604        switch( pa.eventCode )
605                case 2: // mouse up
606                        Variable popNum = pa.popNum
607                        String popStr = pa.popStr
608                        SVAR gStr = root:Packages:NIST:SAS:gFuncStr
609                        gStr = popStr
610                       
611                        break
612        endswitch
613
614        return 0
615End
616
617Function MC_DoItButtonProc(ba) : ButtonControl
618        STRUCT WMButtonAction &ba
619
620        switch( ba.eventCode )
621                case 2: // mouse up
622                        // click code here
623                        ReCalculateInten(1)
624                        break
625        endswitch
626
627        return 0
628End
629
630
631Function MC_Display2DButtonProc(ba) : ButtonControl
632        STRUCT WMButtonAction &ba
633
634        switch( ba.eventCode )
635                case 2: // mouse up
636                        // click code here
637                        Execute "ChangeDisplay(\"SAS\")"
638                        break
639        endswitch
640
641        return 0
642End
643
644/////UNUSED, testing routines that have note been updated to work with SASCALC
645//
646//Macro Simulate2D_MonteCarlo(imon,r1,r2,xCtr,yCtr,sdd,thick,wavelength,sig_incoh,funcStr)
647//      Variable imon=100000,r1=0.6,r2=0.8,xCtr=100,yCtr=64,sdd=400,thick=0.1,wavelength=6,sig_incoh=0.1
648//      String funcStr
649//      Prompt funcStr, "Pick the model function", popup,       MC_FunctionPopupList()
650//     
651//      String coefStr = MC_getFunctionCoef(funcStr)
652//      Variable pixSize = 0.5          // can't have 11 parameters in macro!
653//     
654//      if(!MC_CheckFunctionAndCoef(funcStr,coefStr))
655//              Abort "The coefficients and function type do not match. Please correct the selections in the popup menus."
656//      endif
657//     
658//      Make/O/D/N=10 root:Packages:NIST:SAS:results
659//      Make/O/T/N=10 root:Packages:NIST:SAS:results_desc
660//     
661//      RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr,results)
662//     
663//End
664//
665//Function RunMonte(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,funcStr,coefStr)
666//      Variable imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh
667//      String funcStr,coefStr
668//      WAVE results
669//     
670//      FUNCREF SANSModelAAO_MCproto func=$funcStr
671//     
672//      Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,sig_sas,ran_dev,linear_data,results)
673//End
674//
675////// END UNUSED BLOCK
Note: See TracBrowser for help on using the repository browser.