source: sasmodels/sasmodels/models/hollow_rectangular_prism_infinitely_thin_walls.c @ deb7ee0

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since deb7ee0 was deb7ee0, checked in by gonzalezm, 8 years ago

Added four parallelepiped-like models

  • Property mode set to 100644
File size: 3.0 KB
Line 
1double form_volume(double a_side, double b2a_ratio, double c2a_ratio);
2double Iq(double q, double sld, double solvent_sld, double a_side, 
3          double b2a_ratio, double c2a_ratio);
4double Iqxy(double qx, double qy, double sld, double solvent_sld, 
5            double a_side, double b2a_ratio, double c2a_ratio);
6
7double form_volume(double a_side, double b2a_ratio, double c2a_ratio)
8{
9    double b_side = a_side * b2a_ratio;
10    double c_side = a_side * c2a_ratio;
11    double vol_shell = 2.0 * (a_side*b_side + a_side*c_side + b_side*c_side);
12    return vol_shell;
13}
14
15double Iq(double q,
16    double sld,
17    double solvent_sld,
18    double a_side,
19    double b2a_ratio,
20    double c2a_ratio)
21{
22    double b_side = a_side * b2a_ratio;
23    double c_side = a_side * c2a_ratio;
24    double a_half = 0.5 * a_side;
25    double b_half = 0.5 * b_side;
26    double c_half = 0.5 * c_side;
27
28   //Integration limits to use in Gaussian quadrature
29    double v1a = 0.0;
30    double v1b = 0.5 * M_PI;  //theta integration limits
31    double v2a = 0.0;
32    double v2b = 0.5 * M_PI;  //phi integration limits
33   
34    //Order of integration
35    int nordi=76;                               
36    int nordj=76;
37
38    double sumi = 0.0;
39   
40    for(int i=0; i<nordi; i++) {
41
42            double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
43       
44        // To check potential problems if denominator goes to zero here !!!
45        double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));
46        double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));
47
48            double sumj = 0.0;
49       
50            for(int j=0; j<nordj; j++) {
51
52            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
53           
54            // Amplitude AL from eqn. (7c)
55            double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) * 
56                sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
57
58            // Amplitude AT from eqn. (9)
59            double AT = termAT_theta * (  (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi)) 
60                + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)) );
61
62            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
63
64            }
65
66            sumj = 0.5 * (v2b-v2a) * sumj;
67            sumi += Gauss76Wt[i] * sumj * sin(theta);
68
69    }
70
71    double answer = 0.5*(v1b-v1a)*sumi;
72
73    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization)
74    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
75    answer *= (2.0/M_PI);
76
77    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization.
78    answer *= (sld-solvent_sld)*(sld-solvent_sld);
79
80    // Convert from [1e-12 A-1] to [cm-1]
81    answer *= 1.0e-4;
82
83    return answer;
84   
85}
86
87double Iqxy(double qx, double qy,
88    double sld,
89    double solvent_sld,
90    double a_side,
91    double b2a_ratio,
92    double c2a_ratio)
93{
94    double q = sqrt(qx*qx + qy*qy);
95    double intensity = Iq(q, sld, solvent_sld, a_side, b2a_ratio, c2a_ratio); 
96    return intensity;   
97}
Note: See TracBrowser for help on using the repository browser.