# Changeset 444 for sans/Dev/trunk/NCNR_User_Procedures

Ignore:
Timestamp:
Nov 12, 2008 4:45:15 PM (14 years ago)
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

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

 r396 //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: " 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: " 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: " 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 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: " //////////////////////////////////////////////////// // - 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) End Macro PlotSmearedTwoShell(str) Proc PlotSmearedTwoShell(str) String str Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) Macro PlotSmearedThreeShell(str) Proc PlotSmearedThreeShell(str) String str Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) End Macro PlotSmearedFourShell(str) Proc PlotSmearedFourShell(str) String str Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 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

 r396 //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: " 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 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: " 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: " 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: " //////////////////////////////////////////////////// // - 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) End Macro PlotSmearedPolyTwoShell(str) Proc PlotSmearedPolyTwoShell(str) String str Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) Macro PlotSmearedPolyThreeShell(str) Proc PlotSmearedPolyThreeShell(str) String str Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) 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

 r406 ////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" 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

 r434 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

 r426 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
• ## sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf

 r432 //              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
Note: See TracChangeset for help on using the changeset viewer.