Changeset ab2aea8 in sasmodels for sasmodels/models/rectangular_prism.c
- Timestamp:
- Oct 15, 2016 3:49:53 PM (7 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- e30d645
- Parents:
- ed0827a
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/rectangular_prism.c
r3a48772 rab2aea8 15 15 double c2a_ratio) 16 16 { 17 double termA, termB, termC; 18 19 double b_side = length_a * b2a_ratio; 20 double c_side = length_a * c2a_ratio; 21 double volume = length_a * b_side * c_side; 22 double a_half = 0.5 * length_a; 23 double b_half = 0.5 * b_side; 24 double c_half = 0.5 * c_side; 17 const double length_b = length_a * b2a_ratio; 18 const double length_c = length_a * c2a_ratio; 19 const double a_half = 0.5 * length_a; 20 const double b_half = 0.5 * length_b; 21 const double c_half = 0.5 * length_c; 25 22 26 23 //Integration limits to use in Gaussian quadrature 27 double v1a = 0.0;28 double v1b = M_PI_2; //theta integration limits29 double v2a = 0.0;30 double v2b = M_PI_2; //phi integration limits24 const double v1a = 0.0; 25 const double v1b = M_PI_2; //theta integration limits 26 const double v2a = 0.0; 27 const double v2b = M_PI_2; //phi integration limits 31 28 32 //Order of integration 33 int nordi=76; 34 int nordj=76; 29 double outer_sum = 0.0; 30 for(int i=0; i<76; i++) { 31 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 32 double sin_theta, cos_theta; 33 SINCOS(theta, sin_theta, cos_theta); 35 34 36 double sumi = 0.0; 37 38 for(int i=0; i<nordi; i++) { 35 const double termC = sinc(q * c_half * cos_theta); 39 36 40 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 37 double inner_sum = 0.0; 38 for(int j=0; j<76; j++) { 39 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 40 double sin_phi, cos_phi; 41 SINCOS(phi, sin_phi, cos_phi); 41 42 42 double arg = q * c_half * cos(theta); 43 if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;} 44 45 double sumj = 0.0; 46 47 for(int j=0; j<nordj; j++) { 48 49 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 50 51 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 52 53 arg = q * a_half * sin(theta) * sin(phi); 54 if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;} 55 56 arg = q * b_half * sin(theta) * cos(phi); 57 if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;} 58 59 double AP = termA * termB * termC; 60 61 sumj += Gauss76Wt[j] * (AP*AP); 62 63 } 64 65 sumj = 0.5 * (v2b-v2a) * sumj; 66 sumi += Gauss76Wt[i] * sumj * sin(theta); 67 43 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 44 const double termA = sinc(q * a_half * sin_theta * sin_phi); 45 const double termB = sinc(q * b_half * sin_theta * cos_phi); 46 const double AP = termA * termB * termC; 47 inner_sum += Gauss76Wt[j] * AP * AP; 48 } 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 68 51 } 69 52 70 double answer = 0.5*(v1b-v1a)* sumi;53 double answer = 0.5*(v1b-v1a)*outer_sum; 71 54 72 55 // Normalize by Pi (Eqn. 16). … … 78 61 79 62 // Multiply by contrast^2 and volume^2 80 answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume; 63 const double volume = length_a * length_b * length_c; 64 answer *= square((sld-solvent_sld)*volume); 81 65 82 66 // Convert from [1e-12 A-1] to [cm-1] … … 84 68 85 69 return answer; 86 87 70 }
Note: See TracChangeset
for help on using the changeset viewer.