Changeset ab2aea8 in sasmodels for sasmodels/models/hollow_rectangular_prism_thin_walls.c
- Timestamp:
- Oct 15, 2016 5: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/hollow_rectangular_prism_thin_walls.c
r3a48772 rab2aea8 5 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 6 6 { 7 double b_side= length_a * b2a_ratio;8 double c_side= length_a * c2a_ratio;9 double vol_shell = 2.0 * (length_a* b_side + length_a*c_side + b_side*c_side);7 double length_b = length_a * b2a_ratio; 8 double length_c = length_a * c2a_ratio; 9 double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 10 10 return vol_shell; 11 11 } … … 18 18 double c2a_ratio) 19 19 { 20 double b_side= length_a * b2a_ratio;21 double c_side= length_a * c2a_ratio;22 double a_half = 0.5 * length_a;23 double b_half = 0.5 * b_side;24 double c_half = 0.5 * c_side;20 const double length_b = length_a * b2a_ratio; 21 const double length_c = length_a * c2a_ratio; 22 const double a_half = 0.5 * length_a; 23 const double b_half = 0.5 * length_b; 24 const double c_half = 0.5 * length_c; 25 25 26 26 //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 limits27 const double v1a = 0.0; 28 const double v1b = M_PI_2; //theta integration limits 29 const double v2a = 0.0; 30 const double v2b = M_PI_2; //phi integration limits 31 31 32 //Order of integration33 int nordi=76;34 int nordj=76;32 double outer_sum = 0.0; 33 for(int i=0; i<76; i++) { 34 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 35 35 36 double sumi = 0.0; 37 38 for(int i=0; i<nordi; i++) { 36 double sin_theta, cos_theta; 37 double sin_c, cos_c; 38 SINCOS(theta, sin_theta, cos_theta); 39 SINCOS(q*c_half*cos_theta, sin_c, cos_c); 39 40 40 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );41 42 41 // To check potential problems if denominator goes to zero here !!! 43 double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));44 double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));42 const double termAL_theta = 8.0 * cos_c / (q*q*sin_theta*sin_theta); 43 const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 45 44 46 double sumj= 0.0;47 48 for(int j=0; j<nordj; j++) { 45 double inner_sum = 0.0; 46 for(int j=0; j<76; j++) { 47 const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 49 48 50 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 51 49 double sin_phi, cos_phi; 50 double sin_a, cos_a; 51 double sin_b, cos_b; 52 SINCOS(phi, sin_phi, cos_phi); 53 SINCOS(q*a_half*sin_theta*sin_phi, sin_a, cos_a); 54 SINCOS(q*b_half*sin_theta*cos_phi, sin_b, cos_b); 55 52 56 // Amplitude AL from eqn. (7c) 53 double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) *54 sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));57 const double AL = termAL_theta 58 * sin_a*sin_b / (sin_phi*cos_phi); 55 59 56 60 // Amplitude AT from eqn. (9) 57 double AT = termAT_theta * ( (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi))58 + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)));61 const double AT = termAT_theta 62 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 59 63 60 sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 64 inner_sum += Gauss76Wt[j] * square(AL+AT); 65 } 61 66 62 } 63 64 sumj = 0.5 * (v2b-v2a) * sumj; 65 sumi += Gauss76Wt[i] * sumj * sin(theta); 66 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 67 69 } 68 70 69 double answer = 0.5*(v1b-v1a)*sumi;71 outer_sum *= 0.5*(v1b-v1a); 70 72 71 73 // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 72 74 // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 73 answer /=M_PI_2;75 double answer = outer_sum/M_PI_2; 74 76 75 77 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 76 answer *= (sld-solvent_sld)*(sld-solvent_sld);78 answer *= square(sld-solvent_sld); 77 79 78 80 // Convert from [1e-12 A-1] to [cm-1] … … 80 82 81 83 return answer; 82 83 84 }
Note: See TracChangeset
for help on using the changeset viewer.