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

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

Lots of little changes:

changed debye function to a limiting value at low qRg to avoid numerical errors. There is no XOP, so no changes needed there.

Added MonteCarlo? simulation routine (1st draft) to the SANS Reduction, including changes to SASCALC to implement. Still very alpha with lots of bits to add. See the top of the monte carlo file for the list.

added MonteCarlo? file to SANS includes
added SAS folder to 2D display list to view MonteCarlo? results
removed more "non-functions" from the function popup in analysis. these were a result of motofit/genFit, and simultaneous loading of SANS/USANS reduction. More will follow.

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