source: sasmodels/sasmodels/models/hollow_cylinder.c @ 58210db

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 58210db was 58210db, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

hollow cylinder performance improvements

  • Property mode set to 100644
File size: 3.4 KB
Line 
1double form_volume(double radius, double core_radius, double length);
2
3double Iq(double q, double radius, double core_radius, double length, double sld,
4        double solvent_sld);
5double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld,
6        double solvent_sld, double theta, double phi);
7
8#define INVALID(v) (v.core_radius >= v.radius || v.radius >= v.length)
9
10// From Igor library
11static double hollow_cylinder_scaling(
12    double integrand, double delrho, double volume)
13{
14    double answer;
15    // Multiply by contrast^2
16    answer = integrand*delrho*delrho;
17
18    //normalize by cylinder volume
19    answer *= volume*volume;
20
21    //convert to [cm-1]
22    answer *= 1.0e-4;
23
24    return answer;
25}
26
27
28static double _hollow_cylinder_kernel(
29    double q, double core_radius, double radius, double length, double dum)
30{
31    //Note: lim_{r -> r_c} psi = J0(core_radius*qs)
32    const double qs = q*sqrt(1.0-dum*dum);
33    const double lam1 = sas_J1c(radius*qs);
34    const double lam2 = sas_J1c(core_radius*qs);
35    const double gamma_sq = square(core_radius/radius);
36    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00
37    const double t2 = sinc(q*length*dum/2.0);
38
39    return square(psi*t2);
40}
41
42
43static double hollow_cylinder_analytical_2D_scaled(
44    double q, double q_x, double q_y, double radius, double core_radius,
45    double length, double sld, double solvent_sld, double theta, double phi)
46{
47    double cyl_x, cyl_y; //, cyl_z
48    //double q_z;
49    double vol, cos_val, delrho;
50    double answer;
51    //convert angle degree to radian
52    theta = theta * M_PI_180;
53    phi = phi * M_PI_180;
54    delrho = solvent_sld - sld;
55
56    // Cylinder orientation
57    cyl_x = cos(theta) * cos(phi);
58    cyl_y = sin(theta);
59    //cyl_z = -cos(theta) * sin(phi);
60
61    // q vector
62    //q_z = 0;
63
64    // Compute the angle btw vector q and the
65    // axis of the cylinder
66    cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
67
68    answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val);
69
70    vol = form_volume(radius, core_radius, length);
71    answer = hollow_cylinder_scaling(answer, delrho, vol);
72
73    return answer;
74}
75
76
77double form_volume(double radius, double core_radius, double length)
78{
79    double v_shell = M_PI*length*(radius*radius-core_radius*core_radius);
80    return(v_shell);
81}
82
83
84double Iq(double q, double radius, double core_radius, double length,
85    double sld, double solvent_sld)
86{
87    int i;
88    double lower,upper,zi, inter;               //upper and lower integration limits
89    double summ,answer,delrho;                  //running tally of integration
90    double norm,volume; //final calculation variables
91
92    lower = 0.0;
93    upper = 1.0;                //limits of numerical integral
94
95    summ = 0.0;                 //initialize intergral
96    for (i=0;i<76;i++) {
97        zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0;
98        inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi);
99        summ += inter;
100    }
101
102    norm = summ*(upper-lower)/2.0;
103    volume = form_volume(radius, core_radius, length);
104    delrho = solvent_sld - sld;
105    answer = hollow_cylinder_scaling(norm, delrho, volume);
106
107    return(answer);
108}
109
110
111double Iqxy(double qx, double qy, double radius, double core_radius,
112    double length, double sld, double solvent_sld, double theta, double phi)
113{
114    const double q = sqrt(qx*qx+qy*qy);
115    return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi);
116}
Note: See TracBrowser for help on using the repository browser.