Changeset 632 for sans/XOP_Dev/SANSAnalysis

Timestamp:

Mar 1, 2010 1:57:12 PM (13 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:

• libCylinder.c (modified) (60 diffs)
• libSphere.c (modified) (30 diffs)
• libStructureFactor.c (modified) (5 diffs)
• libTwoPhase.c (modified) (2 diffs)

sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c

@@ -28 +28 @@
         
         Pi = 4.0*atan(1.0);
-        lolim = 0;
+        lolim = 0.0;
         uplim = Pi/2.0;
         
@@ -84 +84 @@
         double va,vb;           //upper and lower integration limits
         double summ,zi,yyy,answer,vell;                 //running tally of integration
-        double summj,vaj,vbj,zij,arg;                   //for the inner integration
-        
+        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration
+
         Pi = 4.0*atan(1.0);
-        va = 0;
-        vb = 1;         //orintational average, outer integral
-        vaj=0;
+        va = 0.0;
+        vb = 1.0;               //orintational average, outer integral
+        vaj=0.0;
         vbj=Pi;         //endpoints of inner integral
         
@@ -107 +107 @@
                 summj=0;
                 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
-                arg = ra*sqrt(1-zi*zi);
+                arg = ra*sqrt(1.0-zi*zi);
                 for(j=0;j<nord;j++) {
                         //76 gauss points for the inner integral as well
@@ -120 +120 @@
                 
                 //now calculate outer integral
-                arg = q*length*zi/2;
-                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+                arg = q*length*zi/2.0;
+                if (arg == 0.0){
+                        si = 1.0;
+                }else{
+                        si = sin(arg) * sin(arg) / arg / arg;
+                }
+                yyy = Gauss76Wt[i] * answer * si;
                 summ += yyy;
         }
@@ -162 +167 @@
         double va,vb;           //upper and lower integration limits
         double summ,zi,yyy,answer,vell;                 //running tally of integration
-        double summj,vaj,vbj,zij,arg;                   //for the inner integration
-        
+        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration
+
         Pi = 4.0*atan(1.0);
-        va = 0;
-        vb = 1;         //orintational average, outer integral
-        vaj=0;
+        va = 0.0;
+        vb = 1.0;               //orintational average, outer integral
+        vaj=0.0;
         vbj=Pi;         //endpoints of inner integral
         
@@ -185 +190 @@
                 summj=0;
                 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
-                arg = ra*sqrt(1-zi*zi);
+                arg = ra*sqrt(1.0-zi*zi);
                 for(j=0;j<nordj;j++) {
                         //20 gauss points for the inner integral
@@ -198 +203 @@
                 
                 //now calculate outer integral
-                arg = q*length*zi/2;
-                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+                arg = q*length*zi/2.0;
+                if (arg == 0.0){
+                        si = 1.0;
+                }else{
+                        si = sin(arg) * sin(arg) / arg / arg;
+                }
+                yyy = Gauss76Wt[i] * answer * si;
                 summ += yyy;
         }
@@ -215 +225 @@
         answer *= scale;
         // add in the background
-        answer += bkg;  
-        
+        answer += bkg;
+
         return answer;
 }
@@ -316 +326 @@
         
         //      Pi = 4.0*atan(1.0);
-        va = 0;
-        vb = 1;         //orintational average, outer integral
-        vaj = 0;
-        vbj = 1;                //endpoints of inner integral
-        
+        va = 0.0;
+        vb = 1.0;               //orintational average, outer integral
+        vaj = 0.0;
+        vbj = 1.0;              //endpoints of inner integral
+
         summ = 0.0;                     //initialize intergral
         
@@ -340 +350 @@
         for(i=0;i<nordi;i++) {
                 //setup inner integral over the ellipsoidal cross-section
-                summj=0;
+                summj=0.0;
                 sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
                 
@@ -346 +356 @@
                         //76 gauss points for the inner integral
                         uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
-                        mudum = mu*sqrt(1-sigma*sigma);
+                        mudum = mu*sqrt(1.0-sigma*sigma);
                         yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
                         summj += yyy;
@@ -352 +362 @@
                 //now calculate the value of the inner integral
                 answer = (vbj-vaj)/2.0*summj;
-                
-                arg = mu*cc*sigma/2;
-                if ( arg == 0 ) {
-                        answer *= 1;
+
+                arg = mu*cc*sigma/2.0;
+                if ( arg == 0.0 ) {
+                        answer *= 1.0;
                 } else {
                         answer *= sin(arg)*sin(arg)/arg/arg;
@@ -401 +411 @@
         
         pi = 4.0*atan(1.0);
-        va = 0;
-        vb = 1;         //limits of numerical integral
-        
+        va = 0.0;
+        vb = 1.0;               //limits of numerical integral
+
         summ = 0.0;                     //initialize intergral
         
@@ -456 +466 @@
         
         pi = 4.0*atan(1.0);
-        va = 0;
-        vb = 1;         //limits of numerical integral
-        
+        va = 0.0;
+        vb = 1.0;               //limits of numerical integral
+
         summ = 0.0;                     //initialize intergral
         
@@ -478 +488 @@
         answer *= delrho*delrho;
         //normalize by volume
-        answer *= 4*pi/3*a*a*nua;
+        answer *= 4.0*pi/3.0*a*a*nua;
         //convert to [cm-1]
         answer *= 1.0e8;
@@ -521 +531 @@
         
         lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
-        if(lolim<0) {
-                lolim = 0;
+        if(lolim<0.0) {
+                lolim = 0.0;
         }
         if(pd>0.3) {
@@ -582 +592 @@
         
         lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
-        if(lolim<0) {
-                lolim = 0;
+        if(lolim<0.0) {
+                lolim = 0.0;
         }
         if(pd>0.3) {
@@ -697 +707 @@
         
         lolim = exp(log(radius)-(4.*sigma));
-        if (lolim<0) {
-                lolim=0;                //to avoid numerical error when  va<0 (-ve r value)
+        if (lolim<0.0) {
+                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
         }
         uplim = exp(log(radius)+(4.*sigma));
@@ -758 +768 @@
         delpc = sldc - slds;    //core - shell
         delps = slds - sld;             //shell - solvent
-        bkg = dp[8];                        
-        
+        bkg = dp[8];
+
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -768 +778 @@
         answer = (uplim-lolim)/2.0*summ;
         // normalize by particle volume
-        oblatevol = 4*Pi/3*trmaj*trmaj*trmin;
+        oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin;
         answer /= oblatevol;
         
@@ -812 +822 @@
         delpc = sldc - slds;            //core - shell
         delps = slds - sld;             //shell  - sovent
-        bkg = dp[8];                        
-        
+        bkg = dp[8];
+
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
@@ -822 +832 @@
         answer = (uplim-lolim)/2.0*summ;
         // normalize by particle volume
-        prolatevol = 4*Pi/3*trmaj*trmin*trmin;
+        prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin;
         answer /= prolatevol;
         
@@ -1232 +1242 @@
         
         lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
-        if(lolim<0) {
-                lolim = 0;
+        if(lolim<0.0) {
+                lolim = 0.0;
         }
         if(pd>0.3) {
@@ -1296 +1306 @@
         
         lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
-        if(lolim<0) {
-                lolim = 0;
+        if(lolim<0.0) {
+                lolim = 0.0;
         }
         if(pd>0.3) {
@@ -1570 +1580 @@
     q02 = q0*q0;
     q0p = pow(q0,(-4.0 + p1short) );
-    
-    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)))))); 
-        
+
+    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))))));
+
     return(yy);
 }
@@ -1610 +1620 @@
     
     t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
-    
-    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;
-    
-    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);
-    
-    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);
-    
-    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);
-    
-    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);
-    
-    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));
-    
-    t8 = ((1 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-    
-    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;
-    
-    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);
-        
-    
-    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))))))))));
-        
+
+    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;
+
+    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);
+
+    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);
+
+    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);
+
+    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);
+
+    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));
+
+    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
+
+    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;
+
+    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);
+
+
+    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))))))))));
+
     return (yy);
 }
@@ -1728 +1738 @@
 {
         // local variables
-        double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi;
-        
+        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;
+
         Pi = 4.0*atan(1.0);
         
@@ -1741 +1751 @@
         vc = pi43*aa*bb*bb;
         vt = pi43*trmaj*trmin*trmin;
-        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
-        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
+        if (uq == 0.0){
+                siq = 1.0/3.0;
+        }else{
+                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+        }
+        if (ut == 0.0){
+                sit = 1.0/3.0;
+        }else{
+                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
+        }
+        gfnc = 3.0*siq*vc*delpc;
+        gfnt = 3.0*sit*vt*delps;
         gfn = gfnc+gfnt;
         gfn2 = gfn*gfn;
@@ -1760 +1780 @@
 {
         // local variables
-        double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
-        
+        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;
+
         Pi = 4.0*atan(1.0);
         pi43=4.0/3.0*Pi;
@@ -1772 +1792 @@
         vc = pi43*aa*aa*bb;
         vt = pi43*trmaj*trmaj*trmin;
-        gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc;
-        gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps;
+        if (uq == 0.0){
+                siq = 1.0/3.0;
+        }else{
+                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+        }
+        if (ut == 0.0){
+                sit = 1.0/3.0;
+        }else{
+                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
+        }
+        gfnc = 3.0*siq*vc*delpc;
+        gfnt = 3.0*sit*vt*delps;
         tgfn = gfnc+gfnt;
         gfn4 = tgfn*tgfn;
@@ -1790 +1820 @@
     
     Pq = Sk_WR(q,zi,lb);                //does not have cross section term
-    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
-    
+    if (qr !=0){
+        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
+    } 
     vcyl=Pi*radius*radius*zi;
     Pq *= vcyl*vcyl;
@@ -1810 +1841 @@
     
     Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
-    Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
-    
+    if (qr !=0){
+        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
+    }
+
     vcyl=Pi*zi*zi*Lc;
     Pq *= vcyl*vcyl;
@@ -1830 +1863 @@
         Pi = 4.0*atan(1.0);
         nord = 76;
-        lolim = 0;
-        uplim = Pi/2;
+        lolim = 0.0;
+        uplim = Pi/2.0;
         halfheight = length/2.0;
         
@@ -1856 +1889 @@
         // length is the *Half* CORE-LENGTH of the cylinder 
         // dum is the dummy variable for the integration (theta)
-        
-        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
+
+        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
         double Pi;
         
@@ -1871 +1904 @@
         sinarg1 = qq*length*cos(dum);
         sinarg2 = qq*(length+facthick)*cos(dum);
-        
-        t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
-        t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
-        
+        if (besarg1 == 0.0){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0.0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0.0){
+                si1 = 1.0;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0.0){
+                si2 = 1.0;
+        }else{
+                si2 = sin(sinarg2)/sinarg2;
+        }
+
+        t1 = 2.0*vol1*dr1*si1*be1;
+        t2 = 2.0*vol2*dr2*si2*be2;
+
         retval = ((t1+t2)*(t1+t2))*sin(dum);
         return (retval);
@@ -1884 +1937 @@
 CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
 {
-        
-    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval;
+
+    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
     double Pi;
     
@@ -1899 +1952 @@
     sinarg1 = qq*length*cos(dum);
     sinarg2 = qq*(length+thick)*cos(dum);
-    
-    t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
-    t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
-    
+
+        if (besarg1 == 0.0){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0.0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0.0){
+                si1 = 1.0;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0.0){
+                si2 = 1.0;
+        }else{
+                si2 = sin(sinarg2)/sinarg2;
+        }
+
+    t1 = 2.0*vol1*dr1*si1*be1;
+    t2 = 2.0*vol2*dr2*si2*be2;
+
     retval = ((t1+t2)*(t1+t2))*sin(dum);
         
@@ -1917 +1991 @@
     
     Pi = 4.0*atan(1.0);
-    lolim = 0;
+    lolim = 0.0;
     uplim = Pi/2.0;
     halfheight = dumLen/2.0;
@@ -1950 +2024 @@
         // length is the *Half* CORE-LENGTH of the cylinder = L (A)
         // dum is the dummy variable for the integration (x in Feigin's notation)
-        
-        //Local variables 
-        double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
+
+        //Local variables
+        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
         double Pi;
         int kk;
@@ -1968 +2042 @@
         sinarg1 = qq*length*cos(dum);
         sinarg2 = qq*(length+thick)*cos(dum);
-        
-        t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1);
-        t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2);
-        
+
+        if (besarg1 == 0.0){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0.0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0.0){
+                si1 = 1.0;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0.0){
+                si2 = 1.0;
+        }else{
+                si2 = sin(sinarg2)/sinarg2;
+        }
+
+        t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1);
+        t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2);
+
         retval =((t1+t2)*(t1+t2))*sin(dum);
         
@@ -1998 +2093 @@
     
     Pi = 4.0*atan(1.0);
-    lolim = 0;
+    lolim = 0.0;
     uplim = Pi/2.0;
     halfheight = length/2.0;
@@ -2061 +2156 @@
         
     nu = nua/a;
-    arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
-    
-    retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0));
+    if (arg == 0.0){
+        retval =1.0/3.0;
+    }else{
+        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    }
     retval *= retval;
-    retval *= 9;
-    
+    retval *= 9.0;
+
     return(retval);
 }//Function EllipsoidKernel()
@@ -2076 +2174 @@
     
     gamma = rcore/rshell;
-    arg1 = qq*rshell*sqrt(1-dum*dum);           //1=shell (outer radius)
-    arg2 = qq*rcore*sqrt(1-dum*dum);                    //2=core (inner radius)
-    lam1 = 2*NR_BessJ1(arg1)/arg1;
-    lam2 = 2*NR_BessJ1(arg2)/arg2;
-    psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);         //SRK 10/19/00
-    
-    sinarg = qq*length*dum/2;
-    t2 = sin(sinarg)/sinarg;
-    
+    arg1 = qq*rshell*sqrt(1.0-dum*dum);         //1=shell (outer radius)
+    arg2 = qq*rcore*sqrt(1.0-dum*dum);                  //2=core (inner radius)
+    if (arg1 == 0.0){
+        lam1 = 1.0;
+    }else{
+        lam1 = 2.0*NR_BessJ1(arg1)/arg1;
+    }
+    if (arg2 == 0.0){
+        lam2 = 1.0;
+    }else{
+        lam2 = 2.0*NR_BessJ1(arg2)/arg2;
+    }
+    //Todo: Need to check psi behavior as gamma goes to 1.
+    psi = 1.0/(1.0-gamma*gamma)*(lam1 -  gamma*gamma*lam2);             //SRK 10/19/00
+    sinarg = qq*length*dum/2.0;
+    if (sinarg == 0.0){
+        t2 = 1.0;
+    }else{
+        t2 = sin(sinarg)/sinarg;
+    }
+
     retval = psi*psi*t2*t2;
     
@@ -2099 +2209 @@
     
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
-    arg1 = (mu/2)*cos(Pi*uu/2);
-    arg2 = (mu*aa/2)*sin(Pi*uu/2);
-    if(arg1==0) {
-                tmp1 = 1;
+    arg1 = (mu/2.0)*cos(Pi*uu/2.0);
+    arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0);
+    if(arg1==0.0) {
+                tmp1 = 1.0;
         } else {
                 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
     }
-    
-    if (arg2==0) {
-                tmp2 = 1;
+
+    if (arg2==0.0) {
+                tmp2 = 1.0;
         } else {
                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
@@ -2130 +2240 @@
     arg += cc*cc*dy*dy;
     arg = q*sqrt(arg);
-    val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
-    
+    if (arg == 0.0){
+        val = 1.0; // as arg --> 0, val should go to 1.0
+    }else{
+        val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
+    }
     return (val);
         
@@ -2144 +2257 @@
         // rr is the radius of the cylinder (A)
         // h is the HALF-LENGTH of the cylinder = L/2 (A)
-        
-    double besarg,bj,retval,d1,t1,b1,t2,b2;              //Local variables 
-        
-        
+
+    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables
+
+
     besarg = qq*rr*sin(theta);
-        
+    siarg = qq * h * cos(theta);
     bj =NR_BessJ1(besarg);
         
         //* Computing 2nd power */
-    d1 = sin(qq * h * cos(theta));
+    d1 = sin(siarg);
     t1 = d1 * d1;
         //* Computing 2nd power */
@@ -2159 +2272 @@
     t2 = d1 * d1 * 4.0 * sin(theta);
         //* Computing 2nd power */
-    d1 = qq * h * cos(theta);
+    d1 = siarg;
     b1 = d1 * d1;
         //* Computing 2nd power */
     d1 = qq * rr * sin(theta);
     b2 = d1 * d1;
-    retval = t1 * t2 / b1 / b2;
-        
+    if (besarg == 0.0){
+        be = sin(theta);
+    }else{
+        be = t2 / b2;
+    }
+    if (siarg == 0.0){
+        si = 1.0;
+    }else{
+        si = t1 / b1;
+    }
+    retval = be * si;
+
     return (retval);
         
@@ -2181 +2304 @@
         double retval,arg;               //Local variables 
         
-        arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
-        
-        retval = 2*NR_BessJ1(arg)/arg;
-        
+        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2);
+        if (arg == 0.0){
+                retval = 1.0;
+        }else{
+                retval = 2.0*NR_BessJ1(arg)/arg;
+        }
+
         //square it
         retval *= retval;
@@ -2381 +2507 @@
                 yyy = inner;
                 
-                if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h
+                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
                 } else {
@@ -2500 +2626 @@
                 yyy = inner;
                 
-                if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h
+                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
                 } else {
@@ -2584 +2710 @@
                 yyy = inner;
                 
-                if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h
+                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
                 } else {
@@ -2705 +2831 @@
                 yyy = inner;
                 
-                if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h
+                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
                 } else {
@@ -2789 +2915 @@
                 yyy = inner;
                 
-                if(arg1 == 0) {         //limiting value of sinc(0) is 1; sinc is not defined in math.h
+                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2;
                 } else {
@@ -2865 +2991 @@
         
         lolim = exp(log(radius)-(4.*sigma));
-        if (lolim<0) {
-                lolim=0;                //to avoid numerical error when  va<0 (-ve r value)
+        if (lolim<0.0) {
+                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
         }
         uplim = exp(log(radius)+(4.*sigma));
@@ -2908 +3034 @@
         Pi = 4.0*atan(1.0);
         nord = 76;
-        lolim = 0;
+        lolim = 0.0;
         uplim = Pi/2;
         halfheight = length/2.0;
@@ -2940 +3066 @@
         dr2 = rhor-rhosolv;
         dr3=  rhoh-rhor;
-        vol1 = Pi*rad*rad*(2*length);
-        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2*length+2*facthick);
-        vol3= Pi*(rad)*(rad)*(2*length+2*facthick);
+        vol1 = Pi*rad*rad*(2.0*length);
+        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
+        vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick);
         besarg1 = qq*rad*sin(dum);
         besarg2 = qq*(rad+radthick)*sin(dum);
@@ -2948 +3074 @@
         sinarg2 = qq*(length+facthick)*cos(dum);
         
-        t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
-        t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
-        t3 = 2*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1;
+        t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1;
+        t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2;
+        t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1;
         
         retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
@@ -2998 +3124 @@
         arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0);
                          
-        if(arg1==0){
+        if(arg1==0.0){
                 t1 = 1.0;
         } else {
                 t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
         }
-        if(arg2==0){
+        if(arg2==0.0){
                 t2 = 1.0;
         } else {
                 t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)    
         }       
-        if(arg3==0){
+        if(arg3==0.0){
                 t3 = 1.0;
         } else {
                 t3 = sin(arg3)/arg3;
         }
-        if(arg4==0){
+        if(arg4==0.0){
                 t4 = 1.0;
         } else {
@@ -3073 +3199 @@
         for(i=0;i<nordi;i++) {
                 //setup inner integral over the ellipsoidal cross-section
-                summj=0;
+                summj=0.0;
                 sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
                 
@@ -3088 +3214 @@
                 //finish the outer integral cc already scaled
                 arg = mu*cc*sigma/2.0;
-                if ( arg == 0 ) {
+                if ( arg == 0.0 ) {
                         answer *= 1.0;
                 } else {

sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

@@ -14 +14 @@
 {
         double scale,radius,delrho,bkg,sldSph,sldSolv;          //my local names
-        double bes,f,vol,f2,pi;
-        
+        double bes,f,vol,f2,pi,qr;
+
         pi = 4.0*atan(1.0);
         scale = dp[0];
@@ -24 +24 @@
         
         delrho = sldSph - sldSolv;
-        //handle q==0 separately
-        if(q==0){
-                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
-                return(f);
+        //handle qr==0 separately
+        qr = q*radius;
+        if(qr == 0.0){
+                bes = 1.0;
+        }else{
+                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
         }
         
-        bes = 3.0*(sin(q*radius)-q*radius*cos(q*radius))/(q*q*q)/(radius*radius*radius);
         vol = 4.0*pi/3.0*radius*radius*radius;
-        f = vol*bes*delrho;             // [=] 
+        f = vol*bes*delrho;             // [=] A-1
                                                         // normalize to single particle volume, convert to 1/cm
         f2 = f * f / vol * 1.0e8;               // [=] 1/cm
@@ -60 +61 @@
         qr=x*rcore;
         contr = rhocore-rhoshel;
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        if(qr == 0.0){
+                bes = 1.0;
+        }else{
+                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        }
         vol = 4.0*pi/3.0*rcore*rcore*rcore;
         f = vol*bes*contr;
@@ -66 +71 @@
         qr=x*(rcore+thick);
         contr = rhoshel-rhosolv;
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        if(qr == 0.0){
+                bes = 1.0;
+        }else{
+                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        }
         vol = 4.0*pi/3.0*pow((rcore+thick),3);
         f += vol*bes*contr;
-        
-        // normalize to particle volume and rescale from [-1] to [cm-1]
+
+        // normalize to particle volume and rescale from [A-1] to [cm-1]
         f2 = f*f/vol*1.0e8;
         
@@ -102 +111 @@
         qr=x*rcore;
         contr = rhocore-rhoshel;
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        if(qr == 0){
+                bes = 1.0;
+        }else{
+                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        }
         vol = 4.0*pi/3.0*rcore*rcore*rcore;
         f = vol*bes*contr;
@@ -108 +121 @@
         qr=x*(rcore+thick);
         contr = rhoshel-rhosolv;
-        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        if(qr == 0.0){
+                bes = 1.0;
+        }else{
+                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
+        }
         vol = 4.0*pi/3.0*pow((rcore+thick),3);
         f += vol*bes*contr;
-        
-        // normalize to the particle volume and rescale from [-1] to [cm-1]
+
+        // normalize to the particle volume and rescale from [A-1] to [cm-1]
         //note that for the vesicle model, the volume is ONLY the shell volume
         vol = 4.0*pi/3.0*(pow((rcore+thick),3)-pow(rcore,3));
@@ -388 +405 @@
         
         //small QR limit - use Guinier approx
-        zp8 = zz+8;
+        zp8 = zz+8.0;
         if(x*ravg < 0.1) {
                 i_zero = scale*delrho*delrho*1.e8*4.*pi/3.*pow(ravg,3);
@@ -399 +416 @@
         //
 
-        aa = (zz+1)/x/ravg;
+        aa = (zz+1.0)/x/ravg;
         
         at1 = atan(1.0/aa);
@@ -420 +437 @@
         
         pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
-        pq = pow(10,pq)*8*pi*pi*delrho*delrho;
+        pq = pow(10,pq)*8.0*pi*pi*delrho*delrho;
         
         //
@@ -483 +500 @@
         // so just use the limiting Guiner value
         if(qr<0.1) {
-                h1 = scale*cont*cont*1.e8*4.*pi/3*pow(rad,3);
+                h1 = scale*cont*cont*1.e8*4.*pi/3.0*pow(rad,3);
                 h1 *=   (1. + 15.*pow(pd,2) + 27.*pow(pd,4) +27./7.*pow(pd,6) );                                //6th moment
                 h1 /= (1.+3.*pd*pd);    //3rd moment
@@ -495 +512 @@
         // normal calculation
         h1 = -0.5*qw + qr*qr*qw + (qw*qw*qw)/3.0;
-        h1 -= 5.0/2.0*cos(2*qr)*sin(qw)*cos(qw);
-        h1 += 0.5*qr*qr*cos(2*qr)*sin(2*qw);
-        h1 += 0.5*qw*qw*cos(2*qr)*sin(2*qw);
-        h1 += qw*qr*sin(2*qr)*cos(2*qw);
+        h1 -= 5.0/2.0*cos(2.0*qr)*sin(qw)*cos(qw);
+        h1 += 0.5*qr*qr*cos(2.0*qr)*sin(2.0*qw);
+        h1 += 0.5*qw*qw*cos(2.0*qr)*sin(2.0*qw);
+        h1 += qw*qr*sin(2.0*qr)*cos(2.0*qw);
         h1 += 3.0*qw*(cos(qr)*cos(qw))*(cos(qr)*cos(qw));
         h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
@@ -550 +567 @@
         
         va = -4.0*sig + rad;
-        if (va<0) {
-                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
+        if (va<0.0) {
+                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value)
         }
         vb = 4.0*sig +rad;
@@ -604 +621 @@
         va = -3.5*sig + mu;
         va = exp(va);
-        if (va<0) {
-                va=0;           //to avoid numerical error when  va<0 (-ve q-value)
+        if (va<0.0) {
+                va=0.0;         //to avoid numerical error when  va<0 (-ve q-value)
         }
         vb = 3.5*sig*(1.0+sig) +mu;
@@ -642 +659 @@
         
         pi = 4.0*atan(1.0);
-        retval = (1/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
+        retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig );
         return(retval);
 }
@@ -791 +808 @@
         double psf11,psf12,psf22;
         double phi1,phi2,phr,a3;
-        double v1,v2,n1,n2,qr1,qr2,b1,b2;
+        double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2;
         int err;
         
@@ -825 +842 @@
         qr1 = r1*x;
         qr2 = r2*x;
-        
-        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*(sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
-        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*(sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+
+        if (qr1 == 0){
+                sc1 = 1.0/3.0;
+        }else{
+                sc1 = (sin(qr1)-qr1*cos(qr1))/qr1/qr1/qr1;
+        }
+        if (qr2 == 0){
+                sc2 = 1.0/3.0;
+        }else{
+                sc2 = (sin(qr2)-qr2*cos(qr2))/qr2/qr2/qr2;
+        }
+        b1 = r1*r1*r1*(rho1-rhos)*4.0*pi*sc1;
+        b2 = r2*r2*r2*(rho2-rhos)*4.0*pi*sc2;
         inten = n1*b1*b1*psf11;
         inten += sqrt(n1*n2)*2.0*b1*b2*psf12;
@@ -941 +968 @@
         double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
         double a1,a2i,a2,b1,b2,b12,gm1,gm12;
-        double err=0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
+        double err=0.0,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3;
         double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
         double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
@@ -1059 +1086 @@
         
         ii=0;
-        fval=0;
+        fval=0.0;
         
         do {
                 ri = rcore + (double)ii*ts + (double)ii*tw;
-                voli = 4*pi/3*ri*ri*ri;
+                voli = 4.0*pi/3.0*ri*ri*ri;
                 sldi = rhocore-rhoshel;
                 fval += voli*sldi*F_func(ri*x);
                 ri += ts;
-                voli = 4*pi/3*ri*ri*ri;
+                voli = 4.0*pi/3.0*ri*ri*ri;
                 sldi = rhoshel-rhocore;
                 fval += voli*sldi*F_func(ri*x);
@@ -1075 +1102 @@
         fval *= fval;           //square it
         fval /= voli;           //normalize by the overall volume
-        fval *= scale*1e8;
+        fval *= scale*1.0e8;
         fval += bkg;
         
@@ -1132 +1159 @@
         
         ii=minPairs;
-        fval=0;
-        d1 = 0;
-        d2 = 0;
-        r1 = 0;
-        r2 = 0;
-        distr = 0;
-        first = 1;
+        fval=0.0;
+        d1 = 0.0;
+        d2 = 0.0;
+        r1 = 0.0;
+        r2 = 0.0;
+        distr = 0.0;
+        first = 1.0;
         do {
                 //make the current values old
@@ -1156 +1183 @@
         } while(ii<=maxPairs);
         
-        fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume
+        fval /= 4.0*pi/3.0*avg*avg*avg;         //normalize by the overall volume
         fval /= distr;
         fval *= scale;
@@ -1172 +1199 @@
         pi = 4.0*atan(1.0);
     ii=0;
-    fval=0;
+    fval=0.0;
     
     do {
         ri = rcore + (double)ii*ts + (double)ii*tw;
-        voli = 4*pi/3*ri*ri*ri;
+        voli = 4.0*pi/3.0*ri*ri*ri;
         sldi = rhocore-rhoshel;
         fval += voli*sldi*F_func(ri*x);
         ri += ts;
-        voli = 4*pi/3*ri*ri*ri;
+        voli = 4.0*pi/3.0*ri*ri*ri;
         sldi = rhoshel-rhocore;
         fval += voli*sldi*F_func(ri*x);
@@ -1188 +1215 @@
     fval *= fval;
     fval /= voli;
-    fval *= 1e8;
+    fval *= 1.0e8;
     
     return(fval);       // this result still needs to be multiplied by scale and have background added
@@ -1198 +1225 @@
     double dr;
     
-    dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
+    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0));
     return (exp(dr));
 }
@@ -1221 +1248 @@
 double
 F_func(double qr) {
-        return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
+        double sc;
+        if (qr == 0.0){
+                sc = 1.0;
+        }else{
+                sc=(3.0*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
+        }
+        return sc;
 }
 
@@ -1310 +1343 @@
         contr = rhocore-rhoshel1;
         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-        vol = 4.0*pi/3*rcore*rcore*rcore;
+        vol = 4.0*pi/3.0*rcore*rcore*rcore;
         f = vol*bes*contr;
         //now the shell (1)
@@ -1377 +1410 @@
         contr = rhocore-rhoshel1;
         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-        vol = 4.0*pi/3*rcore*rcore*rcore;
+        vol = 4.0*pi/3.0*rcore*rcore*rcore;
         f = vol*bes*contr;
         //now the shell (1)
@@ -1453 +1486 @@
         contr = rhocore-rhoshel1;
         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr);
-        vol = 4.0*pi/3*rcore*rcore*rcore;
+        vol = 4.0*pi/3.0*rcore*rcore*rcore;
         f = vol*bes*contr;
         //now the shell (1)
@@ -1515 +1548 @@
         range = 8.0;                    //std deviations for the integration
         va = rcore*(1.0-range*pd);
-        if (va<0) {
-                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
+        if (va<0.0) {
+                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
         }
         if (pd>0.3) {
@@ -1586 +1619 @@
         range = 8.0;                    //std deviations for the integration
         va = rcore*(1.0-range*pd);
-        if (va<0) {
-                va=0;           //otherwise numerical error when pd >= 0.3, making a<0
+        if (va<0.0) {
+                va=0.0;         //otherwise numerical error when pd >= 0.3, making a<0
         }
         if (pd>0.3) {

sans/XOP_Dev/SANSAnalysis/lib/libStructureFactor.c

@@ -146 +146 @@
         eta2 = eta*eta;
         eta3 = eta*eta2;
-        eta4 = eta*eta3;      
-        etam1 = 1. - eta; 
+        eta4 = eta*eta3;
+        etam1 = 1. - eta;
         etam14 = etam1*etam1*etam1*etam1;
         alpha = (  pow((1. + 2.*eta),2) + eta3*( eta-4.0 )  )/etam14;
@@ -187 +187 @@
         pi = 4.0*atan(1.);
         QQ= q;
-        
-        diam=dp[0];             //in  
+
+        diam=dp[0];             //in A
         zz = dp[1];             //# of charges
-        VolFrac=dp[2];  
+        VolFrac=dp[2];
         Temp=dp[3];             //in degrees Kelvin
         csalt=dp[4];            //in molarity
@@ -391 +391 @@
         d = 1.0-reta;
         d2 = d*d;
-        dak = d/rak;                                                  
+        dak = d/rak;
         dd2 = 1.0/d2;
         dd4 = dd2*dd2;
@@ -659 +659 @@
         
         pi = 4.0*atan(1.0);
-        
+        if (rcyl == 0 || hcyl == 0) {
+                return 0.0;
+        }
         a = rcyl;
         b = hcyl/2.0;
@@ -682 +684 @@
         
         double ee,e1,bd,b1,bL,b2,del,ddd,diam;
-        
+
+        if (aa == 0 || bb == 0) {
+                return 0.0;
+        }
+        if (aa == bb) {
+                return 2.0*aa;
+        }
         if(aa>bb) {
                 ee = (aa*aa - bb*bb)/(aa*aa);

sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c

@@ -417 +417 @@
         
         uval = Mw_Mn - 1.0;
-        if(uval == 0) {
+        if(uval == 0.0) {
                 uval = 1e-6;            //avoid divide by zero error
         }
@@ -453 +453 @@
         bgd = dp[4]; 
         
-        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1 + (x*ll)*(x*ll)) + bgd;
+        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1.0 + (x*ll)*(x*ll)) + bgd;
         
         return(inten);