Changeset 14838a3 in sasmodels


Ignore:
Timestamp:
Oct 15, 2016 12:37:09 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
2941abf
Parents:
4962519
Message:

update parallelepiped and core_shell_parallelepiped to use ORIENT macros

Location:
sasmodels
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/conversion_table.py

    r6831fa0 r14838a3  
    173173        "CSParallelepipedModel", 
    174174        { 
     175            "sld_core": "sld_pcore", 
     176            "sld_a": "sld_rimA", 
     177            "sld_b": "sld_rimB", 
     178            "sld_c": "sld_rimC", 
     179            "sld_solvent": "sld_solv", 
     180            "length_a": "shortA", 
     181            "length_b": "midB", 
     182            "length_c": "longC", 
     183            "thick_rim_a": "rimA", 
     184            "thick_rim_c": "rimC", 
     185            "thick_rim_b": "rimB", 
     186            "theta": "parallel_theta", 
    175187            "phi": "parallel_phi", 
    176188            "psi": "parallel_psi", 
    177             "sld_core": "sld_pcore", 
    178             "sld_c": "sld_rimC", 
    179             "sld_b": "sld_rimB", 
    180             "sld_solvent": "sld_solv", 
    181             "length_a": "shortA", 
    182             "sld_a": "sld_rimA", 
    183             "length_b": "midB", 
    184             "thick_rimc": "rimC", 
    185             "theta": "parallel_theta", 
    186             "thick_rim_a": "rimA", 
    187             "length_c": "longC", 
    188             "thick_rim_b": "rimB" 
    189189        } 
    190190    ], 

  • sasmodels/models/core_shell_parallelepiped.c

    r3a48772 r14838a3  
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    3636     
    37     double t1, t2, t3, t4, tmp, answer;    
    38     double mu = q * length_b; 
     37    const double mu = 0.5 * q * length_b; 
    3938     
    4039    //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);             
     40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    4241    //double vol = length_a * length_b * length_c +  
    4342    //       2.0 * thick_rim_a * length_b * length_c +  
    … …  
    4645     
    4746    // Scale sides by B 
    48     double a_scaled = length_a / length_b; 
    49     double c_scaled = length_c / length_b; 
     47    const double a_scaled = length_a / length_b; 
     48    const double c_scaled = length_c / length_b; 
    5049 
    51     // DelRho values (note that drC is not used later)        
    52         double dr0 = core_sld-solvent_sld; 
    53         double drA = arim_sld-solvent_sld; 
    54         double drB = brim_sld-solvent_sld; 
    55         //double drC = crim_sld-solvent_sld; 
     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; 
    5656 
    57     //Order of integration 
    58     int nordi=76;                                
    59     int nordj=76; 
    60      
    61     // outer integral (with nordi gauss points), integration limits = 0, 1 
    62     double summ = 0; //initialize integral 
     57    double Vin = length_a * length_b * length_c; 
     58    //double Vot = (length_a * length_b * length_c + 
     59    //            2.0 * thick_rim_a * length_b * length_c + 
     60    //            2.0 * length_a * thick_rim_b * length_c + 
     61    //            2.0 * length_a * length_b * thick_rim_c); 
     62    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
     63    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    6364 
    64     for( int i=0; i<nordi; i++) { 
    65                  
    66         // inner integral (with nordj gauss points), integration limits = 0, 1 
    67          
    68         double summj = 0.0; 
    69             double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    70                  
    71             for(int j=0; j<nordj; j++) { 
     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 drC_Vot = (crim_sld-solvent_sld)*Vot; 
    7270 
    73             double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    74             double mudum = mu * sqrt(1.0-sigma*sigma); 
     71    // Precompute scale factors for combining cross terms from the shape 
     72    const double scale23 = drhoA*V1; 
     73    const double scale14 = drhoB*V2; 
     74    const double scale12 = drho0*Vin - scale23 - scale14; 
    7575 
    76                 double Vin = length_a * length_b * length_c; 
    77                 //double Vot = (length_a * length_b * length_c + 
    78             //            2.0 * thick_rim_a * length_b * length_c + 
    79             //            2.0 * length_a * thick_rim_b * length_c + 
    80             //            2.0 * length_a * length_b * thick_rim_c); 
    81                 double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    82                 double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    83          
    84             // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
    85             // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
    86             // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
    87             // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)             
    88             double ta = (a_scaled+2.0*thick_rim_a)/length_b;  
    89             double tb = (a_scaled+2.0*thick_rim_b)/length_b; 
    90      
    91                 double arg1 = (0.5*mudum*a_scaled) * sin(M_PI_2*uu); 
    92                 double arg2 = (0.5*mudum) * cos(M_PI_2*uu); 
    93                 double arg3=  (0.5*mudum*ta) * sin(M_PI_2*uu); 
    94                 double arg4=  (0.5*mudum*tb) * cos(M_PI_2*uu); 
     76    // outer integral (with gauss points), integration limits = 0, 1 
     77    double outer_total = 0; //initialize integral 
    9578 
    96                 if(arg1==0.0){ 
    97                         t1 = 1.0; 
    98                 } else { 
    99                         t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)    
    100                 } 
    101                 if(arg2==0.0){ 
    102                         t2 = 1.0; 
    103                 } else { 
    104                         t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)     
    105                 }        
    106                 if(arg3==0.0){ 
    107                         t3 = 1.0; 
    108                 } else { 
    109                         t3 = sin(arg3)/arg3; 
    110                 } 
    111                 if(arg4==0.0){ 
    112                         t4 = 1.0; 
    113                 } else { 
    114                         t4 = sin(arg4)/arg4; 
    115                 } 
    116              
     79    for( int i=0; i<76; i++) { 
     80        double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     81        double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     82 
     83        // inner integral (with gauss points), integration limits = 0, 1 
     84        double inner_total = 0.0; 
     85        for(int j=0; j<76; j++) { 
     86            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     87            double sin_uu, cos_uu; 
     88            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     89            const double si1 = sinc(mu_proj * sin_uu * a_scaled); 
     90            const double si2 = sinc(mu_proj * cos_uu); 
     91            const double si3 = sinc(mu_proj * sin_uu * ta); 
     92            const double si4 = sinc(mu_proj * cos_uu * tb); 
     93 
    11794            // Expression in libCylinder.c (neither drC nor Vot are used) 
    118                 tmp =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors             
    119              
     95            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     96 
    12097            // To note also that in csparallelepiped.cpp, there is a function called 
    12198            // cspkernel, which seems to make more sense and has the following comment: 
    … …  
    126103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    127104             
    128             summj += Gauss76Wt[j] * tmp; 
     105            //  correct FF : sum of square of phase factors 
     106            inner_total += Gauss76Wt[j] * form * form; 
     107        } 
     108        inner_total *= 0.5; 
    129109 
    130         } 
    131                  
    132         // value of the inner integral 
    133         answer = 0.5 * summj; 
     110        // now sum up the outer integral 
     111        const double si = sinc(mu * c_scaled * sigma); 
     112        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     113    } 
     114    outer_total *= 0.5; 
    134115 
    135                 // finish the outer integral  
    136                 double arg = 0.5 * mu* c_scaled *sigma; 
    137                 if ( arg == 0.0 ) { 
    138                         answer *= 1.0; 
    139                 } else { 
    140                 answer *= sin(arg)*sin(arg)/arg/arg; 
    141                 } 
    142                  
    143                 // now sum up the outer integral 
    144                 summ += Gauss76Wt[i] * answer; 
    145  
    146     } 
    147      
    148         answer = 0.5 * summ; 
    149  
    150         //convert from [1e-12 A-1] to [cm-1] 
    151         answer *= 1.0e-4; 
    152          
    153         return answer; 
     116    //convert from [1e-12 A-1] to [cm-1] 
     117    return 1.0e-4 * outer_total; 
    154118} 
    155119 
    … …  
    170134    double psi) 
    171135{ 
    172     double tmp1, tmp2, tmp3, tmpt1, tmpt2, tmpt3;    
     136    double q, cos_val_a, cos_val_b, cos_val_c; 
     137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    173138 
    174     double q = sqrt(qx*qx+qy*qy); 
    175     double qx_scaled = qx/q; 
    176     double qy_scaled = qy/q; 
     139    // cspkernel in csparallelepiped recoded here 
     140    const double dr0 = core_sld-solvent_sld; 
     141    const double drA = arim_sld-solvent_sld; 
     142    const double drB = brim_sld-solvent_sld; 
     143    const double drC = crim_sld-solvent_sld; 
    177144 
    178     // Convert angles given in degrees to radians 
    179     theta *= M_PI_180; 
    180     phi   *= M_PI_180; 
    181     psi   *= M_PI_180; 
    182      
    183     // Parallelepiped c axis orientation 
    184     double cparallel_x = cos(theta) * cos(phi); 
    185     double cparallel_y = sin(theta); 
    186      
    187     // Compute angle between q and parallelepiped axis 
    188     double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 
    189  
    190     // Parallelepiped a axis orientation 
    191     double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    192     double parallel_y = sin(psi)*cos(theta); 
    193     double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 
    194  
    195     // Parallelepiped b axis orientation 
    196     double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    197     double bparallel_y = cos(theta)*cos(psi); 
    198     double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 
    199  
    200     // The following tests should always pass 
    201     if (fabs(cos_val_c)>1.0) { 
    202       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    203       cos_val_c = 1.0; 
    204     } 
    205     if (fabs(cos_val_a)>1.0) { 
    206       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    207       cos_val_a = 1.0; 
    208     } 
    209     if (fabs(cos_val_b)>1.0) { 
    210       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    211       cos_val_b = 1.0; 
    212     } 
    213      
    214     // cspkernel in csparallelepiped recoded here  
    215     double dr0 = core_sld-solvent_sld; 
    216         double drA = arim_sld-solvent_sld; 
    217         double drB = brim_sld-solvent_sld; 
    218         double drC = crim_sld-solvent_sld; 
    219         double Vin = length_a * length_b * length_c; 
     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; 
    220149    // As for the 1D case, Vot is not used 
    221         //double Vot = (length_a * length_b * length_c + 
     150    //double Vot = (length_a * length_b * length_c + 
    222151    //              2.0 * thick_rim_a * length_b * length_c + 
    223152    //              2.0 * length_a * thick_rim_b * length_c + 
    224153    //              2.0 * length_a * length_b * thick_rim_c); 
    225         double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    226         double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    227     double V3 = (2.0 * length_a * length_b * thick_rim_c); 
     154 
    228155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    229156    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    230157    // but for the moment I let it like this until understanding better the code. 
    231         double ta = length_a + 2.0*thick_rim_a; 
     158    double ta = length_a + 2.0*thick_rim_a; 
    232159    double tb = length_a + 2.0*thick_rim_b; 
    233160    double tc = length_a + 2.0*thick_rim_c; 
    234161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    235     double argA = 0.5*q*length_a*cos_val_a; 
    236     double argB = 0.5*q*length_b*cos_val_b; 
    237     double argC = 0.5*q*length_c*cos_val_c; 
    238     double argtA = 0.5*q*ta*cos_val_a; 
    239     double argtB = 0.5*q*tb*cos_val_b; 
    240     double argtC = 0.5*q*tc*cos_val_c; 
     162    double siA = sinc(0.5*q*length_a*cos_val_a); 
     163    double siB = sinc(0.5*q*length_b*cos_val_b); 
     164    double siC = sinc(0.5*q*length_c*cos_val_c); 
     165    double siAt = sinc(0.5*q*ta*cos_val_a); 
     166    double siBt = sinc(0.5*q*tb*cos_val_b); 
     167    double siCt = sinc(0.5*q*tc*cos_val_c); 
    241168     
    242     if(argA==0.0) { 
    243         tmp1 = 1.0; 
    244     } else { 
    245         tmp1 = sin(argA)/argA; 
    246     } 
    247     if (argB==0.0) { 
    248         tmp2 = 1.0; 
    249     } else { 
    250         tmp2 = sin(argB)/argB; 
    251     } 
    252     if (argC==0.0) { 
    253         tmp3 = 1.0; 
    254     } else { 
    255         tmp3 = sin(argC)/argC; 
    256     } 
    257     if(argtA==0.0) { 
    258         tmpt1 = 1.0; 
    259     } else { 
    260         tmpt1 = sin(argtA)/argtA; 
    261     } 
    262     if (argtB==0.0) { 
    263         tmpt2 = 1.0; 
    264     } else { 
    265         tmpt2 = sin(argtB)/argtB; 
    266     } 
    267     if (argtC==0.0) { 
    268         tmpt3 = 1.0; 
    269     } else { 
    270         tmpt3 = sin(argtC)*sin(argtC)/argtC/argtC; 
    271     } 
    272169 
    273170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    274171    // in the 1D code, but should be checked! 
    275     double f = ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1 +  
    276                drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3); 
     172    double f = ( dr0*siA*siB*siC*Vin 
     173               + drA*(siAt-siA)*siB*siC*V1 
     174               + drB*siA*(siBt-siB)*siC*V2 
     175               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    277176    
    278177    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2222134 r14838a3  
    164164# parameters for demo 
    165165demo = dict(scale=1, background=0.0, 
    166             sld_core=1e-6, sld_a=2e-6, sld_b=4e-6, 
    167             sld_c=2e-6, sld_solvent=6e-6, 
     166            sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, 
    168167            length_a=35, length_b=75, length_c=400, 
    169168            thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, 
    … …  
    177176            theta_pd=10, theta_pd_n=1, 
    178177            phi_pd=10, phi_pd_n=1, 
    179             psi_pd=10, psi_pd_n=10) 
     178            psi_pd=10, psi_pd_n=1) 
    180179 
    181180qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 

  • sasmodels/models/parallelepiped.c

    r3a48772 r14838a3  
    11double form_volume(double length_a, double length_b, double length_c); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, double length_b, double length_c); 
     2double Iq(double q, double sld, double solvent_sld, 
     3    double length_a, double length_b, double length_c); 
    34double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    4     double length_a, double length_b, double length_c, double theta, double phi, double psi); 
    5  
    6 // From Igor library 
    7 double _pkernel(double a, double b,double c, double ala, double alb, double alc); 
    8 double _pkernel(double a, double b,double c, double ala, double alb, double alc){ 
    9     double argA,argB,argC,tmp1,tmp2,tmp3; 
    10     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    11     argA = 0.5*a*ala; 
    12     argB = 0.5*b*alb; 
    13     argC = 0.5*c*alc; 
    14     if(argA==0.0) { 
    15         tmp1 = 1.0; 
    16     } else { 
    17         tmp1 = sin(argA)*sin(argA)/argA/argA; 
    18     } 
    19     if (argB==0.0) { 
    20         tmp2 = 1.0; 
    21     } else { 
    22         tmp2 = sin(argB)*sin(argB)/argB/argB; 
    23     } 
    24     if (argC==0.0) { 
    25         tmp3 = 1.0; 
    26     } else { 
    27         tmp3 = sin(argC)*sin(argC)/argC/argC; 
    28     } 
    29     return (tmp1*tmp2*tmp3); 
    30  
    31 } 
    32  
     5    double length_a, double length_b, double length_c, 
     6    double theta, double phi, double psi); 
    337 
    348double form_volume(double length_a, double length_b, double length_c) 
    … …  
    4519    double length_c) 
    4620{ 
    47     double tmp1, tmp2; 
    48      
    49     double mu = q * length_b; 
     21    const double mu = 0.5 * q * length_b; 
    5022     
    5123    // Scale sides by B 
    52     double a_scaled = length_a / length_b; 
    53     double c_scaled = length_c / length_b; 
     24    const double a_scaled = length_a / length_b; 
     25    const double c_scaled = length_c / length_b; 
    5426         
    55     //Order of integration 
    56     int nordi=76;                                
    57     int nordj=76; 
     27    // outer integral (with gauss points), integration limits = 0, 1 
     28    double outer_total = 0; //initialize integral 
    5829 
    59     // outer integral (with nordi gauss points), integration limits = 0, 1 
    60     double summ = 0; //initialize integral 
     30    for( int i=0; i<76; i++) { 
     31        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     32        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    6133 
    62     for( int i=0; i<nordi; i++) { 
    63                  
    64         // inner integral (with nordj gauss points), integration limits = 0, 1 
    65          
    66         double summj = 0.0; 
    67             double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );          
    68                  
    69             for(int j=0; j<nordj; j++) { 
     34        // inner integral (with gauss points), integration limits = 0, 1 
     35        // corresponding to angles from 0 to pi/2. 
     36        double inner_total = 0.0; 
     37        for(int j=0; j<76; j++) { 
     38            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     39            double sin_uu, cos_uu; 
     40            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
     41            const double si1 = sinc(mu_proj * sin_uu * a_scaled); 
     42            const double si2 = sinc(mu_proj * cos_uu); 
     43            inner_total += Gauss76Wt[j] * square(si1 * si2); 
     44        } 
     45        inner_total *= 0.5; 
    7046 
    71             double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    72             double mudum = mu * sqrt(1.0-sigma*sigma); 
    73                 double arg1 = 0.5 * mudum * cos(M_PI_2*uu); 
    74                 double arg2 = 0.5 * mudum * a_scaled * sin(M_PI_2*uu); 
    75             if(arg1==0.0) { 
    76                 tmp1 = 1.0; 
    77             } else { 
    78                 tmp1 = sin(arg1)*sin(arg1)/arg1/arg1; 
    79             } 
    80             if (arg2==0.0) { 
    81                 tmp2 = 1.0; 
    82             } else { 
    83                 tmp2 = sin(arg2)*sin(arg2)/arg2/arg2; 
    84             } 
     47        const double si = sinc(mu * c_scaled * sigma); 
     48        outer_total += Gauss76Wt[i] * inner_total * si * si; 
     49    } 
     50    outer_total *= 0.5; 
    8551 
    86             summj += Gauss76Wt[j] * tmp1 * tmp2; 
    87         } 
    88                  
    89         // value of the inner integral 
    90         double answer = 0.5 * summj; 
    91  
    92         double arg = 0.5 * mu * c_scaled * sigma; 
    93         if ( arg == 0.0 ) { 
    94             answer *= 1.0; 
    95         } else { 
    96             answer *= sin(arg)*sin(arg)/arg/arg; 
    97         } 
    98                  
    99             // sum of outer integral 
    100         summ += Gauss76Wt[i] * answer; 
    101          
    102     }    
    103     
    104     const double vd = (sld-solvent_sld) * form_volume(length_a, length_b, length_c); 
    105      
    106     // convert from [1e-12 A-1] to [cm-1] and 0.5 factor for outer integral 
    107     return 1.0e-4 * 0.5 * vd * vd * summ; 
    108      
     52    // Multiply by contrast^2 and convert from [1e-12 A-1] to [cm-1] 
     53    const double V = form_volume(length_a, length_b, length_c); 
     54    const double drho = (sld-solvent_sld); 
     55    return 1.0e-4 * square(drho * V) * outer_total; 
    10956} 
    11057 
    … …  
    12067    double psi) 
    12168{ 
    122     double q = sqrt(qx*qx+qy*qy); 
    123     double qx_scaled = qx/q; 
    124     double qy_scaled = qy/q; 
     69    double q, cos_val_a, cos_val_b, cos_val_c; 
     70    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 
    12571 
    126     // Convert angles given in degrees to radians 
    127     theta *= M_PI_180; 
    128     phi   *= M_PI_180; 
    129     psi   *= M_PI_180; 
    130      
    131     // Parallelepiped c axis orientation 
    132     double cparallel_x = cos(theta) * cos(phi); 
    133     double cparallel_y = sin(theta); 
    134      
    135     // Compute angle between q and parallelepiped axis 
    136     double cos_val_c = cparallel_x*qx_scaled + cparallel_y*qy_scaled;// + cparallel_z*qz; 
    137  
    138     // Parallelepiped a axis orientation 
    139     double parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    140     double parallel_y = sin(psi)*cos(theta); 
    141     double cos_val_a = parallel_x*qx_scaled + parallel_y*qy_scaled; 
    142  
    143     // Parallelepiped b axis orientation 
    144     double bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    145     double bparallel_y = cos(theta)*cos(psi); 
    146     double cos_val_b = bparallel_x*qx_scaled + bparallel_y*qy_scaled; 
    147  
    148     // The following tests should always pass 
    149     if (fabs(cos_val_c)>1.0) { 
    150       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    151       cos_val_c = 1.0; 
    152     } 
    153     if (fabs(cos_val_a)>1.0) { 
    154       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    155       cos_val_a = 1.0; 
    156     } 
    157     if (fabs(cos_val_b)>1.0) { 
    158       //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    159       cos_val_b = 1.0; 
    160     } 
    161      
    162     // Call the IGOR library function to get the kernel 
    163     double form = _pkernel( q*length_a, q*length_b, q*length_c, cos_val_a, cos_val_b, cos_val_c); 
    164    
    165     // Multiply by contrast^2 
    166     const double vd = (sld - solvent_sld) * form_volume(length_a, length_b, length_c); 
    167     return 1.0e-4 * vd * vd * form; 
     72    const double siA = sinc(0.5*q*length_a*cos_val_a); 
     73    const double siB = sinc(0.5*q*length_b*cos_val_b); 
     74    const double siC = sinc(0.5*q*length_c*cos_val_c); 
     75    const double V = form_volume(length_a, length_b, length_c); 
     76    const double drho = (sld - solvent_sld); 
     77    const double form = V * drho * siA * siB * siC; 
     78    // Square and convert from [1e-12 A-1] to [cm-1] 
     79    return 1.0e-4 * form * form; 
    16880} 
Note: See TracChangeset for help on using the changeset viewer.