Changeset 14838a3 in sasmodels for sasmodels/models/parallelepiped.c


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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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.