Changeset 975ec8e in sasview for sansmodels/src/libigor

Timestamp:

Aug 14, 2009 5:58:58 PM (15 years ago)

Author:

Jae Cho <jhjcho@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

d5bd424

Parents:

cad821b

Message:

working on 2D models. Still need smore corrections and unit tests.

Location:

sansmodels/src/libigor

Files:

• libCylinder.c (modified) (24 diffs)
• libCylinder.h (modified) (1 diff)
• libSphere.c (modified) (79 diffs)

sansmodels/src/libigor/libCylinder.c

@@ -80 +80 @@
         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);
@@ -114 +114 @@
                 //now calculate outer integral
                 arg = q*length*zi/2;
-                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+                if (arg == 0){
+                        si = 1;
+                }else{
+                        si = sin(arg) * sin(arg) / arg / arg;
+                }
+                yyy = Gauss76Wt[i] * answer * si;
                 summ += yyy;
         }
@@ -155 +160 @@
         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);
@@ -189 +194 @@
                 //now calculate outer integral
                 arg = q*length*zi/2;
-                yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg;
+                if (arg == 0){
+                        si = 1;
+                }else{
+                        si = sin(arg) * sin(arg) / arg / arg;
+                }
+                yyy = Gauss76Wt[i] * answer * si;
                 summ += yyy;
         }
@@ -893 +903 @@
         return(inten+bkg);
 }
-
+/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
+                                                -NO polydispersion included
+*/
+double lamellar_kernel(double dp[], double q){
+        double scale,del,sld_bi,sld_sol,contr,bkg;              //local variables of coefficient wave
+        double inten, qval,Pq;
+        double Pi;
+
+
+        Pi = 4.0*atan(1.0);
+        scale = dp[0];
+        del = dp[1];
+        sld_bi = dp[2];
+        sld_sol = dp[3];
+        bkg = dp[4];
+        qval = q;
+        contr = sld_bi -sld_sol;
+
+        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless...
+
+        inten /= del;                   //normalize by the thickness (in A)
+
+        inten *= 1.0e8;         // 1/A to 1/cm
+
+        return(inten+bkg);
+}
 /*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
 -------
@@ -1682 +1719 @@
 {
         // 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);
@@ -1695 +1732 @@
         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){
+                siq = 1/3;
+        }else{
+                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+        }
+        if (ut == 0){
+                sit = 1/3;
+        }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;
@@ -1714 +1761 @@
 {
         // 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);
@@ -1726 +1773 @@
         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){
+                siq = 1/3;
+        }else{
+                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
+        }
+        if (ut == 0){
+                sit = 1/3;
+        }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;
@@ -1744 +1801 @@
 
     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);
+    } //else Pk *=1;
     vcyl=Pi*radius*radius*zi;
     Pq *= vcyl*vcyl;
@@ -1764 +1822 @@
 
     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;
@@ -1811 +1871 @@
         // 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;
 
@@ -1825 +1885 @@
         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){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0){
+                si1 = 1;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0){
+                si2 = 1;
+        }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);
@@ -1839 +1919 @@
 {
 
-    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;
 
@@ -1854 +1934 @@
     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){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0){
+                si1 = 1;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0){
+                si2 = 1;
+        }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);
@@ -1906 +2007 @@
 
         //Local variables
-        double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
+        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
         double Pi;
         int kk;
@@ -1923 +2024 @@
         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){
+                be1 = 0.5;
+        }else{
+                be1 = NR_BessJ1(besarg1)/besarg1;
+        }
+        if (besarg2 == 0){
+                be2 = 0.5;
+        }else{
+                be2 = NR_BessJ1(besarg2)/besarg2;
+        }
+        if (sinarg1 == 0){
+                si1 = 1;
+        }else{
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if (besarg2 == 0){
+                si2 = 1;
+        }else{
+                si2 = sin(sinarg2)/sinarg2;
+        }
+
+        t1 = 2*area*(2*length)*dr1*(si1)*(be1);
+        t2 = 2*area*dr2*(totald*si2-2*length*si1)*(be2);
 
         retval =((t1+t2)*(t1+t2))*sin(dum);
@@ -1937 +2059 @@
         // end of loop for S(q)
         sqq=1.0+2.0*sqq/N;
+
         retval *= sqq;
 
@@ -2016 +2139 @@
     nu = nua/a;
     arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
-
-    retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    if (arg == 0){
+        retval =1/3;
+    }else{
+        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
+    }
     retval *= retval;
     retval *= 9;
@@ -2032 +2158 @@
     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;
+    if (arg1 == 0){
+        lam1 = 1;
+    }else{
+        lam1 = 2*NR_BessJ1(arg1)/arg1;
+    }
+    if (arg2 == 0){
+        lam2 = 1;
+    }else{
+        lam2 = 2*NR_BessJ1(arg2)/arg2;
+    }
+    //Todo: Need to check psi behavior as gamma goes to 1.
     psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);         //SRK 10/19/00
-
     sinarg = qq*length*dum/2;
-    t2 = sin(sinarg)/sinarg;
+    if (sinarg == 0){
+        t2 = 1;
+    }else{
+        t2 = sin(sinarg)/sinarg;
+    }
 
     retval = psi*psi*t2*t2;
@@ -2084 +2222 @@
     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){
+        val = 1; // as arg --> 0, val should go to 1.0
+    }else{
+        val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
+    }
     return (val);
 
@@ -2099 +2240 @@
         // 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 */
@@ -2113 +2254 @@
     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){
+        be = sin(theta);
+    }else{
+        be = t2 / b2;
+    }
+    if (siarg == 0){
+        si = 1;
+    }else{
+        si = t1 / b1;
+    }
+    retval = be * si;
 
     return (retval);
@@ -2136 +2287 @@
 
         arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
-
-        retval = 2*NR_BessJ1(arg)/arg;
+        if (arg == 0){
+                retval = 1;
+        }else{
+                retval = 2*NR_BessJ1(arg)/arg;
+        }
 
         //square it

sansmodels/src/libigor/libCylinder.h

@@ -34 +34 @@
 double EllipCylKernel(double qq, double ra,double nu, double theta);
 double TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy);
+double lamellar_kernel(double dp[], double q);
 double PPKernel(double aa, double mu, double uu);
 double HolCylKernel(double qq, double rcore, double rshell, double length, double dum);

sansmodels/src/libigor/libSphere.c

@@ -14 +14 @@
 {
         double scale,radius,delrho,bkg;         //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];
@@ -21 +21 @@
         delrho = dp[2];
         bkg = dp[3];
-        //handle q==0 separately
-        if(q==0){
-                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*scale*1.0e8 + bkg;
-                return(f);
+
+        qr = q*radius;
+        if(qr == 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;             // [=] A-1
                                                         // normalize to single particle volume, convert to 1/cm
         f2 = f * f / vol * 1.0e8;               // [=] 1/cm
-        
+
         return(scale*f2+bkg);   //scale, and add in the background
 }
@@ -43 +44 @@
         double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names
         double bes,f,vol,qr,contr,f2;
-        
+
         pi = 4.0*atan(1.0);
         x=q;
-        
+
         scale = dp[0];
         rcore = dp[1];
@@ -57 +58 @@
         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;
@@ -63 +68 @@
         qr=x*(rcore+thick);
         contr = rhoshel-rhosolv;
-        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*pow((rcore+thick),3);
         f += vol*bes*contr;
-        
+
         // normalize to particle volume and rescale from [A-1] to [cm-1]
         f2 = f*f/vol*1.0e8;
-        
+
         //scale if desired
         f2 *= scale;
         // then add in the background
         f2 += bkg;
-        
+
         return(f2);
 }
@@ -88 +97 @@
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale = dp[0];
         rcore = dp[1];
@@ -99 +108 @@
         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;
@@ -105 +118 @@
         qr=x*(rcore+thick);
         contr = rhoshel-rhosolv;
-        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*pow((rcore+thick),3);
         f += vol*bes*contr;
-        
+
         // 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));
         f2 = f*f/vol*1.0e8;
-        
+
         //scale if desired
         f2 *= scale;
         // then add in the background
         f2 += bkg;
-        
+
         return(f2);
 }
@@ -135 +152 @@
         double zp3,vpoly;
         double s2y, arg1, arg2, arg3, drh1, drh2;
-        
+
         pi = 4.0*atan(1.0);
         qq= q;
@@ -145 +162 @@
         drho2 = dp[5]-dp[6];            //shell-solvent
         bkg = dp[7];
-        
-        zz = (1.0/sig)*(1.0/sig) - 1.0;   
-        
+
+        zz = (1.0/sig)*(1.0/sig) - 1.0;
+
         h=qq;
-        
+
         drh1 = drho1;
         drh2 = drho2;
@@ -157 +174 @@
         zp3 = zz + 3.;
         vpoly = 4*pi/3*zp3*zp2/zp1/zp1*pow((corrad+del),3);
-        
-        
+
+
         // remember that h is the passed in value of q for the calculation
         y = h *del;
@@ -179 +196 @@
         c8 = c4 - .5 + g * cy - g * .5 * g * (y * y + 1.) * c2y;
         c9 = g * sy * (1. - g * cy);
-        
+
         tb = log(zp1 * zp1 / (zp1 * zp1 + x * 4. * x));
         tb1 = exp(zp1 * .5 * tb);
         tb2 = exp(zp2 * .5 * tb);
         tb3 = exp(zp3 * .5 * tb);
-        
+
         t1 = c1 + c2 * x + c3 * x * x * zp2 / zp1;
         t2 = tb1 * (c4 * cos(arg1) + c7 * sin(arg1));
@@ -199 +216 @@
         // then add in the background
         form += bkg;
-        
+
         return(form);
 }
@@ -218 +235 @@
         double ka2,r1,heff;
         double hh,k;
-        
+
         pi = 4.0*atan(1.0);
         k= q;
-        
+
         rad = dp[0];            // radius (A)
         z2 = dp[1];             //polydispersity (0<z2<1)
@@ -228 +245 @@
         bkg = dp[4];
         sigma = 2*rad;
-        
+
         zz=1.0/(z2*z2)-1.0;
         bb = sigma/(zz+1.0);
         cc = zz+1.0;
-        
+
         //*c   Compute the number density by <r-cubed>, not <r> cubed*/
         r1 = sigma/2.0;
@@ -264 +281 @@
         d5=pow((pi/capd),2)*((rho/k)*(chi-1.0)+0.5*k*t2);
         d6=pow((pi/capd),2)*(rho/k)*(chi2-s2);
-        
+
         e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1.0)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+pow((k*s2-p1),2));
         ee=1.0-(2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2.0*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
         y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-pp)*(chi2-s2)-2.0*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1.0));
-        yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;    
-        
+        yy = (2.0*pi/capd)*(1.0+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1.0)-(2.0*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1;
+
         capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1));
         capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2));
         capmu=2.0*pi*cont*rho/k/k/k*(1.0-chi-0.5*k*p1);
         capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2);
-        
+
         h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4));
         h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5));
         h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2));
         h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3));
-        
+
         //*  This part computes the second integral in equation (1) of the paper.*/
-        
+
         hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy));
-        
+
         //*  This part computes the first integral in equation (1).  It also
         // generates the KC approximated effective structure factor.*/
-        
+
         pq=4.0*pi*cont*(sin(k*sigma/2.0)-0.5*k*sigma*cos(k*sigma/2.0));
         hint2=8.0*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1.0-chi-k*p1+0.25*k*k*(s2+chi2));
-        
+
         ka=k*(sigma/2.0);
         //
@@ -297 +314 @@
         ka2=ka*ka;
         //*
-        //  heff is PY analytical solution for intensity divided by the 
+        //  heff is PY analytical solution for intensity divided by the
         //   form factor.  happ is the KC approximated effective S(q)
-         
+
          //*******************
          //  add in the background then return the intensity value
-         
+
          return(hh+bkg);        //scale, and add in the background
 }
@@ -315 +332 @@
         double scale,ravg,pd,bkg,rho,rhos;              //my local names
         double scale2,ravg2,pd2,rho2;           //my local names
-        
+
         x= q;
-        
+
         scale = dp[0];
         ravg = dp[1];
@@ -328 +345 @@
         rhos = dp[8];
         bkg = dp[9];
-        
+
         pq = SchulzSphere_Fn( scale,  ravg,  pd,  rho,  rhos, x);
         pq += SchulzSphere_Fn( scale2,  ravg2,  pd2,  rho2,  rhos, x);
         // add in the background
         pq += bkg;
-        
+
         return (pq);
 }
@@ -344 +361 @@
         double x,pq;
         double scale,ravg,pd,bkg,rho,rhos;              //my local names
-        
+
         x= q;
-        
+
         scale = dp[0];
         ravg = dp[1];
@@ -356 +373 @@
         // add in the background
         pq += bkg;
-        
+
         return(pq);
 }
@@ -367 +384 @@
         double aa,at1,at2,rt1,rt2,rt3,t1,t2,t3;
         double v1,v2,v3,g1,pq,pi,delrho,zz;
-        
+
         pi = 4.0*atan(1.0);
         delrho = rho-rhos;
-        zz = (1.0/pd)*(1.0/pd) - 1.0;   
-        
+        zz = (1.0/pd)*(1.0/pd) - 1.0;
+
         zp1 = zz + 1.0;
         zp2 = zz + 2.0;
@@ -381 +398 @@
         //
         aa = (zz+1)/x/ravg;
-        
+
         at1 = atan(1.0/aa);
         at2 = atan(2.0/aa);
@@ -387 +404 @@
         //  calculations are performed to avoid  large # errors
         // - trick is to propogate the a^(z+7) term through the g1
-        // 
+        //
         t1 = zp7*log10(aa) - zp1/2.0*log10(aa*aa+4.0);
         t2 = zp7*log10(aa) - zp3/2.0*log10(aa*aa+4.0);
@@ -399 +416 @@
         v3 = -2.0*zp1*rt3*sin(zp2*at2);
         g1 = (v1+v2+v3);
-        
+
         pq = log10(g1) - 6.0*log10(zp1) + 6.0*log10(ravg);
         pq = pow(10,pq)*8*pi*pi*delrho*delrho;
-        
+
         //
-        // beta factor is not used here, but could be for the 
+        // beta factor is not used here, but could be for the
         // decoupling approximation
-        // 
+        //
         //      g11 = g1
         //      gd = -zp7*log(aa)
         //      g1 = log(g11) + gd
-        //                       
+        //
         //      t1 = zp1*at1
         //      t2 = zp2*at1
         //      g2 = sin( t1 ) - zp1/sqrt(aa*aa+1)*cos( t2 )
         //      g22 = g2*g2
-        //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22) 
+        //      beta = zp1*log(aa) - zp1*log(aa*aa+1) - g1 + log(g22)
         //      beta = 2*alog(beta)
-        
+
         //re-normalize by the average volume
         vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*ravg*ravg*ravg;
@@ -423 +440 @@
         //scale, convert to cm^-1
         pq *= scale * 1.0e8;
-    
+
     return(pq);
 }
@@ -435 +452 @@
         double scale,rad,pd,cont,bkg;           //my local names
         double inten,h1,qw,qr,width,sig,averad3;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale = dp[0];
         rad = dp[1];            // radius (A)
@@ -444 +461 @@
         cont = dp[3];           // contrast (A^-2)
         bkg = dp[4];
-        
+
         // as usual, poly = sig/ravg
         // for the rectangular distribution, sig = width/sqrt(3)
         // width is the HALF- WIDTH of the rectangular distrubution
-        
+
         sig = pd*rad;
         width = sqrt(3.0)*sig;
-        
+
         //x is the q-value
         qw = x*width;
@@ -463 +480 @@
         h1+= 3.0*qw*(sin(qr)*sin(qw))*(sin(qr)*sin(qw));
         h1 -= 6.0*qr*cos(qr)*sin(qr)*cos(qw)*sin(qw);
-        
+
         // calculate P(q) = <f^2>
         inten = 8.0*pi*pi*cont*cont/width/pow(x,7)*h1;
-        
+
         // beta(q) would be calculated as 2/width/x/h1*h2*h2
-        // with 
+        // with
         // h2 = 2*sin(x*rad)*sin(x*width)-x*rad*cos(x*rad)*sin(x*width)-x*width*sin(x*rad)*cos(x*width)
-        
+
         // normalize to the average volume
         // <R^3> = ravg^3*(1+3*pd^2)
         // or... "zf"  = (1 + 3*p^2), which will be greater than one
-        
+
         averad3 =  rad*rad*rad*(1.0+3.0*pd*pd);
         inten /= 4.0*pi/3.0*averad3;
@@ -483 +500 @@
         // then add in the background
         inten += bkg;
-        
+
         return(inten);
 }
@@ -497 +514 @@
         double va,vb,zi,yy,summ,inten;
         int nord=20,ii;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale=dp[0];
         rad=dp[1];
@@ -509 +526 @@
         delrho=rho-rhos;
         bkg=dp[5];
-        
+
         va = -4.0*sig + rad;
         if (va<0) {
@@ -515 +532 @@
         }
         vb = 4.0*sig +rad;
-        
+
         summ = 0.0;             // initialize integral
         for(ii=0;ii<nord;ii+=1) {
@@ -525 +542 @@
                 yy *= yy;
                 yy *= Gauss20Wt[ii] *  Gauss_distr(sig,rad,zi);
-                
+
                 summ += yy;             //add to the running total of the quadrature
         }
         // calculate value of integral to return
         inten = (vb-va)/2.0*summ;
-        
+
         //re-normalize by polydisperse sphere volume
         inten /= (4.0*pi/3.0*rad*rad*rad)*(1.0+3.0*pd*pd);
-        
+
         inten *= 1.0e8;
         inten *= scale;
         inten += bkg;
-        
+
     return(inten);      //scale, and add in the background
 }
@@ -550 +567 @@
         double va,vb,zi,yy,summ,inten;
         int nord=76,ii;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale=dp[0];
         rad=dp[1];              //rad is the median radius
@@ -562 +579 @@
         delrho=rho-rhos;
         bkg=dp[5];
-        
+
         va = -3.5*sig + mu;
         va = exp(va);
@@ -570 +587 @@
         vb = 3.5*sig*(1.0+sig) +mu;
         vb = exp(vb);
-        
+
         summ = 0.0;             // initialize integral
         for(ii=0;ii<nord;ii+=1) {
@@ -580 +597 @@
                 yy *= yy;
                 yy *= Gauss76Wt[ii] *  LogNormal_distr(sig,mu,zi);
-                
+
                 summ += yy;             //add to the running total of the quadrature
         }
         // calculate value of integral to return
         inten = (vb-va)/2.0*summ;
-        
+
         //re-normalize by polydisperse sphere volume
         r3 = exp(3.0*mu + 9.0/2.0*sig*sig);             // <R^3> directly
         inten /= (4.0*pi/3.0*r3);               //polydisperse volume
-        
+
         inten *= 1.0e8;
         inten *= scale;
         inten += bkg;
-        
+
         return(inten);
 }
@@ -599 +616 @@
 static double
 LogNormal_distr(double sig, double mu, double pt)
-{       
+{
         double retval,pi;
-        
+
         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 );
@@ -609 +626 @@
 static double
 Gauss_distr(double sig, double avg, double pt)
-{       
+{
         double retval,Pi;
-        
+
         Pi = 4.0*atan(1.0);
         retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0);
@@ -627 +644 @@
         double sld1,sld2,sld3,zp1,zp2,zp3,vpoly;
         double pi43,c1,c2,form,volume,arg1,arg2;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale = dp[0];
         corrad = dp[1];
@@ -639 +656 @@
         sld3 = dp[6];
         bkg = dp[7];
-        
-        zz = (1.0/sig)*(1.0/sig) - 1.0;   
+
+        zz = (1.0/sig)*(1.0/sig) - 1.0;
         shlrad = corrad + thick;
         drho1 = sld1-sld2;              //core-shell
@@ -648 +665 @@
         zp3 = zz + 3.;
         vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((corrad+thick),3);
-        
+
         // the beta factor is not calculated
         // the calculated form factor <f^2> has units [length^2]
         // and must be multiplied by number density [l^-3] and the correct unit
         // conversion to get to absolute scale
-        
+
         pi43=4.0/3.0*pi;
         pp=corrad/shlrad;
@@ -659 +676 @@
         c1=drho1*volume;
         c2=drho2*volume;
-        
+
         arg1 = x*shlrad*pp;
         arg2 = x*shlrad;
-        
+
         form=pow(pp,6)*c1*c1*fnt2(arg1,zz);
         form += c2*c2*fnt2(arg2,zz);
         form += 2.0*c1*c2*fnt3(arg2,pp,zz);
-        
+
         //convert the result to [cm^-1]
-        
+
         //scale the result
         // - divide by the polydisperse volume, mult by 10^8
@@ -674 +691 @@
         form *= 1.0e8;
         form *= scale;
-        
+
         //add in the background
         form += bkg;
-        
+
         return(form);
 }
@@ -689 +706 @@
 {
         double z1,z2,z3,u,ww,term1,term2,term3,ans;
-        
+
         z1=zz+1.0;
         z2=zz+2.0;
@@ -699 +716 @@
         term3=1.0+cos(z3*ww)/pow((1.0+4.0*u*u),(z3/2.0));
         ans=(4.50/z1/pow(yy,6))*(z1*(1.0-term1-term2)+yy*yy*z2*term3);
-        
+
         return(ans);
 }
@@ -709 +726 @@
 double
 fnt3(double yy, double pp, double zz)
-{       
+{
         double z1,z2,z3,yp,yn,up,un,vp,vn,term1,term2,term3,term4,term5,term6,ans;
-        
+
         z1=zz+1.0;
         z2=zz+2.0;
@@ -729 +746 @@
         ans=4.5/z1/pow(yy,6);
         ans *=(z1*(term1-term2)+yy*yy*pp*z2*(term3+term4)+z1*(term5-term6));
-        
+
         return(ans);
 }
@@ -752 +769 @@
         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;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
@@ -765 +782 @@
         rhos = dp[6];
         bgd = dp[7];
-        
-        
+
+
         phi = phi1 + phi2;
         aa = r1/r2;
@@ -775 +792 @@
         // calculate the PSF's here
         err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-        
+
         // /* do form factor calculations  */
-        
+
         v1 = 4.0*pi/3.0*r1*r1*r1;
         v2 = 4.0*pi/3.0*r2*r2*r2;
-        
+
         n1 = phi1/v1;
         n2 = phi2/v2;
-        
+
         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;
@@ -794 +821 @@
         ///* convert I(1/A) to (1/cm)  */
         inten *= 1.0e8;
-        
+
         inten += bgd;
-        
+
         return(inten);
 }
@@ -808 +835 @@
         double phi1,phi2,phr,a3;
         int err;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
@@ -827 +854 @@
         // calculate the PSF's here
         err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-        
+
     return(psf11);      //scale, and add in the background
 }
@@ -839 +866 @@
         double phi1,phi2,phr,a3;
         int err;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
@@ -858 +885 @@
         // calculate the PSF's here
         err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-        
+
     return(psf12);      //scale, and add in the background
 }
@@ -870 +897 @@
         double phi1,phi2,phr,a3;
         int err;
-        
+
         pi = 4.0*atan(1.0);
         x= q;
-        
+
         r2 = dp[0];
         r1 = dp[1];
@@ -890 +917 @@
         // calculate the PSF's here
         err = ashcroft(x,r2,nf2,aa,phi,&psf11,&psf22,&psf12);
-        
+
     return(psf22);      //scale, and add in the background
 }
@@ -898 +925 @@
 {
         //      variable qval,r2,nf2,aa,phi,&s11,&s22,&s12
-        
+
         //   calculate constant terms
         double s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4;
@@ -905 +932 @@
         double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13;
         double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1;
-        
+
         s2 = 2.0*r2;
         s1 = aa*s2;
@@ -926 +953 @@
         gm1=(v1*a1+a3*v2*a2)*.5;
         gm12=2.*gm1*(1.-aa)/aa;
-        //c  
+        //c
         //c   calculate the direct correlation functions and print results
         //c
         //      do 20 j=1,npts
-        
+
         yy=qval*s2;
         //c   calculate direct correlation functions
@@ -941 +968 @@
         t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3;
         f11=24.*v1*(t1+t2+t3)/ay3;
-        
+
         //c ------c22
         y2=yy*yy;
@@ -949 +976 @@
         tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3;
         f22=24.*v2*(tt1+tt2+tt3)/y3;
-        
+
         //c   -----c12
         yl=.5*yy*(1.-aa);
@@ -969 +996 @@
         ttt4=a1*(t41+t42)/y1;
         f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3);
-        
+
         c11=f11;
         c22=f22;
         c12=f12;
         *s11=1./(1.+c11-(c12)*c12/(1.+c22));
-        *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
-        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));   
-        
+        *s22=1./(1.+c22-(c12)*c12/(1.+c11));
+        *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));
+
         return(err);
 }
@@ -989 +1016 @@
  // bragg peaks arise naturally from the periodicity of the sample
  // resolution smeared version gives he most appropriate view of the model
- 
+
         Warning:
  The call to WaveData() below returns a pointer to the middle
@@ -1004 +1031 @@
         double fval,voli,ri,sldi;
         double pi;
-        
-        pi = 4.0*atan(1.0);
-        
+
+        pi = 4.0*atan(1.0);
+
         x= q;
         scale = dp[0];
@@ -1016 +1043 @@
         num = dp[6];
         bkg = dp[7];
-        
+
         //calculate with a loop, two shells at a time
-        
+
         ii=0;
         fval=0;
-        
+
         do {
                 ri = rcore + (double)ii*ts + (double)ii*tw;
@@ -1033 +1060 @@
                 ii+=1;          //do 2 layers at a time
         } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
-        
+
         fval *= fval;           //square it
         fval /= voli;           //normalize by the overall volume
         fval *= scale*1e8;
         fval += bkg;
-        
+
         return(fval);
 }
@@ -1068 +1095 @@
         double fval,ri,pi;
         double avg,pd,zz,lo,hi,r1,r2,d1,d2,distr;
-        
-        pi = 4.0*atan(1.0);     
+
+        pi = 4.0*atan(1.0);
         x= q;
-        
+
         scale = dp[0];
         avg = dp[1];            // average (total) outer radius
@@ -1081 +1108 @@
         rhoshel = dp[7];
         bkg = dp[8];
-        
+
         zz = (1.0/pd)*(1.0/pd)-1.0;
-        
+
         //max radius set to be 5 std deviations past mean
         hi = avg + pd*avg*5.0;
         lo = avg - pd*avg*5.0;
-        
+
         maxPairs = trunc( (hi-rcore+tw)/(ts+tw) );
         minPairs = trunc( (lo-rcore+tw)/(ts+tw) );
         minPairs = (minPairs < 1) ? 1 : minPairs;       // need a minimum of one
-        
+
         ii=minPairs;
         fval=0;
@@ -1104 +1131 @@
                 r1 = r2;
                 d1 = d2;
-                
+
                 ri = (double)ii*(ts+tw) - tw + rcore;
                 fval += SchulzPoint(ri,avg,zz) * MultiShellGuts(x, rcore, ts, tw, rhocore, rhoshel, ii) * (4*pi/3*ri*ri*ri);
@@ -1116 +1143 @@
                 first = 0;
         } while(ii<=maxPairs);
-        
+
         fval /= 4*pi/3*avg*avg*avg;             //normalize by the overall volume
         fval /= distr;
         fval *= scale;
         fval += bkg;
-        
+
         return(fval);
 }
@@ -1127 +1154 @@
 double
 MultiShellGuts(double x,double rcore,double ts,double tw,double rhocore,double rhoshel,int num) {
-        
+
     double ri,sldi,fval,voli,pi;
     int ii;
-    
+
         pi = 4.0*atan(1.0);
     ii=0;
     fval=0;
-    
+
     do {
         ri = rcore + (double)ii*ts + (double)ii*tw;
@@ -1146 +1173 @@
         ii+=1;          //do 2 layers at a time
     } while(ii<=num-1);  //change to make 0 < num < 2 correspond to unilamellar vesicles (C. Glinka, 11/24/03)
-    
+
     fval *= fval;
     fval /= voli;
     fval *= 1e8;
-    
+
     return(fval);       // this result still needs to be multiplied by scale and have background added
 }
@@ -1156 +1183 @@
 static double
 SchulzPoint(double x, double avg, double zz) {
-        
+
     double dr;
-    
+
     dr = zz*log(x) - gammln(zz+1)+(zz+1)*log((zz+1)/avg)-(x/avg*(zz+1));
     return (exp(dr));
@@ -1165 +1192 @@
 static double
 gammln(double xx) {
-        
+
     double x,y,tmp,ser;
     static double cof[6]={76.18009172947146,-86.50532032941677,
@@ -1171 +1198 @@
                 0.1208650973866179e-2,-0.5395239384953e-5};
     int j;
-        
+
     y=x=xx;
     tmp=x+5.5;
@@ -1182 +1209 @@
 double
 F_func(double qr) {
-        return(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
-}
-
+        double sc;
+        if (qr == 0){
+                sc = 1.0;
+        }else{
+                sc=(3*(sin(qr) - qr*cos(qr))/(qr*qr*qr));
+        }
+        return sc;
+}
+