-                      rdd71228
+                      r54bcd4a
 static double form_volume(
     int n_shells,
+    double radius_core,
+    double thick_inter0,
+    double thick_flat[],
+    double thick_inter[])
+    double thickness[],
+    double interface[])
+{
+    int i;
+    double r = radius_core;
+    r += thick_inter0;
+    for (i=0; i < n_shells; i++) {
+        r += thick_inter[i];
+        r += thick_flat[i];
+    double r = 0.0;
+    for (int i=0; i < n_shells; i++) {
+        r += thickness[i] + interface[i];
+    }
     return M_4PI_3*cube(r);
+}
+static double blend(int shape, double nu, double z)
+{
+    if (shape==0) {
+        const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0));
+        const double denom = 2.0 * sas_erf(nu * M_SQRT1_2);
+        return num/denom + 0.5;
+    } else if (shape==1) {
+        return pow(z, nu);
+    } else if (shape==2) {
+        return 1.0 - pow(1. - z, nu);
+    } else if (shape==3) {
+        return expm1(-nu*z)/expm1(-nu);
+    } else if (shape==4) {
+        return expm1(nu*z)/expm1(nu);
+    } else {
+        return NAN;
+    }
+}
+static double sphere_sld_kernel(
+static double f_linear(double q, double r, double contrast, double slope)
+{
+    const double qr = q * r;
+    const double qrsq = qr * qr;
+    const double bes = sph_j1c(qr);
+    double sinqr, cosqr;
+    SINCOS(qr, sinqr, cosqr);
+    const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
+    const double vol = M_4PI_3 * cube(r);
+    return vol*(bes*contrast + fun*slope);
+}
+static double Iq(
     double q,
     int n_shells,
-    int npts_inter,
-    double radius_core,
-    double sld_core,
     double sld_solvent,
+    double func_inter_core,
+    double thick_inter_core,
+    double nu_inter_core,
+    double sld_flat[],
+    double thick_flat[],
+    double func_inter[],
+    double thick_inter[],
+    double nu_inter[] ) {
+    double sld[],
+    double thickness[],
+    double interface[],
+    double shape[],
+    double nu[],
+    int n_steps)
+{
+    // iteration for # of shells + core + solvent
+    double f=0.0;
+    double r=0.0;
+    for (int shell=0; shell<n_shells; shell++){
+        const double sld_l = sld[shell];
+    int i,j,k;
+    int n_s;
+        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
+        f -= M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
+        r += thickness[shell];
+        f += M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
+    double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2;
+    double sign,slope=0.0;
+    double pi;
+        // iterate over sub_shells in the interface
+        const double dr = interface[shell]/n_steps;
+        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l;
+        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14);
+        const int shape_shell = (int)(shape[shell]);
+    double total_thick=0.0;
+        // if there is no interface the equations don't work
+        if (dr == 0.) continue;
+    //TODO: This part can be further simplified
+    int fun_type[12];
+    double sld[12];
+    double thick_internal[12];
+    double thick[12];
+    double fun_coef[12];
+        double sld_in = sld_l;
+        for (int step=1; step <= n_steps; step++) {
+            // find sld_i at the outer boundary of sub-shell step
+            //const double z = (double)step/(double)n_steps;
+            const double z = (double)step/(double)n_steps;
+            const double fraction = blend(shape_shell, nu_shell, z);
+            const double sld_out = fraction*delta + sld_l;
+            // calculate slope
+            const double slope = (sld_out - sld_in)/dr;
+            const double contrast = sld_in - slope*r;
+    fun_type[0] = func_inter_core;
+    fun_coef[0] = fabs(nu_inter_core);
+    sld[0] = sld_core;
+    thick[0] = radius_core;
+    thick_internal[0] = thick_inter_core;
+    for (i =1; i<=n_shells; i++){
+        sld[i] = sld_flat[i-1];
+        thick_internal[i]= thick_inter[i-1];
+        thick[i] = thick_flat[i-1];
+        fun_type[i] = func_inter[i-1];
+        fun_coef[i] = fabs(nu_inter[i-1]);
+        total_thick += thick[i];
+        total_thick += thick_internal[i]; //doesn't account for core layer
+    }
+    sld[n_shells+1] = sld_solvent;
+    thick[n_shells+1] = total_thick/5.0;
+    thick_internal[n_shells+1] = 0.0;
+    fun_coef[n_shells+1] = 0.0;
+    fun_type[n_shells+1] = 0;
+    pi = 4.0*atan(1.0);
+    f = 0.0;
+    r = 0.0;
+    vol = 0.0;
+    sld_f = sld_core;
+    dz = 0.0;
+    // iteration for # of shells + core + solvent
+    for (i=0;i<=n_shells+1; i++){
+        // iteration for flat and interface
+        for (j=0;j<2;j++){
+            // iteration for sub_shells in the interface
+            // starts from #1 sub-layer
+            for (n_s=1;n_s<=npts_inter; n_s++){
+                // for solvent, it doesn't have an interface
+                if (i==n_shells+1 && j==1)
+                    break;
+                // for flat layers
+                if (j==0){
+                    dz = thick[i];
+                    sld_i = sld[i];
+                    slope = 0.0;
+                }
+                // for interfacial sub_shells
+                else{
+                    dz = thick_internal[i]/npts_inter;
+                    // find sld_i at the outer boundary of sub-layer #n_s
+                    sld_i = intersldfunc(fun_type[i], npts_inter, n_s,
+                            fun_coef[i], sld[i], sld[i+1]);
+                    // calculate slope
+                    slope= (sld_i -sld_f)/dz;
+                }
+                contr = sld_f-slope*r;
+                // iteration for the left and right boundary of the shells
+                for (k=0; k<2; k++){
+                    // At r=0, the contribution to I is zero, so skip it.
+                    if ( i == 0 && j == 0 && k == 0){
+                        continue;
+                    }
+                    // On the top of slovent there is no interface; skip it.
+                    if (i == n_shells+1 && k == 1){
+                        continue;
+                    }
+                    // At the right side (outer) boundary
+                    if ( k == 1){
+                        sign = 1.0;
+                        r += dz;
+                    }
+                    // At the left side (inner) boundary
+                    else{
+                        sign = -1.0;
+                    }
+                    qr = q * r;
+                    fun = 0.0;
+                    if(qr == 0.0){
+                        bes = sign * 1.0;
+                    }
+                    else{
+                        // for flat sub-layer
+                        //TODO: Single precision calculation fails here
+                        bes = sign *  sph_j1c(qr);
+                        if (fabs(slope) > 0.0 ){
+                            const double qrsq = qr*qr;
+                            double sinqr, cosqr;
+                            SINCOS(qr, sinqr, cosqr);
+                            fun = sign * 3.0 * r *
+                            (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq);
+                            // In the onion model Jae-He's formula is rederived
+                            // and gives following:
+                            //fun = 3.0 * sign * r *
+                            //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
+                            //But this seems not to be working in this case...
+                        }
+                    }
+                    // update total volume
+                    vol = M_4PI_3 * cube(r);
+                    f += vol * (bes * contr + fun * slope);
+                }
+                sld_f = sld_i;
+                // no sub-layer iteration (n_s loop) for the flat layer
+                if (j==0)
+                    break;
+            }
+            // iteration for the left and right boundary of the shells
+            f -= f_linear(q, r, contrast, slope);
+            r += dr;
+            f += f_linear(q, r, contrast, slope);
+            sld_in = sld_out;
+        }
+    }
+    // add in solvent effect
+    f -= M_4PI_3 * cube(r) * sld_solvent * sph_j1c(q*r);
     f2 = f * f * 1.0e-4;
     return (f2);
+    const double f2 = f * f * 1.0e-4;
+    return f2;
+}
-/**
- * Function to evaluate 1D SphereSLD function
- * @param q: q-value
- * @return: function value
- */
-static double Iq(double q,
-    int n_shells,
-    int npts_inter,
-    double radius_core,
-    double sld_core,
-    double sld_solvent,
-    double func_inter0,
-    double thick_inter0,
-    double nu_inter0,
-    double sld_flat[],
-    double thick_flat[],
-    double func_inter[],
-    double thick_inter[],
-    double nu_inter[] ) {
-    double intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core,
-                sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0,
-                sld_flat, thick_flat, func_inter, thick_inter, nu_inter);
-    return intensity;
+}
-/**
- * Function to evaluate 2D SphereSLD function
- * @param q_x: value of Q along x
- * @param q_y: value of Q along y
- * @return: function value
- */
-/*static double Iqxy(double qx, double qy,
-    int n_shells,
-    int npts_inter,
-    double radius_core
-    double sld_core,
-    double sld_solvent,
-    double sld_flat[],
-    double thick_flat[],
-    double func_inter[],
-    double thick_inter[],
-    double nu_inter[],
-    ) {
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent,
-    sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[])
-}*/

sasmodels/models/spherical_sld.py

-                      r4fd2c63
+                      r54bcd4a
 import numpy as np
+from numpy import inf
+from numpy import inf, expm1, sqrt
+from scipy.special import erf
 name = "spherical_sld"
 …
 category = "shape:sphere"
+SHAPES = ["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)",
+          "Rexp(-|nu|z)", "Lexp(-|nu|z)"],
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
+parameters = [["n_shells",          "",               1,      [0, 10],        "volume", "number of shells"],
+              ["npts_inter",        "",               35,     [0, inf],       "", "number of points in each sublayer Must be odd number"],
+              ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
+              ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "sld", "sld function flat"],
+              ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "sld", "sld function solvent"],
+              ["func_inter0",       "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"],
+              ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"],
+              ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"],
+              ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "sld", "sld function flat"],
+              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"],
+              ["func_inter[n_shells]",    "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"],
+              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
+              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
+parameters = [["n_shells",             "",           1,      [1, 11],        "volume", "number of shells"],
+              ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"],
+              ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"],
+              ["thickness[n_shells]",  "Ang",        100.0,  [0, inf],       "volume", "thickness shell"],
+              ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"],
+              ["shape[n_shells]",      "",           0,      SHAPES,         "", "interface shape"],
+              ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"],
+              ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"],
+              ]
 # pylint: enable=bad-whitespace, line-too-long
 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"]
+source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sph_j1c.c", "spherical_sld.c"]
 single = False  # TODO: fix low q behaviour
 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
+def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat,
+            thick_flat, func_inter, thick_inter, nu_inter, npts_inter):
+SHAPE_FUNCTIONS = [
+    lambda z, nu: erf(nu/sqrt(2)*(2*z-1))/(2*erf(nu/sqrt(2))) + 0.5,  # erf
+    lambda z, nu: z**nu,                    # Rpow
+    lambda z, nu: 1 - (1-z)**nu,            # Lpow
+    lambda z, nu: expm1(-nu*z)/expm1(-nu),  # Rexp
+    lambda z, nu: expm1(nu*z)/expm1(nu),    # Lexp
+]
+def profile(n_shells, sld_solvent, sld, thickness,
+            interface, shape, nu, n_steps):
     """
     Returns shape profile with x=radius, y=SLD.
 …
     z = []
     beta = []
+    rho = []
     z0 = 0
     # two sld points for core
     z.append(0)
+    beta.append(sld_core)
+    z.append(radius_core)
+    beta.append(sld_core)
+    z0 += radius_core
+    for i in range(1, n_shells+2):
+        dz = thick_inter[i-1]/npts_inter
+        # j=0 for interface, j=1 for flat layer
+        for j in range(0, 2):
+            # interation for sub-layers
+            for n_s in range(0, npts_inter+1):
+                if j == 1:
+                    if i == n_shells+1:
+                        break
+                    # shift half sub thickness for the first point
+                    z0 -= dz#/2.0
+                    z.append(z0)
+                    #z0 -= dz/2.0
+                    z0 += thick_flat[i]
+                    sld_i = sld_flat[i]
+                    beta.append(sld_flat[i])
+                    dz = 0
+                else:
+                    nu = nu_inter[i-1]
+                    # decide which sld is which, sld_r or sld_l
+                    if i == 1:
+                        sld_l = sld_core
+                    else:
+                        sld_l = sld_flat[i-1]
+                    if i == n_shells+1:
+                        sld_r = sld_solvent
+                    else:
+                        sld_r = sld_flat[i]
+                    # get function type
+                    func_idx = func_inter[i-1]
+                    # calculate the sld
+                    sld_i = intersldfunc(func_idx, npts_inter, n_s, nu,
+                                            sld_l, sld_r)
+                # append to the list
+                z.append(z0)
+                beta.append(sld_i)
+                z0 += dz
+                if j == 1:
+                    break
+    z.append(z0)
+    beta.append(sld_solvent)
+    z_ext = z0/5.0
+    z.append(z0+z_ext)
+    beta.append(sld_solvent)
+    rho.append(sld[0])
+    for i in range(0, n_shells):
+        z0 += thickness[i]
+        z.append(z0)
+        rho.append(sld[i])
+        dz = interface[i]/n_steps
+        sld_l = sld[i]
+        sld_r = sld[i+1] if i < n_shells-1 else sld_solvent
+        interface = SHAPE_FUNCTIONS[int(np.clip(shape[i], 0, len(SHAPES)-1))]
+        for step in range(1, n_steps+1):
+            portion = interface(float(step)/n_steps, max(abs(nu[i]), 1e-14))
+            z0 += dz
+            z.append(z0)
+            rho.append((sld_r - sld_l)*portion + sld_l)
+    z.append(z0*1.2)
+    rho.append(sld_solvent)
     # return sld profile (r, beta)
+    return np.asarray(z), np.asarray(beta)*1e-6
+def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat):
+    return np.asarray(z), np.asarray(rho)*1e-6
+def ER(n_shells, thickness, interface):
     n_shells = int(n_shells)
+    total_thickness = thick_inter0
+    total_thickness += np.sum(thick_inter[:n_shells], axis=0)
+    total_thickness += np.sum(thick_flat[:n_shells], axis=0)
+    return total_thickness + radius_core
+    total = (np.sum(thickness[:n_shells], axis=1)
+             + np.sum(interface[:n_shells], axis=1))
+    return total
 demo = {
+    "n_shells": 4,
+    "npts_inter": 35.0,
+    "radius_core": 50.0,
+    "sld_core": 2.07,
+    "n_shells": 5,
+    "n_steps": 35.0,
     "sld_solvent": 1.0,
+    "thick_inter0": 50.0,
+    "func_inter0": 0,
+    "nu_inter0": 2.5,
+    "sld_flat":[4.0,3.5,4.0,3.5],
+    "thick_flat":[100.0,100.0,100.0,100.0],
+    "func_inter":[0,0,0,0],
+    "thick_inter":[50.0,50.0,50.0,50.0],
+    "nu_inter":[2.5,2.5,2.5,2.5],
+    "sld":[2.07,4.0,3.5,4.0,3.5],
+    "thickness":[50.0,100.0,100.0,100.0,100.0],
+    "interface":[50.0,50.0,50.0,50.0],
+    "shape": [0,0,0,0,0],
+    "nu":[2.5,2.5,2.5,2.5,2.5],
+    }
 …
 tests = [
     # Accuracy tests based on content in test/utest_extra_models.py
+    [{"n_shells":4,
+        'npts_inter':35,
+        "radius_core":50.0,
+        "sld_core":2.07,
+    [{"n_shells": 5,
+        "n_steps": 35,
         "sld_solvent": 1.0,
         "sld_flat":[4.0,3.5,4.0,3.5],
         "thick_flat":[100.0,100.0,100.0,100.0],
         "func_inter":[0,0,0,0],
         "thick_inter":[50.0,50.0,50.0,50.0],
         "nu_inter":[2.5,2.5,2.5,2.5]
+        "sld": [2.07, 4.0, 3.5, 4.0, 3.5],
+        "thickness": [50.0, 100.0, 100.0, 100.0, 100.0],
+        "interface": [50]*5,
+        "shape": [0]*5,
+        "nu": [2.5]*5,
     }, 0.001, 0.001],
+]

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

SasView

Changeset 54bcd4a in sasmodels for sasmodels/models

Legend:

sasmodels/models/spherical_sld.c

sasmodels/models/spherical_sld.py

Download in other formats: