Three-layer displaced core-shell ellipsoid - absolute scale version - ABS_core_shell_extra_inner_shell.c

    /*
 * Three-layer displaced core-shell ellipsoid — ABSOLUTE SCALE version
 *
 * Same geometry as core_shell_extra_inner_shell.c but contrasts are
 * derived from explicit SLD values for core, inner shell, outer shell,
 * and solvent.  The particle volume fraction is an explicit parameter and
 * the kernel returns absolute intensity in cm^-1.
 *
 * Scattering length densities (SLDs) are in 10^-6 Ang^-2.
 * Contrasts:  delta rho_X = sld_X − sld_solvent
 *
 * Form factor amplitude at a given orientation:
 *
 *   F = FFS + FFC * phase
 *
 *   FFS = delta rho_out * V_RO * phi(q*r_RO) * DW_outer
 *
 *   FFC = [ (delta rho_in − delta rho_out) * V_R  * phi(q*r_R )
 *         + (delta rho_core− delta rho_in) * V_RI * phi(q*r_RI) ] * DW_core
 *
 * Absolute intensity:
 *
 *   I(q) = 1e-4 * (Phi/<V>) * [<F^2> + <F>^2*(S(q)−1)] * B(q) + I_poly(q)
 *
 * The 1e-4 factor converts (10^-6 Ang^-2)^2 * Ang^3 * Ang^-3 = 10^-12 Ang^-1
 * to cm^-1.  Keep SASView scale=1 for true absolute units.
 *
 * 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,
   double sld_core,
   double sld_inner_shell,
   double sld_outer_shell,
   double sld_solvent,
   double thickness_shell,
   double thickness_inner_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_out  = sld_outer_shell - sld_solvent;
    const double drho_in   = sld_inner_shell - sld_solvent;
    const double drho_core = sld_core        - sld_solvent;

    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;
    double sum1  = 0.0;
    double sum_w = 0.0;
    double sum_V = 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));

        const double ro   = r + dshell;
        const double epso = (r*eps + dshell) / ro;
        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;

        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 delta rho_out */
            const double ffs = drho_out * vol_ro
                             * sas_3j1x_x(q*rao)
                             * exp(-0.5*(q*sigma_outer)*(q*sigma_outer));

            /* Inner structure: onion-peel decomposition.
             * Each term adds the *extra* contrast of one layer over the layer outside it. */
            const double ffc = ((drho_in   - drho_out)  * vol_r  * sas_3j1x_x(q*ra)
                               +(drho_core  - drho_in)   * vol_ri * sas_3j1x_x(q*rai))
                             * exp(-0.5*(q*sigma_core)*(q*sigma_core));

            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_ro*f_exp;
    }

    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;

    const double n_part = volfraction / V_avg;
    const double sq     = hard_sphere_sf(q, volfraction_hs, radius_hs);

    /* 1e-4: (10^-6 Ang^-2)^2 * Ang^6 * Ang^-3 * 10^8 (Ang^-1 -> cm^-1) = 10^-4 */
    const double i_particle = 1.0e-4 * n_part
                            * (F2_avg + F_avg*F_avg*(sq - 1.0));

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

    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;
    default: return radius;
    }
}
Back to Model Download