static double sc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) { // Equations from Matsuoka 9-10-11, multiplied by |q| const double a1 = qa; const double a2 = qb; const double a3 = qc; // Matsuoka 13-14-15 // Z_k numerator: 1 - exp(a)^2 // Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) // Rewriting numerator // => -(exp(2a) - 1) // => -expm1(2a) // Rewriting denominator // => exp(a)^2 - 2 cos(d ak) exp(a) + 1) // => (exp(a) - 2 cos(d ak)) * exp(a) + 1 const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); const double exp_arg = exp(arg); const double Zq = -cube(expm1(2.0*arg)) / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); return Zq; } // occupied volume fraction calculated from lattice symmetry and sphere radius static double sc_volume_fraction(double radius, double lattice_spacing) { return sphere_volume(radius/lattice_spacing); } static double form_volume(double radius) { return sphere_volume(radius); } static double Iq(double q, double lattice_spacing, double lattice_distortion, double radius, double sld, double solvent_sld) { // translate a point in [-1,1] to a point in [0, 2 pi] const double phi_m = M_PI_4; const double phi_b = M_PI_4; // translate a point in [-1,1] to a point in [0, pi] const double theta_m = M_PI_4; const double theta_b = M_PI_4; double outer_sum = 0.0; for(int i=0; i