source: sasmodels/sasmodels/models/parallelepiped.c @ d507c3a

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

Added parallelepiped model

  • Property mode set to 100644
File size: 4.6 KB
Line 
1double form_volume(double a_side, double b_side, double c_side);
2double Iq(double q, double sld, double solvent_sld, double a_side, double b_side, double c_side);
3double Iqxy(double qx, double qy, double sld, double solvent_sld,
4    double a_side, double b_side, double c_side, double theta, double phi, double psi);
5
6// From Igor library
7double _pkernel(double a, double b,double c, double ala, double alb, double alc);
8double _pkernel(double a, double b,double c, double ala, double alb, double alc){
9    double argA,argB,argC,tmp1,tmp2,tmp3;
10    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
11    argA = a*ala/2.0;
12    argB = b*alb/2.0;
13    argC = c*alc/2.0;
14    if(argA==0.0) {
15        tmp1 = 1.0;
16    } else {
17        tmp1 = sin(argA)*sin(argA)/argA/argA;
18    }
19    if (argB==0.0) {
20        tmp2 = 1.0;
21    } else {
22        tmp2 = sin(argB)*sin(argB)/argB/argB;
23    }
24    if (argC==0.0) {
25        tmp3 = 1.0;
26    } else {
27        tmp3 = sin(argC)*sin(argC)/argC/argC;
28    }
29    return (tmp1*tmp2*tmp3);
30
31}
32
33double form_volume(double a_side, double b_side, double c_side)
34{
35    return a_side * b_side * c_side;
36}
37
38double Iq(double q,
39    double sld,
40    double solvent_sld,
41    double a_side,
42    double b_side,
43    double c_side)
44{
45    double tmp1, tmp2, yyy;
46   
47    double mu = q * b_side;
48   
49    // Scale sides by B
50    double a_scaled = a_side / b_side;
51    double c_scaled = c_side / b_side;
52   
53    // outer integral (with 76 gauss points), integration limits = 0, 1
54    double summ = 0; //initialize integral
55
56    for( int i=0; i<76; i++) {
57               
58        // inner integral (with 76 gauss points), integration limits = 0, 1
59       
60        double summj = 0.0;
61        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );             
62               
63        for(int j=0; j<76; j++) {
64
65            double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
66            double mudum = mu * sqrt(1.0-sigma*sigma);
67            double arg1 = 0.5 * mudum * cos(0.5*M_PI*uu);
68            double arg2 = 0.5 * mudum * a_scaled * sin(0.5*M_PI*uu);
69            if(arg1==0.0) {
70                tmp1 = 1.0;
71            } else {
72                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
73            }
74            if (arg2==0.0) {
75                tmp2 = 1.0;
76            } else {
77                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
78            }
79            yyy = Gauss76Wt[j] * tmp1 * tmp2;
80
81            summj += yyy;
82        }
83               
84        // value of the inner integral
85        double answer = 0.5 * summj;
86
87        double arg = 0.5 * mu * c_scaled * sigma;
88        if ( arg == 0.0 ) {
89            answer *= 1.0;
90        } else {
91            answer *= sin(arg)*sin(arg)/arg/arg;
92        }
93               
94        // sum of outer integral
95        yyy = Gauss76Wt[i] * answer;
96        summ += yyy;
97       
98    }   
99   
100    const double vd = (sld-solvent_sld) * form_volume(a_side, b_side, c_side);
101    return 1.0e-4 * 0.5 * vd * vd * summ;
102   
103}
104
105
106double Iqxy(double qx, double qy,
107    double sld,
108    double solvent_sld,
109    double a_side,
110    double b_side,
111    double c_side,
112    double theta,
113    double phi,
114    double psi)
115{
116    double q = sqrt(qx*qx+qy*qy);
117    double qx_scaled = qx/q;
118    double qy_scaled = qy/q;
119
120    // Convert angles given in degrees to radians
121    theta *= M_PI_180;
122    phi   *= M_PI_180;
123    psi   *= M_PI_180;
124   
125    // Parallelepiped c axis orientation
126    double cparallel_x = cos(theta) * cos(phi);
127    double cparallel_y = sin(theta);
128   
129    // Compute angle between q and parallelepiped axis
130    double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz;
131
132    // Parallelepiped a axis orientation
133    double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
134    double parallel_y = sin(psi)*cos(theta);
135    double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled;
136
137    // Parallelepiped b axis orientation
138    double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
139    double bparallel_y = cos(theta)*cos(psi);
140    double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled;
141
142    // The following tests should always pass
143    if (fabs(cos_val_c)>1.0) {
144      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
145      cos_val_c = 1.0;
146    }
147    if (fabs(cos_val_a)>1.0) {
148      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
149      cos_val_a = 1.0;
150    }
151    if (fabs(cos_val_b)>1.0) {
152      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
153      cos_val_b = 1.0;
154    }
155   
156    // Call the IGOR library function to get the kernel
157    double form = _pkernel( q*a_side, q*b_side, q*c_side, cos_val_a, cos_val_b, cos_val_c);
158 
159    // Multiply by contrast^2
160    const double vd = (sld - solvent_sld) * form_volume(a_side, b_side, c_side);
161    return 1.0e-4 * vd * vd * form;
162}
Note: See TracBrowser for help on using the repository browser.