Changeset 444


Ignore:
Timestamp:
Nov 12, 2008 4:45:15 PM (14 years ago)
Author:
srkline
Message:

Changes to the analysis package to add a few more model functions. Documentation and XOPs are to follow later.

General n-point gaussian quadrature has been added to GaussUtils? by including a Gauss-Laguere point generator from Numerical Recipes. See the Paracrystal models for an example, since they needed more than the 76 point quadrature.

Location:
sans/Dev/trunk/NCNR_User_Procedures
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Core_and_NShells_v40.ipf

    r396 r444  
    2323//needed to calculate the model function. 
    2424// 
    25 Macro PlotOneShell(num,qmin,qmax) 
     25Proc PlotOneShell(num,qmin,qmax) 
    2626        Variable num=200, qmin=0.001, qmax=0.7 
    2727        Prompt num "Enter number of data points for model: " 
     
    4848End 
    4949 
    50 Macro PlotTwoShell(num,qmin,qmax) 
     50Proc PlotTwoShell(num,qmin,qmax) 
    5151        Variable num=200, qmin=0.001, qmax=0.7 
    5252        Prompt num "Enter number of data points for model: " 
     
    7373End 
    7474 
    75 Macro PlotThreeShell(num,qmin,qmax) 
     75Proc PlotThreeShell(num,qmin,qmax) 
    7676        Variable num=200, qmin=0.001, qmax=0.7 
    7777        Prompt num "Enter number of data points for model: " 
     
    8181        Make/O/D/n=(num) xwave_ThreeShell, ywave_ThreeShell 
    8282        xwave_ThreeShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) 
    83         Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001} 
    84         make/o/t parameters_ThreeShell = {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"} 
     83        Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001} 
     84        make/o/t parameters_ThreeShell = {"scale","core radius (A)","Core SLD (A-2)","Shell 1 thickness","Shell 1 SLD (A-2)","Shell 2 thickness","Shell 2 SLD (A-2)","Shell 3 thickness","Shell 3 SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"} 
    8585        Edit parameters_ThreeShell, coef_ThreeShell 
    8686         
     
    9898End 
    9999 
    100 Macro PlotFourShell(num,qmin,qmax) 
     100Proc PlotFourShell(num,qmin,qmax) 
    101101        Variable num=200, qmin=0.001, qmax=0.7 
    102102        Prompt num "Enter number of data points for model: " 
     
    134134//////////////////////////////////////////////////// 
    135135// - sets up a dependency to a wrapper, not the actual SmearedModelFunction 
    136 Macro PlotSmearedOneShell(str)                                                           
     136Proc PlotSmearedOneShell(str)                                                            
    137137        String str 
    138138        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    168168End 
    169169 
    170 Macro PlotSmearedTwoShell(str)                                                           
     170Proc PlotSmearedTwoShell(str)                                                            
    171171        String str 
    172172        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    203203 
    204204 
    205 Macro PlotSmearedThreeShell(str)                                                                 
     205Proc PlotSmearedThreeShell(str)                                                          
    206206        String str 
    207207        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    237237End 
    238238 
    239 Macro PlotSmearedFourShell(str)                                                          
     239Proc PlotSmearedFourShell(str)                                                           
    240240        String str 
    241241        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    676676        return (0) 
    677677End 
     678 
     679//wrapper to calculate the smeared model as an AAO-Struct 
     680// fills the struct and calls the ususal function with the STRUCT parameter 
     681// 
     682// used only for the dependency, not for fitting 
     683// 
     684Function fSmearedTwoShell(coefW,yW,xW) 
     685        Wave coefW,yW,xW 
     686         
     687        String str = getWavesDataFolder(yW,0) 
     688        String DF="root:"+str+":" 
     689         
     690        WAVE resW = $(DF+str+"_res") 
     691         
     692        STRUCT ResSmearAAOStruct fs 
     693        WAVE fs.coefW = coefW    
     694        WAVE fs.yW = yW 
     695        WAVE fs.xW = xW 
     696        WAVE fs.resW = resW 
     697         
     698        Variable err 
     699        err = SmearedTwoShell(fs) 
     700         
     701        return (0) 
     702End 
     703 
     704//wrapper to calculate the smeared model as an AAO-Struct 
     705// fills the struct and calls the ususal function with the STRUCT parameter 
     706// 
     707// used only for the dependency, not for fitting 
     708// 
     709Function fSmearedThreeShell(coefW,yW,xW) 
     710        Wave coefW,yW,xW 
     711         
     712        String str = getWavesDataFolder(yW,0) 
     713        String DF="root:"+str+":" 
     714         
     715        WAVE resW = $(DF+str+"_res") 
     716         
     717        STRUCT ResSmearAAOStruct fs 
     718        WAVE fs.coefW = coefW    
     719        WAVE fs.yW = yW 
     720        WAVE fs.xW = xW 
     721        WAVE fs.resW = resW 
     722         
     723        Variable err 
     724        err = SmearedThreeShell(fs) 
     725         
     726        return (0) 
     727End 
     728 
     729//wrapper to calculate the smeared model as an AAO-Struct 
     730// fills the struct and calls the ususal function with the STRUCT parameter 
     731// 
     732// used only for the dependency, not for fitting 
     733// 
     734Function fSmearedFourShell(coefW,yW,xW) 
     735        Wave coefW,yW,xW 
     736         
     737        String str = getWavesDataFolder(yW,0) 
     738        String DF="root:"+str+":" 
     739         
     740        WAVE resW = $(DF+str+"_res") 
     741         
     742        STRUCT ResSmearAAOStruct fs 
     743        WAVE fs.coefW = coefW    
     744        WAVE fs.yW = yW 
     745        WAVE fs.xW = xW 
     746        WAVE fs.resW = resW 
     747         
     748        Variable err 
     749        err = SmearedFourShell(fs) 
     750         
     751        return (0) 
     752End 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/PolyCore_and_NShells_v40.ipf

    r396 r444  
    2727//needed to calculate the model function. 
    2828// 
    29 Macro PlotPolyPolyOneShell(num,qmin,qmax) 
     29Proc PlotPolyOneShell(num,qmin,qmax) 
    3030        Variable num=200, qmin=0.001, qmax=0.7 
    3131        Prompt num "Enter number of data points for model: " 
     
    3333        Prompt qmax "Enter maximum q-value (^-1) for model: " 
    3434// 
    35         Make/O/D/n=(num) xwave_PolyPolyOneShell, ywave_PolyPolyOneShell 
    36         xwave_PolyPolyOneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) 
    37         Make/O/D coef_PolyPolyOneShell = {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001} 
    38         make/o/t parameters_PolyPolyOneShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"} 
    39         Edit parameters_PolyPolyOneShell, coef_PolyPolyOneShell 
    40          
    41         Variable/G root:g_PolyPolyOneShell 
     35        Make/O/D/n=(num) xwave_PolyOneShell, ywave_PolyOneShell 
     36        xwave_PolyOneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) 
     37        Make/O/D coef_PolyOneShell = {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001} 
     38        make/o/t parameters_PolyOneShell = {"scale","core radius (A)","Core Polydispersity(0,1)","Core SLD (A-2)","Shell thickness (A)","Shell SLD (A-2)","Solvent SLD (A-2)","bkg (cm-1)"} 
     39        Edit parameters_PolyOneShell, coef_PolyOneShell 
     40         
     41        Variable/G root:g_PolyOneShell 
    4242        g_PolyOneShell := PolyOneShell(coef_PolyOneShell, ywave_PolyOneShell, xwave_PolyOneShell) 
    4343        Display ywave_PolyOneShell vs xwave_PolyOneShell 
     
    5252End 
    5353 
    54 Macro PlotPolyPolyTwoShell(num,qmin,qmax) 
     54Proc PlotPolyTwoShell(num,qmin,qmax) 
    5555        Variable num=200, qmin=0.001, qmax=0.7 
    5656        Prompt num "Enter number of data points for model: " 
     
    7777End 
    7878 
    79 Macro PlotPolyThreeShell(num,qmin,qmax) 
     79Proc PlotPolyThreeShell(num,qmin,qmax) 
    8080        Variable num=200, qmin=0.001, qmax=0.7 
    8181        Prompt num "Enter number of data points for model: " 
     
    102102End 
    103103 
    104 Macro PlotPolyFourShell(num,qmin,qmax) 
     104Proc PlotPolyFourShell(num,qmin,qmax) 
    105105        Variable num=200, qmin=0.001, qmax=0.7 
    106106        Prompt num "Enter number of data points for model: " 
     
    138138//////////////////////////////////////////////////// 
    139139// - sets up a dependency to a wrapper, not the actual SmearedModelFunction 
    140 Macro PlotSmearedPolyOneShell(str)                                                               
     140Proc PlotSmearedPolyOneShell(str)                                                                
    141141        String str 
    142142        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    172172End 
    173173 
    174 Macro PlotSmearedPolyTwoShell(str)                                                               
     174Proc PlotSmearedPolyTwoShell(str)                                                                
    175175        String str 
    176176        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    207207 
    208208 
    209 Macro PlotSmearedPolyThreeShell(str)                                                             
     209Proc PlotSmearedPolyThreeShell(str)                                                              
    210210        String str 
    211211        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
     
    241241End 
    242242 
    243 Macro PlotSmearedPolyFourShell(str)                                                              
     243Proc PlotSmearedPolyFourShell(str)                                                               
    244244        String str 
    245245        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/ModelPicker/SANSModelPicker_v40.ipf

    r406 r444  
    127127////paste here... after deleting the old make statement and list 
    128128         
    129   Make/O/T/N=69  SANS_Model_List 
     129  Make/O/T/N=80  SANS_Model_List 
    130130 
    131131  SANS_Model_List[0] = "BE_Polyelectrolyte.ipf" 
     
    195195  SANS_Model_List[64] = "Vesicle_UL_and_Struct.ipf" 
    196196  SANS_Model_List[65] = "Vesicle_UL.ipf" 
    197   //Beta Models 
     197  //2008 Models 
    198198  SANS_Model_List[66] = "Core_and_NShells.ipf" 
    199199  SANS_Model_List[67] = "PolyCore_and_NShells.ipf" 
    200200  SANS_Model_List[68] = "Fractal_Polysphere.ipf" 
     201  SANS_Model_List[69] = "GaussLorentzGel.ipf" 
     202  SANS_Model_List[70] = "PolyGaussCoil.ipf" 
     203  SANS_Model_List[71] = "Two_Power_Law.ipf" 
     204  SANS_Model_List[72] = "BroadPeak.ipf" 
     205  SANS_Model_List[73] = "CorrelationLengthModel.ipf" 
     206  SANS_Model_List[74] = "TwoLorentzian.ipf" 
     207  SANS_Model_List[75] = "PolyGaussShell.ipf" 
     208  SANS_Model_List[76] = "LamellarParacrystal.ipf" 
     209  SANS_Model_List[77] = "SC_ParaCrystal.ipf" 
     210  SANS_Model_List[78] = "BCC_ParaCrystal.ipf" 
     211  SANS_Model_List[79] = "FCC_ParaCrystal.ipf" 
    201212 
    202213 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf

    r434 r444  
    807807        Variable/G V_FitError=0                         //0=no err, 1=error,(2^1+2^0)=3=singular matrix 
    808808        Variable/G V_FitQuitReason=0            //0=ok,1=maxiter,2=user stop,3=no chisq decrease 
     809 
    809810         
    810811// don't use the auto-destination with no flag, it doesn't appear to work correctly 
    811812// dispatch the fit 
     813 
     814// currently, none of the fit functions are defined as threadsafe, so I don't think that the /NTHR flag really 
     815// does anything. The functions themselves can be threaded since they are AAO, and that is probably enough, 
     816// since it doesn't make much sense to thread threads. In addition, there is a little-publicized warning 
     817// in the WM help file that /C=texWave cannot be used to specify constraints for threadsafe functions! 
     818// The textwave would have to be parsed into a constraint matrix first, then passed as /C={cMat,cVec}. 
     819// -- just something to watch out for. 
     820 
    812821        do 
    813822//              Variable t0 = stopMStimer(-2)           // corresponding print is at the end of the do-while loop (outside) 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf

    r426 r444  
    323323End             //Make76GaussPoints() 
    324324 
     325// !!!!! reduces the length of qt and zi by one !!!!! 
     326// 
     327Function Make_N_GaussPoints(wt,zi) 
     328        Wave wt,zi 
     329         
     330        Variable num 
     331        num = numpnts(wt) - 1 
     332         
     333        gauleg(-1,1,zi,wt,num) 
     334         
     335        DeletePoints 0,1,wt,zi 
     336         
     337        return(0) 
     338End 
     339 
     340/// gauleg subroutine from NR to calculate weights and abscissae for  
     341// Gauss-Legendre quadrature 
     342// 
     343// 
     344// arrays are indexed from 1 
     345// 
     346Function gauleg( x1, x2, x, w, n) 
     347        Variable x1, x2 
     348        Wave x, w 
     349        Variable n 
     350         
     351        variable m,j,i 
     352        variable z1,z,xm,xl,pp,p3,p2,p1 
     353        Variable eps = 3e-11 
     354 
     355        m=(n+1)/2 
     356        xm=0.5*(x2+x1) 
     357        xl=0.5*(x2-x1) 
     358        for (i=1;i<=m;i+=1)  
     359                z=cos(pi*(i-0.25)/(n+0.5)) 
     360                do  
     361                        p1=1.0 
     362                        p2=0.0 
     363                        for (j=1;j<=n;j+=1)  
     364                                p3=p2 
     365                                p2=p1 
     366                                p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j 
     367                        endfor 
     368                        pp=n*(z*p1-p2)/(z*z-1.0) 
     369                        z1=z 
     370                        z=z1-p1/pp 
     371                while (abs(z-z1) > EPS) 
     372                x[i]=xm-xl*z 
     373                x[n+1-i]=xm+xl*z 
     374                w[i]=2.0*xl/((1.0-z*z)*pp*pp) 
     375                w[n+1-i]=w[i] 
     376        Endfor 
     377End 
     378 
     379 
     380/// uses a user-supplied number of Gauss points, and generates them on-the-fly as needed 
     381// using a Numerical Recipes routine 
     382// 
     383// - note that this has an extra input parameter, nord 
     384// 
     385//////////// 
     386Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                                 
     387        FUNCREF GenericQuadrature_proto fcn 
     388        Variable loLim,upLim    //limits of integration 
     389        Wave w                  //coefficients of function fcn(w,x) 
     390        Variable x      //x-value (q) for the calculation 
     391        Variable nord                   //number of quadrature points to used 
     392 
     393// local variables 
     394        Variable ii,va,vb,summ,yyy,zi 
     395        Variable answer,dum 
     396        String weightStr,zStr 
     397         
     398        weightStr = "gauss"+num2iStr(nord)+"wt" 
     399        zStr = "gauss"+num2istr(nord)+"z" 
     400                 
     401        if (WaveExists($weightStr) == 0) // wave reference is not valid,  
     402                Make/D/N=(nord+1) $weightStr,$zStr 
     403                Wave wt = $weightStr 
     404                Wave xx = $zStr         // wave references to pass 
     405                Make_N_GaussPoints(wt,xx)                               //generates the gauss points and removes the extra point 
     406        else 
     407                if(exists(weightStr) > 1)  
     408                         Abort "wave name is already in use"            //executed only if name is in use elsewhere 
     409                endif 
     410                Wave wt = $weightStr 
     411                Wave xx = $zStr         // create the wave references 
     412        endif 
     413 
     414        //limits of integration are input to function 
     415        va = loLim 
     416        vb = upLim 
     417        // Using 5 Gauss points              
     418        // remember to index from 0,size-1 
     419 
     420        summ = 0.0              // initialize integral 
     421        ii=0                    // loop counter 
     422        do 
     423                // calculate Gauss points on integration interval (q-value for evaluation) 
     424                zi = ( xx[ii]*(vb-va) + vb + va )/2.0 
     425                //calculate partial sum for the passed-in model function         
     426                yyy = wt[ii] *  fcn(w,x,zi)                                              
     427                summ += yyy             //add to the running total of the quadrature 
     428       ii+=1             
     429        while (ii<nord)                         // end of loop over quadrature points 
     430    
     431        // calculate value of integral to return 
     432        answer = (vb-va)/2.0*summ 
     433 
     434        Return (answer) 
     435End 
     436 
     437 
     438 
    325439//////////// 
    326440Function IntegrateFn5(fcn,loLim,upLim,w,x)                               
     
    12011315        return(0) 
    12021316end 
     1317 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf

    r432 r444  
    106106//              print "filename :"+filename 
    107107                Variable numCols = V_flag 
     108                 
     109                //changes JIL to allow 2-column data to be read in, "faking" a 3rd column of errors 
     110                if(numCols==2) //no errors 
     111                        n1 = StringFromList(1, S_waveNames ,";" ) 
     112                        Duplicate/O $("root:"+n1), errorTmp 
     113                        errorTmp = 0.01*(errorTmp)+ 0.03*sqrt(errorTmp) 
     114                        S_waveNames+="errorTmp;" 
     115                        numCols=3 
     116                endif 
    108117                 
    109118                if(numCols==3)          //simple 3-column data with no resolution information 
Note: See TracChangeset for help on using the changeset viewer.