Displaced core-shell ellipsoids - core_shell_displ_core.c

    /*
 * Displaced core-shell ellipsoid model with polydispersity
 *
 * Scattering from core-shell ellipsoids where the core centre is displaced
 * from the shell centre. Integrates over a Gaussian size distribution and
 * over all orientations using the decoupling approximation:
 *
 *   I(q) = scale * [ P(q) + <F(q)>^2 * (S(q) - 1) ] * B(q) + background + I_poly(q)
 *
 * where P(q) = <F^2(q)> / F_{\theta}^2 (normalised),  <F(q)>/F_{\theta} is the size- and
 * orientation-averaged amplitude,  S(q) is the Percus-Yevick hard-sphere
 * structure factor, and B(q) is an empirical power-law correction.
 *
 * Translated from Fortran (Spinozzi et al.).
 */



#define NSTP 50
#define NPOI 50

/* Percus-Yevick hard-sphere structure factor */
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_shell,
   double sigma_rel,
   double power_law,
   double exponent_2,
   double sigma_outer,
   double rg_polymer,
   double scale_polymer,
   double displacement,
   double sigma_core)
{
    const double pi  = M_PI;
    const double eps = aspect_ratio;

    /* Gaussian size distribution; clamp width away from zero */
    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;

    /* Effective shell thickness including surface roughness Debye-Waller term */
    const double dshell = fabs(thickness_shell) + 2.0*fabs(sigma_outer);
    const double step   = 0.5*pi / (double)NPOI;

    /* Clamp displacement so core stays inside shell */
    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 geometry (angle-independent) */
        const double ro   = r + dshell;
        const double epso = (r*eps + dshell) / ro;

        /* Forward-scattering amplitude used for normalisation */
        const double f0 = (4.0*pi/3.0)*(ro*ro*ro*epso
                          - r*r*r*eps*(1.0-contrast_shell));

        double sumx = 0.0, sum1x = 0.0;

        for (int ii = 0; ii < NPOI; ii++) {
            double xx = ((double)ii + 0.5)*step;
            xx = sin(xx);  /* integration variable: cos(theta) */

            /* Effective radii along q for this orientation */
            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 ra  = r *sc;
            const double rao = ro*sco;

            /* Outer-shell form factor amplitude */
            const double ffs =  (4.0*pi/3.0)*ro*ro*ro*epso
                              * sas_3j1x_x(q*rao)
                              * exp(-0.5*(q*sigma_outer)*(q*sigma_outer));

            /* Core form factor amplitude (negative: hollow interior) */
            const double ffc = -(4.0*pi/3.0)*r*r*r*eps*(1.0-contrast_shell)
                              * sas_3j1x_x(q*ra)
                              * exp(-0.5*(q*sigma_core)*(q*sigma_core));

            /* Displacement phase factor (perpendicular to symmetry axis) */
            const double phase = cos(q*displ*sqrt(1.0-xx*xx));

            sumx  += (ffs*ffs + ffc*ffc + 2.0*ffs*ffc*phase)*xx;
            sum1x += (ffs + ffc*phase)*xx;
        }

        sca   += sumx *step*f_exp;
        sum1  += sum1x*step*f_exp;   /* NOTE: Fortran used SUM1*STEP (bug); corrected here */
        sum2  += f0*f0*f_exp;
        sum1n += f0*f_exp;
    }

    /* Normalised form factor and mean amplitude squared */
    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;

    /* Hard-sphere structure factor */
    const double sq = hard_sphere_sf(q, volfraction_hs, radius_hs);

    /* Debye (polymer) scattering contribution */
    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;
    }

    /* Empirical power-law correction B(q) = 1 + A10*(q0/q)^m */
    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;  /* outer equatorial radius */
    default: return radius;                    /* core equatorial radius  */
    }
}
Back to Model Download