Changeset f0fb9fe in sasmodels


Ignore:
Timestamp:
Feb 9, 2016 7:03:28 AM (7 years ago)
Author:
richardh
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
093f754
Parents:
ec2ca99
Message:

Taylor series at smallest Q in hardsphere

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/hardsphere.py

    r13ed84c rf0fb9fe  
    3535from numpy import inf 
    3636 
    37 name = "hardsphere" 
    38 title = "Hard sphere structure factor, with Percus-Yevick closure" 
     37name = "hardsphere_fish" 
     38title = "Hard sphere structure factor from FISH, with Percus-Yevick closure" 
    3939description = """\ 
    4040    [Hard sphere structure factor, with Percus-Yevick closure] 
     
    5555               "volume fraction of hard spheres"], 
    5656             ] 
    57 single = False 
    5857 
    5958# No volume normalization despite having a volume parameter 
     
    6463 
    6564Iq = """ 
    66     double denom,dnum,alpha,beta,gamm,a,asq,ath,afor,rca,rsa; 
    67     double calp,cbeta,cgam,prefac,c,vstruc; 
    68     double struc; 
     65      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH; 
    6966 
    70     //  compute constants 
    71     denom = pow((1.0-volfraction),4); 
    72     dnum = pow((1.0 + 2.0*volfraction),2); 
    73     alpha = dnum/denom; 
    74     beta = -6.0*volfraction*pow((1.0 + volfraction/2.0),2)/denom; 
    75     gamm = 0.50*volfraction*dnum/denom; 
    76     // 
    77     //  calculate the structure factor 
    78     // 
    79     a = 2.0*q*effect_radius; 
    80     asq = a*a; 
    81     ath = asq*a; 
    82     afor = ath*a; 
    83     SINCOS(a,rsa,rca); 
    84     //rca = cos(a); 
    85     //rsa = sin(a); 
    86     calp = alpha*(rsa/asq - rca/a); 
    87     cbeta = beta*(2.0*rsa/asq - (asq - 2.0)*rca/ath - 2.0/ath); 
    88     cgam = gamm*(-rca/a + (4.0/a)*((3.0*asq - 6.0)*rca/afor + (asq - 6.0)*rsa/ath + 6.0/afor)); 
    89     prefac = -24.0*volfraction/a; 
    90     c = prefac*(calp + cbeta + cgam); 
    91     vstruc = 1.0/(1.0-c); 
    92     struc = vstruc; 
     67      if(fabs(effect_radius) < 1.E-12) { 
     68               HARDSPH=1.0; 
     69               return(HARDSPH); 
     70      } 
     71      D=pow((1.-volfraction),2); 
     72      A=pow((1.+2.*volfraction)/D, 2); 
     73      X=fabs(q*effect_radius*2.0); 
    9374 
    94     return(struc); 
     75      if(X < 5.E-06) { 
     76                 HARDSPH=1./A; 
     77                 return(HARDSPH); 
     78      } 
     79      X2=pow(X,2); 
     80      X4=pow(X2,2); 
     81      B=-6.*volfraction* pow((1.+0.5*volfraction)/D ,2); 
     82      G=0.5*volfraction*A; 
     83 
     84      if(X < 0.2) { 
     85      // use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
     86      // No obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
     87      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient  
     88      // for 5 or 6 significant figures but I put the X^4 one in anyway  
     89            FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
     90            // combining the terms makes things worse at smallest Q in single precision 
     91            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4; 
     92            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q 
     93            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40)*2*B; 
     94            HARDSPH= 1./(1. + volfraction*FF ); 
     95            return(HARDSPH); 
     96      } 
     97      SINCOS(X,S,C); 
     98 
     99// RKH Feb 2016, use version from FISH code as it is better than original sasview one at small Q in single precision 
     100      FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4; 
     101      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
     102 
     103// rearrange the terms, is now about same as sasmodels 
     104//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 ); 
     105//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
     106// remove 1/X2 from final line, take more powers of X inside the brackets, stil bad 
     107//      FF=A*(S/X3-C/X2) + B*(2.*S/X3 - C/X2 +2.0*(C-1.0)/X4) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 )/X2; 
     108//      HARDSPH= 1./(1. + 24.*volfraction*FF ); 
     109      return(HARDSPH); 
    95110   """ 
    96111 
Note: See TracChangeset for help on using the changeset viewer.