#define INVALID(v) (v.radius_cap < v.radius) // Integral over a convex lens kernel for t in [h/R,1]. See the docs for // the definition of the function being integrated. // q is the magnitude of the q vector. // h is the length of the lens "inside" the cylinder. This negative wrt the // definition of h in the docs. // radius_cap is the radius of the lens // length is the cylinder length, or the separation between the lens halves // theta is the angle of the cylinder wrt q. static double _cap_kernel(double qab, double qc, double h, double radius_cap, double half_length) { // translate a point in [-1,1] to a point in [lower,upper] const double upper = 1.0; const double lower = h/radius_cap; // integral lower bound const double zm = 0.5*(upper-lower); const double zb = 0.5*(upper+lower); // cos term in integral is: // cos (q (R t - h + L/2) cos(theta)) // so turn it into: // cos (m t + b) // where: // m = q R cos(theta) // b = q(L/2-h) cos(theta) const double m = radius_cap*qc; // cos argument slope const double b = (half_length-h)*qc; // cos argument intercept const double qab_r = radius_cap*qab; // Q*R*sin(theta) double total = 0.0; for (int i=0; i h_c^2 // (2R+h) => 2R+ h_c-h_c + h => 2R + (R-h)-h_c + h => 3R-h_c // And so: // 2/3 pi (R-h)^2 (2R + h) => 2/3 pi h_c^2 (3 r_cap - h_c) // which is 2 V_cap, using the second form above. // // In this function we are going to use the first form of V_cap // // V = V_cyl + 2 V_cap // = pi r^2 L + pi hc (r^2 + hc^2/3) // = pi (r^2 (L+hc) + hc^3/3) const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); return M_PI*(radius*radius*(length+hc) + hc*hc*hc/3.0); } static double radius_from_excluded_volume(double radius, double radius_cap, double length) { const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); const double length_tot = length + 2.0*hc; return 0.5*cbrt(0.75*radius*(2.0*radius*length_tot + (radius + length_tot)*(M_PI*radius + length_tot))); } static double radius_from_volume(double radius, double radius_cap, double length) { const double vol_cappedcyl = form_volume(radius,radius_cap,length); return cbrt(vol_cappedcyl/M_4PI_3); } static double radius_from_totallength(double radius, double radius_cap, double length) { const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); return 0.5*length + hc; } static double radius_effective(int mode, double radius, double radius_cap, double length) { switch (mode) { default: case 1: // equivalent cylinder excluded volume return radius_from_excluded_volume(radius, radius_cap, length); case 2: // equivalent volume sphere return radius_from_volume(radius, radius_cap, length); case 3: // radius return radius; case 4: // half length return 0.5*length; case 5: // half total length return radius_from_totallength(radius, radius_cap,length); } } static void Fq(double q,double *F1, double *F2, double sld, double solvent_sld, double radius, double radius_cap, double length) { const double h = sqrt(radius_cap*radius_cap - radius*radius); const double half_length = 0.5*length; // translate a point in [-1,1] to a point in [0, pi/2] const double zm = M_PI_4; const double zb = M_PI_4; double total_F1 = 0.0; double total_F2 = 0.0; for (int i=0; i