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