- Categories
- Ellipsoid
- Magnetically oriented, rotating and precessing anisometric particle (MORP)
- morp_ellipsoid.c
Magnetically oriented, rotating and precessing anisometric particle (MORP) - morp_ellipsoid.c
static double
form_volume(double radius_polar, double radius_equatorial)
{
return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial;
}
static double
radius_from_volume(double radius_polar, double radius_equatorial)
{
return cbrt(radius_polar*radius_equatorial*radius_equatorial);
}
static double
radius_from_curvature(double radius_polar, double radius_equatorial)
{
// Trivial cases
if (radius_polar * radius_equatorial == 0.) return 0.;
// see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449
const double ratio = (radius_polar < radius_equatorial
? radius_polar / radius_equatorial
: radius_equatorial / radius_polar);
const double e1 = sqrt(1.0 - ratio*ratio);
const double b1 = 1.0 + asin(e1) / (e1 * ratio);
const double bL = (1.0 + e1) / (1.0 - e1);
const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL);
const double delta = 0.75 * b1 * b2;
const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial;
return 0.5 * cbrt(ddd);
}
static double
radius_effective(int mode, double radius_polar, double radius_equatorial)
{
switch (mode) {
default:
case 1: // average curvature
return radius_from_curvature(radius_polar, radius_equatorial);
case 2: // equivalent volume sphere
return radius_from_volume(radius_polar, radius_equatorial);
case 3: // min radius
return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial);
case 4: // max radius
return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial);
}
}
void SET_VECTOR(double *vector, double v0, double v1, double v2) {
vector[0] = v0;
vector[1] = v1;
vector[2] = v2;
}
void SCALE_VECTOR(double *vector, double a) {
vector[0] = a*vector[0];
vector[1] = a*vector[1];
vector[2] = a*vector[2];
}
void ADD_VECTOR(double *result_vec, double *vec1, double *vec2) {
result_vec[0] = vec1[0] + vec2[0];
result_vec[1] = vec1[1] + vec2[1];
result_vec[2] = vec1[2] + vec2[2];
}
static void rotated_long_axis(
double gamma_1, double gamma_2,
double psi,
double *result_rotated_axis) {
double cos_psi, sin_psi;
double cos_gamma_1, sin_gamma_1;
double cos_gamma_2, sin_gamma_2;
SINCOS(psi, sin_psi, cos_psi);
SINCOS(gamma_1, sin_gamma_1, cos_gamma_1);
SINCOS(gamma_2, sin_gamma_2, cos_gamma_2);
double vector_long[3];
double easy_cross_long[3];
//double vector_long_rotated[3];
// B field is along (1,0,0), orientation of easy axis
// is precessing in a cone around B
// vector_long: long axis direction,
// vector lies in plane perpendicular to easy axis
SET_VECTOR(vector_long, -sin_psi, cos_psi*sin_gamma_1, cos_psi*cos_gamma_1);
// to obtain all possible vectors of vector_long, integrate gamma_2
// use for this purpose Rodrigues rotation formula, need cross product:
SET_VECTOR(easy_cross_long, 0.0, cos_gamma_1, -sin_gamma_1);
// rotation formula:
// vector_long_rotated = vector_long*cos_gamma_2 + easy_cross_long*sin_gamma_2
SCALE_VECTOR(vector_long, cos_gamma_2);
SCALE_VECTOR(easy_cross_long, sin_gamma_2);
ADD_VECTOR(
result_rotated_axis,
vector_long,
easy_cross_long
);
// last part rotation formula not necessary as r_easy, r_long
// are perpendicular by definition
// see: r_easy*dot_product(r_easy, r_long)*(1d0-cos(gamma_2))
// rotate coordinate frame in case the magnetic field is not
// pointing along (1,0,0)
// SET_VEC(result_rotated_axis,
// cos_beta*vector_long_rotated[0] - sin_beta*vector_long_rotated[2],
// vector_long_rotated[1],
// sin_beta*vector_long_rotated[0] + cos_beta*vector_long_rotated[2]
// );
}
static double oriented_spindle_amplitude(
double qx, double qy, double radius_polar, double radius_equatorial,
double gamma_1, double gamma_2, double psi) {
// get orientation of long axis for given gamma_1, psi, gamma_2 angles:
double vec_long_rotated[3];
rotated_long_axis(gamma_1, gamma_2, psi, vec_long_rotated);
// what is angle between q vector and long axis?
// acos results in NaN if argument is numerically slightly larger than 1
const double q = sqrt(square(qx) + square(qy));
double cos_alpha = 1.0;
if (q > 0) {
cos_alpha = (qx*vec_long_rotated[0] + qy*vec_long_rotated[1])/q;
}
// calculate the effective radius
const double req2 = square(radius_equatorial);
const double r_eff = sqrt( req2 + (square(radius_polar) - req2)*square(cos_alpha) );
// effective radius for ellipsoidal model:
// http://gisaxs.com/index.php/Form_Factor:Ellipsoid_of_revolution
return sas_3j1x_x(q*r_eff);
}
static double boltzmann_statistics(
double xi, double psi) {
// xi * exp(xi * (cos(psi) - 1) ) / ( 1 - exp(-2xi) )
if (xi < 0.0001) {
// avoid dividing by zero
// use first order taylor approximation of exp(x) in this case
// which is accurate to 2e-8
return 0.5*(1.0 + xi * (cos(psi) - 1.0));
} else {
return xi * exp(xi*(cos(psi) - 1.0)) / (1.0 - exp(-2*xi));
}
}
static void
Fq(double q,
double *F1,
double *F2,
double sld,
double sld_solvent,
double radius_polar,
double radius_equatorial,
double xi)
{
// perform integration over gamma_1=0..2pi and gamma_2=0..2pi
// to account for the precession motion of the easy axis around B with
// cone angle psi
// as well as the motion of the long axis of the spindle
// translate a point in [-1,1] to a point in [0, 1]
// const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;
double total_F1 = 0.0;
double total_F2 = 0.0;
for (int i_detangle=0; i_detangle < GAUSS_N; i_detangle++) {
double detangle = M_PI_2 * (GAUSS_Z[i_detangle] + 1.0); // 0 .. pi
double const qx = q*cos(detangle);
double const qy = q*sin(detangle);
double integral_psi_F1 = 0.0;
double integral_psi_F2 = 0.0;
for (int i_psi=0; i_psi < GAUSS_N; i_psi++) {
const double psi = M_PI_2 * (GAUSS_Z[i_psi] + 1.0); // 0 .. pi
double integral_g1_F1 = 0.0;
double integral_g1_F2 = 0.0;
for (int i_g1=0; i_g1 < GAUSS_N; i_g1++) {
const double gamma_1 = M_PI * (GAUSS_Z[i_g1] + 1.0); // 0 .. 2 pi
double integral_g2_F1 = 0.0;
double integral_g2_F2 = 0.0;
for (int i_g2=0; i_g2 < GAUSS_N; i_g2++) {
const double gamma_2 = M_PI * (GAUSS_Z[i_g2] + 1.0); // 0 .. 2 pi
const double f = oriented_spindle_amplitude(
qx, qy, radius_polar, radius_equatorial,
gamma_1, gamma_2, psi);
integral_g2_F1 += GAUSS_W[i_g2] * f;
integral_g2_F2 += GAUSS_W[i_g2] * f * f;
}
integral_g1_F1 += GAUSS_W[i_g1] * integral_g2_F1;
integral_g1_F2 += GAUSS_W[i_g1] * integral_g2_F2;
}
const double sin_psi = sin(psi);
const double boltzmann_weight = boltzmann_statistics(xi, psi);
integral_psi_F1 += GAUSS_W[i_psi] * sin_psi * boltzmann_weight * integral_g1_F1;
integral_psi_F2 += GAUSS_W[i_psi] * sin_psi * boltzmann_weight * integral_g1_F2;
}
total_F1 += GAUSS_W[i_detangle] * integral_psi_F1;
total_F2 += GAUSS_W[i_detangle] * integral_psi_F2;
}
// factors for the translation of the [-1,1] range to
// 0 .. 2pi, 0 .. 2pi, 0 .. pi -> pi^3 / 2
// divide by 4 pi^2 for the integration over the whole space
// divide by pi for integration over detector
// -> factor of 1/8
total_F1 *= 0.125;
total_F2 *= 0.125;
const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
*F1 = 1e-2 * s * total_F1;
*F2 = 1e-4 * s * s * total_F2;
}
static double
Iqxy(double qx, double qy,
double sld,
double sld_solvent,
double radius_polar,
double radius_equatorial,
double xi)
{
// mu_bohr / k_Boltzmann = 0.6717140430498562
// const double xi = 0.6717140430498562 * magnetic_moment * magnetic_field / temperature;
double integral_psi_F2 = 0.0;
for (int i_psi=0; i_psi < GAUSS_N; i_psi++) {
const double psi = M_PI_2 * (GAUSS_Z[i_psi] + 1.0); // 0 .. pi
double integral_g1_F2 = 0.0;
for (int i_g1=0; i_g1 < GAUSS_N; i_g1++) {
const double gamma_1 = M_PI * (GAUSS_Z[i_g1] + 1.0); // 0 .. 2 pi
double integral_g2_F2 = 0.0;
for (int i_g2=0; i_g2 < GAUSS_N; i_g2++) {
const double gamma_2 = M_PI * (GAUSS_Z[i_g2] + 1.0); // 0 .. 2 pi
const double f = oriented_spindle_amplitude(
qx, qy, radius_polar, radius_equatorial,
gamma_1, gamma_2, psi);
integral_g2_F2 += GAUSS_W[i_g2] * f * f;
}
integral_g1_F2 += GAUSS_W[i_g1] * integral_g2_F2;
}
const double sin_psi = sin(psi);
const double boltzmann_weight = boltzmann_statistics(xi, psi);
integral_psi_F2 += GAUSS_W[i_psi] * sin_psi * boltzmann_weight * integral_g1_F2;
}
// factors for the translation of the [-1,1] range to
// 0 .. 2pi, 0 .. 2pi, 0 .. pi -> pi^3 / 2
// divide by 4 pi^2 for the integration over the whole space
// -> factor of pi/8
integral_psi_F2 *= 0.125 * M_PI;
const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial);
return 1.0e-4 * square(s) * integral_psi_F2;
}
Back to Model
Download