Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rc69d6d6  
    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 
    7 static double 
    8 form_volume(double length_a, double length_b, double length_c, 
    9     double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     1double form_volume(double length_a, double length_b, double length_c, 
     2                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
     3double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
     4          double solvent_sld, double length_a, double length_b, double length_c, 
     5          double thick_rim_a, double thick_rim_b, double thick_rim_c); 
     6double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, 
     7            double crim_sld, double solvent_sld, double length_a, double length_b, 
     8            double length_c, double thick_rim_a, double thick_rim_b, 
     9            double thick_rim_c, double theta, double phi, double psi); 
     10 
     11double form_volume(double length_a, double length_b, double length_c, 
     12                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1013{ 
    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 
     14    //return length_a * length_b * length_c; 
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
     17           2.0 * thick_rim_b * length_a * length_c + 
     18           2.0 * thick_rim_c * length_a * length_b; 
    2719} 
    2820 
    29 static double 
    30 Iq(double q, 
     21double Iq(double q, 
    3122    double core_sld, 
    3223    double arim_sld, 
    … …  
    4132    double thick_rim_c) 
    4233{ 
    43     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
     34    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    4435    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    45     // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker) 
    46     // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK) 
     36    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    4737 
    48     const double half_q = 0.5*q; 
     38    const double mu = 0.5 * q * length_b; 
    4939 
    50     const double tA = length_a + 2.0*thick_rim_a; 
    51     const double tB = length_b + 2.0*thick_rim_b; 
    52     const double tC = length_c + 2.0*thick_rim_c; 
     40    //calculate volume before rescaling (in original code, but not used) 
     41    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     42    //double vol = length_a * length_b * length_c + 
     43    //       2.0 * thick_rim_a * length_b * length_c + 
     44    //       2.0 * thick_rim_b * length_a * length_c + 
     45    //       2.0 * thick_rim_c * length_a * length_b; 
    5346 
    54     // Scale factors 
    55     const double dr0 = (core_sld-solvent_sld); 
    56     const double drA = (arim_sld-solvent_sld); 
    57     const double drB = (brim_sld-solvent_sld); 
    58     const double drC = (crim_sld-solvent_sld); 
     47    // Scale sides by B 
     48    const double a_scaled = length_a / length_b; 
     49    const double c_scaled = length_c / length_b; 
     50 
     51    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
     52    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     53    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
     54 
     55    double Vin = length_a * length_b * length_c; 
     56    //double Vot = (length_a * length_b * length_c + 
     57    //            2.0 * thick_rim_a * length_b * length_c + 
     58    //            2.0 * length_a * thick_rim_b * length_c + 
     59    //            2.0 * length_a * length_b * thick_rim_c); 
     60    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
     61    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     62    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
     63    double Vot = Vin + V1 + V2 + V3; 
     64 
     65    // Scale factors (note that drC is not used later) 
     66    const double drho0 = (core_sld-solvent_sld); 
     67    const double drhoA = (arim_sld-solvent_sld); 
     68    const double drhoB = (brim_sld-solvent_sld); 
     69    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     70 
     71 
     72    // Precompute scale factors for combining cross terms from the shape 
     73    const double scale23 = drhoA*V1; 
     74    const double scale14 = drhoB*V2; 
     75    const double scale24 = drhoC*V3; 
     76    const double scale11 = drho0*Vin; 
     77    const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
    5978 
    6079    // outer integral (with gauss points), integration limits = 0, 1 
    61     double outer_sum = 0; //initialize integral 
    62     for( int i=0; i<GAUSS_N; i++) { 
    63         const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    64         const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
     80    double outer_total = 0; //initialize integral 
    6581 
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    67         const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    68         const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
    69         double inner_sum = 0.0; 
    70         for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    72             double sin_beta, cos_beta; 
    73             SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
    74             const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    75             const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    76             const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
    77             const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
     82    for( int i=0; i<76; i++) { 
     83        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     84        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    7885 
    79 #if OVERLAPPING 
    80             const double f = dr0*siA*siB*siC 
    81                 + drA*(siAt-siA)*siB*siC 
    82                 + drB*siAt*(siBt-siB)*siC 
    83                 + drC*siAt*siBt*(siCt-siC); 
    84 #else 
    85             const double f = dr0*siA*siB*siC 
    86                 + drA*(siAt-siA)*siB*siC 
    87                 + drB*siA*(siBt-siB)*siC 
    88                 + drC*siA*siB*(siCt-siC); 
    89 #endif 
     86        // inner integral (with gauss points), integration limits = 0, 1 
     87        double inner_total = 0.0; 
     88        double inner_total_crim = 0.0; 
     89        for(int j=0; j<76; j++) { 
     90            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     91            double sin_uu, cos_uu; 
     92            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     93            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
     94            const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
     95            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
     96            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
    9097 
    91             inner_sum += GAUSS_W[j] * f * f; 
     98            // Expression in libCylinder.c (neither drC nor Vot are used) 
     99            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     100            const double form_crim = scale11*si1*si2; 
     101 
     102 
     103            //  correct FF : sum of square of phase factors 
     104            inner_total += Gauss76Wt[j] * form * form; 
     105            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
    92106        } 
    93         inner_sum *= 0.5; 
     107        inner_total *= 0.5; 
     108        inner_total_crim *= 0.5; 
    94109        // now sum up the outer integral 
    95         outer_sum += GAUSS_W[i] * inner_sum; 
     110        const double si = sas_sinx_x(mu * c_scaled * sigma); 
     111        const double si_crim = sas_sinx_x(mu * tc * sigma); 
     112        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
    96113    } 
    97     outer_sum *= 0.5; 
     114    outer_total *= 0.5; 
    98115 
    99116    //convert from [1e-12 A-1] to [cm-1] 
    100     return 1.0e-4 * outer_sum; 
     117    return 1.0e-4 * outer_total; 
    101118} 
    102119 
    103 static double 
    104 Iqabc(double qa, double qb, double qc, 
     120double Iqxy(double qx, double qy, 
    105121    double core_sld, 
    106122    double arim_sld, 
    … …  
    113129    double thick_rim_a, 
    114130    double thick_rim_b, 
    115     double thick_rim_c) 
     131    double thick_rim_c, 
     132    double theta, 
     133    double phi, 
     134    double psi) 
    116135{ 
     136    double q, zhat, yhat, xhat; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138 
    117139    // cspkernel in csparallelepiped recoded here 
    118140    const double dr0 = core_sld-solvent_sld; 
    … …  
    121143    const double drC = crim_sld-solvent_sld; 
    122144 
    123     const double tA = length_a + 2.0*thick_rim_a; 
    124     const double tB = length_b + 2.0*thick_rim_b; 
    125     const double tC = length_c + 2.0*thick_rim_c; 
    126     const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
    127     const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
    128     const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
    129     const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
    130     const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
    131     const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
     145    double Vin = length_a * length_b * length_c; 
     146    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc) 
     147    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     148    double V3 = 2.0 * length_a * length_b * thick_rim_c; 
     149    // As for the 1D case, Vot is not used 
     150    //double Vot = (length_a * length_b * length_c + 
     151    //              2.0 * thick_rim_a * length_b * length_c + 
     152    //              2.0 * length_a * thick_rim_b * length_c + 
     153    //              2.0 * length_a * length_b * thick_rim_c); 
    132154 
    133 #if OVERLAPPING 
    134     const double f = dr0*siA*siB*siC 
    135         + drA*(siAt-siA)*siB*siC 
    136         + drB*siAt*(siBt-siB)*siC 
    137         + drC*siAt*siBt*(siCt-siC); 
    138 #else 
    139     const double f = dr0*siA*siB*siC 
    140         + drA*(siAt-siA)*siB*siC 
    141         + drB*siA*(siBt-siB)*siC 
    142         + drC*siA*siB*(siCt-siC); 
    143 #endif 
     155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
     156    // the scaling by B. 
     157    double ta = length_a + 2.0*thick_rim_a; 
     158    double tb = length_b + 2.0*thick_rim_b; 
     159    double tc = length_c + 2.0*thick_rim_c; 
     160    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
     161    double siA = sas_sinx_x(0.5*q*length_a*xhat); 
     162    double siB = sas_sinx_x(0.5*q*length_b*yhat); 
     163    double siC = sas_sinx_x(0.5*q*length_c*zhat); 
     164    double siAt = sas_sinx_x(0.5*q*ta*xhat); 
     165    double siBt = sas_sinx_x(0.5*q*tb*yhat); 
     166    double siCt = sas_sinx_x(0.5*q*tc*zhat); 
     167 
     168 
     169    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
     170    // in the 1D code, but should be checked! 
     171    double f = ( dr0*siA*siB*siC*Vin 
     172               + drA*(siAt-siA)*siB*siC*V1 
     173               + drB*siA*(siBt-siB)*siC*V2 
     174               + drC*siA*siB*(siCt-siC)*V3); 
    144175 
    145176    return 1.0e-4 * f * f; 
Note: See TracChangeset for help on using the changeset viewer.