Changeset 6e93a02 in sasview for sansmodels/src/libigor/libTwoPhase.c

Timestamp:

Mar 30, 2010 5:53:47 PM (14 years ago)

Author:

Jae Cho <jhjcho@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

f10063e

Parents:

64017a8

Message:

updated libigor files from NIST svn

File:

• sansmodels/src/libigor/libTwoPhase.c (modified) (8 diffs)

sansmodels/src/libigor/libTwoPhase.c

@@ -116 +116 @@
 }
 
+// 6 JUL 2009 SRK changed definition of Izero scale factor to be uncorrelated with range
+//
 double
 DAB_Model(double dp[], double q)
@@ -127 +129 @@
         incoh = dp[2]; 
         
-        inten = Izero/pow((1.0 + (qval*range)*(qval*range)),2) + incoh;
+        inten = (Izero*range*range*range)/pow((1.0 + (qval*range)*(qval*range)),2) + incoh;
         
         return(inten);
@@ -152 +154 @@
         ans = G1*exp(-x*x*Rg1*Rg1/3.0);
         ans += B1*pow((erf1*erf1*erf1/x),Pow1);
+        
+        if(x == 0) {
+                ans = G1;
+        }
         
         ans *= scale;
@@ -188 +194 @@
         ans += G2*exp(-x*x*Rg2*Rg2/3.0);
         ans += B2*pow((erf2*erf2*erf2/x),Pow2);
+
+        if(x == 0) {
+                ans = G1+G2;
+        }
         
         ans *= scale;
@@ -230 +240 @@
         ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2);
         ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3);
+
+        if(x == 0) {
+                ans = G1+G2+G3;
+        }
         
         ans *= scale;
@@ -278 +292 @@
         ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*exp(-x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3);
         ans += G4*exp(-x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4);
+
+        if(x == 0) {
+                ans = G1+G2+G3+G4;
+        }
         
         ans *= scale;
@@ -283 +301 @@
         
     return(ans);
+}
+
+double
+BroadPeak(double dp[], double q)
+{
+        // variables are:                                                       
+        //[0] Porod term scaling
+        //[1] Porod exponent
+        //[2] Lorentzian term scaling
+        //[3] Lorentzian screening length [A]
+        //[4] peak location [1/A]
+        //[5] Lorentzian exponent
+        //[6] background
+        
+        double aa,nn,cc,LL,Qzero,mm,bgd,inten,qval;
+        qval= q;
+        aa = dp[0];
+        nn = dp[1];
+        cc = dp[2]; 
+        LL = dp[3]; 
+        Qzero = dp[4]; 
+        mm = dp[5]; 
+        bgd = dp[6]; 
+        
+        inten = aa/pow(qval,nn);
+        inten += cc/(1.0 + pow((fabs(qval-Qzero)*LL),mm) );
+        inten += bgd;
+        
+        return(inten);
+}
+
+double
+CorrLength(double dp[], double q)
+{
+        // variables are:                                                       
+        //[0] Porod term scaling
+        //[1] Porod exponent
+        //[2] Lorentzian term scaling
+        //[3] Lorentzian screening length [A]
+        //[4] Lorentzian exponent
+        //[5] background
+        
+        double aa,nn,cc,LL,mm,bgd,inten,qval;
+        qval= q;
+        aa = dp[0];
+        nn = dp[1];
+        cc = dp[2]; 
+        LL = dp[3]; 
+        mm = dp[4]; 
+        bgd = dp[5]; 
+        
+        inten = aa/pow(qval,nn);
+        inten += cc/(1.0 + pow((qval*LL),mm) );
+        inten += bgd;
+        
+        return(inten);
+}
+
+double
+TwoLorentzian(double dp[], double q)
+{
+        // variables are:                                                       
+        //[0] Lorentzian term scaling
+        //[1] Lorentzian screening length [A]
+        //[2] Lorentzian exponent
+        //[3] Lorentzian #2 term scaling
+        //[4] Lorentzian #2 screening length [A]
+        //[5] Lorentzian #2 exponent
+        //[6] background
+                
+        double aa,LL1,nn,cc,LL2,mm,bgd,inten,qval;
+        qval= q;
+        aa = dp[0];
+        LL1 = dp[1];
+        nn = dp[2]; 
+        cc = dp[3]; 
+        LL2 = dp[4]; 
+        mm = dp[5]; 
+        bgd= dp[6];
+        
+        inten = aa/(1.0 + pow((qval*LL1),nn) );
+        inten += cc/(1.0 + pow((qval*LL2),mm) );
+        inten += bgd;
+        
+        return(inten);
+}
+
+double
+TwoPowerLaw(double dp[], double q)
+{
+        //[0] Coefficient
+        //[1] (-) Power @ low Q
+        //[2] (-) Power @ high Q
+        //[3] crossover Q-value
+        //[4] incoherent background
+                
+        double A, m1,m2,qc,bgd,scale,inten,qval;
+        qval= q;
+        A = dp[0];
+        m1 = dp[1];
+        m2 = dp[2]; 
+        qc = dp[3]; 
+        bgd = dp[4]; 
+        
+        if(qval<=qc){
+                inten = A*pow(qval,-1.0*m1);
+        } else {
+                scale = A*pow(qc,-1.0*m1) / pow(qc,-1.0*m2);
+                inten = scale*pow(qval,-1.0*m2);
+        }
+        
+        inten += bgd;
+        
+        return(inten);
+}
+
+double
+PolyGaussCoil(double dp[], double x)
+{
+        //w[0] = scale
+        //w[1] = radius of gyration []
+        //w[2] = polydispersity, ratio of Mw/Mn
+        //w[3] = bkg [cm-1]
+                
+        double scale,bkg,Rg,uval,Mw_Mn,inten,xi;
+
+        scale = dp[0];
+        Rg = dp[1];
+        Mw_Mn = dp[2]; 
+        bkg = dp[3]; 
+        
+        uval = Mw_Mn - 1.0;
+        if(uval == 0.0) {
+                uval = 1e-6;            //avoid divide by zero error
+        }
+        
+        xi = Rg*Rg*x*x/(1.0+2.0*uval);
+        
+        if(xi < 1e-3) {
+                return(scale+bkg);              //limiting value
+        }
+        
+        inten = 2.0*(pow((1.0+uval*xi),(-1.0/uval))+xi-1.0);
+        inten /= (1.0+uval)*xi*xi;
+
+        inten *= scale;
+        //add in the background
+        inten += bkg;   
+        return(inten);
+}
+
+double
+GaussLorentzGel(double dp[], double x)
+{
+        //[0] Gaussian scale factor
+        //[1] Gaussian (static) screening length
+        //[2] Lorentzian (fluctuation) scale factor
+        //[3] Lorentzian screening length
+        //[4] incoherent background
+                
+        double Ig0,gg,Il0,ll,bgd,inten;
+
+        Ig0 = dp[0];
+        gg = dp[1];
+        Il0 = dp[2]; 
+        ll = dp[3];
+        bgd = dp[4]; 
+        
+        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1.0 + (x*ll)*(x*ll)) + bgd;
+        
+        return(inten);
 }
 
@@ -303 +492 @@
 }
 
-
+double
+GaussianShell(double w[], double x)
+{
+        // variables are:                                                       
+        //[0] scale
+        //[1] radius ()
+        //[2] thick () (thickness parameter - this is the std. dev. of the Gaussian width of the shell)
+        //[3] polydispersity of the radius
+        //[4] sld shell (-2)
+        //[5] sld solvent
+        //[6] background (cm-1)
+        
+        double scale,rad,delrho,bkg,del,thick,pd,sig,pi;
+        double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent;
+        scale = w[0];
+        rad = w[1];
+        thick = w[2];
+        pd = w[3];
+        sldShell = w[4];
+        sldSolvent = w[5];
+        bkg = w[6];
+        
+        delrho = w[4] - w[5];
+        sig = pd*rad;
+        
+        pi = 4.0*atan(1.0);
+        
+        ///APPROXIMATION (see eqn 4 - but not a bad approximation)
+        // del is the equivalent shell thickness with sharp boundaries, centered at mean radius
+        del = thick*sqrt(2.0*pi);
+        
+        // calculate the polydisperse shell volume and the excluded volume
+        vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd);
+        vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd);
+        
+        //intensity, eqn 9(a-d)
+        exfact = exp(-2.0*sig*sig*x*x);
+        
+        t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact);
+        t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact;
+        t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact);
+        t4 = 0.5*sig*sig*(1.0+4.0*x*rad*sin(2.0*x*rad)*exfact+cos(2.0*x*rad)*(4.0*sig*sig*x*x-1.0)*exfact);
+        
+        retval = t1+t2+t3+t4;
+        retval *= exp(-1.0*x*x*thick*thick);
+        retval *= (del*del/x/x);
+        retval *= 16.0*pi*pi*delrho*delrho*scale;
+        retval *= 1.0e8;
+        
+        //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material
+//      retval /= vshell
+        retval /= vexcl;
+        //re-normalize by polydisperse sphere volume, Gaussian distribution
+        retval /= (1.0+3.0*pd*pd);
+        
+        retval += bkg;
+        
+        return(retval);
+}
+
+