-                      r222d75c7
+                      r6e93a02
 #include "GaussWeights.h"
 #include "libCylinder.h"
 /*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x
 …
         int i;
         double Pi;
         double scale,radius,length,delrho,bkg,halfheight;               //local variables of coefficient wave
+        double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv;                //local variables of coefficient wave
         int nord=76;                    //order of integration
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,vcyl;                 //running tally of integration
         Pi = 4.0*atan(1.0);
         lolim = 0;
+        lolim = 0.0;
         uplim = Pi/2.0;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         radius = dp[1];
         length = dp[2];
+        delrho = dp[3];
+        bkg = dp[4];
+        sldCyl = dp[3];
+        sldSolv = dp[4];
+        bkg = dp[5];
+        delrho = sldCyl-sldSolv;
         halfheight = length/2.0;
         for(i=0;i<nord;i++) {
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         // Multiply by contrast^2
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
+{
         int i,j;
         double Pi;
+        double Pi,slde,sld;
         double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
         int nord=76;                    //order of integration
 …
         vaj=0.0;
         vbj=Pi;         //endpoints of inner integral
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         ra = dp[1];
         nu = dp[2];
         length = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        slde = dp[4];
+        sld = dp[5];
+        delrho = slde - sld;
+        bkg = dp[6];
         for(i=0;i<nord;i++) {
                 //setup inner integral over the ellipsoidal cross-section
                 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
 …
                 //divide integral by Pi
                 answer /=Pi;
                 //now calculate outer integral
                 arg = q*length*zi/2.0;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
+{
         int i,j;
         double Pi;
+        double Pi,slde,sld;
         double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
         int nordi=76;                   //order of integration
 …
         vaj=0.0;
         vbj=Pi;         //endpoints of inner integral
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         ra = dp[1];
         nu = dp[2];
         length = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        slde = dp[4];
+        sld = dp[5];
+        delrho = slde - sld;
+        bkg = dp[6];
         for(i=0;i<nordi;i++) {
                 //setup inner integral over the ellipsoidal cross-section
                 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
 …
                 //divide integral by Pi
                 answer /=Pi;
                 //now calculate outer integral
                 arg = q*length*zi/2;
+                arg = q*length*zi/2.0;
                 if (arg == 0.0){
                         si = 1.0;
 …
                 summ += yyy;
+        }
         answer = (vb-va)/2.0*summ;
         // Multiply by contrast^2
 …
         double va,vb;           //upper and lower integration limits
         double summ,zi,yyy,answer;                      //running tally of integration
         double summj,vaj,vbj,zij;                       //for the inner integration
+        double summj,vaj,vbj,zij,slde,sld;                      //for the inner integration
         Pi = 4.0*atan(1.0);
         va = 0.0;
 …
         vaj = 0.0;
         vbj = 1.0;              //endpoints of inner integral
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         aa = dp[1];
         bb = dp[2];
         cc = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        slde = dp[4];
+        sld = dp[5];
+        delrho = slde - sld;
+        bkg = dp[6];
         for(i=0;i<nordi;i++) {
                 //setup inner integral over the ellipsoidal cross-section
                 summj=0;
+                summj=0.0;
                 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
                 for(j=0;j<nordj;j++) {
 …
                 //now calculate the value of the inner integral
                 answer = (vbj-vaj)/2.0*summj;
                 //now calculate outer integral
                 yyy = Gauss76Wt[i] * answer;
                 summ += yyy;
         }               //final scaling is done at the end of the function, after the NT_FP64 case
         answer = (vb-va)/2.0*summ;
         // Multiply by contrast^2
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double summ,yyy,answer;                 //running tally of integration
         double summj,vaj,vbj;                   //for the inner integration
         double mu,mudum,arg,sigma,uu,vol;
+        double mu,mudum,arg,sigma,uu,vol,sldp,sld;
         //      Pi = 4.0*atan(1.0);
         va = 0.0;
 …
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         aa = dp[1];
         bb = dp[2];
         cc = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        sldp = dp[4];
+        sld = dp[5];
+        delrho = sldp - sld;
+        bkg = dp[6];
         mu = q*bb;
         vol = aa*bb*cc;
 …
         aa = aa/bb;
         cc = cc/bb;
         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
                 for(j=0;j<nordj;j++) {
                         //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;
 …
                         answer *= sin(arg)*sin(arg)/arg/arg;
+                }
                 //now sum up the outer integral
                 yyy = Gauss76Wt[i] * answer;
                 summ += yyy;
         }               //final scaling is done at the end of the function, after the NT_FP64 case
         answer = (vb-va)/2.0*summ;
         // Multiply by contrast^2
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         int nord=76;                    //order of integration
         double va,vb,zi;                //upper and lower integration limits
         double summ,answer,pi;                  //running tally of integration
+        double summ,answer,pi,sldc,sld;                 //running tally of integration
         pi = 4.0*atan(1.0);
         va = 0.0;
 …
         summ = 0.0;                     //initialize intergral
         scale = dp[0];          //make local copies in case memory moves
         rcore = dp[1];
         rshell = dp[2];
         length = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        sldc = dp[4];
+        sld = dp[5];
+        delrho = sldc - sld;
+        bkg = dp[6];
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
                 summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
+        }
         answer = (vb-va)/2.0*summ;
         // Multiply by contrast^2
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         int nord=76;                    //order of integration
         double va,vb,zi;                //upper and lower integration limits
         double summ,answer,pi;                  //running tally of integration
+        double summ,answer,pi,slde,sld;                 //running tally of integration
         pi = 4.0*atan(1.0);
         va = 0.0;
 …
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         nua = dp[1];
         a = dp[2];
+        delrho = dp[3];
+        bkg = dp[4];
+        slde = dp[3];
+        sld = dp[4];
+        delrho = slde - sld;
+        bkg = dp[5];
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
                 summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
+        }
         answer = (vb-va)/2.0*summ;
         // Multiply by contrast^2
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
         double range,zz,Pi;
+        double range,zz,Pi,sldc,sld;
         Pi = 4.0*atan(1.0);
         range = 3.4;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         radius = dp[1];
         length = dp[2];
         pd = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        sldc = dp[4];
+        sld = dp[5];
+        delrho = sldc - sld;
+        bkg = dp[6];
         zz = (1.0/pd)*(1.0/pd) - 1.0;
         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) {
 …
+        }
         uplim = radius*(1.0+range*pd);
         for(i=0;i<nord;i++) {
                 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         //normalize by average cylinder volume
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
         double range,zz,Pi;
+        double range,zz,Pi,sldc,sld;
         Pi = 4.0*atan(1.0);
         range = 3.4;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         radius = dp[1];
         length = dp[2];
         pd = dp[3];
+        delrho = dp[4];
+        bkg = dp[5];
+        sldc = dp[4];
+        sld = dp[5];
+        delrho = sldc - sld;
+        bkg = dp[6];
         zz = (1.0/pd)*(1.0/pd) - 1.0;
         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) {
 …
+        }
         uplim = length*(1.0+range*pd);
         for(i=0;i<nord;i++) {
                 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         //normalize by average cylinder volume (first moment)
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double summ,zi,yyy,answer,Vcyl;                 //running tally of integration
         double Pi;
         Pi = 4.0*atan(1.0);
         lolim = 0.0;
         uplim = Pi/2.0;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];          //make local copies in case memory moves
         rcore = dp[1];
 …
         rhosolv = dp[6];
         bkg = dp[7];
         halfheight = length/2.0;
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         // length is the total core length
+        // length is the total core length
         Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
         answer /= Vcyl;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double summ,yyy,answer,Vpoly;                   //running tally of integration
         double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
         Pi = 4.0*atan(1.0);
         summ = 0.0;                     //initialize intergral
         Rsqrsumm = 0.0;
         scale = dp[0];
         radius = dp[1];
 …
         rhosolv = dp[8];
         bkg = dp[9];
         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));
         for(i=0;i<nord;i++) {
                 rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 Rsqrsumm += Rsqryyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,oblatevol;                    //running tally of integration
         double Pi;
+        double Pi,sldc,slds,sld;
         Pi = 4.0*atan(1.0);
         lolim = 0.0;
         uplim = 1.0;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];          //make local copies in case memory moves
         crmaj = dp[1];
 …
         trmaj = dp[3];
         trmin = dp[4];
+        delpc = dp[5];
+        delps = dp[6];
+        bkg = dp[7];
+        sldc = dp[5];
+        slds = dp[6];
+        sld = dp[7];
+        delpc = sldc - slds;    //core - shell
+        delps = slds - sld;             //shell - solvent
+        bkg = dp[8];
         for(i=0;i<nord;i++) {
 …
                 summ += yyy;
+        }
         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;
         //convert to [cm-1]
         answer *= 1.0e8;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,prolatevol;                   //running tally of integration
         double Pi;
+        double Pi,sldc,slds,sld;
         Pi = 4.0*atan(1.0);
         lolim = 0.0;
         uplim = 1.0;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];          //make local copies in case memory moves
         crmaj = dp[1];
 …
         trmaj = dp[3];
         trmin = dp[4];
+        delpc = dp[5];
+        delps = dp[6];
+        bkg = dp[7];
+        sldc = dp[5];
+        slds = dp[6];
+        sld = dp[7];
+        delpc = sldc - slds;            //core - shell
+        delps = slds - sld;             //shell  - sovent
+        bkg = dp[8];
         for(i=0;i<nord;i++) {
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         // normalize by particle volume
         prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin;
         answer /= prolatevol;
         //convert to [cm-1]
         answer *= 1.0e8;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         int nord=76;                    //order of integration
         double Pi;
         Pi = 4.0*atan(1.0);
         va = 0.0;
         vb = Pi/2.0;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];
         rcore = dp[1];
 …
         gsd = dp[8];
         bkg = dp[9];
         d=2.0*thick+length;
         halfheight = length/2.0;
         for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
 …
                 summ += yyy;
+        }
         answer = (vb-va)/2.0*summ;
         // length is the total core length
+        // length is the total core length
         vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
         answer /= vcyl;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double scale,del,sig,contr,bkg;         //local variables of coefficient wave
         double inten, qval,Pq;
         double Pi;
+        double Pi,sldb,sld;
         Pi = 4.0*atan(1.0);
         scale = dp[0];
         del = dp[1];
         sig = dp[2]*del;
+        contr = dp[3];
+        bkg = dp[4];
+        qval = q;
+        sldb = dp[3];
+        sld = dp[4];
+        contr = sldb - sld;
+        bkg = dp[5];
+        qval=q;
         Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
         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
+-------
+------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!!
+        */
+--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
+model is smeared just like any other function
+*/
 double
 LamellarPS(double dp[], double q)
+{
         double scale,dd,del,sig,contr,NN,Cp,bkg;                //local variables of coefficient wave
         double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
         double Pi,Euler,dQDefault,fii;
+        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
+        double Pi,Euler,dQDefault,fii,sldb,sld;
         int ii,NNint;
+//      char buf[256];
         Euler = 0.5772156649;           // Euler's constant
+        dQDefault = 0.0;//0.0025;               //[=] 1/A, q-resolution, default value
+//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
+        dQDefault = 0.0;
         dQ = dQDefault;
         Pi = 4.0*atan(1.0);
         qval = q;
         scale = dp[0];
         dd = dp[1];
         del = dp[2];
         sig = dp[3]*del;
+        contr = dp[4];
+        NN = trunc(dp[5]);              //be sure that NN is an integer
+        Cp = dp[6];
+        bkg = dp[7];
+        sldb = dp[4];
+        sld = dp[5];
+        contr = sldb - sld;
+        NN = trunc(dp[6]);              //be sure that NN is an integer
+        Cp = dp[7];
+        bkg = dp[8];
         Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
         NNint = (int)NN;                //cast to an integer for the loop
+//      sprintf(buf, "qval = %g\r", qval);
+//      XOPNotice(buf);
         ii=0;
         Sq = 0.0;
         for(ii=1;ii<(NNint-1);ii+=1) {
                 fii = (double)ii;               //do I really need to do this?
                 temp = 0.0;
                 alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
+                alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler);
                 t1 = 2.0*dQ*dQ*dd*dd*alpha;
                 t2 = 2.0*qval*qval*dd*dd*alpha;
                 t3 = dQ*dQ*dd*dd*ii*ii;
                 temp = 1.0-ii/NN;
                 temp *= cos(dd*qval*ii/(1.0+t1));
+                t3 = dQ*dQ*dd*dd*fii*fii;
+                temp = 1.0-fii/NN;
+                temp *= cos(dd*qval*fii/(1.0+t1));
                 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
                 temp /= sqrt(1.0+t1);
                 Sq += temp;
+        }
         Sq *= 2.0;
         Sq += 1.0;
         inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
         inten *= 1.0e8;         // 1/A to 1/cm
     return(inten+bkg);
+}
 …
 /*      LamellarPS_HGX  :  calculates the form factor of a lamellar structure - with S(q) effects included
+-------
+------- resolution effects ARE included, but only a CONSTANT default value, not the real q-dependent resolution!!
+        */
+--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
+model is smeared just like any other function
+*/
 double
 LamellarPS_HG(double dp[], double q)
 …
         double Pi,Euler,dQDefault,fii;
         int ii,NNint;
         Euler = 0.5772156649;           // Euler's constant
+        dQDefault = 0.0; //0.0025;              //[=] 1/A, q-resolution, default value
+//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
+        dQDefault = 0.0;
         dQ = dQDefault;
         Pi = 4.0*atan(1.0);
         qval= q;
         scale = dp[0];
         dd = dp[1];
 …
         Cp = dp[8];
         bkg = dp[9];
         drh = SLD_H - SLD_S;
         drt = SLD_T - SLD_S;    //correction 13FEB06 by L.Porcar
         Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
         Pq *= Pq;
         Pq *= 4.0/(qval*qval);
         NNint = (int)NN;                //cast to an integer for the loop
         ii=0;
         Sq = 0.0;
         for(ii=1;ii<(NNint-1);ii+=1) {
                 fii = (double)ii;               //do I really need to do this?
                 temp = 0.0;
                 alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
 …
                 t2 = 2.0*qval*qval*dd*dd*alpha;
                 t3 = dQ*dQ*dd*dd*ii*ii;
                 temp = 1.0-ii/NN;
                 temp *= cos(dd*qval*ii/(1.0+t1));
                 temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
                 temp /= sqrt(1.0+t1);
                 Sq += temp;
+        }
         Sq *= 2.0;
         Sq += 1.0;
         inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
         inten *= 1.0e8;         // 1/A to 1/cm
         return(inten+bkg);
+}
 …
         double inten, qval,Pq,drh,drt;
         double Pi;
         Pi = 4.0*atan(1.0);
         qval= q;
 …
         slds = dp[5];
         bkg = dp[6];
         drh = sldh - slds;
         drt = sldt - slds;              //correction 13FEB06 by L.Porcar
         Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
         Pq *= Pq;
         Pq *= 4.0/(qval*qval);
         inten = 2.0*Pi*scale*Pq/(qval*qval);            //dimensionless...
         inten /= 2.0*(delT+delH);                       //normalize by the bilayer thickness
         inten *= 1.0e8;         // 1/A to 1/cm
         return(inten+bkg);
+}
 …
 FlexExclVolCyl(double dp[], double q)
+{
         double scale,L,B,bkg,rad,qr,cont;
+        double scale,L,B,bkg,rad,qr,cont,sldc,slds;
         double Pi,flex,crossSect,answer;
         Pi = 4.0*atan(1.0);
         scale = dp[0];          //make local copies in case memory moves
         L = dp[1];
         B = dp[2];
         rad = dp[3];
+        cont = dp[4];
+        bkg = dp[5];
+        sldc = dp[4];
+        slds = dp[5];
+        cont = sldc-slds;
+        bkg = dp[6];
         qr = q*rad;
         flex = Sk_WR(q,L,B);
         crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
         flex *= crossSect;
 …
         flex *= 1.0e8;
         answer = scale*flex + bkg;
         return answer;
+}
 …
 FlexCyl_Ellip(double dp[], double q)
+{
         double scale,L,B,bkg,rad,qr,cont,ellRatio;
+        double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc;
         double Pi,flex,crossSect,answer;
         Pi = 4.0*atan(1.0);
         scale = dp[0];          //make local copies in case memory moves
 …
         rad = dp[3];
         ellRatio = dp[4];
+        cont = dp[5];
+        bkg = dp[6];
+        sldc = dp[5];
+        slds = dp[6];
+        bkg = dp[7];
+        cont = sldc - slds;
         qr = q*rad;
         flex = Sk_WR(q,L,B);
         crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
         flex *= crossSect;
 …
         flex *= 1.0e8;
         answer = scale*flex + bkg;
         return answer;
+}
 …
     double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
     int i,nord=76;
     Pi = 4.0*atan(1.0);
     lolim=0.0;
     uplim=Pi/2.0;
     summ=0.0;
     for(i=0;i<nord;i++) {
                 zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
     answer *= 2.0/Pi;
     return(answer);
+}
 /*      FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
 …
+{
         int i;
         double scale,radius,length,pd,delrho,bkg,lb;            //local variables of coefficient wave
+        double scale,radius,length,pd,bkg,lb,delrho,sldc,slds;          //local variables of coefficient wave
         int nord=20;                    //order of integration
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
         double range,zz,Pi;
         Pi = 4.0*atan(1.0);
         range = 3.4;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
 …
         lb = dp[3];
         radius = dp[4];
+        delrho = dp[5];
+        bkg = dp[6];
+        sldc = dp[5];
+        slds = dp[6];
+        bkg = dp[7];
+        delrho = sldc - slds;
         zz = (1.0/pd)*(1.0/pd) - 1.0;
         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) {
 …
+        }
         uplim = length*(1.0+range*pd);
         for(i=0;i<nord;i++) {
                 zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         //normalize by average cylinder volume (first moment), using the average length
         Vpoly=Pi*radius*radius*length;
         answer /= Vpoly;
         answer *=delrho*delrho;
         //convert to [cm-1]
         answer *= 1.0e8;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
+{
         int i;
         double scale,radius,length,pd,delrho,bkg,lb;            //local variables of coefficient wave
+        double scale,radius,length,pd,delrho,bkg,lb,sldc,slds;          //local variables of coefficient wave
         int nord=76;                    //order of integration
         double uplim,lolim;             //upper and lower integration limits
         double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
         double range,zz,Pi;
         Pi = 4.0*atan(1.0);
         range = 3.4;
         summ = 0.0;                     //initialize intergral
         scale = dp[0];                  //make local copies in case memory moves
         length = dp[1];                 //radius
 …
         radius = dp[3];
         pd = dp[4];
+        delrho = dp[5];
+        bkg = dp[6];
+        sldc = dp[5];
+        slds = dp[6];
+        bkg = dp[7];
+        delrho = sldc-slds;
         zz = (1.0/pd)*(1.0/pd) - 1.0;
         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) {
 …
+        }
         uplim = radius*(1.0+range*pd);
         for(i=0;i<nord;i++) {
                 //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
 …
                 summ += yyy;
+        }
         answer = (uplim-lolim)/2.0*summ;
         //normalize by average cylinder volume (second moment), using the average radius
         Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
         answer /= Vpoly;
         answer *=delrho*delrho;
         //convert to [cm-1]
         answer *= 1.0e8;
 …
         // add in the background
         answer += bkg;
         return answer;
+}
 …
         double p1,p2,p1short,p2short,q0,qconnect;
         double C,epsilon,ans,q0short,Sexvmodify,pi;
         pi = 4.0*atan(1.0);
         p1 = 4.12;
         p2 = 4.42;
 …
         q0 = 3.1;
         qconnect = q0/b;
         //
+        //
         q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
         //
         if(L/b > 10.0) {
 …
+        }
         //
         if( L > 4*b ) { // Longer Chains
                 if (q*b <= 3.1) {               //Modified by Yun on Oct. 15,
 …
                 } else { //q(i)*b > 3.1
                         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
+                }
+                }
         } else { //L <= 4*b Shorter Chains
                 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
 …
+                }
+        }
         return(ans);
         //return(a2long(q, L, b, p1, p2, q0));
 …
     double yy;
     yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
     return (yy);
+}
 …
+{
     double yy;
     yy = Rgsquareshort(q,L,b)*q*q;
     return (yy);
+}
 …
 static double
 Rgsquarezero(double q, double L, double b)
+{
+{
     double yy;
     yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));
     return (yy);
+}
 …
 static double
 Rgsquareshort(double q, double L, double b)
+{
+{
     double yy;
     yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
     return (yy);
+}
 …
     double yy;
     yy = AlphaSquare(L/b)*L*b/6.0;
     return (yy);
+}
 …
 static double
 AlphaSquare(double x)
+{
+{
     double yy;
     yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
     return (yy);
+}
 …
 static double
 miu(double x)
+{
+{
     double yy;
     yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x)));
     return (yy);
+}
 …
 static double
 Sdebye(double q, double L, double b)
+{
+{
     double yy;
     yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
     return (yy);
+}
 …
 static double
 Sdebye1(double q, double L, double b)
+{
+{
     double yy;
     yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
     return (yy);
+}
 …
 static double
 Sexv(double q, double L, double b)
+{
+{
     double yy,C1,C2,C3,miu,Rg2;
     C1=1.22;
 …
     C3=-1.651;
     miu = 0.585;
     Rg2 = Rgsquare(q,L,b);
     yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
     return (yy);
+}
 …
 static double
 Sexvnew(double q, double L, double b)
+{
+{
     double yy,C1,C2,C3,miu;
     double del=1.05,C_star2,Rg2;
     C1=1.22;
     C2=0.4288;
     C3=-1.651;
     miu = 0.585;
     //calculating the derivative to decide on the corection (cutoff) term?
     // I have modified this from WRs original code
     if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
         C_star2 = 0.0;
 …
         C_star2 = 1.0;
+    }
     Rg2 = Rgsquare(q,L,b);
         yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));
     return (yy);
+}
 …
 static double
 a2short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+{
     double yy,Rg2_sh;
     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
     double pi;
     E = 2.718281828459045091;
         pi = 4.0*atan(1.0);
 …
     t1 = ((q0*q0*Rg2_sh)/(b*b));
     Et1 = pow(E,t1);
     Emt1 =pow(E,-t1);
+    Emt1 =pow(E,-t1);
     q02 = q0*q0;
     q0p = pow(q0,(-4.0 + p2short) );
     //E is the number e
     yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2)))))));
     return (yy);
+}
 …
 static double
 a1short(double q, double L, double b, double p1short, double p2short, double q0)
+{
+{
     double yy,Rg2_sh;
     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
         double pi;
     E = 2.718281828459045091;
         pi = 4.0*atan(1.0);
 …
     t1 = ((q0*q0*Rg2_sh)/(b*b));
     Et1 = pow(E,t1);
     Emt1 =pow(E,-t1);
+    Emt1 =pow(E,-t1);
     q02 = q0*q0;
     q0p = pow(q0,(-4.0 + p1short) );
 …
     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
     double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
     pi = 4.0*atan(1.0);
     E = 2.718281828459045091;
 …
         C = 1.0;
+    }
     C1 = 1.22;
     C2 = 0.4288;
 …
     C5 = 0.1477;
     miu = 0.585;
     Rg2 = Rgsquare(q,L,b);
     Rg22 = Rg2*Rg2;
 …
     q04 = q03*q0;
     q05 = q04*q0;
     t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
 …
     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*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);
+    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);
 …
 static double
 sech_WR(double x)
+{
+{
         return(1/cosh(x));
+}
 …
 static double
 a1long(double q, double L, double b, double p1, double p2, double q0)
+{
+{
     double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
     double E,pi;
     double b2,b3,b4,q02,q03,q04,q05,Rg22;
     pi = 4.0*atan(1.0);
     E = 2.718281828459045091;
     if( L/b > 10.0) {
         C = 3.06/pow((L/b),0.44);
 …
         C = 1.0;
+    }
     C1 = 1.22;
     C2 = 0.4288;
 …
     C5 = 0.1477;
     miu = 0.585;
     Rg2 = Rgsquare(q,L,b);
     Rg22 = Rg2*Rg2;
 …
     q04 = q03*q0;
     q05 = q04*q0;
     t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
     t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
     t3 = ((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))));
     t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
     t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
     t6 = (b*C*(((-((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.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));
     t7 = (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));
     t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
     t9 = (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))))));
     t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
     t11 = (((-((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));
     t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
     t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));
     t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));
     t15 = ((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))));
     yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));
     return (yy);
+}
 …
 //
 //     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
 //                       BY (53) AND (56,57) IN CHEN AND
+//                       BY (53) AND (56,57) IN CHEN AND
 //                       KOTLARCHYK REFERENCE
 //
 …
         Pi = 4.0*atan(1.0);
         pi43=4.0/3.0*Pi;
         aa = crmaj;
 …
         gfn = gfnc+gfnt;
         gfn2 = gfn*gfn;
         return (gfn2);
+}
 …
 //
 //       <OBLATE ELLIPSOID>
 // function gfn4 for oblate ellipsoids
+// function gfn4 for oblate ellipsoids
 double
 gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
 …
         tgfn = gfnc+gfnt;
         gfn4 = tgfn*tgfn;
         return (gfn4);
+}
 …
     double Pq,vcyl,dl;
     double Pi,qr;
     Pi = 4.0*atan(1.0);
     qr = q*radius;
     Pq = Sk_WR(q,zi,lb);                //does not have cross section term
     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;
     dl = SchulzPoint_cpr(zi,length,zz);
     return (Pq*dl);
+    return (Pq*dl);
+}
 …
     double Pq,vcyl,dr;
     double Pi,qr;
     Pi = 4.0*atan(1.0);
     qr = q*zi;
     Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
     if (qr !=0){
 …
     vcyl=Pi*zi*zi*Lc;
     Pq *= vcyl*vcyl;
     dr = SchulzPoint_cpr(zi,ravg,zz);
     return (Pq*dr);
+    return (Pq*dr);
+}
 …
         double lolim,uplim,summ,yyy,zi;
         int nord,i;
         // set up the integration end points
+        // set up the integration end points
         Pi = 4.0*atan(1.0);
         nord = 76;
         lolim = 0;
+        lolim = 0.0;
         uplim = Pi/2.0;
         halfheight = length/2.0;
         summ = 0.0;                             // initialize integral
         i=0;
 …
                 summ += yyy;
+        }
         // calculate value of integral to return
         answer = (uplim-lolim)/2.0*summ;
 …
 double
 CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
+{
+{
         // qq is the q-value for the calculation (1/A)
         // radius is the core radius of the cylinder (A)
         //  radthick and facthick are the radial and face layer thicknesses
         // rho(n) are the respective SLD's
         // length is the *Half* CORE-LENGTH of the cylinder
+        // 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,si1,si2,be1,be2,t1,t2,retval;
         double Pi;
         Pi = 4.0*atan(1.0);
+        Pi = 4.0*atan(1.0);
         dr1 = rhoc-rhos;
         dr2 = rhos-rhosolv;
         vol1 = Pi*rad*rad*(2.0*length);
         vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
         besarg1 = qq*rad*sin(dum);
         besarg2 = qq*(rad+radthick)*sin(dum);
 …
                 si1 = sin(sinarg1)/sinarg1;
+        }
         if (besarg2 == 0.0){
+        if (sinarg2 == 0.0){
                 si2 = 1.0;
         }else{
 …
         retval = ((t1+t2)*(t1+t2))*sin(dum);
         return (retval);
+}
 …
     double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
     double Pi;
     Pi = 4.0*atan(1.0);
     dr1 = rhoc-rhos;
     dr2 = rhos-rhosolv;
     vol1 = Pi*rcore*rcore*(2.0*length);
     vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
     besarg1 = qq*rcore*sin(dum);
     besarg2 = qq*(rcore+thick)*sin(dum);
 …
                 si1 = sin(sinarg1)/sinarg1;
+        }
         if (besarg2 == 0.0){
+        if (sinarg2 == 0.0){
                 si2 = 1.0;
         }else{
 …
     retval = ((t1+t2)*(t1+t2))*sin(dum);
     return (retval);
+}
 …
 Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
+{
     double halfheight,uplim,lolim,zi,summ,yyy,Pi;
     double answer,dr,Vcyl;
     int i,nord;
     Pi = 4.0*atan(1.0);
     lolim = 0;
+    lolim = 0.0;
     uplim = Pi/2.0;
     halfheight = dumLen/2.0;
     nord=20;
     summ = 0.0;
     //do the cylinder orientational average
     for(i=0;i<nord;i++) {
 …
     Vcyl = Pi*radius*radius*dumLen;
     answer *= Vcyl*Vcyl;
     dr = SchulzPoint_cpr(dumLen,len_avg,zz);
     return(dr*answer);
 …
 double
 Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
+{
+{
         // qq is the q-value for the calculation (1/A)
         // rcore is the core radius of the cylinder (A)
 …
         double Pi;
         int kk;
         Pi = 4.0*atan(1.0);
         dr1 = rhoc-rhosolv;
         dr2 = rhol-rhosolv;
         area = Pi*rcore*rcore;
         totald=2.0*(thick+length);
         besarg1 = qq*rcore*sin(dum);
         besarg2 = qq*rcore*sin(dum);
         sinarg1 = qq*length*cos(dum);
         sinarg2 = qq*(length+thick)*cos(dum);
 …
                 si1 = sin(sinarg1)/sinarg1;
+        }
         if (besarg2 == 0.0){
+        if (sinarg2 == 0.0){
                 si2 = 1.0;
         }else{
 …
         retval =((t1+t2)*(t1+t2))*sin(dum);
         // loop for the structure facture S(q)
         sqq=0.0;
 …
                 dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
                 sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
+        }
+        }
         // end of loop for S(q)
         sqq=1.0+2.0*sqq/N;
         retval *= sqq;
         return(retval);
+}
 …
 Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
+{
     double halfheight,uplim,lolim,zi,summ,yyy,Pi;
     double answer,dr,Vcyl;
     int i,nord;
     Pi = 4.0*atan(1.0);
     lolim = 0;
+    lolim = 0.0;
     uplim = Pi/2.0;
     halfheight = length/2.0;
 …
     nord=76;
     summ = 0.0;
     //do the cylinder orientational average
         //    for(i=0;i<nord;i++) {
 …
     Vcyl = Pi*dumRad*dumRad*length;
     answer *= Vcyl*Vcyl;
     dr = SchulzPoint_cpr(dumRad,radius,zz);
     return(dr*answer);
 …
+{
     double dr;
     dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
     return(exp(dr));
 …
 .1208650973866179e-2,-0.5395239384953e-5};
         int j;
         y=x=xx;
         tmp=x+5.5;
 …
+{
     double arg,nu,retval;               //local variables
     nu = nua/a;
     arg = qq*a*sqrt(1+dum*dum*(nu*nu-1));
+    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0));
     if (arg == 0.0){
         retval =1.0/3.0;
 …
+{
     double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
     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)
+    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;
 …
     retval = psi*psi*t2*t2;
     return(retval);
 }//Function HolCylKernel()
 …
     // mu passed in is really mu*sqrt(1-sig^2)
     double arg1,arg2,Pi,tmp1,tmp2;                      //local variables
     Pi = 4.0*atan(1.0);
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     arg1 = (mu/2.0)*cos(Pi*uu/2.0);
 …
                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
+    }
     return (tmp1*tmp2);
 }//Function PPKernel()
 …
 TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
+{
     double arg,val,pi;                  //local variables
     pi = 4.0*atan(1.0);
     arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
     arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
 …
+    }
     return (val);
 }//Function TriaxialKernel()
 …
 CylKernel(double qq, double rr,double h, double theta)
+{
         // qq is the q-value for the calculation (1/A)
         // rr is the radius of the cylinder (A)
 …
     siarg = qq * h * cos(theta);
     bj =NR_BessJ1(besarg);
         //* Computing 2nd power */
     d1 = sin(siarg);
 …
     return (retval);
 }//Function CylKernel()
 …
 EllipCylKernel(double qq, double ra,double nu, double theta)
+{
         //this is the function LAMBDA1^2 in Feigin's notation
+        //this is the function LAMBDA1^2 in Feigin's notation
         // qq is the q-value for the calculation (1/A)
         // ra is the transformed radius"a" in Feigin's notation
         // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
         // theta is the dummy variable of the integration
         double retval,arg;               //Local variables
         arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2);
+        double retval,arg;               //Local variables
+        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2);
         if (arg == 0.0){
                 retval = 1.0;
 …
         //square it
         retval *= retval;
     return(retval);
 }//Function EllipCylKernel()
 …
         double ax,z;
         double xx,y,ans,ans1,ans2;
         if ((ax=fabs(x)) < 8.0) {
                 y=x*x;
 …
                 if (x < 0.0) ans = -ans;
+        }
         return(ans);
+}
+/*      Lamellar_ParaCrystal - Pedersen's model
+*/
+double
+Lamellar_ParaCrystal(double w[], double q)
+{
+//       Input (fitting) variables are:
+        //[0] scale factor
+        //[1]   thickness
+        //[2]   number of layers
+        //[3]   spacing between layers
+        //[4]   polydispersity of spacing
+        //[5] SLD lamellar
+        //[6] SLD solvent
+        //[7] incoherent background
+//      give them nice names
+        double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn;
+        double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi;
+        long n1,n2;
+        pi = 4.0*atan(1.0);
+        scale = w[0];
+        th = w[1];
+        nl = w[2];
+        davg = w[3];
+        pd = w[4];
+        sldLayer = w[5];
+        sldSolvent = w[6];
+        bkg = w[7];
+        contr = w[5] - w[6];
+        qval = q;
+        //get the fractional part of nl, to determine the "mixing" of N's
+        n1 = trunc(nl);         //rounds towards zero
+        n2 = n1 + 1;
+        xn = (double)n2 - nl;                   //fractional contribution of n1
+        ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
+        //calculate the n1 contribution
+        an = paraCryst_an(ww,qval,davg,n1);
+        Snq = paraCryst_sn(ww,qval,davg,n1,an);
+        Znq = xn*Snq;
+        //calculate the n2 contribution
+        an = paraCryst_an(ww,qval,davg,n2);
+        Snq = paraCryst_sn(ww,qval,davg,n2,an);
+        Znq += (1.0-xn)*Snq;
+        //and the independent contribution
+        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
+        //the limit when NL approaches infinity
+//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
+        xi = th/2.0;            //use 1/2 the bilayer thickness
+        Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
+        inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval);
+        inten *= 1.0e8;
+        return(scale*inten+bkg);
+}
+// functions for the lamellar paracrystal model
+double
+paraCryst_sn(double ww, double qval, double davg, long nl, double an) {
+        double Snq;
+        Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
+        return(Snq);
+}
+double
+paraCryst_an(double ww, double qval, double davg, long nl) {
+        double an;
+        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
+        an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg);
+        an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg);
+        an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg);
+        return(an);
+}
+/*      Spherocylinder  :
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+Spherocylinder(double w[], double x)
+{
+        int i,j;
+        double Pi;
+        double scale,contr,bkg,sldc,slds;
+        double len,rad,hDist,endRad;
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,zi,yyy,answer;                      //running tally of integration
+        double summj,vaj,vbj,zij;                       //for the inner integration
+        double SphCyl_tmp[7],arg1,arg2,inner,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        sldc = w[3];
+        slds = w[4];
+        bkg = w[5];
+        SphCyl_tmp[0] = w[0];
+        SphCyl_tmp[1] = w[1];
+        SphCyl_tmp[2] = w[2];
+        SphCyl_tmp[3] = w[1];           //end radius is same as cylinder radius
+        SphCyl_tmp[4] = w[3];
+        SphCyl_tmp[5] = w[4];
+        SphCyl_tmp[6] = w[5];
+        hDist = 0;              //by definition for this model
+        endRad = rad;
+        contr = sldc-slds;
+        Pi = 4.0*atan(1.0);
+        va = 0.0;
+        vb = Pi/2.0;            //orintational average, outer integral
+        vaj = -1.0*hDist/endRad;
+        vbj = 1.0;              //endpoints of inner integral
+        summ = 0.0;                     //initialize intergral
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
+                for(j=0;j<nordj;j++) {
+                        //20 gauss points for the inner integral
+                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
+                        yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi);
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                inner = (vbj-vaj)/2.0*summj;
+                inner *= 4.0*Pi*endRad*endRad*endRad;
+                //now calculate outer integrand
+                arg1 = x*len/2.0*cos(zi);
+                arg2 = x*rad*sin(zi);
+                yyy = inner;
+                if(arg2 == 0) {
+                        be = 0.5;
+                } else {
+                        be = NR_BessJ1(arg2)/arg2;
+                }
+                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*be;
+                } else {
+                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+                }
+                yyy *= yyy;
+                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
+                yyy *= Gauss76Wt[i];
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0;             //divide by volume
+        answer *= 1.0e8;                //convert to cm^-1
+        answer *= contr*contr;
+        answer *= scale;
+        answer += bkg;
+        return answer;
+}
+// inner integral of the sphereocylinder model, special case of hDist = 0
+//
+double
+SphCyl_kernel(double w[], double x, double tt, double theta) {
+        double val,arg1,arg2;
+        double scale,bkg,sldc,slds;
+        double len,rad,hDist,endRad,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        endRad = w[3];
+        sldc = w[4];
+        slds = w[5];
+        bkg = w[6];
+        hDist = 0.0;
+        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
+        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
+        if(arg2 == 0) {
+                be = 0.5;
+        } else {
+                be = NR_BessJ1(arg2)/arg2;
+        }
+        val = cos(arg1)*(1.0-tt*tt)*be;
+        return(val);
+}
+/*      Convex Lens  : special case where L ~ 0 and hDist < 0
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+ConvexLens(double w[], double x)
+{
+        int i,j;
+        double Pi;
+        double scale,contr,bkg,sldc,slds;
+        double len,rad,hDist,endRad;
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,zi,yyy,answer;                      //running tally of integration
+        double summj,vaj,vbj,zij;                       //for the inner integration
+        double CLens_tmp[7],arg1,arg2,inner,hh,be;
+        scale = w[0];
+        rad = w[1];
+//      len = w[2]
+        endRad = w[2];
+        sldc = w[3];
+        slds = w[4];
+        bkg = w[5];
+        len = 0.01;
+        CLens_tmp[0] = w[0];
+        CLens_tmp[1] = w[1];
+        CLens_tmp[2] = len;                     //length is some small number, essentially zero
+        CLens_tmp[3] = w[2];
+        CLens_tmp[4] = w[3];
+        CLens_tmp[5] = w[4];
+        CLens_tmp[6] = w[5];
+        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
+        contr = sldc-slds;
+        Pi = 4.0*atan(1.0);
+        va = 0.0;
+        vb = Pi/2.0;            //orintational average, outer integral
+        vaj = -1.0*hDist/endRad;
+        vbj = 1.0;              //endpoints of inner integral
+        summ = 0.0;                     //initialize intergral
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
+                for(j=0;j<nordj;j++) {
+                        //20 gauss points for the inner integral
+                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
+                        yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi);
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                inner = (vbj-vaj)/2.0*summj;
+                inner *= 4.0*Pi*endRad*endRad*endRad;
+                //now calculate outer integrand
+                arg1 = x*len/2.0*cos(zi);
+                arg2 = x*rad*sin(zi);
+                yyy = inner;
+                if(arg2 == 0) {
+                        be = 0.5;
+                } else {
+                        be = NR_BessJ1(arg2)/arg2;
+                }
+                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*be;
+                } else {
+                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+                }
+                yyy *= yyy;
+                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
+                yyy *= Gauss76Wt[i];
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        hh = fabs(hDist);               //need positive value for spherical cap volume
+        answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));             //divide by volume
+        answer *= 1.0e8;                //convert to cm^-1
+        answer *= contr*contr;
+        answer *= scale;
+        answer += bkg;
+        return answer;
+}
+/*      Capped Cylinder  : special case where L is nonzero and hDist < 0
+ -- uses the same Kernel as the Convex Lens
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+CappedCylinder(double w[], double x)
+{
+        int i,j;
+        double Pi;
+        double scale,contr,bkg,sldc,slds;
+        double len,rad,hDist,endRad;
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,zi,yyy,answer;                      //running tally of integration
+        double summj,vaj,vbj,zij;                       //for the inner integration
+        double arg1,arg2,inner,hh,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        endRad = w[3];
+        sldc = w[4];
+        slds = w[5];
+        bkg = w[6];
+        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
+        contr = sldc-slds;
+        Pi = 4.0*atan(1.0);
+        va = 0.0;
+        vb = Pi/2.0;            //orintational average, outer integral
+        vaj = -1.0*hDist/endRad;
+        vbj = 1.0;              //endpoints of inner integral
+        summ = 0.0;                     //initialize intergral
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
+                for(j=0;j<nordj;j++) {
+                        //20 gauss points for the inner integral
+                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
+                        yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi);               //uses the same kernel as ConvexLens, except here L != 0
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                inner = (vbj-vaj)/2.0*summj;
+                inner *= 4.0*Pi*endRad*endRad*endRad;
+                //now calculate outer integrand
+                arg1 = x*len/2.0*cos(zi);
+                arg2 = x*rad*sin(zi);
+                yyy = inner;
+                if(arg2 == 0) {
+                        be = 0.5;
+                } else {
+                        be = NR_BessJ1(arg2)/arg2;
+                }
+                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*be;
+                } else {
+                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+                }
+                yyy *= yyy;
+                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
+                yyy *= Gauss76Wt[i];
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        hh = fabs(hDist);               //need positive value for spherical cap volume
+        answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));            //divide by volume
+        answer *= 1.0e8;                //convert to cm^-1
+        answer *= contr*contr;
+        answer *= scale;
+        answer += bkg;
+        return answer;
+}
+// inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0
+//
+double
+ConvLens_kernel(double w[], double x, double tt, double theta) {
+        double val,arg1,arg2;
+        double scale,bkg,sldc,slds;
+        double len,rad,hDist,endRad,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        endRad = w[3];
+        sldc = w[4];
+        slds = w[5];
+        bkg = w[6];
+        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));
+        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
+        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
+        if(arg2 == 0) {
+                be = 0.5;
+        } else {
+                be = NR_BessJ1(arg2)/arg2;
+        }
+        val = cos(arg1)*(1.0-tt*tt)*be;
+        return(val);
+}
+/*      Dumbbell  : special case where L ~ 0 and hDist > 0
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+Dumbbell(double w[], double x)
+{
+        int i,j;
+        double Pi;
+        double scale,contr,bkg,sldc,slds;
+        double len,rad,hDist,endRad;
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,zi,yyy,answer;                      //running tally of integration
+        double summj,vaj,vbj,zij;                       //for the inner integration
+        double Dumb_tmp[7],arg1,arg2,inner,be;
+        scale = w[0];
+        rad = w[1];
+//      len = w[2]
+        endRad = w[2];
+        sldc = w[3];
+        slds = w[4];
+        bkg = w[5];
+        len = 0.01;
+        Dumb_tmp[0] = w[0];
+        Dumb_tmp[1] = w[1];
+        Dumb_tmp[2] = len;              //length is some small number, essentially zero
+        Dumb_tmp[3] = w[2];
+        Dumb_tmp[4] = w[3];
+        Dumb_tmp[5] = w[4];
+        Dumb_tmp[6] = w[5];
+        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
+        contr = sldc-slds;
+        Pi = 4.0*atan(1.0);
+        va = 0.0;
+        vb = Pi/2.0;            //orintational average, outer integral
+        vaj = -1.0*hDist/endRad;
+        vbj = 1.0;              //endpoints of inner integral
+        summ = 0.0;                     //initialize intergral
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
+                for(j=0;j<nordj;j++) {
+                        //20 gauss points for the inner integral
+                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
+                        yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi);
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                inner = (vbj-vaj)/2.0*summj;
+                inner *= 4.0*Pi*endRad*endRad*endRad;
+                //now calculate outer integrand
+                arg1 = x*len/2.0*cos(zi);
+                arg2 = x*rad*sin(zi);
+                yyy = inner;
+                if(arg2 == 0) {
+                        be = 0.5;
+                } else {
+                        be = NR_BessJ1(arg2)/arg2;
+                }
+                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*be;
+                } else {
+                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+                }
+                yyy *= yyy;
+                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
+                yyy *= Gauss76Wt[i];
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);              //divide by volume
+        answer *= 1.0e8;                //convert to cm^-1
+        answer *= contr*contr;
+        answer *= scale;
+        answer += bkg;
+        return answer;
+}
+/*      Barbell  : "normal" case where L is nonzero 0 and hDist > 0
+-- uses the same kernel as the Dumbbell case
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+Barbell(double w[], double x)
+{
+        int i,j;
+        double Pi;
+        double scale,contr,bkg,sldc,slds;
+        double len,rad,hDist,endRad;
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,zi,yyy,answer;                      //running tally of integration
+        double summj,vaj,vbj,zij;                       //for the inner integration
+        double arg1,arg2,inner,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        endRad = w[3];
+        sldc = w[4];
+        slds = w[5];
+        bkg = w[6];
+        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
+        contr = sldc-slds;
+        Pi = 4.0*atan(1.0);
+        va = 0.0;
+        vb = Pi/2.0;            //orintational average, outer integral
+        vaj = -1.0*hDist/endRad;
+        vbj = 1.0;              //endpoints of inner integral
+        summ = 0.0;                     //initialize intergral
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
+                for(j=0;j<nordj;j++) {
+                        //20 gauss points for the inner integral
+                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
+                        yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi);           //uses the same Kernel as the Dumbbell, here L>0
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                inner = (vbj-vaj)/2.0*summj;
+                inner *= 4.0*Pi*endRad*endRad*endRad;
+                //now calculate outer integrand
+                arg1 = x*len/2.0*cos(zi);
+                arg2 = x*rad*sin(zi);
+                yyy = inner;
+                if(arg2 == 0) {
+                        be = 0.5;
+                } else {
+                        be = NR_BessJ1(arg2)/arg2;
+                }
+                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*be;
+                } else {
+                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+                }
+                yyy *= yyy;
+                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
+                yyy *= Gauss76Wt[i];
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);             //divide by volume
+        answer *= 1.0e8;                //convert to cm^-1
+        answer *= contr*contr;
+        answer *= scale;
+        answer += bkg;
+        return answer;
+}
+// inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0
+//
+// inner integral of the Barbell model if L is nonzero
+//
+double
+Dumb_kernel(double w[], double x, double tt, double theta) {
+        double val,arg1,arg2;
+        double scale,bkg,sldc,slds;
+        double len,rad,hDist,endRad,be;
+        scale = w[0];
+        rad = w[1];
+        len = w[2];
+        endRad = w[3];
+        sldc = w[4];
+        slds = w[5];
+        bkg = w[6];
+        hDist = sqrt(fabs(endRad*endRad-rad*rad));
+        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
+        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
+        if(arg2 == 0) {
+                be = 0.5;
+        } else {
+                be = NR_BessJ1(arg2)/arg2;
+        }
+        val = cos(arg1)*(1.0-tt*tt)*be;
+        return(val);
+}
+double PolyCoreBicelle(double dp[], double q)
+{
+        int i;
+        int nord = 20;
+        double scale, length, sigma, bkg, radius, radthick, facthick;
+        double rhoc, rhoh, rhor, rhosolv;
+        double answer, Vpoly;
+        double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy;
+        scale = dp[0];
+        radius = dp[1];
+        sigma = dp[2];                          //sigma is the standard mean deviation
+        length = dp[3];
+        radthick = dp[4];
+        facthick= dp[5];
+        rhoc = dp[6];
+        rhoh = dp[7];
+        rhor=dp[8];
+        rhosolv = dp[9];
+        bkg = dp[10];
+        Pi = 4.0*atan(1.0);
+        lolim = exp(log(radius)-(4.*sigma));
+        if (lolim<0.0) {
+                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
+        }
+        uplim = exp(log(radius)+(4.*sigma));
+        summ = 0.0;
+        Rsqrsumm = 0.0;
+        for(i=0;i<nord;i++) {
+                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
+                yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length);
+                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
+                summ += yyy;
+                Rsqrsumm += Rsqryyy;
+        }
+        answer = (uplim-lolim)/2.0*summ;
+        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
+        //normalize by average cylinder volume
+        Vpoly = Pi*Rsqr*(length+2*facthick);
+        answer /= Vpoly;
+        //convert to [cm-1]
+        answer *= 1.0e8;
+        //Scale
+        answer *= scale;
+        // add in the background
+        answer += bkg;
+        return answer;
+}
+double
+BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){
+        double answer,halfheight,Pi;
+        double lolim,uplim,summ,yyy,zi;
+        int nord,i;
+        // set up the integration end points
+        Pi = 4.0*atan(1.0);
+        nord = 76;
+        lolim = 0.0;
+        uplim = Pi/2;
+        halfheight = length/2.0;
+        summ = 0.0;                             // initialize integral
+        i=0;
+        for(i=0;i<nord;i++) {
+                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
+                yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi);
+                summ += yyy;
+        }
+        // calculate value of integral to return
+        answer = (uplim-lolim)/2.0*summ;
+        return(answer);
+}
+double
+BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum)
+{
+        double dr1,dr2,dr3;
+        double besarg1,besarg2;
+        double vol1,vol2,vol3;
+        double sinarg1,sinarg2;
+        double t1,t2,t3;
+        double retval,si1,si2,be1,be2;
+        double Pi = 4.0*atan(1.0);
+        dr1 = rhoc-rhoh;
+        dr2 = rhor-rhosolv;
+        dr3=  rhoh-rhor;
+        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);
+        sinarg1 = qq*length*cos(dum);
+        sinarg2 = qq*(length+facthick)*cos(dum);
+        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.0;
+        } else {
+                si1 = sin(sinarg1)/sinarg1;
+        }
+        if(sinarg2 == 0) {
+                si2 = 1.0;
+        } else {
+                si2 = sin(sinarg2)/sinarg2;
+        }
+        t1 = 2.0*vol1*dr1*si1*be1;
+        t2 = 2.0*vol2*dr2*si2*be2;
+        t3 = 2.0*vol3*dr3*si2*be1;
+        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
+        return(retval);
+}
+double
+CSPPKernel(double dp[], double mu, double uu)
+{
+        double aa,bb,cc, ta,tb,tc;
+        double Vin,Vot,V1,V2;
+        double rhoA,rhoB,rhoC, rhoP, rhosolv;
+        double dr0, drA,drB, drC;
+        double arg1,arg2,arg3,arg4,t1,t2, t3, t4;
+        double Pi,retVal;
+        aa = dp[1];
+        bb = dp[2];
+        cc = dp[3];
+        ta = dp[4];
+        tb = dp[5];
+        tc = dp[6];
+        rhoA=dp[7];
+        rhoB=dp[8];
+        rhoC=dp[9];
+        rhoP=dp[10];
+        rhosolv=dp[11];
+        dr0=rhoP-rhosolv;
+        drA=rhoA-rhosolv;
+        drB=rhoB-rhosolv;
+        drC=rhoC-rhosolv;
+        Vin=(aa*bb*cc);
+        Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
+        V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
+        V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+        aa = aa/bb;
+        ta=(aa+2.0*ta)/bb;
+        tb=(aa+2.0*tb)/bb;
+        Pi = 4.0*atan(1.0);
+        arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0);
+        arg2 = (mu/2.0)*cos(Pi*uu/2.0);
+        arg3=  (mu*ta/2.0)*sin(Pi*uu/2.0);
+        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.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.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.0){
+                t3 = 1.0;
+        } else {
+                t3 = sin(arg3)/arg3;
+        }
+        if(arg4==0.0){
+                t4 = 1.0;
+        } else {
+                t4 = sin(arg4)/arg4;
+        }
+        retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors
+        return(retVal);
+}
+/*      CSParallelepiped  :  calculates the form factor of a Parallelepiped with a core-shell structure
+ -- different SLDs can be used for the face and rim
+Uses 76 pt Gaussian quadrature for both integrals
+*/
+double
+CSParallelepiped(double dp[], double q)
+{
+        int i,j;
+        double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg;         //local variables of coefficient wave
+        int nordi=76;                   //order of integration
+        int nordj=76;
+        double va,vb;           //upper and lower integration limits
+        double summ,yyy,answer;                 //running tally of integration
+        double summj,vaj,vbj;                   //for the inner integration
+        double mu,mudum,arg,sigma,uu,vol;
+        //      Pi = 4.0*atan(1.0);
+        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
+        scale = dp[0];
+        aa = dp[1];
+        bb = dp[2];
+        cc = dp[3];
+        ta  = dp[4];
+        tb  = dp[5];
+        tc  = dp[6];   // is 0 at the moment
+        rhoA=dp[7];   //rim A SLD
+        rhoB=dp[8];   //rim B SLD
+        rhoC=dp[9];    //rim C SLD
+        rhoP = dp[10];   //Parallelpiped core SLD
+        rhosolv=dp[11];  // Solvent SLD
+        bkg = dp[12];
+        mu = q*bb;
+        vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc;          //calculate volume before rescaling
+        // do the rescaling here, not in the kernel
+        // normalize all WRT bb
+        aa = aa/bb;
+        cc = cc/bb;
+        for(i=0;i<nordi;i++) {
+                //setup inner integral over the ellipsoidal cross-section
+                summj=0.0;
+                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
+                for(j=0;j<nordj;j++) {
+                        //76 gauss points for the inner integral
+                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
+                        mudum = mu*sqrt(1.0-sigma*sigma);
+                        yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu);
+                        summj += yyy;
+                }
+                //now calculate the value of the inner integral
+                answer = (vbj-vaj)/2.0*summj;
+                //finish the outer integral cc already scaled
+                arg = mu*cc*sigma/2.0;
+                if ( arg == 0.0 ) {
+                        answer *= 1.0;
+                } else {
+                        answer *= sin(arg)*sin(arg)/arg/arg;
+                }
+                //now sum up the outer integral
+                yyy = Gauss76Wt[i] * answer;
+                summ += yyy;
+        }               //final scaling is done at the end of the function, after the NT_FP64 case
+        answer = (vb-va)/2.0*summ;
+        //normalize by volume
+        answer /= vol;
+        //convert to [cm-1]
+        answer *= 1.0e8;
+        //Scale
+        answer *= scale;
+        // add in the background
+        answer += bkg;
+        return answer;
+}

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 6e93a02 in sasview for sansmodels/src/libigor/libCylinder.c

Legend:

sansmodels/src/libigor/libCylinder.c

Download in other formats: