Changeset fc0b7aa in sasmodels
- Timestamp:
- Nov 5, 2017 12:33:25 PM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 8de1477
- Parents:
- 510277d
- Location:
- sasmodels/models
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
TabularUnified sasmodels/models/core_shell_parallelepiped.c ¶
r3a1fc7d rfc0b7aa 1 // Set OVERLAPPING to 1 in order to fill in the edges of the box, with 2 // c endcaps and b overlapping a. With the proper choice of parameters, 3 // (setting rim slds to sld, core sld to solvent, rim thickness to thickness 4 // and subtracting 2*thickness from length, this should match the hollow 5 // rectangular prism.) Set it to 0 for the documented behaviour. 6 #define OVERLAPPING 0 1 7 static double 2 8 form_volume(double length_a, double length_b, double length_c, 3 9 double thick_rim_a, double thick_rim_b, double thick_rim_c) 4 10 { 5 return length_a * length_b * length_c + 6 2.0 * thick_rim_a * length_b * length_c + 7 2.0 * thick_rim_b * length_a * length_c + 8 2.0 * thick_rim_c * length_a * length_b; 11 return 12 #if OVERLAPPING 13 // Hollow rectangular prism only includes the volume of the shell 14 // so uncomment the next line when comparing. Solid rectangular 15 // prism, or parallelepiped want filled cores, so comment when 16 // comparing. 17 //-length_a * length_b * length_c + 18 (length_a + 2.0*thick_rim_a) * 19 (length_b + 2.0*thick_rim_b) * 20 (length_c + 2.0*thick_rim_c); 21 #else 22 length_a * length_b * length_c + 23 2.0 * thick_rim_a * length_b * length_c + 24 2.0 * length_a * thick_rim_b * length_c + 25 2.0 * length_a * length_b * thick_rim_c; 26 #endif 9 27 } 10 28 … … 38 56 39 57 double Vin = length_a * length_b * length_c; 40 double VtA = (2.0 * thick_rim_a * length_b * length_c); 41 double VtB = (2.0 * length_a * thick_rim_b * length_c); 42 double VtC = (2.0 * length_a * length_b * thick_rim_c); 58 #if OVERLAPPING 59 const double capA_area = length_b*length_c; 60 const double capB_area = (length_a+2.*thick_rim_a)*length_c; 61 const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 62 #else 63 const double capA_area = length_b*length_c; 64 const double capB_area = length_a*length_c; 65 const double capC_area = length_a*length_b; 66 #endif 67 const double Va = length_a * capA_area; 68 const double Vb = length_b * capB_area; 69 const double Vc = length_c * capC_area; 70 const double Vat = Va + 2.0 * thick_rim_a * capA_area; 71 const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 72 const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 43 73 44 74 // Scale factors (note that drC is not used later) … … 47 77 const double drB = (brim_sld-solvent_sld); 48 78 const double drC = (crim_sld-solvent_sld); 49 50 // Precompute scale factors for combining cross terms from the shape51 const double dr0_Vin = dr0*Vin;52 const double drA_VtA = drA*VtA;53 const double drB_VtB = drB*VtB;54 const double drC_VtC = drC*VtC;55 const double drV_delta = dr0_Vin - drA_VtA - drB_VtB - drC_VtC;56 57 /* *************** algorithm description ******************58 59 // Rewrite f as x*siC + y*siCt to move the siC/siCt calculation out60 // of the inner loop. That is:61 62 f = (di-da-db-dc) sa sb sc + da sa' sb sc + db sa sb' sc + dc sa sb sc'63 = [ (di-da-db-dc) sa sb + da sa' sb + db sa sb' ] sc + [dc sa sb] sc'64 = x sc + y sc'65 66 // where:67 di = delta rho_core V_core68 da = delta rho_rimA V_rimA69 db = delta rho_rimB V_rimB70 dc = delta rho_rimC V_rimC71 sa = j_0 (q_a a/2) // siA the code72 sb = j_0 (q_b b/2)73 sc = j_0 (q_c c/2)74 sa' = j_0(q_a a_rim/2) // siAt the code75 sb' = j_0(q_b b_rim/2)76 sc' = j_0(q_c c_rim/2)77 78 // qa, qb, and qc are generated using polar coordinates, with the79 // outer loop integrating over [0,1] after the u-substitution80 // sigma = cos(theta), sqrt(1-sigma^2) = sin(theta)81 // and inner loop integrating over [0, pi/2] as82 // uu = phi83 84 ************************************************************ */85 79 86 80 // outer integral (with gauss points), integration limits = 0, 1 … … 103 97 const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); 104 98 105 const double x = drV_delta*siA*siB + drA_VtA*siB*siAt + drB_VtB*siA*siBt; 106 const double form = x*siC + drC_VtC*siA*siB*siCt; 99 #if OVERLAPPING 100 const double f = dr0*Vin*siA*siB*siC 101 + drA*(Vat*siAt-Va*siA)*siB*siC 102 + drB*siAt*(Vbt*siBt-Vb*siB)*siC 103 + drC*siAt*siBt*(Vct*siCt-Vc*siC); 104 #else 105 const double f = dr0*Vin*siA*siB*siC 106 + drA*(Vat*siAt-Va*siA)*siB*siC 107 + drB*siA*(Vbt*siBt-Vb*siB)*siC 108 + drC*siA*siB*(Vct*siCt-Vc*siC); 109 #endif 107 110 108 inner_sum += Gauss76Wt[j] * f orm * form;111 inner_sum += Gauss76Wt[j] * f * f; 109 112 } 110 113 inner_sum *= 0.5; … … 139 142 140 143 double Vin = length_a * length_b * length_c; 141 double VtA = 2.0 * thick_rim_a * length_b * length_c; 142 double VtB = 2.0 * length_a * thick_rim_b * length_c; 143 double VtC = 2.0 * length_a * length_b * thick_rim_c; 144 #if OVERLAPPING 145 const double capA_area = length_b*length_c; 146 const double capB_area = (length_a+2.*thick_rim_a)*length_c; 147 const double capC_area = (length_a+2.*thick_rim_a)*(length_b+2.*thick_rim_b); 148 #else 149 const double capA_area = length_b*length_c; 150 const double capB_area = length_a*length_c; 151 const double capC_area = length_a*length_b; 152 #endif 153 const double Va = length_a * capA_area; 154 const double Vb = length_b * capB_area; 155 const double Vc = length_c * capC_area; 156 const double Vat = Va + 2.0 * thick_rim_a * capA_area; 157 const double Vbt = Vb + 2.0 * thick_rim_b * capB_area; 158 const double Vct = Vc + 2.0 * thick_rim_c * capC_area; 144 159 145 160 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 146 161 // the scaling by B. 147 double tA = length_a + 2.0*thick_rim_a; 148 double tB = length_b + 2.0*thick_rim_b; 149 double tC = length_c + 2.0*thick_rim_c; 150 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 151 double siA = sas_sinx_x(0.5*length_a*qa); 152 double siB = sas_sinx_x(0.5*length_b*qb); 153 double siC = sas_sinx_x(0.5*length_c*qc); 154 double siAt = sas_sinx_x(0.5*tA*qa); 155 double siBt = sas_sinx_x(0.5*tB*qb); 156 double siCt = sas_sinx_x(0.5*tC*qc); 162 const double tA = length_a + 2.0*thick_rim_a; 163 const double tB = length_b + 2.0*thick_rim_b; 164 const double tC = length_c + 2.0*thick_rim_c; 165 const double siA = sas_sinx_x(0.5*length_a*qa); 166 const double siB = sas_sinx_x(0.5*length_b*qb); 167 const double siC = sas_sinx_x(0.5*length_c*qc); 168 const double siAt = sas_sinx_x(0.5*tA*qa); 169 const double siBt = sas_sinx_x(0.5*tB*qb); 170 const double siCt = sas_sinx_x(0.5*tC*qc); 157 171 158 159 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 160 // in the 1D code, but should be checked! 161 double f = ( dr0*Vin*siA*siB*siC 162 + drA*VtA*(siAt-siA)*siB*siC 163 + drB*VtB*siA*(siBt-siB)*siC 164 + drC*VtC*siA*siB*(siCt-siC)); 172 #if OVERLAPPING 173 const double f = dr0*Vin*siA*siB*siC 174 + drA*(Vat*siAt-Va*siA)*siB*siC 175 + drB*siAt*(Vbt*siBt-Vb*siB)*siC 176 + drC*siAt*siBt*(Vct*siCt-Vc*siC); 177 #else 178 const double f = dr0*Vin*siA*siB*siC 179 + drA*(Vat*siAt-Va*siA)*siB*siC 180 + drB*siA*(Vbt*siBt-Vb*siB)*siC 181 + drC*siA*siB*(Vct*siCt-Vc*siC); 182 #endif 165 183 166 184 return 1.0e-4 * f * f; -
TabularUnified sasmodels/models/core_shell_parallelepiped.py ¶
r904cd9c rfc0b7aa 109 109 Equations (1), (13-14). (in German) 110 110 .. [#] D Singh (2009). *Small angle scattering studies of self assembly in 111 lipid mixtures*, John 's Hopkins University Thesis (2009) 223-225. `Available111 lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 112 112 from Proquest <http://search.proquest.com/docview/304915826?accountid 113 113 =26379>`_
Note: See TracChangeset
for help on using the changeset viewer.