Ignore:
Timestamp:
Nov 18, 2008 3:26:32 PM (14 years ago)
Author:
srkline
Message:

Additions to the library of the 2008 model functions. Direct proting of the Igor code, duplicated by these XOPs

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c

    r97 r453  
    285285} 
    286286 
     287double 
     288BroadPeak(double dp[], double q) 
     289{ 
     290        // variables are:                                                        
     291        //[0] Porod term scaling 
     292        //[1] Porod exponent 
     293        //[2] Lorentzian term scaling 
     294        //[3] Lorentzian screening length [A] 
     295        //[4] peak location [1/A] 
     296        //[5] Lorentzian exponent 
     297        //[6] background 
     298         
     299        double aa,nn,cc,LL,Qzero,mm,bgd,inten,qval; 
     300        qval= q; 
     301        aa = dp[0]; 
     302        nn = dp[1]; 
     303        cc = dp[2];  
     304        LL = dp[3];  
     305        Qzero = dp[4];  
     306        mm = dp[5];  
     307        bgd = dp[6];  
     308         
     309        inten = aa/pow(qval,nn); 
     310        inten += cc/(1.0 + pow((fabs(qval-Qzero)*LL),mm) ); 
     311        inten += bgd; 
     312         
     313        return(inten); 
     314} 
     315 
     316double 
     317CorrLength(double dp[], double q) 
     318{ 
     319        // variables are:                                                        
     320        //[0] Porod term scaling 
     321        //[1] Porod exponent 
     322        //[2] Lorentzian term scaling 
     323        //[3] Lorentzian screening length [A] 
     324        //[4] Lorentzian exponent 
     325        //[5] background 
     326         
     327        double aa,nn,cc,LL,mm,bgd,inten,qval; 
     328        qval= q; 
     329        aa = dp[0]; 
     330        nn = dp[1]; 
     331        cc = dp[2];  
     332        LL = dp[3];  
     333        mm = dp[4];  
     334        bgd = dp[5];  
     335         
     336        inten = aa/pow(qval,nn); 
     337        inten += cc/(1.0 + pow((qval*LL),mm) ); 
     338        inten += bgd; 
     339         
     340        return(inten); 
     341} 
     342 
     343double 
     344TwoLorentzian(double dp[], double q) 
     345{ 
     346        // variables are:                                                        
     347        //[0] Lorentzian term scaling 
     348        //[1] Lorentzian screening length [A] 
     349        //[2] Lorentzian exponent 
     350        //[3] Lorentzian #2 term scaling 
     351        //[4] Lorentzian #2 screening length [A] 
     352        //[5] Lorentzian #2 exponent 
     353        //[6] background 
     354                 
     355        double aa,LL1,nn,cc,LL2,mm,bgd,inten,qval; 
     356        qval= q; 
     357        aa = dp[0]; 
     358        LL1 = dp[1]; 
     359        nn = dp[2];  
     360        cc = dp[3];  
     361        LL2 = dp[4];  
     362        mm = dp[5];  
     363        bgd= dp[6]; 
     364         
     365        inten = aa/(1.0 + pow((qval*LL1),nn) ); 
     366        inten += cc/(1.0 + pow((qval*LL2),mm) ); 
     367        inten += bgd; 
     368         
     369        return(inten); 
     370} 
     371 
     372double 
     373TwoPowerLaw(double dp[], double q) 
     374{ 
     375        //[0] Coefficient 
     376        //[1] (-) Power @ low Q 
     377        //[2] (-) Power @ high Q 
     378        //[3] crossover Q-value 
     379        //[4] incoherent background 
     380                 
     381        double A, m1,m2,qc,bgd,scale,inten,qval; 
     382        qval= q; 
     383        A = dp[0]; 
     384        m1 = dp[1]; 
     385        m2 = dp[2];  
     386        qc = dp[3];  
     387        bgd = dp[4];  
     388         
     389        if(qval<=qc){ 
     390                inten = A*pow(qval,-1.0*m1); 
     391        } else { 
     392                scale = A*pow(qc,-1.0*m1) / pow(qc,-1.0*m2); 
     393                inten = scale*pow(qval,-1.0*m2); 
     394        } 
     395         
     396        inten += bgd; 
     397         
     398        return(inten); 
     399} 
     400 
     401double 
     402PolyGaussCoil(double dp[], double x) 
     403{ 
     404        //w[0] = scale 
     405        //w[1] = radius of gyration [] 
     406        //w[2] = polydispersity, ratio of Mw/Mn 
     407        //w[3] = bkg [cm-1] 
     408                 
     409        double scale,bkg,Rg,uval,Mw_Mn,inten,xi; 
     410 
     411        scale = dp[0]; 
     412        Rg = dp[1]; 
     413        Mw_Mn = dp[2];  
     414        bkg = dp[3];  
     415         
     416        uval = Mw_Mn - 1.0; 
     417        if(uval == 0) { 
     418                uval = 1e-6;            //avoid divide by zero error 
     419        } 
     420         
     421        xi = Rg*Rg*x*x/(1.0+2.0*uval); 
     422         
     423        if(xi < 1e-3) { 
     424                return(scale+bkg);              //limiting value 
     425        } 
     426         
     427        inten = 2.0*(pow((1.0+uval*xi),(-1.0/uval))+xi-1.0); 
     428        inten /= (1.0+uval)*xi*xi; 
     429 
     430        inten *= scale; 
     431        //add in the background 
     432        inten += bkg;    
     433        return(inten); 
     434} 
     435 
     436double 
     437GaussLorentzGel(double dp[], double x) 
     438{ 
     439        //[0] Gaussian scale factor 
     440        //[1] Gaussian (static) screening length 
     441        //[2] Lorentzian (fluctuation) scale factor 
     442        //[3] Lorentzian screening length 
     443        //[4] incoherent background 
     444                 
     445        double Ig0,gg,Il0,ll,bgd,inten; 
     446 
     447        Ig0 = dp[0]; 
     448        gg = dp[1]; 
     449        Il0 = dp[2];  
     450        ll = dp[3]; 
     451        bgd = dp[4];  
     452         
     453        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1 + (x*ll)*(x*ll)) + bgd; 
     454         
     455        return(inten); 
     456} 
     457 
    287458 
    288459static double 
     
    303474} 
    304475 
    305  
     476double 
     477GaussianShell(double w[], double x) 
     478{ 
     479        // variables are:                                                        
     480        //[0] scale 
     481        //[1] radius () 
     482        //[2] thick () (thickness parameter - this is the std. dev. of the Gaussian width of the shell) 
     483        //[3] polydispersity of the radius 
     484        //[4] sld shell (-2) 
     485        //[5] sld solvent 
     486        //[6] background (cm-1) 
     487         
     488        double scale,rad,delrho,bkg,del,thick,pd,sig,pi; 
     489        double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent; 
     490        scale = w[0]; 
     491        rad = w[1]; 
     492        thick = w[2]; 
     493        pd = w[3]; 
     494        sldShell = w[4]; 
     495        sldSolvent = w[5]; 
     496        bkg = w[6]; 
     497         
     498        delrho = w[4] - w[5]; 
     499        sig = pd*rad; 
     500         
     501        pi = 4.0*atan(1.0); 
     502         
     503        ///APPROXIMATION (see eqn 4 - but not a bad approximation) 
     504        // del is the equivalent shell thickness with sharp boundaries, centered at mean radius 
     505        del = thick*sqrt(2.0*pi); 
     506         
     507        // calculate the polydisperse shell volume and the excluded volume 
     508        vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd); 
     509        vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd); 
     510         
     511        //intensity, eqn 9(a-d) 
     512        exfact = exp(-2.0*sig*sig*x*x); 
     513         
     514        t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact); 
     515        t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact; 
     516        t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact); 
     517        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); 
     518         
     519        retval = t1+t2+t3+t4; 
     520        retval *= exp(-1.0*x*x*thick*thick); 
     521        retval *= (del*del/x/x); 
     522        retval *= 16.0*pi*pi*delrho*delrho*scale; 
     523        retval *= 1.0e8; 
     524         
     525        //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material 
     526//      retval /= vshell 
     527        retval /= vexcl; 
     528        //re-normalize by polydisperse sphere volume, Gaussian distribution 
     529        retval /= (1.0+3.0*pd*pd); 
     530         
     531        retval += bkg; 
     532         
     533        return(retval); 
     534} 
     535 
     536 
Note: See TracChangeset for help on using the changeset viewer.