Displaced core-shell ellipsoid model - ABS_core_shell_displ_core.c

    /*
 * Displaced core-shell ellipsoid model — ABSOLUTE SCALE version
 *
 * Same geometry and physics as core_shell_displ_core.c, but the scattering
 * contrast is specified via explicit SLD values (sld_core, sld_shell,
 * sld_solvent) rather than a dimensionless contrast ratio, and the particle
 * volume fraction is an explicit model parameter.  The kernel returns the
 * scattering intensity directly in cm^-1 (absolute scale).
 *
 * The form factor amplitudes are decomposed as:
 *
 *   F(q) = delta rho_shell  V_outer  phi(q·r_outer)
 *         + (delta rho_core − delta rho_shell)  V_core  phi(q r_core)  exp(iq delta)
 *
 * where delta rho = SLD − sld_solvent (in 10^-6 Ang^-2).
 *
 * The absolute intensity is:
 *
 *   I(q) = 1e-4  (phi/<V>)  [<F^2(q)> + <F(q)>^2  (S(q)−1)]  B(q)
 *         + I_poly(q)
 *
 * The factor 1e-4 converts (10^-6 Ang^-2)^2  Ang^6  Ang^-3 = 10^-12 Ang^-1
 * to cm^-1  (1 Ang^-1 = 10^8 cm^-1, times 10^-12 gives 10^-4).
 *
 * 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,
   double sld_core,
   double sld_shell,
   double sld_solvent,
   double thickness_shell,
   double volfraction_hs,
   double radius_hs,
   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;

    /* SLD contrasts in 10^-6 Ang^-2 */
    const double drho_shell    = sld_shell - sld_solvent;
    const double drho_core_rel = sld_core  - sld_shell;   /* extra contrast of core vs shell */

    /* 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 Debye-Waller surface roughness */
    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;   /* Sigma w  <F^2>  (size+orientation averaged F^2) */
    double sum1  = 0.0;   /* Sigma w · <F>   (size+orientation averaged F)  */
    double sum_w = 0.0;   /* Sigma w         (Gaussian weight sum)           */
    double sum_V = 0.0;   /* Sigma w · V     (weight-averaged particle volume)*/

    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));

        const double ro          = r + dshell;
        const double epso        = (r*eps + dshell) / ro;
        const double vol_particle= (4.0*pi/3.0)*ro*ro*ro*epso;

        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 ra  = r *sc;
            const double rao = ro*sco;

            const double vol_outer = (4.0*pi/3.0)*ro*ro*ro*epso;
            const double vol_core  = (4.0*pi/3.0)*r *r *r *eps;

            /* Shell contribution: delta rho_shell fills the entire outer ellipsoid */
            const double ffs = drho_shell * vol_outer
                             * sas_3j1x_x(q*rao)
                             * exp(-0.5*(q*sigma_outer)*(q*sigma_outer));

            /* Core contribution: (delta rho_core − delta rho_shell) fills the displaced core */
            const double ffc = drho_core_rel * vol_core
                             * sas_3j1x_x(q*ra)
                             * exp(-0.5*(q*sigma_core)*(q*sigma_core));

            /* Phase factor from core-centre displacement (perpendicular to axis) */
            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;
        sum_w += f_exp;
        sum_V += vol_particle*f_exp;
    }

    /* Weight-averaged quantities */
    const double F2_avg = (sum_w > 0.0) ? sca  / sum_w : 0.0;
    const double F_avg  = (sum_w > 0.0) ? sum1 / sum_w : 0.0;
    const double V_avg  = (sum_w > 0.0) ? sum_V/ sum_w : 1.0;

    /* Number density of particles */
    const double n_part = volfraction / V_avg;

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

    /* Absolute intensity in cm^-1
     * Units: 10^-4 * [Ang^-3] * [(10^-6 Ang^-2)^2 * Ang^6] = cm^-1 */
    const double i_particle = 1.0e-4 * n_part
                            * (F2_avg + F_avg*F_avg*(sq - 1.0));

    /* 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;

    /* Debye polymer contribution (already in cm^-1) */
    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;
    }

    return i_particle*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