static double shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) { const double length_b = length_a * b2a_ratio; const double length_c = length_a * c2a_ratio; const double form_volume = length_a * length_b * length_c; const double a_core = length_a - 2.0*thickness; const double b_core = length_b - 2.0*thickness; const double c_core = length_c - 2.0*thickness; const double core_volume = a_core * b_core * c_core; return form_volume - core_volume; } static double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) { const double length_b = length_a * b2a_ratio; const double length_c = length_a * c2a_ratio; const double form_volume = length_a * length_b * length_c; return form_volume; } static double radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio) { const double r_equiv = sqrt(length_a*length_a*b2a_ratio/M_PI); const double length_c = length_a*c2a_ratio; return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c))); } static double radius_effective(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) // NOTE length_a is external dimension! { switch (mode) { default: case 1: // equivalent cylinder excluded volume return radius_from_excluded_volume(length_a, b2a_ratio, c2a_ratio); case 2: // equivalent outer volume sphere return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); case 3: // half length_a return 0.5 * length_a; case 4: // half length_b return 0.5 * length_a*b2a_ratio; case 5: // half length_c return 0.5 * length_a*c2a_ratio; case 6: // equivalent outer circular cross-section return length_a*sqrt(b2a_ratio/M_PI); case 7: // half ab diagonal return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); case 8: // half diagonal return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); } } static void Fq(double q, double *F1, double *F2, double sld, double solvent_sld, double length_a, double b2a_ratio, double c2a_ratio, double thickness) { const double length_b = length_a * b2a_ratio; const double length_c = length_a * c2a_ratio; const double a_half = 0.5 * length_a; const double b_half = 0.5 * length_b; const double c_half = 0.5 * length_c; const double vol_total = length_a * length_b * length_c; const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); //Integration limits to use in Gaussian quadrature const double v1a = 0.0; const double v1b = M_PI_2; //theta integration limits const double v2a = 0.0; const double v2b = M_PI_2; //phi integration limits double outer_sum_F1 = 0.0; double outer_sum_F2 = 0.0; for(int i=0; i