Changeset 71b751d in sasmodels for sasmodels/models/rectangular_prism.c
- Timestamp:
- Aug 14, 2018 10:09:34 AM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 86aa992
- Parents:
- 2f8cbb9
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/rectangular_prism.c
rd86f0fc r71b751d 68 68 } 69 69 70 static void 71 Fq(double q, 72 double *F1, 73 double *F2, 74 double sld, 75 double solvent_sld, 76 double length_a, 77 double b2a_ratio, 78 double c2a_ratio) 79 { 80 const double length_b = length_a * b2a_ratio; 81 const double length_c = length_a * c2a_ratio; 82 const double a_half = 0.5 * length_a; 83 const double b_half = 0.5 * length_b; 84 const double c_half = 0.5 * length_c; 85 86 //Integration limits to use in Gaussian quadrature 87 const double v1a = 0.0; 88 const double v1b = M_PI_2; //theta integration limits 89 const double v2a = 0.0; 90 const double v2b = M_PI_2; //phi integration limits 91 92 double outer_sum_F1 = 0.0; 93 double outer_sum_F2 = 0.0; 94 for(int i=0; i<GAUSS_N; i++) { 95 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 96 double sin_theta, cos_theta; 97 SINCOS(theta, sin_theta, cos_theta); 98 99 const double termC = sas_sinx_x(q * c_half * cos_theta); 100 101 double inner_sum_F1 = 0.0; 102 double inner_sum_F2 = 0.0; 103 for(int j=0; j<GAUSS_N; j++) { 104 double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 105 double sin_phi, cos_phi; 106 SINCOS(phi, sin_phi, cos_phi); 107 108 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 109 const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi); 110 const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 111 const double AP = termA * termB * termC; 112 inner_sum_F1 += GAUSS_W[j] * AP; 113 inner_sum_F2 += GAUSS_W[j] * AP * AP; 114 } 115 inner_sum_F1 = 0.5 * (v2b-v2a) * inner_sum_F1; 116 inner_sum_F2 = 0.5 * (v2b-v2a) * inner_sum_F2; 117 outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin_theta; 118 outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin_theta; 119 } 120 121 outer_sum_F1 *= 0.5*(v1b-v1a); 122 outer_sum_F2 *= 0.5*(v1b-v1a); 123 124 // Normalize by Pi (Eqn. 16). 125 // The term (ABC)^2 does not appear because it was introduced before on 126 // the definitions of termA, termB, termC. 127 // The factor 2 appears because the theta integral has been defined between 128 // 0 and pi/2, instead of 0 to pi. 129 outer_sum_F1 /= M_PI_2; 130 outer_sum_F2 /= M_PI_2; 131 132 // Multiply by contrast and volume 133 const double s = (sld-solvent_sld) * (length_a * length_b * length_c); 134 135 // Convert from [1e-12 A-1] to [cm-1] 136 *F1 = 1e-2 * s * outer_sum_F1; 137 *F2 = 1e-4 * s * s * outer_sum_F2; 138 } 139 70 140 71 141 static double … … 82 152 const double b_half = 0.5 * length_b; 83 153 const double c_half = 0.5 * length_c; 84 const double volume = length_a * length_b * length_c;85 154 86 155 // Amplitude AP from eqn. (13) 87 88 156 const double termA = sas_sinx_x(qa * a_half); 89 157 const double termB = sas_sinx_x(qb * b_half); 90 158 const double termC = sas_sinx_x(qc * c_half); 91 92 159 const double AP = termA * termB * termC; 93 160 94 // Multiply by contrast ^2. Factor corresponding to volume^2 cancels with previous normalization.95 const double delrho = sld - solvent_sld;161 // Multiply by contrast and volume 162 const double s = (sld-solvent_sld) * (length_a * length_b * length_c); 96 163 97 164 // Convert from [1e-12 A-1] to [cm-1] 98 return 1.0e-4 * square( volume * delrho* AP);165 return 1.0e-4 * square(s * AP); 99 166 }
Note: See TracChangeset
for help on using the changeset viewer.