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

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since e5cb3df was 99658f6, checked in by grethevj, 5 years ago

updated ER functions including cylinder excluded volume, to match 4.x

  • Property mode set to 100644
File size: 5.0 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
8radius_from_excluded_volume(double radius_minor, double r_ratio, double length)
9{
10    const double r_equiv = sqrt(radius_minor*radius_minor*r_ratio);
11    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length)));
12}
13
14static double
15radius_from_volume(double radius_minor, double r_ratio, double length)
16{
17    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length);
18    return cbrt(volume_ellcyl/M_4PI_3);
19}
20
21static double
22radius_from_min_dimension(double radius_minor, double r_ratio, double hlength)
23{
24    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
25    return (rad_min < hlength ? rad_min : hlength);
26}
27
28static double
29radius_from_max_dimension(double radius_minor, double r_ratio, double hlength)
30{
31    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
32    return (rad_max > hlength ? rad_max : hlength);
33}
34
35static double
36radius_from_diagonal(double radius_minor, double r_ratio, double length)
37{
38    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor);
39    return sqrt(radius_max*radius_max + 0.25*length*length);
40}
41
42static double
43effective_radius(int mode, double radius_minor, double r_ratio, double length)
44{
45    switch (mode) {
46    default:
47    case 1: // equivalent cylinder excluded volume
48        return radius_from_excluded_volume(radius_minor, r_ratio, length);
49    case 2: // equivalent volume sphere
50        return radius_from_volume(radius_minor, r_ratio, length);
51    case 3: // average radius
52        return 0.5*radius_minor*(1.0 + r_ratio);
53    case 4: // min radius
54        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor);
55    case 5: // max radius
56        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor);
57    case 6: // equivalent circular cross-section
58        return sqrt(radius_minor*radius_minor*r_ratio);
59    case 7: // half length
60        return 0.5*length;
61    case 8: // half min dimension
62        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length);
63    case 9: // half max dimension
64        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length);
65    case 10: // half diagonal
66        return radius_from_diagonal(radius_minor,r_ratio,length);
67    }
68}
69
70static void
71Fq(double q, double *F1, double *F2, double radius_minor, double r_ratio, double length,
72   double sld, double solvent_sld)
73{
74    // orientational average limits
75    const double va = 0.0;
76    const double vb = 1.0;
77    // inner integral limits
78    const double vaj=0.0;
79    const double vbj=M_PI;
80
81    const double radius_major = r_ratio * radius_minor;
82    const double rA = 0.5*(square(radius_major) + square(radius_minor));
83    const double rB = 0.5*(square(radius_major) - square(radius_minor));
84
85    //initialize integral
86    double outer_sum_F1 = 0.0;
87    double outer_sum_F2 = 0.0;
88    for(int i=0;i<GAUSS_N;i++) {
89        //setup inner integral over the ellipsoidal cross-section
90        const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
91        const double sin_val = sqrt(1.0 - cos_val*cos_val);
92        //const double arg = radius_minor*sin_val;
93        double inner_sum_F1 = 0.0;
94        double inner_sum_F2 = 0.0;
95        for(int j=0;j<GAUSS_N;j++) {
96            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
97            const double r = sin_val*sqrt(rA - rB*cos(theta));
98            const double be = sas_2J1x_x(q*r);
99            inner_sum_F1 += GAUSS_W[j] * be;
100            inner_sum_F2 += GAUSS_W[j] * be * be;
101        }
102        //now calculate the value of the inner integral
103        inner_sum_F1 *= 0.5*(vbj-vaj);
104        inner_sum_F2 *= 0.5*(vbj-vaj);
105
106        //now calculate outer integral
107        const double si = sas_sinx_x(q*0.5*length*cos_val);
108        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * si;
109        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * si * si;
110    }
111    // correct limits and divide integral by pi
112    outer_sum_F1 *= 0.5*(vb-va)/M_PI;
113    outer_sum_F2 *= 0.5*(vb-va)/M_PI;
114
115    // scale by contrast and volume, and convert to to 1/cm units
116    const double volume = form_volume(radius_minor, r_ratio, length);
117    const double contrast = sld - solvent_sld;
118    const double s = contrast*volume;
119    *F1 = 1.0e-2*s*outer_sum_F1;
120    *F2 = 1.0e-4*s*s*outer_sum_F2;
121}
122
123
124static double
125Iqabc(double qa, double qb, double qc,
126     double radius_minor, double r_ratio, double length,
127     double sld, double solvent_sld)
128{
129    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
130    // Given:    radius_major = r_ratio * radius_minor
131    const double qr = radius_minor*sqrt(square(r_ratio*qb) + square(qa));
132    const double be = sas_2J1x_x(qr);
133    const double si = sas_sinx_x(qc*0.5*length);
134    const double fq = be * si;
135    const double contrast = sld - solvent_sld;
136    const double volume = form_volume(radius_minor, r_ratio, length);
137    return 1.0e-4 * square(contrast * volume * fq);
138}
Note: See TracBrowser for help on using the repository browser.