Changeset 632


Ignore:
Timestamp:
Mar 1, 2010 1:57:12 PM (12 years ago)
Author:
srkline
Message:

incorporating Jae-Hie's changes to the library functions to explicitly handle the cases at qr==0. Also a lot of changes to explicitly declare constant values as floating point, rather than chance the interpretation as integer. (probably not necessary)

Location:
sans/XOP_Dev/SANSAnalysis/lib
Files:
4 edited

Legend:

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

    r594 r632  
    2828         
    2929        Pi = 4.0*atan(1.0); 
    30         lolim = 0; 
     30        lolim = 0.0; 
    3131        uplim = Pi/2.0; 
    3232         
     
    8484        double va,vb;           //upper and lower integration limits 
    8585        double summ,zi,yyy,answer,vell;                 //running tally of integration 
    86         double summj,vaj,vbj,zij,arg;                   //for the inner integration 
    87          
     86        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration 
     87 
    8888        Pi = 4.0*atan(1.0); 
    89         va = 0; 
    90         vb = 1;         //orintational average, outer integral 
    91         vaj=0; 
     89        va = 0.0; 
     90        vb = 1.0;               //orintational average, outer integral 
     91        vaj=0.0; 
    9292        vbj=Pi;         //endpoints of inner integral 
    9393         
     
    107107                summj=0; 
    108108                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy 
    109                 arg = ra*sqrt(1-zi*zi); 
     109                arg = ra*sqrt(1.0-zi*zi); 
    110110                for(j=0;j<nord;j++) { 
    111111                        //76 gauss points for the inner integral as well 
     
    120120                 
    121121                //now calculate outer integral 
    122                 arg = q*length*zi/2; 
    123                 yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg; 
     122                arg = q*length*zi/2.0; 
     123                if (arg == 0.0){ 
     124                        si = 1.0; 
     125                }else{ 
     126                        si = sin(arg) * sin(arg) / arg / arg; 
     127                } 
     128                yyy = Gauss76Wt[i] * answer * si; 
    124129                summ += yyy; 
    125130        } 
     
    162167        double va,vb;           //upper and lower integration limits 
    163168        double summ,zi,yyy,answer,vell;                 //running tally of integration 
    164         double summj,vaj,vbj,zij,arg;                   //for the inner integration 
    165          
     169        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration 
     170 
    166171        Pi = 4.0*atan(1.0); 
    167         va = 0; 
    168         vb = 1;         //orintational average, outer integral 
    169         vaj=0; 
     172        va = 0.0; 
     173        vb = 1.0;               //orintational average, outer integral 
     174        vaj=0.0; 
    170175        vbj=Pi;         //endpoints of inner integral 
    171176         
     
    185190                summj=0; 
    186191                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy 
    187                 arg = ra*sqrt(1-zi*zi); 
     192                arg = ra*sqrt(1.0-zi*zi); 
    188193                for(j=0;j<nordj;j++) { 
    189194                        //20 gauss points for the inner integral 
     
    198203                 
    199204                //now calculate outer integral 
    200                 arg = q*length*zi/2; 
    201                 yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg; 
     205                arg = q*length*zi/2.0; 
     206                if (arg == 0.0){ 
     207                        si = 1.0; 
     208                }else{ 
     209                        si = sin(arg) * sin(arg) / arg / arg; 
     210                } 
     211                yyy = Gauss76Wt[i] * answer * si; 
    202212                summ += yyy; 
    203213        } 
     
    215225        answer *= scale; 
    216226        // add in the background 
    217         answer += bkg;   
    218          
     227        answer += bkg; 
     228 
    219229        return answer; 
    220230} 
     
    316326         
    317327        //      Pi = 4.0*atan(1.0); 
    318         va = 0; 
    319         vb = 1;         //orintational average, outer integral 
    320         vaj = 0; 
    321         vbj = 1;                //endpoints of inner integral 
    322          
     328        va = 0.0; 
     329        vb = 1.0;               //orintational average, outer integral 
     330        vaj = 0.0; 
     331        vbj = 1.0;              //endpoints of inner integral 
     332 
    323333        summ = 0.0;                     //initialize intergral 
    324334         
     
    340350        for(i=0;i<nordi;i++) { 
    341351                //setup inner integral over the ellipsoidal cross-section 
    342                 summj=0; 
     352                summj=0.0; 
    343353                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy 
    344354                 
     
    346356                        //76 gauss points for the inner integral 
    347357                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy 
    348                         mudum = mu*sqrt(1-sigma*sigma); 
     358                        mudum = mu*sqrt(1.0-sigma*sigma); 
    349359                        yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu); 
    350360                        summj += yyy; 
     
    352362                //now calculate the value of the inner integral 
    353363                answer = (vbj-vaj)/2.0*summj; 
    354                  
    355                 arg = mu*cc*sigma/2; 
    356                 if ( arg == 0 ) { 
    357                         answer *= 1; 
     364 
     365                arg = mu*cc*sigma/2.0; 
     366                if ( arg == 0.0 ) { 
     367                        answer *= 1.0; 
    358368                } else { 
    359369                        answer *= sin(arg)*sin(arg)/arg/arg; 
     
    401411         
    402412        pi = 4.0*atan(1.0); 
    403         va = 0; 
    404         vb = 1;         //limits of numerical integral 
    405          
     413        va = 0.0; 
     414        vb = 1.0;               //limits of numerical integral 
     415 
    406416        summ = 0.0;                     //initialize intergral 
    407417         
     
    456466         
    457467        pi = 4.0*atan(1.0); 
    458         va = 0; 
    459         vb = 1;         //limits of numerical integral 
    460          
     468        va = 0.0; 
     469        vb = 1.0;               //limits of numerical integral 
     470 
    461471        summ = 0.0;                     //initialize intergral 
    462472         
     
    478488        answer *= delrho*delrho; 
    479489        //normalize by volume 
    480         answer *= 4*pi/3*a*a*nua; 
     490        answer *= 4.0*pi/3.0*a*a*nua; 
    481491        //convert to [cm-1] 
    482492        answer *= 1.0e8; 
     
    521531         
    522532        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution 
    523         if(lolim<0) { 
    524                 lolim = 0; 
     533        if(lolim<0.0) { 
     534                lolim = 0.0; 
    525535        } 
    526536        if(pd>0.3) { 
     
    582592         
    583593        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution 
    584         if(lolim<0) { 
    585                 lolim = 0; 
     594        if(lolim<0.0) { 
     595                lolim = 0.0; 
    586596        } 
    587597        if(pd>0.3) { 
     
    697707         
    698708        lolim = exp(log(radius)-(4.*sigma)); 
    699         if (lolim<0) { 
    700                 lolim=0;                //to avoid numerical error when  va<0 (-ve r value) 
     709        if (lolim<0.0) { 
     710                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value) 
    701711        } 
    702712        uplim = exp(log(radius)+(4.*sigma)); 
     
    758768        delpc = sldc - slds;    //core - shell 
    759769        delps = slds - sld;             //shell - solvent 
    760         bkg = dp[8];                         
    761          
     770        bkg = dp[8]; 
     771 
    762772        for(i=0;i<nord;i++) { 
    763773                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     
    768778        answer = (uplim-lolim)/2.0*summ; 
    769779        // normalize by particle volume 
    770         oblatevol = 4*Pi/3*trmaj*trmaj*trmin; 
     780        oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin; 
    771781        answer /= oblatevol; 
    772782         
     
    812822        delpc = sldc - slds;            //core - shell 
    813823        delps = slds - sld;             //shell  - sovent 
    814         bkg = dp[8];                         
    815          
     824        bkg = dp[8]; 
     825 
    816826        for(i=0;i<nord;i++) { 
    817827                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     
    822832        answer = (uplim-lolim)/2.0*summ; 
    823833        // normalize by particle volume 
    824         prolatevol = 4*Pi/3*trmaj*trmin*trmin; 
     834        prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin; 
    825835        answer /= prolatevol; 
    826836         
     
    12321242         
    12331243        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution 
    1234         if(lolim<0) { 
    1235                 lolim = 0; 
     1244        if(lolim<0.0) { 
     1245                lolim = 0.0; 
    12361246        } 
    12371247        if(pd>0.3) { 
     
    12961306         
    12971307        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution 
    1298         if(lolim<0) { 
    1299                 lolim = 0; 
     1308        if(lolim<0.0) { 
     1309                lolim = 0.0; 
    13001310        } 
    13011311        if(pd>0.3) { 
     
    15701580    q02 = q0*q0; 
    15711581    q0p = pow(q0,(-4.0 + p1short) ); 
    1572      
    1573     yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2))))));  
    1574          
     1582 
     1583    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); 
     1584 
    15751585    return(yy); 
    15761586} 
     
    16101620     
    16111621    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
    1612      
    1613     t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14*b3*pow(E,(-((q02*Rg2)/b2))))/(15*q03*Rg2) + (2*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15*q02*Rg2)))*Rg2)/b)))/L; 
    1614      
    1615     t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(2*C5); 
    1616      
    1617     t4 = (b4*sqrt(Rg2)*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
    1618      
    1619     t5 = (2*b4*(((2*q0*Rg2)/b - (2*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    1620      
    1621     t6 = (8*b4*b*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
    1622      
    1623     t7 = (((-((3*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 3/miu)))/miu)) - (2*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 2/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1) - 1/miu)))/miu)); 
    1624      
    1625     t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    1626      
    1627     t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7*b2)/(15*q02*Rg2))) + (7*b2)/(15*q02*Rg2))))/L; 
    1628      
    1629     t10 = (2*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    1630          
    1631      
    1632     yy = ((-1*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1)/miu))))*((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
    1633          
     1622 
     1623    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 
     1624 
     1625    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5); 
     1626 
     1627    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
     1628 
     1629    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     1630 
     1631    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
     1632 
     1633    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
     1634 
     1635    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     1636 
     1637    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 
     1638 
     1639    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     1640 
     1641 
     1642    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
     1643 
    16341644    return (yy); 
    16351645} 
     
    17281738{ 
    17291739        // local variables 
    1730         double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi; 
    1731          
     1740        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi; 
     1741 
    17321742        Pi = 4.0*atan(1.0); 
    17331743         
     
    17411751        vc = pi43*aa*bb*bb; 
    17421752        vt = pi43*trmaj*trmin*trmin; 
    1743         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc; 
    1744         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps; 
     1753        if (uq == 0.0){ 
     1754                siq = 1.0/3.0; 
     1755        }else{ 
     1756                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; 
     1757        } 
     1758        if (ut == 0.0){ 
     1759                sit = 1.0/3.0; 
     1760        }else{ 
     1761                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; 
     1762        } 
     1763        gfnc = 3.0*siq*vc*delpc; 
     1764        gfnt = 3.0*sit*vt*delps; 
    17451765        gfn = gfnc+gfnt; 
    17461766        gfn2 = gfn*gfn; 
     
    17601780{ 
    17611781        // local variables 
    1762         double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 
    1763          
     1782        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 
     1783 
    17641784        Pi = 4.0*atan(1.0); 
    17651785        pi43=4.0/3.0*Pi; 
     
    17721792        vc = pi43*aa*aa*bb; 
    17731793        vt = pi43*trmaj*trmaj*trmin; 
    1774         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc; 
    1775         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps; 
     1794        if (uq == 0.0){ 
     1795                siq = 1.0/3.0; 
     1796        }else{ 
     1797                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; 
     1798        } 
     1799        if (ut == 0.0){ 
     1800                sit = 1.0/3.0; 
     1801        }else{ 
     1802                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; 
     1803        } 
     1804        gfnc = 3.0*siq*vc*delpc; 
     1805        gfnt = 3.0*sit*vt*delps; 
    17761806        tgfn = gfnc+gfnt; 
    17771807        gfn4 = tgfn*tgfn; 
     
    17901820     
    17911821    Pq = Sk_WR(q,zi,lb);                //does not have cross section term 
    1792     Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
    1793      
     1822    if (qr !=0){ 
     1823        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
     1824    }  
    17941825    vcyl=Pi*radius*radius*zi; 
    17951826    Pq *= vcyl*vcyl; 
     
    18101841     
    18111842    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term 
    1812     Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
    1813      
     1843    if (qr !=0){ 
     1844        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
     1845    } 
     1846 
    18141847    vcyl=Pi*zi*zi*Lc; 
    18151848    Pq *= vcyl*vcyl; 
     
    18301863        Pi = 4.0*atan(1.0); 
    18311864        nord = 76; 
    1832         lolim = 0; 
    1833         uplim = Pi/2; 
     1865        lolim = 0.0; 
     1866        uplim = Pi/2.0; 
    18341867        halfheight = length/2.0; 
    18351868         
     
    18561889        // length is the *Half* CORE-LENGTH of the cylinder  
    18571890        // dum is the dummy variable for the integration (theta) 
    1858          
    1859         double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 
     1891 
     1892        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 
    18601893        double Pi; 
    18611894         
     
    18711904        sinarg1 = qq*length*cos(dum); 
    18721905        sinarg2 = qq*(length+facthick)*cos(dum); 
    1873          
    1874         t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    1875         t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
    1876          
     1906        if (besarg1 == 0.0){ 
     1907                be1 = 0.5; 
     1908        }else{ 
     1909                be1 = NR_BessJ1(besarg1)/besarg1; 
     1910        } 
     1911        if (besarg2 == 0.0){ 
     1912                be2 = 0.5; 
     1913        }else{ 
     1914                be2 = NR_BessJ1(besarg2)/besarg2; 
     1915        } 
     1916        if (sinarg1 == 0.0){ 
     1917                si1 = 1.0; 
     1918        }else{ 
     1919                si1 = sin(sinarg1)/sinarg1; 
     1920        } 
     1921        if (besarg2 == 0.0){ 
     1922                si2 = 1.0; 
     1923        }else{ 
     1924                si2 = sin(sinarg2)/sinarg2; 
     1925        } 
     1926 
     1927        t1 = 2.0*vol1*dr1*si1*be1; 
     1928        t2 = 2.0*vol2*dr2*si2*be2; 
     1929 
    18771930        retval = ((t1+t2)*(t1+t2))*sin(dum); 
    18781931        return (retval); 
     
    18841937CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum) 
    18851938{ 
    1886          
    1887     double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 
     1939 
     1940    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 
    18881941    double Pi; 
    18891942     
     
    18991952    sinarg1 = qq*length*cos(dum); 
    19001953    sinarg2 = qq*(length+thick)*cos(dum); 
    1901      
    1902     t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    1903     t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
    1904      
     1954 
     1955        if (besarg1 == 0.0){ 
     1956                be1 = 0.5; 
     1957        }else{ 
     1958                be1 = NR_BessJ1(besarg1)/besarg1; 
     1959        } 
     1960        if (besarg2 == 0.0){ 
     1961                be2 = 0.5; 
     1962        }else{ 
     1963                be2 = NR_BessJ1(besarg2)/besarg2; 
     1964        } 
     1965        if (sinarg1 == 0.0){ 
     1966                si1 = 1.0; 
     1967        }else{ 
     1968                si1 = sin(sinarg1)/sinarg1; 
     1969        } 
     1970        if (besarg2 == 0.0){ 
     1971                si2 = 1.0; 
     1972        }else{ 
     1973                si2 = sin(sinarg2)/sinarg2; 
     1974        } 
     1975 
     1976    t1 = 2.0*vol1*dr1*si1*be1; 
     1977    t2 = 2.0*vol2*dr2*si2*be2; 
     1978 
    19051979    retval = ((t1+t2)*(t1+t2))*sin(dum); 
    19061980         
     
    19171991     
    19181992    Pi = 4.0*atan(1.0); 
    1919     lolim = 0; 
     1993    lolim = 0.0; 
    19201994    uplim = Pi/2.0; 
    19211995    halfheight = dumLen/2.0; 
     
    19502024        // length is the *Half* CORE-LENGTH of the cylinder = L (A) 
    19512025        // dum is the dummy variable for the integration (x in Feigin's notation) 
    1952          
    1953         //Local variables  
    1954         double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; 
     2026 
     2027        //Local variables 
     2028        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; 
    19552029        double Pi; 
    19562030        int kk; 
     
    19682042        sinarg1 = qq*length*cos(dum); 
    19692043        sinarg2 = qq*(length+thick)*cos(dum); 
    1970          
    1971         t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1); 
    1972         t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2); 
    1973          
     2044 
     2045        if (besarg1 == 0.0){ 
     2046                be1 = 0.5; 
     2047        }else{ 
     2048                be1 = NR_BessJ1(besarg1)/besarg1; 
     2049        } 
     2050        if (besarg2 == 0.0){ 
     2051                be2 = 0.5; 
     2052        }else{ 
     2053                be2 = NR_BessJ1(besarg2)/besarg2; 
     2054        } 
     2055        if (sinarg1 == 0.0){ 
     2056                si1 = 1.0; 
     2057        }else{ 
     2058                si1 = sin(sinarg1)/sinarg1; 
     2059        } 
     2060        if (besarg2 == 0.0){ 
     2061                si2 = 1.0; 
     2062        }else{ 
     2063                si2 = sin(sinarg2)/sinarg2; 
     2064        } 
     2065 
     2066        t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1); 
     2067        t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2); 
     2068 
    19742069        retval =((t1+t2)*(t1+t2))*sin(dum); 
    19752070         
     
    19982093     
    19992094    Pi = 4.0*atan(1.0); 
    2000     lolim = 0; 
     2095    lolim = 0.0; 
    20012096    uplim = Pi/2.0; 
    20022097    halfheight = length/2.0; 
     
    20612156         
    20622157    nu = nua/a; 
    2063     arg = qq*a*sqrt(1+dum*dum*(nu*nu-1)); 
    2064      
    2065     retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); 
     2158    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0)); 
     2159    if (arg == 0.0){ 
     2160        retval =1.0/3.0; 
     2161    }else{ 
     2162        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); 
     2163    } 
    20662164    retval *= retval; 
    2067     retval *= 9; 
    2068      
     2165    retval *= 9.0; 
     2166 
    20692167    return(retval); 
    20702168}//Function EllipsoidKernel() 
     
    20762174     
    20772175    gamma = rcore/rshell; 
    2078     arg1 = qq*rshell*sqrt(1-dum*dum);           //1=shell (outer radius) 
    2079     arg2 = qq*rcore*sqrt(1-dum*dum);                    //2=core (inner radius) 
    2080     lam1 = 2*NR_BessJ1(arg1)/arg1; 
    2081     lam2 = 2*NR_BessJ1(arg2)/arg2; 
    2082     psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);         //SRK 10/19/00 
    2083      
    2084     sinarg = qq*length*dum/2; 
    2085     t2 = sin(sinarg)/sinarg; 
    2086      
     2176    arg1 = qq*rshell*sqrt(1.0-dum*dum);         //1=shell (outer radius) 
     2177    arg2 = qq*rcore*sqrt(1.0-dum*dum);                  //2=core (inner radius) 
     2178    if (arg1 == 0.0){ 
     2179        lam1 = 1.0; 
     2180    }else{ 
     2181        lam1 = 2.0*NR_BessJ1(arg1)/arg1; 
     2182    } 
     2183    if (arg2 == 0.0){ 
     2184        lam2 = 1.0; 
     2185    }else{ 
     2186        lam2 = 2.0*NR_BessJ1(arg2)/arg2; 
     2187    } 
     2188    //Todo: Need to check psi behavior as gamma goes to 1. 
     2189    psi = 1.0/(1.0-gamma*gamma)*(lam1 -  gamma*gamma*lam2);             //SRK 10/19/00 
     2190    sinarg = qq*length*dum/2.0; 
     2191    if (sinarg == 0.0){ 
     2192        t2 = 1.0; 
     2193    }else{ 
     2194        t2 = sin(sinarg)/sinarg; 
     2195    } 
     2196 
    20872197    retval = psi*psi*t2*t2; 
    20882198     
     
    20992209     
    21002210    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    2101     arg1 = (mu/2)*cos(Pi*uu/2); 
    2102     arg2 = (mu*aa/2)*sin(Pi*uu/2); 
    2103     if(arg1==0) { 
    2104                 tmp1 = 1; 
     2211    arg1 = (mu/2.0)*cos(Pi*uu/2.0); 
     2212    arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0); 
     2213    if(arg1==0.0) { 
     2214                tmp1 = 1.0; 
    21052215        } else { 
    21062216                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 
    21072217    } 
    2108      
    2109     if (arg2==0) { 
    2110                 tmp2 = 1; 
     2218 
     2219    if (arg2==0.0) { 
     2220                tmp2 = 1.0; 
    21112221        } else { 
    21122222                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
     
    21302240    arg += cc*cc*dy*dy; 
    21312241    arg = q*sqrt(arg); 
    2132     val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); 
    2133      
     2242    if (arg == 0.0){ 
     2243        val = 1.0; // as arg --> 0, val should go to 1.0 
     2244    }else{ 
     2245        val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); 
     2246    } 
    21342247    return (val); 
    21352248         
     
    21442257        // rr is the radius of the cylinder (A) 
    21452258        // h is the HALF-LENGTH of the cylinder = L/2 (A) 
    2146          
    2147     double besarg,bj,retval,d1,t1,b1,t2,b2;              //Local variables  
    2148          
    2149          
     2259 
     2260    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables 
     2261 
     2262 
    21502263    besarg = qq*rr*sin(theta); 
    2151          
     2264    siarg = qq * h * cos(theta); 
    21522265    bj =NR_BessJ1(besarg); 
    21532266         
    21542267        //* Computing 2nd power */ 
    2155     d1 = sin(qq * h * cos(theta)); 
     2268    d1 = sin(siarg); 
    21562269    t1 = d1 * d1; 
    21572270        //* Computing 2nd power */ 
     
    21592272    t2 = d1 * d1 * 4.0 * sin(theta); 
    21602273        //* Computing 2nd power */ 
    2161     d1 = qq * h * cos(theta); 
     2274    d1 = siarg; 
    21622275    b1 = d1 * d1; 
    21632276        //* Computing 2nd power */ 
    21642277    d1 = qq * rr * sin(theta); 
    21652278    b2 = d1 * d1; 
    2166     retval = t1 * t2 / b1 / b2; 
    2167          
     2279    if (besarg == 0.0){ 
     2280        be = sin(theta); 
     2281    }else{ 
     2282        be = t2 / b2; 
     2283    } 
     2284    if (siarg == 0.0){ 
     2285        si = 1.0; 
     2286    }else{ 
     2287        si = t1 / b1; 
     2288    } 
     2289    retval = be * si; 
     2290 
    21682291    return (retval); 
    21692292         
     
    21812304        double retval,arg;               //Local variables  
    21822305         
    2183         arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2); 
    2184          
    2185         retval = 2*NR_BessJ1(arg)/arg; 
    2186          
     2306        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2); 
     2307        if (arg == 0.0){ 
     2308                retval = 1.0; 
     2309        }else{ 
     2310                retval = 2.0*NR_BessJ1(arg)/arg; 
     2311        } 
     2312 
    21872313        //square it 
    21882314        retval *= retval; 
     
    23812507                yyy = inner; 
    23822508                 
    2383                 if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h 
     2509                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    23842510                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
    23852511                } else { 
     
    25002626                yyy = inner; 
    25012627                 
    2502                 if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h 
     2628                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    25032629                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
    25042630                } else { 
     
    25842710                yyy = inner; 
    25852711                 
    2586                 if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h 
     2712                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    25872713                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
    25882714                } else { 
     
    27052831                yyy = inner; 
    27062832                 
    2707                 if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h 
     2833                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    27082834                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
    27092835                } else { 
     
    27892915                yyy = inner; 
    27902916                 
    2791                 if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h 
     2917                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    27922918                        yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
    27932919                } else { 
     
    28652991         
    28662992        lolim = exp(log(radius)-(4.*sigma)); 
    2867         if (lolim<0) { 
    2868                 lolim=0;                //to avoid numerical error when  va<0 (-ve r value) 
     2993        if (lolim<0.0) { 
     2994                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value) 
    28692995        } 
    28702996        uplim = exp(log(radius)+(4.*sigma)); 
     
    29083034        Pi = 4.0*atan(1.0); 
    29093035        nord = 76; 
    2910         lolim = 0; 
     3036        lolim = 0.0; 
    29113037        uplim = Pi/2; 
    29123038        halfheight = length/2.0; 
     
    29403066        dr2 = rhor-rhosolv; 
    29413067        dr3=  rhoh-rhor; 
    2942         vol1 = Pi*rad*rad*(2*length); 
    2943         vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick); 
    2944         vol3= Pi*(rad)*(rad)*(2*length+2*facthick); 
     3068        vol1 = Pi*rad*rad*(2.0*length); 
     3069        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick); 
     3070        vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick); 
    29453071        besarg1 = qq*rad*sin(dum); 
    29463072        besarg2 = qq*(rad+radthick)*sin(dum); 
     
    29483074        sinarg2 = qq*(length+facthick)*cos(dum); 
    29493075         
    2950         t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    2951         t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
    2952         t3 = 2*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; 
     3076        t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
     3077        t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
     3078        t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; 
    29533079         
    29543080        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 
     
    29983124        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0); 
    29993125                          
    3000         if(arg1==0){ 
     3126        if(arg1==0.0){ 
    30013127                t1 = 1.0; 
    30023128        } else { 
    30033129                t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)    
    30043130        } 
    3005         if(arg2==0){ 
     3131        if(arg2==0.0){ 
    30063132                t2 = 1.0; 
    30073133        } else { 
    30083134                t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)     
    30093135        }        
    3010         if(arg3==0){ 
     3136        if(arg3==0.0){ 
    30113137                t3 = 1.0; 
    30123138        } else { 
    30133139                t3 = sin(arg3)/arg3; 
    30143140        } 
    3015         if(arg4==0){ 
     3141        if(arg4==0.0){ 
    30163142                t4 = 1.0; 
    30173143        } else { 
     
    30733199        for(i=0;i<nordi;i++) { 
    30743200                //setup inner integral over the ellipsoidal cross-section 
    3075                 summj=0; 
     3201                summj=0.0; 
    30763202                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy 
    30773203                 
     
    30883214                //finish the outer integral cc already scaled 
    30893215                arg = mu*cc*sigma/2.0; 
    3090                 if ( arg == 0 ) { 
     3216                if ( arg == 0.0 ) { 
    30913217                        answer *= 1.0; 
    30923218                } else { 
  • sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

    r593 r632  
    1414{ 
    1515        double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names 
    16         double bes,f,vol,f2,pi; 
    17          
     16        double bes,f,vol,f2,pi,qr; 
     17 
    1818        pi = 4.0*atan(1.0); 
    1919        scale = dp[0]; 
     
    2424         
    2525        delrho = sldSph - sldSolv; 
    26         //handle q==0 separately 
    27         if(q==0){ 
    28                 f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg; 
    29                 return(f); 
     26        //handle qr==0 separately 
     27        qr = q*radius; 
     28        if(qr == 0.0){ 
     29                bes = 1.0; 
     30        }else{ 
     31                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    3032        } 
    3133         
    32         bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius); 
    3334        vol = 4.0*pi/3.0*radius*radius*radius; 
    34         f = vol*bes*delrho;             // [=]  
     35        f = vol*bes*delrho;             // [=] A-1 
    3536                                                        // normalize to single particle volume, convert to 1/cm 
    3637        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 
     
    6061        qr=x*rcore; 
    6162        contr = rhocore-rhoshel; 
    62         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     63        if(qr == 0.0){ 
     64                bes = 1.0; 
     65        }else{ 
     66                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     67        } 
    6368        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    6469        f = vol*bes*contr; 
     
    6671        qr=x*(rcore+thick); 
    6772        contr = rhoshel-rhosolv; 
    68         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     73        if(qr == 0.0){ 
     74                bes = 1.0; 
     75        }else{ 
     76                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     77        } 
    6978        vol = 4.0*pi/3.0*pow((rcore+thick),3); 
    7079        f += vol*bes*contr; 
    71          
    72         // normalize to particle volume and rescale from [-1] to [cm-1] 
     80 
     81        // normalize to particle volume and rescale from [A-1] to [cm-1] 
    7382        f2 = f*f/vol*1.0e8; 
    7483         
     
    102111        qr=x*rcore; 
    103112        contr = rhocore-rhoshel; 
    104         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     113        if(qr == 0){ 
     114                bes = 1.0; 
     115        }else{ 
     116                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     117        } 
    105118        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    106119        f = vol*bes*contr; 
     
    108121        qr=x*(rcore+thick); 
    109122        contr = rhoshel-rhosolv; 
    110         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     123        if(qr == 0.0){ 
     124                bes = 1.0; 
     125        }else{ 
     126                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     127        } 
    111128        vol = 4.0*pi/3.0*pow((rcore+thick),3); 
    112129        f += vol*bes*contr; 
    113          
    114         // normalize to the particle volume and rescale from [-1] to [cm-1] 
     130 
     131        // normalize to the particle volume and rescale from [A-1] to [cm-1] 
    115132        //note that for the vesicle model, the volume is ONLY the shell volume 
    116133        vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3)); 
     
    388405         
    389406        //small QR limit - use Guinier approx 
    390         zp8 = zz+8; 
     407        zp8 = zz+8.0; 
    391408        if(x*ravg < 0.1) { 
    392409                i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3); 
     
    399416        // 
    400417 
    401         aa = (zz+1)/x/ravg; 
     418        aa = (zz+1.0)/x/ravg; 
    402419         
    403420        at1 = atan(1.0/aa); 
     
    420437         
    421438        pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg); 
    422         pq = pow(10,pq)*8*pi*pi*delrho*delrho; 
     439        pq = pow(10,pq)*8.0*pi*pi*delrho*delrho; 
    423440         
    424441        // 
     
    483500        // so just use the limiting Guiner value 
    484501        if(qr<0.1) { 
    485                 h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3); 
     502                h1 = scale*cont*cont*1.e8*4.*pi/3.0*pow(rad,3); 
    486503                h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment 
    487504                h1 /= (1.+3.*pd*pd);    //3rd moment 
     
    495512        // normal calculation 
    496513        h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0; 
    497         h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw); 
    498         h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw); 
    499         h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw); 
    500         h1 += qw*qr*sin(2*qr)*cos(2*qw); 
     514        h1 -= 5.0/2.0*cos(2.0*qr)*sin(qw)*cos(qw); 
     515        h1 += 0.5*qr*qr*cos(2.0*qr)*sin(2.0*qw); 
     516        h1 += 0.5*qw*qw*cos(2.0*qr)*sin(2.0*qw); 
     517        h1 += qw*qr*sin(2.0*qr)*cos(2.0*qw); 
    501518        h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw)); 
    502519        h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw)); 
     
    550567         
    551568        va = -4.0*sig + rad; 
    552         if (va<0) { 
    553                 va=0;           //to avoid numerical error when  va<0 (-ve q-value) 
     569        if (va<0.0) { 
     570                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) 
    554571        } 
    555572        vb = 4.0*sig +rad; 
     
    604621        va = -3.5*sig + mu; 
    605622        va = exp(va); 
    606         if (va<0) { 
    607                 va=0;           //to avoid numerical error when  va<0 (-ve q-value) 
     623        if (va<0.0) { 
     624                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value) 
    608625        } 
    609626        vb = 3.5*sig*(1.0+sig) +mu; 
     
    642659         
    643660        pi = 4.0*atan(1.0); 
    644         retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
     661        retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
    645662        return(retval); 
    646663} 
     
    791808        double psf11,psf12,psf22; 
    792809        double phi1,phi2,phr,a3; 
    793         double v1,v2,n1,n2,qr1,qr2,b1,b2; 
     810        double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; 
    794811        int err; 
    795812         
     
    825842        qr1 = r1*x; 
    826843        qr2 = r2*x; 
    827          
    828         b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1; 
    829         b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2; 
     844 
     845        if (qr1 == 0){ 
     846                sc1 = 1.0/3.0; 
     847        }else{ 
     848                sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1; 
     849        } 
     850        if (qr2 == 0){ 
     851                sc2 = 1.0/3.0; 
     852        }else{ 
     853                sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2; 
     854        } 
     855        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*sc1; 
     856        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*sc2; 
    830857        inten = n1*b1*b1*psf11; 
    831858        inten += sqrt(n1*n2)*2.0*b1*b2*psf12; 
     
    941968        double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; 
    942969        double a1,a2i,a2,b1,b2,b12,gm1,gm12; 
    943         double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; 
     970        double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3; 
    944971        double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; 
    945972        double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; 
     
    10591086         
    10601087        ii=0; 
    1061         fval=0; 
     1088        fval=0.0; 
    10621089         
    10631090        do { 
    10641091                ri = rcore + (double)ii*ts + (double)ii*tw; 
    1065                 voli = 4*pi/3*ri*ri*ri; 
     1092                voli = 4.0*pi/3.0*ri*ri*ri; 
    10661093                sldi = rhocore-rhoshel; 
    10671094                fval += voli*sldi*F_func(ri*x); 
    10681095                ri += ts; 
    1069                 voli = 4*pi/3*ri*ri*ri; 
     1096                voli = 4.0*pi/3.0*ri*ri*ri; 
    10701097                sldi = rhoshel-rhocore; 
    10711098                fval += voli*sldi*F_func(ri*x); 
     
    10751102        fval *= fval;           //square it 
    10761103        fval /= voli;           //normalize by the overall volume 
    1077         fval *= scale*1e8; 
     1104        fval *= scale*1.0e8; 
    10781105        fval += bkg; 
    10791106         
     
    11321159         
    11331160        ii=minPairs; 
    1134         fval=0; 
    1135         d1 = 0; 
    1136         d2 = 0; 
    1137         r1 = 0; 
    1138         r2 = 0; 
    1139         distr = 0; 
    1140         first = 1; 
     1161        fval=0.0; 
     1162        d1 = 0.0; 
     1163        d2 = 0.0; 
     1164        r1 = 0.0; 
     1165        r2 = 0.0; 
     1166        distr = 0.0; 
     1167        first = 1.0; 
    11411168        do { 
    11421169                //make the current values old 
     
    11561183        } while(ii<=maxPairs); 
    11571184         
    1158         fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume 
     1185        fval /= 4.0*pi/3.0*avg*avg*avg;         //normalize by the overall volume 
    11591186        fval /= distr; 
    11601187        fval *= scale; 
     
    11721199        pi = 4.0*atan(1.0); 
    11731200    ii=0; 
    1174     fval=0; 
     1201    fval=0.0; 
    11751202     
    11761203    do { 
    11771204        ri = rcore + (double)ii*ts + (double)ii*tw; 
    1178         voli = 4*pi/3*ri*ri*ri; 
     1205        voli = 4.0*pi/3.0*ri*ri*ri; 
    11791206        sldi = rhocore-rhoshel; 
    11801207        fval += voli*sldi*F_func(ri*x); 
    11811208        ri += ts; 
    1182         voli = 4*pi/3*ri*ri*ri; 
     1209        voli = 4.0*pi/3.0*ri*ri*ri; 
    11831210        sldi = rhoshel-rhocore; 
    11841211        fval += voli*sldi*F_func(ri*x); 
     
    11881215    fval *= fval; 
    11891216    fval /= voli; 
    1190     fval *= 1e8; 
     1217    fval *= 1.0e8; 
    11911218     
    11921219    return(fval);       // this result still needs to be multiplied by scale and have background added 
     
    11981225    double dr; 
    11991226     
    1200     dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1)); 
     1227    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 
    12011228    return (exp(dr)); 
    12021229} 
     
    12211248double 
    12221249F_func(double qr) { 
    1223         return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); 
     1250        double sc; 
     1251        if (qr == 0.0){ 
     1252                sc = 1.0; 
     1253        }else{ 
     1254                sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr)); 
     1255        } 
     1256        return sc; 
    12241257} 
    12251258 
     
    13101343        contr = rhocore-rhoshel1; 
    13111344        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    1312         vol = 4.0*pi/3*rcore*rcore*rcore; 
     1345        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    13131346        f = vol*bes*contr; 
    13141347        //now the shell (1) 
     
    13771410        contr = rhocore-rhoshel1; 
    13781411        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    1379         vol = 4.0*pi/3*rcore*rcore*rcore; 
     1412        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    13801413        f = vol*bes*contr; 
    13811414        //now the shell (1) 
     
    14531486        contr = rhocore-rhoshel1; 
    14541487        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
    1455         vol = 4.0*pi/3*rcore*rcore*rcore; 
     1488        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    14561489        f = vol*bes*contr; 
    14571490        //now the shell (1) 
     
    15151548        range = 8.0;                    //std deviations for the integration 
    15161549        va = rcore*(1.0-range*pd); 
    1517         if (va<0) { 
    1518                 va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1550        if (va<0.0) { 
     1551                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 
    15191552        } 
    15201553        if (pd>0.3) { 
     
    15861619        range = 8.0;                    //std deviations for the integration 
    15871620        va = rcore*(1.0-range*pd); 
    1588         if (va<0) { 
    1589                 va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1621        if (va<0.0) { 
     1622                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0 
    15901623        } 
    15911624        if (pd>0.3) { 
  • sans/XOP_Dev/SANSAnalysis/lib/libStructureFactor.c

    r356 r632  
    146146        eta2 = eta*eta; 
    147147        eta3 = eta*eta2; 
    148         eta4 = eta*eta3;       
    149         etam1 = 1. - eta;  
     148        eta4 = eta*eta3; 
     149        etam1 = 1. - eta; 
    150150        etam14 = etam1*etam1*etam1*etam1; 
    151151        alpha = (  pow((1. + 2.*eta),2) + eta3*( eta-4.0 )  )/etam14; 
     
    187187        pi = 4.0*atan(1.); 
    188188        QQ= q; 
    189          
    190         diam=dp[0];             //in   
     189 
     190        diam=dp[0];             //in A 
    191191        zz = dp[1];             //# of charges 
    192         VolFrac=dp[2];   
     192        VolFrac=dp[2]; 
    193193        Temp=dp[3];             //in degrees Kelvin 
    194194        csalt=dp[4];            //in molarity 
     
    391391        d = 1.0-reta; 
    392392        d2 = d*d; 
    393         dak = d/rak;                                                   
     393        dak = d/rak; 
    394394        dd2 = 1.0/d2; 
    395395        dd4 = dd2*dd2; 
     
    659659         
    660660        pi = 4.0*atan(1.0); 
    661          
     661        if (rcyl == 0 || hcyl == 0) { 
     662                return 0.0; 
     663        } 
    662664        a = rcyl; 
    663665        b = hcyl/2.0; 
     
    682684         
    683685        double ee,e1,bd,b1,bL,b2,del,ddd,diam; 
    684          
     686 
     687        if (aa == 0 || bb == 0) { 
     688                return 0.0; 
     689        } 
     690        if (aa == bb) { 
     691                return 2.0*aa; 
     692        } 
    685693        if(aa>bb) { 
    686694                ee = (aa*aa - bb*bb)/(aa*aa); 
  • sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c

    r554 r632  
    417417         
    418418        uval = Mw_Mn - 1.0; 
    419         if(uval == 0) { 
     419        if(uval == 0.0) { 
    420420                uval = 1e-6;            //avoid divide by zero error 
    421421        } 
     
    453453        bgd = dp[4];  
    454454         
    455         inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1 + (x*ll)*(x*ll)) + bgd; 
     455        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1.0 + (x*ll)*(x*ll)) + bgd; 
    456456         
    457457        return(inten); 
Note: See TracChangeset for help on using the changeset viewer.