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)

File:
1 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 { 
Note: See TracChangeset for help on using the changeset viewer.