Changeset c5ac2b2 in sasmodels

Timestamp:

Jul 18, 2016 2:05:12 AM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

4a21b121

Parents:

256dfe1

Message:

working core multi shell

Location:

sasmodels

Files:

• convert.py (modified) (3 diffs)
• modelinfo.py (modified) (1 diff)
• models/core_multi_shell.c (modified) (2 diffs)
• models/core_multi_shell.py (modified) (6 diffs)

sasmodels/convert.py

@@ -5 +5 @@
 
 from os.path import join as joinpath, abspath, dirname
+import math
 import warnings
 import json
@@ -34 +35 @@
     'vesicle',
     'multilayer_vesicle',
-    'core_multi_shell',
 ]
 
@@ -287 +287 @@
         elif name == 'poly_gauss_coil':
             pars['i_zero'] = 1
-
+        elif name == 'core_multi_shell':
+            pars['n'] = max(math.ceil(pars['n']), 4)
+

sasmodels/modelinfo.py

@@ -526 +526 @@
         Return the list of parameters for the given data type.
 
-        Vector parameters are expanded as in place.  If multiple parameters
+        Vector parameters are expanded in place.  If multiple parameters
         share the same vector length, then the parameters will be interleaved
         in the result.  The control parameters come first.  For example,

sasmodels/models/core_multi_shell.c

@@ -1 @@
-FourShell(double dp[], double q)
-{
-    // variables are:
-    //[0] scale factor
-    //[1] radius of core [ￜ]
-    //[2] SLD of the core   [ￜ-2]
-    //[3] thickness of shell 1 [ￜ]
-    //[4] SLD of shell 1
-    //[5] thickness of shell 2 [ￜ]
-    //[6] SLD of shell 2
-    //[7] thickness of shell 3
-    //[8] SLD of shell 3
-    //[9] thickness of shell 3
-    //[10] SLD of shell 3
-    //[11] SLD of solvent
-    //[12] background   [cm-1]
-    
-    double x,pi;
-    double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;     //my local names
-    double bes,f,vol,qr,contr,f2;
-    double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4;
-    
-    pi = 4.0*atan(1.0);
-    x=q;
-    
-    scale = dp[0];
-    rcore = dp[1];
-    rhocore = dp[2];
-    thick1 = dp[3];
-    rhoshel1 = dp[4];
-    thick2 = dp[5];
-    rhoshel2 = dp[6];   
-    thick3 = dp[7];
-    rhoshel3 = dp[8];
-    thick4 = dp[9];
-    rhoshel4 = dp[10];  
-    rhosolv = dp[11];
-    bkg = dp[12];
-    
-        // core first, then add in shells
-    qr=x*rcore;
-    contr = rhocore-rhoshel1;
-    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;
-    //now the shell (1)
-    qr=x*(rcore+thick1);
-    contr = rhoshel1-rhoshel2;
-    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+thick1)*(rcore+thick1)*(rcore+thick1);
-    f += vol*bes*contr;
-    //now the shell (2)
-    qr=x*(rcore+thick1+thick2);
-    contr = rhoshel2-rhoshel3;
-    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+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2);
-    f += vol*bes*contr;
-    //now the shell (3)
-    qr=x*(rcore+thick1+thick2+thick3);
-    contr = rhoshel3-rhoshel4;
-    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+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3);
-    f += vol*bes*contr;
-    //now the shell (4)
-    qr=x*(rcore+thick1+thick2+thick3+thick4);
-    contr = rhoshel4-rhosolv;
-    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+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4);
-    f += vol*bes*contr;
-    
-        
-    // normalize to particle volume and rescale from [ￜ-1] to [cm-1]
-    f2 = f*f/vol*1.0e8;
-    
-    //scale if desired
-    f2 *= scale;
-    // then add in the background
-    f2 += bkg;
-    
-    return(f2);
-}
-
-// above is cut and past with no changes from libigor four shell
-static double
-f_exp(double q, double r, double sld_in, double sld_out,
-    double thickness, double A)
-{
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
-  const double qr = q * r;
-  const double alpha = A * r/thickness;
-  const double bes = sph_j1c(qr);
-  const double B = (sld_out - sld_in)/expm1(A);
-  const double C = sld_in - B;
-  double fun;
-  if (qr == 0.0) {
-    fun = 1.0;
-  } else if (fabs(A) > 0.0) {
-    const double qrsq = qr*qr;
-    const double alphasq = alpha*alpha;
-    const double sumsq = alphasq + qrsq;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    fun = -3.0(
-            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq
-                - (alpha*sinqr/qr - cosqr)
-        ) / sumsq;
-  } else {
-    fun = bes;
-  }
-  return vol * (B*fun + C*bes);
-}
-
-static double
-f_linear(double q, double r, double sld, double slope)
-{
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
-  const double qr = q * r;
-  const double bes = sph_j1c(qr);
-  double fun = 0.0;
-  if (qr > 0.0) {
-    const double qrsq = qr*qr;
-    double sinqr, cosqr;
-    SINCOS(qr, sinqr, cosqr);
-    // Jae-He's code seems to simplify to this
-    //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
-    // Rederiving the math, we get the following instead:
-    fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
-  }
-  return vol * (sld*bes + fun);
-}
 
 static double
@@ -154 +4 @@
 {
   const double bes = sph_j1c(q * r);
-  const double vol = 4.0/3.0 * M_PI * r * r * r;
+  const double vol = M_4PI_3 * cube(r);
   return sld * vol * bes;
 }
 
-static double
+double
+form_volume(double core_radius, double n, double thickness[]);
+double
 form_volume(double core_radius, double n, double thickness[])
 {
-  int i;
   double r = core_radius;
-  for (i=0; i < n; i++) {
+  for (int i=0; i < n; i++) {
     r += thickness[i];
   }
-  return 4.0/3.0 * M_PI * r * r * r;
+  return M_4PI_3 * cube(r);
 }
 
-static double
-Iq(double q, double core_sld, double core_radius, double solvent_sld,
-    double n, double in_sld[], double out_sld[], double thickness[],
-    double A[])
+double
+Iq(double q, double core_sld, double core_radius,
+   double solvent_sld, double num_shells, double sld[], double thickness[]);
+double
+Iq(double q, double core_sld, double core_radius,
+   double solvent_sld, double num_shells, double sld[], double thickness[])
 {
-  int i;
+  const int n = (int)ceil(num_shells);
+  double f, r, last_sld;
   r = core_radius;
-  f = f_constant(q, r, core_sld);
-  for (i=0; i<n; i++){
-    const double r0 = r;
+  last_sld = core_sld;
+  f = 0.;
+  for (int i=0; i<n; i++) {
+    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sph_j1c(q*r);
+    last_sld = sld[i];
     r += thickness[i];
-    if (r == r0) {
-      // no thickness, so nothing to add
-    } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) {
-      f -= f_constant(q, r0, sld_in[i]);
-      f += f_constant(q, r, sld_in[i]);
-    } else if (fabs(A[i]) < 1.0e-4) {
-      const double slope = (sld_out[i] - sld_in[i])/thickness[i];
-      f -= f_linear(q, r0, sld_in[i], slope);
-      f += f_linear(q, r, sld_out[i], slope);
-    } else {
-      f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]);
-      f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]);
-    }
   }
-  f -= f_constant(q, r, solvent_sld);
-  f2 = f * f * 1.0e-4;
-
-  return f2;
+  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sph_j1c(q*r);
+  return f * f * 1.0e-4;
 }

sasmodels/models/core_multi_shell.py

@@ -87 +87 @@
 
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["volfraction", "", 0.05, [0,1],"",
-               "volume fraction of spheres"],
-              ["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
+parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "",
                "Core scattering length density"],
               ["radius", "Ang", 200., [0, inf], "volume",
@@ -95 +93 @@
               ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "",
                "Solvent scattering length density"],
-              ["n", "", 1, [0, 4], "volume",
+              ["n", "", 1, [0, 10], "volume",
                "number of shells"],
               ["sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "",
@@ -103 +101 @@
               ]
 
-#source = ["lib/sph_j1c.c", "onion.c"]
+source = ["lib/sph_j1c.c", "core_multi_shell.c"]
 
-def Iq(q, *args, **kw):
-    return q
-
-def Iqxy(qx, *args, **kw):
-    return qx
-
-
-def profile(volfraction, sld_core, radius, sld_solvent, n, sld, thickness):
+def profile(sld_core, radius, sld_solvent, n, sld, thickness):
     """
     SLD profile
@@ -120 +111 @@
 
     """
-#        r = []
-#        beta = []
-#        # for core at r=0
-#        r.append(0)
-#        beta.append(self.params['sld_core0'])
-#        # for core at r=rad_core
-#        r.append(self.params['rad_core0'])
-#        beta.append(self.params['sld_core0'])
-#
-#        # for shells
-#        for n in range(1, self.n_shells+1):
-#            # Left side of each shells
-#            r0 = r[len(r)-1]
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#            # Right side of each shells
-#            exec "r0 += self.params['thick_shell%s'% str(n)]"
-#            r.append(r0)
-#            exec "beta.append(self.params['sld_shell%s'% str(n)])"
-#
-#        # for solvent
-#        r0 = r[len(r)-1]            
-#        r.append(r0)
-#        beta.append(self.params['sld_solv'])
-#        r_solv = 5*r0/4
-#        r.append(r_solv)
-#        beta.append(self.params['sld_solv'])
-#
-#        return r, beta
-# above is directly from old code -- below is alotered from Woitek's first
-# cut an the onion.  Seems like we should be consistant?
-
     total_radius = 1.25*(sum(thickness[:n]) + radius + 1)
 
@@ -172 +130 @@
         r.append(r0 + thickness[k])
         beta.append(sld[k])
-   # add in the solvent
+    # add in the solvent
     r.append(r[-1])
     beta.append(sld_solvent)
@@ -187 +145 @@
     return 1.0, 1.0
 
-demo = dict(volfraction = 1.0,
-            sld_core = 6.4,
+demo = dict(sld_core = 6.4,
             radius = 60,
             sld_solvent = 6.4,
-            n = 1,
-            sld = [2.0],
-            thickness = [10])
+            n = 2,
+            sld = [2.0, 3.0],
+            thickness = 20,
+            thickness1_pd = 0.3,
+            thickness2_pd = 0.3,
+            thickness1_pd_n = 10,
+            thickness2_pd_n = 10,
+            )