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/Common
Files:
2 edited

Legend:

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