 Timestamp:
 Nov 12, 2008 4:45:15 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Common/GaussUtils_v40.ipf
r426 r444 323 323 End //Make76GaussPoints() 324 324 325 // !!!!! reduces the length of qt and zi by one !!!!! 326 // 327 Function 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) 338 End 339 340 /// gauleg subroutine from NR to calculate weights and abscissae for 341 // GaussLegendre quadrature 342 // 343 // 344 // arrays are indexed from 1 345 // 346 Function 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 = 3e11 354 355 m=(n+1)/2 356 xm=0.5*(x2+x1) 357 xl=0.5*(x2x1) 358 for (i=1;i<=m;i+=1) 359 z=cos(pi*(i0.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*j1.0)*z*p2(j1.0)*p3)/j 367 endfor 368 pp=n*(z*p1p2)/(z*z1.0) 369 z1=z 370 z=z1p1/pp 371 while (abs(zz1) > EPS) 372 x[i]=xmxl*z 373 x[n+1i]=xm+xl*z 374 w[i]=2.0*xl/((1.0z*z)*pp*pp) 375 w[n+1i]=w[i] 376 Endfor 377 End 378 379 380 /// uses a usersupplied number of Gauss points, and generates them onthefly as needed 381 // using a Numerical Recipes routine 382 // 383 //  note that this has an extra input parameter, nord 384 // 385 //////////// 386 Function 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 //xvalue (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,size1 419 420 summ = 0.0 // initialize integral 421 ii=0 // loop counter 422 do 423 // calculate Gauss points on integration interval (qvalue for evaluation) 424 zi = ( xx[ii]*(vbva) + vb + va )/2.0 425 //calculate partial sum for the passedin 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 = (vbva)/2.0*summ 433 434 Return (answer) 435 End 436 437 438 325 439 //////////// 326 440 Function IntegrateFn5(fcn,loLim,upLim,w,x) … … 1201 1315 return(0) 1202 1316 end 1317
Note: See TracChangeset
for help on using the changeset viewer.