Ticket #817: micelles_ellipsoid.2.c

static double effective_radius(double radius_core, double epsilon, double alpha)
{
        return radius_core
                 * sqrt(square(sin(alpha)) + square(epsilon * cos(alpha)));
}

static double bes_core(double q, double or_radius)
{
        return sas_3j1x_x(q*or_radius);
}

static double bes_corona(double q, double or_radius, double radius_gyr, double d_penetration)
{
        return sas_sinx_x(q * (or_radius + d_penetration * radius_gyr));
}

static double Fs(double q, double radius_core, double epsilon)
{
        const int bins = 100;
        const double dalpha = M_PI_2 / bins;
        double total = 0.0;
        for(int i = 0; i != bins; ++i) {
                const double alpha = i * dalpha;
                const double or_radius=effective_radius(radius_core,epsilon,alpha);
                total += square(bes_core(q, or_radius)) * sin(alpha);
        }
        return total * dalpha;
}

static double Ssc(double q, double radius_core, double radius_gyr, double d_penetration, double epsilon, double chain_ampl)
{
        const int bins = 100;
        const double dalpha = M_PI_2 / bins;
        double total = 0.0;
        for(int i = 0; i != bins; ++i) {
                const double alpha = i * dalpha;
                const double or_radius=effective_radius(radius_core, epsilon, alpha);
                total += bes_core(q, or_radius) * bes_corona(q, or_radius, radius_gyr, d_penetration) * sin(alpha);
        }
        return total * dalpha * chain_ampl;
}

static double Scc(double q, double radius_core, double radius_gyr, double d_penetration, double epsilon, double chain_ampl)
{
        const int bins = 100;
        const double dalpha = M_PI_2 / bins;
        double total = 0.0;
        for(int i = 0; i != bins; ++i) {
                const double alpha = i * dalpha;
                const double or_radius=effective_radius(radius_core,epsilon,alpha);
                total += square(bes_corona(q, or_radius, radius_gyr, d_penetration)) * sin(alpha);
        }
        return total * dalpha * chain_ampl * chain_ampl;
}

static double Iq(
        double q,
        double ndensity,
        double v_core,
        double v_corona,
        double solvent_sld,
        double core_sld,
        double corona_sld,
        double radius_core,
        double radius_gyr,
        double d_penetration,
        double n_aggreg,
        double epsilon)
{
        const double beta_core = v_core * (core_sld - solvent_sld);
        const double beta_corona = v_corona * (corona_sld - solvent_sld);

        const double term1 = square(n_aggreg*beta_core)*Fs(q, radius_core, epsilon);

        // Self-correlation term of the chains
        const double qrg2 = square(q*radius_gyr);
        const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2);
        const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain;


        const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2;

        // Interference cross-term between core and chains
        const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona * Ssc(q, radius_core, radius_gyr, d_penetration, epsilon, chain_ampl);

        // Interference cross-term between chains
        const double term4 = n_aggreg * (n_aggreg - 1.0) * beta_corona * beta_corona * Scc(q, radius_core, radius_gyr, d_penetration, epsilon, chain_ampl);

        // I(q)_micelle : Sum of 4 terms computed above
        double i_micelle = term1 + term2 + term3 + term4;

        // rescale from [A^2] to [cm^2]
        i_micelle *= 1.0e-13;

        // "normalize" by number density --> intensity in [cm-1]
        i_micelle *= ndensity;

        return(i_micelle);
}