static double stacked_disks_kernel( double qab, double qc, double halfheight, double thick_layer, double radius, int n_stacking, double sigma_dnn, double core_sld, double layer_sld, double solvent_sld, double d) { // q is the q-value for the calculation (1/A) // radius is the core radius of the cylinder (A) // *_sld are the respective SLD's // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A) // zi is the dummy variable for the integration (x in Feigin's notation) const double besarg1 = radius*qab; //const double besarg2 = radius*qab; const double sinarg1 = halfheight*qc; const double sinarg2 = (halfheight+thick_layer)*qc; const double be1 = sas_2J1x_x(besarg1); //const double be2 = sas_2J1x_x(besarg2); const double be2 = be1; const double si1 = sas_sinx_x(sinarg1); const double si2 = sas_sinx_x(sinarg2); const double dr1 = core_sld - solvent_sld; const double dr2 = layer_sld - solvent_sld; const double area = M_PI*radius*radius; const double totald = 2.0*(thick_layer + halfheight); const double t1 = area * (2.0*halfheight) * dr1 * si1 * be1; const double t2 = area * dr2 * (totald*si2 - 2.0*halfheight*si1) * be2; double pq = square(t1 + t2); // loop for the structure factor S(q) double qd_cos_alpha = d*qc; //d*cos_alpha is the projection of d onto q (in other words the component //of d that is parallel to q. double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn); double sq=0.0; for (int kk=1; kk