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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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.