Ignore:
Timestamp:
Nov 5, 2017 12:33:25 PM (5 years ago)
Author:
Paul Kienzle <pkienzle@…>
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
Message:

core_shell_parallelepiped: modify calculation to be consistent with hollow rectangular prism

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 
    17static double 
    28form_volume(double length_a, double length_b, double length_c, 
    39    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    410{ 
    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 
    927} 
    1028 
     
    3856 
    3957    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; 
    4373 
    4474    // Scale factors (note that drC is not used later) 
     
    4777    const double drB = (brim_sld-solvent_sld); 
    4878    const double drC = (crim_sld-solvent_sld); 
    49  
    50     // Precompute scale factors for combining cross terms from the shape 
    51     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 out 
    60     // 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_core 
    68     da = delta rho_rimA V_rimA 
    69     db = delta rho_rimB V_rimB 
    70     dc = delta rho_rimC V_rimC 
    71     sa = j_0 (q_a a/2)    // siA the code 
    72     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 code 
    75     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 the 
    79     // outer loop integrating over [0,1] after the u-substitution 
    80     //    sigma = cos(theta), sqrt(1-sigma^2) = sin(theta) 
    81     // and inner loop integrating over [0, pi/2] as 
    82     //    uu = phi 
    83  
    84     ************************************************************  */ 
    8579 
    8680    // outer integral (with gauss points), integration limits = 0, 1 
     
    10397            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b); 
    10498 
    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 
    107110 
    108             inner_sum += Gauss76Wt[j] * form * form; 
     111            inner_sum += Gauss76Wt[j] * f * f; 
    109112        } 
    110113        inner_sum *= 0.5; 
     
    139142 
    140143    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; 
    144159 
    145160    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    146161    // 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); 
    157171 
    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 
    165183 
    166184    return 1.0e-4 * f * f; 
Note: See TracChangeset for help on using the changeset viewer.