Changeset 430


Ignore:
Timestamp:
Oct 17, 2008 3:12:18 PM (14 years ago)
Author:
srkline
Message:

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.

Location:
sans/Dev/trunk
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf

    r429 r430  
    219219        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;" 
    220220        list = RemoveFromList(tmp, list  ,";") 
    221          
     221        list = RemoveFromList("Monte_SANS", list) 
     222 
     223        // USANS Reduction bits 
     224        tmp = "DSM_Guinier_Fit;RemoveMaskedPoints;" 
     225        list = RemoveFromList(tmp, list  ,";") 
     226 
     227 
    222228        tmp = FunctionList("f*",";","NPARAMS:2")                //point calculations 
    223229        list = RemoveFromList(tmp, list  ,";") 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r429 r430  
    77// This code simulates the scattering for a selected model based on the instrument configuration 
    88// 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() 
    910// 
    1011// 
    1112// - 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// X how can I get the "data" on absolute scale? This would be a great comparison vs. the ideal model calculation 
    1314// - why does my integrated tau not match up with John's analytical calculations? where are the assumptions? 
    1415// - get rid of all small angle assumptions - to make sure that the calculation is correct at all angles 
    1516// - 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()) 
    1717 
    1818// - at the larger angles, is the "flat" detector being properly accounted for - in terms of 
    1919//   the solid angle and how many counts fall in that pixel. Am I implicitly defining a spherical detector 
    2020//   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. 
     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. 
    2324// - the background parameter for the model MUST be zero, or the integration for scattering 
    2425//    power will be incorrect. 
    2526// - 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 
     27// X if no MC desired, still use the selected model 
     28// X better display of MC results on panel 
    2829// - settings for "count for X seconds" or "how long to 1E6 cts on detector" (run short sim, then multiply) 
    2930// - add quartz window scattering to the simulation somehow 
    3031// - do smeared models make any sense?? 
    31 // 
    32  
    33 Macro 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          
    50 End 
    51  
    52 Function 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) 
    60 End 
     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 
    6153 
    6254////////// 
    … …  
    8072// 
    8173 
    82 Function 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  
     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 
    8778        Variable NUM_BINS,N_INDEX 
    8879        Variable RHO,SIGSAS,SIGABS_0 
    … …  
    9485        // at 0, making things much easier for me... 
    9586        //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          
     87 
    10388        Variable G0,E_NT,E_NN,TRANS_th,Trans_exp,rat 
    10489        Variable GG,GG_ED,dS_dW,ds_dw_double,ds_dw_single 
    10590        Variable DONE,FIND_THETA,err            //used as logicals 
    10691 
    107         Variable Vx,Vy,Vz,Theta_z,qq,SIG_SAS 
     92        Variable Vx,Vy,Vz,Theta_z,qq 
    10893        Variable Sig_scat,Sig_abs,Ratio,Sig_total 
    10994        Variable isOn=0,testQ,testPhi,xPixel,yPixel 
    11095         
    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  
     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         
    118108//      SetRandomSeed 0.1               //to get a reproduceable sequence 
    119109 
    120110//scattering power and maximum qvalue to bin 
    121111//      zpow = .1               //scattering power, calculated below 
    122         qmax = 0.1      //maximum Q to bin 1D data. (A-1) (not really used) 
     112        qmax = 4*pi/wavelength          //maximum Q to bin 1D data. (A-1) (not really used, so set to a big value) 
    123113        sigabs_0 = 0.0          // ignore absorption cross section/wavelength [1/(cm A)] 
    124114        N_INDEX = 50            // maximum number of scattering events per neutron 
    … …  
    128118// and calculate the scattering power from the model function 
    129119// 
    130         CalculateRandomDeviate(func,coef,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS) 
    131         WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev" 
    132120        Variable left = leftx(ran_dev) 
    133121        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 
    141122         
    142123//c       total SAS cross-section 
    … …  
    156137//      RATIO = SIG_ABS / SIG_TOTAL 
    157138        RATIO = sig_incoh / SIG_TOTAL 
    158 //!!!! the ratio is not yer porperly weighted for the volume fractions of each component!!! 
     139//!!!! the ratio is not yet properly weighted for the volume fractions of each component!!! 
    159140         
    160141//c       assuming theta = sin(theta)...OK 
    … …  
    174155        nt = 0 
    175156        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 
    180157         
    181158//C     MONITOR LOOP - looping over the number of incedent neutrons 
    … …  
    187164                Vy = 0.0 
    188165                Vz = 1.0 
    189                 d_vector[0] = 0 
    190                 d_vector[1] = 0 
    191                 d_vector[2] = 1 
    192166                 
    193167                Theta = 0.0             //      Initialize scattering angle. 
    … …  
    197171                INDEX = 0                       //      Set counter for number of scattering events. 
    198172                zz = 0.0                        //      Set entering dimension of sample. 
    199                 //RR = 2.0*R1           //      Make sure next loop runs at least once. 
    200173                 
    201174                do                                      //      Makes sure position is within circle. 
    … …  
    212185                        ll = PATH_len(ran,Sig_total) 
    213186                        //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] 
     187                        err = NewDirection(vx,vy,vz,Theta,Phi)          //vx,vy,vz is updated, theta, phi unchanged by function 
     188 
    218189                        //X,Y,Z-POSITION OF SCATTERING EVENT. 
    219190                        xx += ll*vx 
    … …  
    223194 
    224195                        //Check whether interaction occurred within sample volume. 
    225                         IF (((zz > 0.0) %& (zz < THICK)) %& (rr < r2)) 
     196                        IF (((zz > 0.0) && (zz < THICK)) && (rr < r2)) 
    226197                                //NEUTRON INTERACTED. 
    227198                                INDEX += 1                      //Increment counter of scattering events. 
    … …  
    234205                                        FIND_THETA = 0                  //false 
    235206                                        DO 
    236                                                 ran = abs(enoise(1))            //[0,1] 
    237                                                  
     207                                                //ran = abs(enoise(1))          //[0,1] 
    238208                                                //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 
    240209                                                //Q0 = 2*PI*THETA/WAVElength                                    // John chose theta, calculated Q 
    241210 
    242211                                                // pick a q-value from the deviate function 
    243212                                                // pnt2x truncates the point to an integer before returning the x 
    244                                                 //Q0 = pnt2x(ran_dev,binarysearchinterp(ran_dev,abs(enoise(1)))) 
    245213                                                // so get it from the wave scaling instead 
    246214                                                Q0 =left + binarysearchinterp(ran_dev,abs(enoise(1)))*delta 
    247215                                                theta = Q0/2/Pi*wavelength              //SAS approximation 
    248216                                                 
    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 //                                               
     217                                                //Print "q0, theta = ",q0,theta 
     218                                                 
     219                                                FIND_THETA = 1          //always accept 
     220 
    260221                                        while(!find_theta) 
    261222                                        ran = abs(enoise(1))            //[0,1] 
    … …  
    289250                                // since the scattering is isotropic, I can safely pick a new, random value 
    290251                                // this would not be true if simulating anisotropic scattering. 
    291                                 testPhi = abs(enoise(1,2))*2*Pi 
     252                                testPhi = abs(enoise(1))*2*Pi 
    292253                                // is it on the detector?        
    293254                                FindPixel(testQ,testPhi,wavelength,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
     255                                 
     256//                              if(xPixel != xCtr && yPixel != yCtr) 
     257//                                      Print "testQ,testPhi,xPixel,yPixel",testQ,testPhi,xPixel,yPixel 
     258//                              endif 
     259                                 
    294260                                if(xPixel != -1 && yPixel != -1) 
    295261                                        isOn += 1 
    … …  
    328294        results[3] = isOn/iMon 
    329295         
    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                                  
    337296                                 
    338297// OUTPUT of the 1D data, not necessary now since I want the 2D 
    … …  
    413372        WAVE Gq = root:Packages:NIST:SAS:gQ 
    414373        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 
     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 
    417375        xw=x                                                                                            //for the AAO 
    418376        func(coef,Gq,xw)                                                                        //call as AAO 
    … …  
    428386         
    429387        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] 
    442388 
    443389        return(0) 
    … …  
    539485        tmp = "ASStandardFunction;Ann_1D_Graph;Avg_1D_Graph;BStandardFunction;CStandardFunction;Draw_Plot1D;MyMat2XYZ;NewDirection;SANSModelAAO_MCproto;" 
    540486        list = RemoveFromList(tmp, list  ,";") 
     487        list = RemoveFromList("Monte_SANS", list) 
    541488 
    542489        tmp = FunctionList("f*",";","NPARAMS:2")                //point calculations 
    … …  
    560507        list = SortList(list) 
    561508        return(list) 
    562 End 
    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                     
     509End               
    589510 
    590511 
    … …  
    605526//calculates new direction (xyz) from an old direction 
    606527//theta and phi don't change 
    607 Function NewDirection(d_vector,theta,phi) 
    608         Wave d_vector 
     528Function NewDirection(vx,vy,vz,theta,phi) 
     529        Variable &vx,&vy,&vz 
    609530        Variable theta,phi 
    610531         
    611532        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] 
     533        Variable nx,ny,mag_xy,tx,ty,tz 
    616534         
    617535        //store old direction vector 
    … …  
    641559        Vz = cos(phi)*sin(theta)*Tz + cos(theta)*Vz0 
    642560         
    643         //reassign d_vector before exiting 
    644         d_vector[0] = vx 
    645         d_vector[1] = vy 
    646         d_vector[2] = vz 
    647          
    648561        Return(err) 
    649562End 
    … …  
    654567        Variable retval 
    655568         
    656         retval = -1*log(1-aval)/sig_tot 
     569        retval = -1*ln(1-aval)/sig_tot 
    657570         
    658571        return(retval) 
    … …  
    661574Window MC_SASCALC() : Panel 
    662575        PauseUpdate; Silent 1           // building window... 
    663         NewPanel /W=(787,44,1088,563) /K=1 /N=MC_SASCALC as "SANS Simulator" 
     576        NewPanel /W=(787,44,1088,563) /N=MC_SASCALC as "SANS Simulator" 
    664577        CheckBox MC_check0,pos={11,11},size={98,14},title="Use MC Simulation" 
    665578        CheckBox MC_check0,variable= root:Packages:NIST:SAS:gDoMonteCarlo 
    … …  
    674587        PopupMenu MC_popup0,pos={11,63},size={162,20},proc=MC_ModelPopMenuProc,title="Model Function" 
    675588        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 
     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" 
    676591         
    677592        SetDataFolder root:Packages:NIST:SAS: 
    … …  
    698613        return 0 
    699614End 
     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 

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/SASCALC.ipf

    r429 r430  
    647647        // do the simulation here 
    648648        Variable r1,xCtr,yCtr,sdd,pixSize,wavelength 
    649         String coefStr 
    650          
    651 //      Variable imon,thick,r2,sig_incoh 
    652 //      String funcStr 
    653 //      imon = 10000 
    654 //      thick = 0.1 
    655 //      sig_incoh = 0.1 
    656 //      funcStr = "SphereForm" 
    657 //      r2 = 2.54/2                     //typical 1" diameter sample, convert to radius in cm 
     649        String coefStr,abortStr 
    658650 
    659651        NVAR doMonteCarlo = root:Packages:NIST:SAS:gDoMonteCarlo                // == 1 if MC, 0 if other 
    … …  
    680672                endif 
    681673                 
     674                Variable sig_sas 
    682675                FUNCREF SANSModelAAO_MCproto func=$funcStr 
    683676                WAVE results = root:Packages:NIST:SAS:results 
     677                WAVE linear_data = root:Packages:NIST:SAS:linear_data 
     678                WAVE data = root:Packages:NIST:SAS:data 
     679 
    684680                results = 0 
     681                linear_data = 0 
    685682                 
    686                 Monte_SANS(imon,r1,r2,xCtr,yCtr,sdd,pixSize,thick,wavelength,sig_incoh,func,$coefStr,results) 
     683                CalculateRandomDeviate(func,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",SIG_SAS) 
     684                if(sig_sas > 100) 
     685                        sprintf abortStr,"sig_sas = %g. Please check that the model coefficients have a zero background, or the low q is well-behaved.",sig_sas 
     686                        Abort abortStr 
     687                endif 
     688                 
     689                WAVE ran_dev=$"root:Packages:NIST:SAS:ran_dev" 
     690                 
     691                Make/O/D/N=5000 root:Packages:NIST:SAS:nt=0,root:Packages:NIST:SAS:j1=0,root:Packages:NIST:SAS:j2=0 
     692                Make/O/D/N=100 root:Packages:NIST:SAS:nn=0 
     693                Make/O/D/N=11 root:Packages:NIST:SAS:inputWave=0 
     694                 
     695                WAVE nt = root:Packages:NIST:SAS:nt 
     696                WAVE j1 = root:Packages:NIST:SAS:j1 
     697                WAVE j2 = root:Packages:NIST:SAS:j2 
     698                WAVE nn = root:Packages:NIST:SAS:nn 
     699                WAVE inputWave = root:Packages:NIST:SAS:inputWave 
     700 
     701                inputWave[0] = imon 
     702                inputWave[1] = r1 
     703                inputWave[2] = r2 
     704                inputWave[3] = xCtr 
     705                inputWave[4] = yCtr 
     706                inputWave[5] = sdd 
     707                inputWave[6] = pixSize 
     708                inputWave[7] = thick 
     709                inputWave[8] = wavelength 
     710                inputWave[9] = sig_incoh 
     711                inputWave[10] = sig_sas 
     712 
     713                //initialize ran1 in the XOP by passing a negative integer 
     714                // does nothing in the Igor code 
     715                results[0] = -1*trunc(datetime)/10 
     716 
     717//              Variable t0 = stopMStimer(-2) 
     718         
     719#if exists("Monte_SANSX") 
     720        Monte_SANSX(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     721#else 
     722        Monte_SANS(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
     723#endif 
     724//              t0 = (stopMSTimer(-2) - t0)*1e-6 
     725//              Printf  "Mc sim time = %g seconds\r\r",t0        
    687726                 
    688727                // convert to absolute scale 
    … …  
    690729                // results[6] is the fraction transmitted 
    691730//              kappa = beamInten*pi*r1*r1*thick*(pixSize/sdd)^2*results[6]*(iMon/beaminten) 
    692                 kappa = thick*(pixSize/sdd)^2*results[6]*iMon *2        //why the factor of 2? 
    693                  
    694                 WAVE linear_data = root:Packages:NIST:SAS:linear_data 
    695                 WAVE data = root:Packages:NIST:SAS:data 
     731                kappa = thick*(pixSize/sdd)^2*results[6]*iMon 
     732 
    696733                linear_data = linear_data / kappa 
     734                linear_data[xCtr][yCtr] = 0                     //snip out the transmitted spike 
    697735                data = linear_data 
    698736         
    … …  
    725763                        aveint = S_Debye(1000,100,0.0,qval) 
    726764                endif 
     765                WAVE sigave=root:Packages:NIST:SAS:sigave 
     766                sigave = 0              //reset for model calculation 
    727767        endif 
    728768         
    … …  
    894934 
    895935// fake mask that uses all of the detector 
    896         Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask") 
     936        Make/D/O/N=(pixelsX,pixelsY) $(destPath + ":mask") 
    897937        Wave mask = $(destPath + ":mask") 
    898938        mask = 0 
    … …  
    914954        // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good) 
    915955        Variable defWavePts=500 
    916         Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint") 
    917         Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave") 
    918         Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar") 
     956        Make/O/D/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint") 
     957        Make/O/D/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave") 
     958        Make/O/D/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar") 
    919959 
    920960        WAVE qval = $(destPath + ":qval") 
Note: See TracChangeset for help on using the changeset viewer.