Changeset 217


Ignore:
Timestamp:
Dec 3, 2007 12:11:32 PM (15 years ago)
Author:
srkline
Message:

Fixed the low QR behavior of the Schulz and Rectangular distributed spheres to use the (polydisperse) Guinier approximation at QR<0.1 where machine precision is an issue. The XOP code now duplicates the Igor code.

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  
    125125                h1 *=   (1 + 15*pd^2 + 27*pd^4 +27/7*pd^6)                              //6th moment 
    126126                h1 /= (1+3*pd^2)                //3rd moment 
     127                 
    127128                Rg2 = 3/5*rad*rad*(1+28*pd^2+126*pd^4+108*pd^6+27*pd^8) 
    128129                Rg2 /= (1+15*pd^2+27*pd^4+27/7*pd^6) 
    129                  
     130                                 
    130131                h1 *= exp(-1/3*Rg2*x*x) 
    131132                h1 += bkg 
  • sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/lib/libSphere.c

    r97 r217  
    367367        double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3; 
    368368        double v1,v2,v3,g1,pq,pi,delrho,zz; 
     369        double i_zero,Rg2,zp8; 
    369370         
    370371        pi = 4.0*atan(1.0); 
     
    380381        zp7 = zz + 7.0; 
    381382        // 
     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 
    382396        aa = (zz+1)/x/ravg; 
    383397         
     
    434448        double pi,x; 
    435449        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; 
    437451         
    438452        pi = 4.0*atan(1.0); 
     
    455469        qw = x*width; 
    456470        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 
    457489        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0; 
    458490        h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw); 
Note: See TracChangeset for help on using the changeset viewer.