Changeset 444

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:

• Analysis/Models/NewModels_2008/Core_and_NShells_v40.ipf (modified) (10 diffs)
• Analysis/Models/NewModels_2008/PolyCore_and_NShells_v40.ipf (modified) (9 diffs)
• Analysis/Packages/ModelPicker/SANSModelPicker_v40.ipf (modified) (2 diffs)
• Analysis/Packages/Wrapper_v40.ipf (modified) (1 diff)
• Common/GaussUtils_v40.ipf (modified) (2 diffs)
• Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf (modified) (1 diff)

sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Core_and_NShells_v40.ipf

@@ -23 +23 @@
 //needed to calculate the model function.
 //
-Macro PlotOneShell(num,qmin,qmax)
+Proc PlotOneShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -48 +48 @@
 End
 
-Macro PlotTwoShell(num,qmin,qmax)
+Proc PlotTwoShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -73 +73 @@
 End
 
-Macro PlotThreeShell(num,qmin,qmax)
+Proc PlotThreeShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -81 +81 @@
         Make/O/D/n=(num) xwave_ThreeShell, ywave_ThreeShell
         xwave_ThreeShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
-        Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,6.4e-6,0.001}
-        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)"}
+        Make/O/D coef_ThreeShell ={1.,60,6.4e-6,10,1e-6,10,2e-6,10,3e-6,6.4e-6,0.001}
+        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)"}
         Edit parameters_ThreeShell, coef_ThreeShell
         
@@ -98 +98 @@
 End
 
-Macro PlotFourShell(num,qmin,qmax)
+Proc PlotFourShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -134 +134 @@
 ////////////////////////////////////////////////////
 // - sets up a dependency to a wrapper, not the actual SmearedModelFunction
-Macro PlotSmearedOneShell(str)                                                          
+Proc PlotSmearedOneShell(str)                                                           
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -168 +168 @@
 End
 
-Macro PlotSmearedTwoShell(str)                                                          
+Proc PlotSmearedTwoShell(str)                                                           
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -203 +203 @@
 
 
-Macro PlotSmearedThreeShell(str)                                                                
+Proc PlotSmearedThreeShell(str)                                                         
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -237 +237 @@
 End
 
-Macro PlotSmearedFourShell(str)                                                         
+Proc PlotSmearedFourShell(str)                                                          
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -676 +676 @@
         return (0)
 End
+
+//wrapper to calculate the smeared model as an AAO-Struct
+// fills the struct and calls the ususal function with the STRUCT parameter
+//
+// used only for the dependency, not for fitting
+//
+Function fSmearedTwoShell(coefW,yW,xW)
+        Wave coefW,yW,xW
+        
+        String str = getWavesDataFolder(yW,0)
+        String DF="root:"+str+":"
+        
+        WAVE resW = $(DF+str+"_res")
+        
+        STRUCT ResSmearAAOStruct fs
+        WAVE fs.coefW = coefW   
+        WAVE fs.yW = yW
+        WAVE fs.xW = xW
+        WAVE fs.resW = resW
+        
+        Variable err
+        err = SmearedTwoShell(fs)
+        
+        return (0)
+End
+
+//wrapper to calculate the smeared model as an AAO-Struct
+// fills the struct and calls the ususal function with the STRUCT parameter
+//
+// used only for the dependency, not for fitting
+//
+Function fSmearedThreeShell(coefW,yW,xW)
+        Wave coefW,yW,xW
+        
+        String str = getWavesDataFolder(yW,0)
+        String DF="root:"+str+":"
+        
+        WAVE resW = $(DF+str+"_res")
+        
+        STRUCT ResSmearAAOStruct fs
+        WAVE fs.coefW = coefW   
+        WAVE fs.yW = yW
+        WAVE fs.xW = xW
+        WAVE fs.resW = resW
+        
+        Variable err
+        err = SmearedThreeShell(fs)
+        
+        return (0)
+End
+
+//wrapper to calculate the smeared model as an AAO-Struct
+// fills the struct and calls the ususal function with the STRUCT parameter
+//
+// used only for the dependency, not for fitting
+//
+Function fSmearedFourShell(coefW,yW,xW)
+        Wave coefW,yW,xW
+        
+        String str = getWavesDataFolder(yW,0)
+        String DF="root:"+str+":"
+        
+        WAVE resW = $(DF+str+"_res")
+        
+        STRUCT ResSmearAAOStruct fs
+        WAVE fs.coefW = coefW   
+        WAVE fs.yW = yW
+        WAVE fs.xW = xW
+        WAVE fs.resW = resW
+        
+        Variable err
+        err = SmearedFourShell(fs)
+        
+        return (0)
+End

sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/PolyCore_and_NShells_v40.ipf

@@ -27 +27 @@
 //needed to calculate the model function.
 //
-Macro PlotPolyPolyOneShell(num,qmin,qmax)
+Proc PlotPolyOneShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -33 +33 @@
         Prompt qmax "Enter maximum q-value (^-1) for model: "
 //
-        Make/O/D/n=(num) xwave_PolyPolyOneShell, ywave_PolyPolyOneShell
-        xwave_PolyPolyOneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
-        Make/O/D coef_PolyPolyOneShell = {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001}
-        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)"}
-        Edit parameters_PolyPolyOneShell, coef_PolyPolyOneShell
-        
-        Variable/G root:g_PolyPolyOneShell
+        Make/O/D/n=(num) xwave_PolyOneShell, ywave_PolyOneShell
+        xwave_PolyOneShell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
+        Make/O/D coef_PolyOneShell = {1.,60,0.1,6.4e-6,10,1e-6,6.4e-6,0.001}
+        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)"}
+        Edit parameters_PolyOneShell, coef_PolyOneShell
+        
+        Variable/G root:g_PolyOneShell
         g_PolyOneShell := PolyOneShell(coef_PolyOneShell, ywave_PolyOneShell, xwave_PolyOneShell)
         Display ywave_PolyOneShell vs xwave_PolyOneShell
@@ -52 +52 @@
 End
 
-Macro PlotPolyPolyTwoShell(num,qmin,qmax)
+Proc PlotPolyTwoShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -77 +77 @@
 End
 
-Macro PlotPolyThreeShell(num,qmin,qmax)
+Proc PlotPolyThreeShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -102 +102 @@
 End
 
-Macro PlotPolyFourShell(num,qmin,qmax)
+Proc PlotPolyFourShell(num,qmin,qmax)
         Variable num=200, qmin=0.001, qmax=0.7
         Prompt num "Enter number of data points for model: "
@@ -138 +138 @@
 ////////////////////////////////////////////////////
 // - sets up a dependency to a wrapper, not the actual SmearedModelFunction
-Macro PlotSmearedPolyOneShell(str)                                                              
+Proc PlotSmearedPolyOneShell(str)                                                               
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -172 +172 @@
 End
 
-Macro PlotSmearedPolyTwoShell(str)                                                              
+Proc PlotSmearedPolyTwoShell(str)                                                               
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -207 +207 @@
 
 
-Macro PlotSmearedPolyThreeShell(str)                                                            
+Proc PlotSmearedPolyThreeShell(str)                                                             
         String str
         Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
@@ -241 +241 @@
 End
 
-Macro PlotSmearedPolyFourShell(str)                                                             
+Proc PlotSmearedPolyFourShell(str)                                                              
         String str
         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

@@ -127 +127 @@
 ////paste here... after deleting the old make statement and list
         
-  Make/O/T/N=69  SANS_Model_List
+  Make/O/T/N=80  SANS_Model_List
 
   SANS_Model_List[0] = "BE_Polyelectrolyte.ipf"
@@ -195 +195 @@
   SANS_Model_List[64] = "Vesicle_UL_and_Struct.ipf"
   SANS_Model_List[65] = "Vesicle_UL.ipf"
-  //Beta Models
+  //2008 Models
   SANS_Model_List[66] = "Core_and_NShells.ipf"
   SANS_Model_List[67] = "PolyCore_and_NShells.ipf"
   SANS_Model_List[68] = "Fractal_Polysphere.ipf"
+  SANS_Model_List[69] = "GaussLorentzGel.ipf"
+  SANS_Model_List[70] = "PolyGaussCoil.ipf"
+  SANS_Model_List[71] = "Two_Power_Law.ipf"
+  SANS_Model_List[72] = "BroadPeak.ipf"
+  SANS_Model_List[73] = "CorrelationLengthModel.ipf"
+  SANS_Model_List[74] = "TwoLorentzian.ipf"
+  SANS_Model_List[75] = "PolyGaussShell.ipf"
+  SANS_Model_List[76] = "LamellarParacrystal.ipf"
+  SANS_Model_List[77] = "SC_ParaCrystal.ipf"
+  SANS_Model_List[78] = "BCC_ParaCrystal.ipf"
+  SANS_Model_List[79] = "FCC_ParaCrystal.ipf"
 
 

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

@@ -807 +807 @@
         Variable/G V_FitError=0                         //0=no err, 1=error,(2^1+2^0)=3=singular matrix
         Variable/G V_FitQuitReason=0            //0=ok,1=maxiter,2=user stop,3=no chisq decrease
+
         
 // don't use the auto-destination with no flag, it doesn't appear to work correctly
 // dispatch the fit
+
+// currently, none of the fit functions are defined as threadsafe, so I don't think that the /NTHR flag really
+// does anything. The functions themselves can be threaded since they are AAO, and that is probably enough,
+// since it doesn't make much sense to thread threads. In addition, there is a little-publicized warning
+// in the WM help file that /C=texWave cannot be used to specify constraints for threadsafe functions!
+// The textwave would have to be parsed into a constraint matrix first, then passed as /C={cMat,cVec}.
+// -- just something to watch out for.
+
         do
 //              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

@@ -323 +323 @@
 End             //Make76GaussPoints()
 
+// !!!!! reduces the length of qt and zi by one !!!!!
+//
+Function Make_N_GaussPoints(wt,zi)
+        Wave wt,zi
+        
+        Variable num
+        num = numpnts(wt) - 1
+        
+        gauleg(-1,1,zi,wt,num)
+        
+        DeletePoints 0,1,wt,zi
+        
+        return(0)
+End
+
+/// gauleg subroutine from NR to calculate weights and abscissae for 
+// Gauss-Legendre quadrature
+//
+//
+// arrays are indexed from 1
+//
+Function gauleg( x1, x2, x, w, n)
+        Variable x1, x2
+        Wave x, w
+        Variable n
+        
+        variable m,j,i
+        variable z1,z,xm,xl,pp,p3,p2,p1
+        Variable eps = 3e-11
+
+        m=(n+1)/2
+        xm=0.5*(x2+x1)
+        xl=0.5*(x2-x1)
+        for (i=1;i<=m;i+=1) 
+                z=cos(pi*(i-0.25)/(n+0.5))
+                do 
+                        p1=1.0
+                        p2=0.0
+                        for (j=1;j<=n;j+=1) 
+                                p3=p2
+                                p2=p1
+                                p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j
+                        endfor
+                        pp=n*(z*p1-p2)/(z*z-1.0)
+                        z1=z
+                        z=z1-p1/pp
+                while (abs(z-z1) > EPS)
+                x[i]=xm-xl*z
+                x[n+1-i]=xm+xl*z
+                w[i]=2.0*xl/((1.0-z*z)*pp*pp)
+                w[n+1-i]=w[i]
+        Endfor
+End
+
+
+/// uses a user-supplied number of Gauss points, and generates them on-the-fly as needed
+// using a Numerical Recipes routine
+//
+// - note that this has an extra input parameter, nord
+//
+////////////
+Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                                
+        FUNCREF GenericQuadrature_proto fcn
+        Variable loLim,upLim    //limits of integration
+        Wave w                  //coefficients of function fcn(w,x)
+        Variable x      //x-value (q) for the calculation
+        Variable nord                   //number of quadrature points to used
+
+// local variables
+        Variable ii,va,vb,summ,yyy,zi
+        Variable answer,dum
+        String weightStr,zStr
+        
+        weightStr = "gauss"+num2iStr(nord)+"wt"
+        zStr = "gauss"+num2istr(nord)+"z"
+                
+        if (WaveExists($weightStr) == 0) // wave reference is not valid, 
+                Make/D/N=(nord+1) $weightStr,$zStr
+                Wave wt = $weightStr
+                Wave xx = $zStr         // wave references to pass
+                Make_N_GaussPoints(wt,xx)                               //generates the gauss points and removes the extra point
+        else
+                if(exists(weightStr) > 1) 
+                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
+                endif
+                Wave wt = $weightStr
+                Wave xx = $zStr         // create the wave references
+        endif
+
+        //limits of integration are input to function
+        va = loLim
+        vb = upLim
+        // Using 5 Gauss points             
+        // remember to index from 0,size-1
+
+        summ = 0.0              // initialize integral
+        ii=0                    // loop counter
+        do
+                // calculate Gauss points on integration interval (q-value for evaluation)
+                zi = ( xx[ii]*(vb-va) + vb + va )/2.0
+                //calculate partial sum for the passed-in model function        
+                yyy = wt[ii] *  fcn(w,x,zi)                                             
+                summ += yyy             //add to the running total of the quadrature
+       ii+=1            
+        while (ii<nord)                         // end of loop over quadrature points
+   
+        // calculate value of integral to return
+        answer = (vb-va)/2.0*summ
+
+        Return (answer)
+End
+
+
+
 ////////////
 Function IntegrateFn5(fcn,loLim,upLim,w,x)                              
@@ -1201 +1315 @@
         return(0)
 end
+

sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf

@@ -106 +106 @@
 //              print "filename :"+filename
                 Variable numCols = V_flag
+                
+                //changes JIL to allow 2-column data to be read in, "faking" a 3rd column of errors
+                if(numCols==2) //no errors
+                        n1 = StringFromList(1, S_waveNames ,";" )
+                        Duplicate/O $("root:"+n1), errorTmp
+                        errorTmp = 0.01*(errorTmp)+ 0.03*sqrt(errorTmp)
+                        S_waveNames+="errorTmp;"
+                        numCols=3
+                endif
                 
                 if(numCols==3)          //simple 3-column data with no resolution information