Changeset 217
 Timestamp:
 Dec 3, 2007 12:11:32 PM (15 years ago)
 Location:
 sans/Analysis/branches/ajj_23APR07
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/RectPolySpheres.ipf
r166 r217 125 125 h1 *= (1 + 15*pd^2 + 27*pd^4 +27/7*pd^6) //6th moment 126 126 h1 /= (1+3*pd^2) //3rd moment 127 127 128 Rg2 = 3/5*rad*rad*(1+28*pd^2+126*pd^4+108*pd^6+27*pd^8) 128 129 Rg2 /= (1+15*pd^2+27*pd^4+27/7*pd^6) 129 130 130 131 h1 *= exp(1/3*Rg2*x*x) 131 132 h1 += bkg 
sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/lib/libSphere.c
r97 r217 367 367 double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3; 368 368 double v1,v2,v3,g1,pq,pi,delrho,zz; 369 double i_zero,Rg2,zp8; 369 370 370 371 pi = 4.0*atan(1.0); … … 380 381 zp7 = zz + 7.0; 381 382 // 383 384 //small QR limit  use Guinier approx 385 zp8 = zz+8; 386 if(x*ravg < 0.1) { 387 i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3); 388 i_zero *= zp6*zp5*zp4/zp1/zp1/zp1; //6th moment / 3rd moment 389 Rg2 = 3.*zp8*zp7/5./(zp1*zp1)*ravg*ravg; 390 pq = i_zero*exp(x*x*Rg2/3.); 391 //pq += bkg; //unlike the Igor code, the backgorund is added in the wrapper (above) 392 return(pq); 393 } 394 // 395 382 396 aa = (zz+1)/x/ravg; 383 397 … … 434 448 double pi,x; 435 449 double scale,rad,pd,cont,bkg; //my local names 436 double inten,h1,qw,qr,width,sig,averad3 ;450 double inten,h1,qw,qr,width,sig,averad3,Rg2; 437 451 438 452 pi = 4.0*atan(1.0); … … 455 469 qw = x*width; 456 470 qr = x*rad; 471 472 // as for the numerical inatabilities at low QR, the function is calculating the sines and cosines 473 // just fine  the problem seems to be that the 474 // leading terms nearly cancel with the last term (the 6*qr... term), to within machine 475 // precision  the difference is on the order of 10^20 476 // so just use the limiting Guiner value 477 if(qr<0.1) { 478 h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3); 479 h1 *= (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) ); //6th moment 480 h1 /= (1.+3.*pd*pd); //3rd moment 481 Rg2 = 3.0/5.0*rad*rad*( 1.+28.*pow(pd,2)+126.*pow(pd,4)+108.*pow(pd,6)+27.*pow(pd,8) ); 482 Rg2 /= (1.+15.*pow(pd,2)+27.*pow(pd,4)+27./7.*pow(pd,6)); 483 h1 *= exp(1./3.*Rg2*x*x); 484 h1 += bkg; 485 return(h1); 486 } 487 488 // normal calculation 457 489 h1 = 0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0; 458 490 h1 = 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
Note: See TracChangeset
for help on using the changeset viewer.