Three-layer displaced core-shell ellipsoid - core_shell_extra_inner_shell.c

    /*
 * Core-shell ellipsoid with extra inner shell and polydispersity
 *
 * Three-layer ellipsoidal particle (normalised contrast, scale set by SASView):
 *
 *   Layer         Radial extent        Contrast (normalised to outer shell = 1)
 *   -----------   ------------------   ----------------------------------------
 *   Outer shell   R < r < R+d_shell    1
 *   Inner shell   R-d_in < r < R       contrast_inner_shell
 *   Core          r < R-d_in           contrast_core
 *
 * The effective outer shell extent is:
 *   d_shell = |thickness_shell| + 2*|sigma_outer|
 *
 * The inner shell occupies the radial band [R - |thickness_inner_shell|, R].
 * The displacement phase applies to the entire inner structure (inner shell
 * + core) rigidly displaced perpendicular to the symmetry axis.
 *
 * The scattering amplitude at a given size and orientation is decomposed as:
 *
 *   F = FFS + FFC * phase
 *
 *   FFS = V_RO * phi(q*r_RO) * DW_outer          [outer ellipsoid, contrast 1]
 *   FFC = [-V_R  * (1 - c_in) * phi(q*r_R)
 *          -V_RI * (c_in - c_core) * phi(q*r_RI)] * DW_core
 *
 * where c_in = contrast_inner_shell, c_core = contrast_core.
 *
 * Same normalization and decoupling approximation as core_shell_displ_core.
 * Fortran bugs corrected: SUM1X (not SUM1) used for amplitude accumulation;
 * F0 computed from angle-independent RO (not RAO).
 *
 * Translated from Fortran (Spinozzi et al., NXS=31).
 */


#define NSTP 50
#define NPOI 50

static double
hard_sphere_sf(double q, double eta, double r_hs)
{
    if (eta <= 0.0 || q <= 0.0) return 1.0;
    double aln   = (1.0-eta)*(1.0-eta)*(1.0-eta)*(1.0-eta);
    double alpha = (1.0+2.0*eta)*(1.0+2.0*eta) / aln;
    double beta  = -6.0*eta*(1.0+0.5*eta)*(1.0+0.5*eta) / aln;
    double gamma = 0.5*eta*alpha;
    double ar    = 2.0*r_hs*q + 1.0e-4;
    double sa = sin(ar), ca = cos(ar);
    double ar2=ar*ar, ar3=ar2*ar, ar4=ar3*ar, ar5=ar4*ar;
    double gg = alpha*(sa-ar*ca)/ar2
              + beta *(2.0*ar*sa+(2.0-ar2)*ca-2.0)/ar3
              + gamma*(-ar4*ca+4.0*((3.0*ar2-6.0)*ca
                       +(ar3-6.0*ar)*sa+6.0))/ar5;
    return 1.0 / (1.0 + 24.0*eta*gg/ar);
}

double
Iq(double q,
   double radius,
   double aspect_ratio,
   double volfraction_hs,
   double radius_hs,
   double thickness_shell,
   double contrast_core,
   double sigma_rel,
   double power_law,
   double exponent_2,
   double sigma_outer,
   double rg_polymer,
   double scale_polymer,
   double displacement,
   double sigma_core,
   double thickness_inner_shell,
   double contrast_inner_shell)
{
    const double pi  = M_PI;
    const double eps = aspect_ratio;

    double sw = fabs(sigma_rel) * radius;
    if (sw < 1.0e-4 * radius) sw = 1.0e-4 * radius;
    const double dw   = 6.0*sw / (double)NSTP;
    double rbeg = radius - 3.0*sw;
    if (rbeg < 0.0) rbeg = 0.0;

    const double dshell  = fabs(thickness_shell) + 2.0*fabs(sigma_outer);
    const double d_inner = fabs(thickness_inner_shell);
    const double step    = 0.5*pi / (double)NPOI;

    double displ = displacement;
    if (displ > dshell) displ = dshell - fabs(dshell - displ);

    double sca = 0.0, sum2 = 0.0, sum1 = 0.0, sum1n = 0.0;

    for (int jj = 0; jj < NSTP; jj++) {
        const double r     = rbeg + ((double)jj + 0.5)*dw;
        const double dr    = r - radius;
        const double f_exp = exp(-0.5*(dr/sw)*(dr/sw));

        /* Outer ellipsoid */
        const double ro   = r + dshell;
        const double epso = (r*eps + dshell) / ro;
        /* Inner ellipsoid (inner boundary of inner shell = core surface) */
        const double ri   = r - d_inner;
        const double epsi = (r*eps - d_inner) / ri;

        const double vol_ro = (4.0*pi/3.0)*ro*ro*ro*epso;
        const double vol_r  = (4.0*pi/3.0)*r *r *r *eps;
        const double vol_ri = (4.0*pi/3.0)*ri*ri*ri*epsi;

        /* Forward-scattering amplitude (q=0 limit, angle-independent) */
        const double f0 = vol_ro
                        - vol_r  * (1.0 - contrast_inner_shell)
                        - vol_ri * (contrast_inner_shell - contrast_core);

        double sumx = 0.0, sum1x = 0.0;

        for (int ii = 0; ii < NPOI; ii++) {
            double xx = ((double)ii + 0.5)*step;
            xx = sin(xx);

            const double sc  = sqrt(xx*xx + eps *eps *(1.0-xx*xx));
            const double sco = sqrt(xx*xx + epso*epso*(1.0-xx*xx));
            const double sci = sqrt(xx*xx + epsi*epsi*(1.0-xx*xx));
            const double ra  = r *sc;
            const double rao = ro*sco;
            const double rai = ri*sci;

            /* Outer shell: fills entire outer ellipsoid with unit contrast */
            const double ffs = vol_ro * sas_3j1x_x(q*rao)
                             * exp(-0.5*(q*sigma_outer)*(q*sigma_outer));

            /* Inner structure: inner shell + core, with shared DW roughness.
             * The two-term decomposition peels off each layer's extra contrast. */
            const double ffc = -(vol_r  * (1.0 - contrast_inner_shell) * sas_3j1x_x(q*ra)
                                + vol_ri * (contrast_inner_shell - contrast_core) * sas_3j1x_x(q*rai))
                             * exp(-0.5*(q*sigma_core)*(q*sigma_core));

            /* Displacement phase for the entire inner structure */
            const double phase   = cos(q*displ*sqrt(1.0-xx*xx));
            const double F_total = ffs + ffc*phase;

            sumx  += F_total*F_total*xx;
            sum1x += F_total*xx;
        }

        sca   += sumx *step*f_exp;
        sum1  += sum1x*step*f_exp;
        sum2  += f0*f0*f_exp;
        sum1n += f0*f_exp;
    }

    const double p_form   = (sum2  > 0.0) ? sca /sum2  : 0.0;
    const double f_rel_sq = (sum1n > 0.0) ? (sum1/sum1n)*(sum1/sum1n) : 0.0;

    const double sq = hard_sphere_sf(q, volfraction_hs, radius_hs);

    double i_poly = 0.0;
    if (scale_polymer != 0.0 && rg_polymer > 0.0 && q > 0.0) {
        const double u    = q*q*rg_polymer*rg_polymer;
        const double poly = (u > 0.1) ? 2.0*(exp(-u)-1.0+u)/(u*u) : 1.0-u/3.0;
        i_poly = scale_polymer*poly;
    }

    const double b_q = (q > 0.0)
                     ? 1.0 + power_law*pow(0.001/q, fabs(exponent_2)+2.0)
                     : 1.0;

    return (p_form + f_rel_sq*(sq-1.0))*b_q + i_poly;
}

double
form_volume(double radius, double aspect_ratio, double thickness_shell)
{
    const double r_outer   = radius + thickness_shell;
    const double eps_outer = (radius*aspect_ratio + thickness_shell) / r_outer;
    return 4.0*M_PI/3.0 * r_outer*r_outer*r_outer * eps_outer;
}

double
radius_effective(int mode, double radius, double aspect_ratio, double thickness_shell)
{
    (void)aspect_ratio;
    switch (mode) {
    case 1:  return radius + thickness_shell;
    default: return radius;
    }
}
Back to Model Download