source: sasmodels/sasmodels/models/elliptical_cylinder.c @ ef07e95

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ef07e95 was 82592da, checked in by richardh, 6 years ago

swapped 2d qa with qb in three elliptical particles, edits to orientation doc

  • Property mode set to 100644
File size: 2.5 KB
Line 
1static double
2form_volume(double radius_minor, double r_ratio, double length)
3{
4    return M_PI * radius_minor * radius_minor * r_ratio * length;
5}
6
7static double
8Iq(double q, double radius_minor, double r_ratio, double length,
9   double sld, double solvent_sld)
10{
11    // orientational average limits
12    const double va = 0.0;
13    const double vb = 1.0;
14    // inner integral limits
15    const double vaj=0.0;
16    const double vbj=M_PI;
17
18    const double radius_major = r_ratio * radius_minor;
19    const double rA = 0.5*(square(radius_major) + square(radius_minor));
20    const double rB = 0.5*(square(radius_major) - square(radius_minor));
21
22    //initialize integral
23    double outer_sum = 0.0;
24    for(int i=0;i<76;i++) {
25        //setup inner integral over the ellipsoidal cross-section
26        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
27        const double sin_val = sqrt(1.0 - cos_val*cos_val);
28        //const double arg = radius_minor*sin_val;
29        double inner_sum=0;
30        for(int j=0;j<76;j++) {
31            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017
32            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
33            const double r = sin_val*sqrt(rA - rB*cos(theta));
34            const double be = sas_2J1x_x(q*r);
35            inner_sum += Gauss76Wt[j] * be * be;
36        }
37        //now calculate the value of the inner integral
38        inner_sum *= 0.5*(vbj-vaj);
39
40        //now calculate outer integral
41        const double si = sas_sinx_x(q*0.5*length*cos_val);
42        outer_sum += Gauss76Wt[i] * inner_sum * si * si;
43    }
44    outer_sum *= 0.5*(vb-va);
45
46    //divide integral by Pi
47    const double form = outer_sum/M_PI;
48
49    // scale by contrast and volume, and convert to to 1/cm units
50    const double vol = form_volume(radius_minor, r_ratio, length);
51    const double delrho = sld - solvent_sld;
52    return 1.0e-4*square(delrho*vol)*form;
53}
54
55
56static double
57Iqxy(double qa, double qb, double qc,
58     double radius_minor, double r_ratio, double length,
59     double sld, double solvent_sld)
60{
61    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
62    // Given:    radius_major = r_ratio * radius_minor
63    const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa));
64    const double be = sas_2J1x_x(qr);
65    const double si = sas_sinx_x(qc*0.5*length);
66    const double Aq = be * si;
67    const double delrho = sld - solvent_sld;
68    const double vol = form_volume(radius_minor, r_ratio, length);
69    return 1.0e-4 * square(delrho * vol * Aq);
70}
Note: See TracBrowser for help on using the repository browser.