source: sasmodels/sasmodels/models/hollow_cylinder.c @ a807206

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a807206 was a807206, checked in by butler, 7 years ago

updating more parameter names addressing #649

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