Changes in / [edb0f85:bd36af0] in sasmodels


Ignore:
Location:
sasmodels
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/generate.py

    r30b60d2 r6db17bd  
    370370    """ 
    371371    # Note: need 0xffffffff&val to force an unsigned 32-bit number 
     372    try: 
     373        source = source.encode('utf8') 
     374    except AttributeError: # bytes has no encode attribute in python 3 
     375        pass 
    372376    return "%08X"%(0xffffffff&crc32(source)) 
    373377 
  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c rc69d6d6  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
     
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //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 +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
     
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     37 
    3738    const double mu = 0.5 * q * length_b; 
    38      
     39 
    3940    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     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 + 
    4344    //       2.0 * thick_rim_b * length_a * length_c + 
    4445    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     46 
    4647    // Scale sides by B 
    4748    const double a_scaled = length_a / length_b; 
    4849    const double c_scaled = length_c / length_b; 
    4950 
    50     // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
    51     // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
    52     // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
    53     // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 
    54     double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    55     double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     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 
    5654 
    5755    double Vin = length_a * length_b * length_c; 
     
    6260    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    6361    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; 
    6464 
    6565    // Scale factors (note that drC is not used later) 
     
    6767    const double drhoA = (arim_sld-solvent_sld); 
    6868    const double drhoB = (brim_sld-solvent_sld); 
    69     //const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     69    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     70 
    7071 
    7172    // Precompute scale factors for combining cross terms from the shape 
    7273    const double scale23 = drhoA*V1; 
    7374    const double scale14 = drhoB*V2; 
    74     const double scale12 = drho0*Vin - scale23 - scale14; 
     75    const double scale24 = drhoC*V3; 
     76    const double scale11 = drho0*Vin; 
     77    const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
    7578 
    7679    // outer integral (with gauss points), integration limits = 0, 1 
     
    8386        // inner integral (with gauss points), integration limits = 0, 1 
    8487        double inner_total = 0.0; 
     88        double inner_total_crim = 0.0; 
    8589        for(int j=0; j<76; j++) { 
    8690            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     
    8892            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    8993            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    90             const double si2 = sas_sinx_x(mu_proj * cos_uu); 
     94            const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    9195            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    9296            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     
    9498            // Expression in libCylinder.c (neither drC nor Vot are used) 
    9599            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     100            const double form_crim = scale11*si1*si2; 
    96101 
    97             // To note also that in csparallelepiped.cpp, there is a function called 
    98             // cspkernel, which seems to make more sense and has the following comment: 
    99             //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010. 
    100             //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)* 
    101             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102             // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     102 
    105103            //  correct FF : sum of square of phase factors 
    106104            inner_total += Gauss76Wt[j] * form * form; 
     105            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
    107106        } 
    108107        inner_total *= 0.5; 
    109  
     108        inner_total_crim *= 0.5; 
    110109        // now sum up the outer integral 
    111110        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    112         outer_total += Gauss76Wt[i] * inner_total * si * si; 
     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); 
    113113    } 
    114114    outer_total *= 0.5; 
     
    154154 
    155155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    156     // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    157     // but for the moment I let it like this until understanding better the code. 
     156    // the scaling by B. 
    158157    double ta = length_a + 2.0*thick_rim_a; 
    159     double tb = length_a + 2.0*thick_rim_b; 
    160     double tc = length_a + 2.0*thick_rim_c; 
     158    double tb = length_b + 2.0*thick_rim_b; 
     159    double tc = length_c + 2.0*thick_rim_c; 
    161160    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162161    double siA = sas_sinx_x(0.5*q*length_a*xhat); 
     
    166165    double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167166    double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     167 
    169168 
    170169    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
     
    173172               + drA*(siAt-siA)*siB*siC*V1 
    174173               + drB*siA*(siBt-siB)*siC*V2 
    175                + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     174               + drC*siA*siB*(siCt-siC)*V3); 
     175 
    177176    return 1.0e-4 * f * f; 
    178177} 
Note: See TracChangeset for help on using the changeset viewer.