Changeset 2a0b2b1 in sasmodels for sasmodels/models/fcc_paracrystal.c


Ignore:
Timestamp:
Apr 14, 2017 8:30:29 AM (7 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:
fdd56a1
Parents:
9901384
Message:

restructure all 2D models to work with (qa,qb,qc) = rotate(qx,qy) rather than working with angles directly in preparation for revised jitter algorithm

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r2a0b2b1  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2_sq_fcc(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Rewriting equations for efficiency, accuracy and readability, and so 
     5    // code is reusable between 1D and 2D models. 
     6    const double a1 = qb + qa; 
     7    const double a2 = qa + qc; 
     8    const double a3 = qb + qc; 
    69 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     10    const double half_dnn = 0.5*dnn; 
     11    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     12 
     13    // Numerator: (1 - exp(a)^2)^3 
     14    //         => (-(exp(2a) - 1))^3 
     15    //         => -expm1(2a)^3 
     16    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2) 
     17    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1 
     18    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1 
     19    const double exp_arg = exp(-arg); 
     20    const double Sq = -cube(expm1(-2.0*arg)) 
     21        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0)); 
     24 
     25    return Sq; 
     26} 
    927 
    1028 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     29// occupied volume fraction calculated from lattice symmetry and sphere radius 
     30static double 
     31_fcc_volume_fraction(double radius, double dnn) 
     32{ 
     33    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1934} 
    2035 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     36static double 
     37form_volume(double radius) 
     38{ 
    4339    return sphere_volume(radius); 
    4440} 
    4541 
    4642 
    47 double Iq(double q, double dnn, 
     43static double Iq(double q, double dnn, 
    4844  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     45  double sld, double solvent_sld) 
     46{ 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI; 
     49    const double phi_b = M_PI; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_2; 
     52    const double theta_b = M_PI_2; 
    5053 
    51         //Volume fraction calculated from lattice symmetry and sphere radius 
    52         const double s1 = dnn*sqrt(2.0); 
    53         const double latticescale = 4.0*sphere_volume(radius/s1); 
     54    double outer_sum = 0.0; 
     55    for(int i=0; i<150; i++) { 
     56        double inner_sum = 0.0; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        double sin_theta, cos_theta; 
     59        SINCOS(theta, sin_theta, cos_theta); 
     60        const double qc = q*cos_theta; 
     61        const double qab = q*sin_theta; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     64            double sin_phi, cos_phi; 
     65            SINCOS(phi, sin_phi, cos_phi); 
     66            const double qa = qab*cos_phi; 
     67            const double qb = qab*sin_phi; 
     68            const double fq = _sq_fcc(qa, qb, qc, dnn, d_factor); 
     69            inner_sum += Gauss150Wt[j] * fq; 
     70        } 
     71        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73    } 
     74    outer_sum *= theta_m; 
     75    const double Sq = outer_sum/(4.0*M_PI); 
     76    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5477 
    55     const double va = 0.0; 
    56     const double vb = 2.0*M_PI; 
    57     const double vaj = 0.0; 
    58     const double vbj = M_PI; 
    59  
    60     double summ = 0.0; 
    61     double answer = 0.0; 
    62         for(int i=0; i<150; i++) { 
    63                 //setup inner integral over the ellipsoidal cross-section 
    64                 double summj=0.0; 
    65                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    66                 for(int j=0;j<150;j++) { 
    67                         //20 gauss points for the inner integral 
    68                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    69                         double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    70                         summj += yyy; 
    71                 } 
    72                 //now calculate the value of the inner integral 
    73                 double answer = (vbj-vaj)/2.0*summj; 
    74  
    75                 //now calculate outer integral 
    76                 summ = summ+(Gauss150Wt[i] * answer); 
    77         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    78  
    79         answer = (vb-va)/2.0*summ; 
    80         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    81  
    82     return answer; 
     78    return _fcc_volume_fraction(radius, dnn) * Pq * Sq; 
     79} 
    8380 
    8481 
    85 } 
    86  
    87 double Iqxy(double qx, double qy, 
     82static double Iqxy(double qx, double qy, 
    8883    double dnn, double d_factor, double radius, 
    8984    double sld, double solvent_sld, 
    … …  
    9287    double q, zhat, yhat, xhat; 
    9388    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
    9492 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
    98     const double qd = 0.5*q*dnn; 
    99     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    100     const double tanh_qd = tanh(arg); 
    101     const double cosh_qd = cosh(arg); 
    102     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    103                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    105  
    106     //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 
    107  
    108     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    109     //the occupied volume of the lattice 
    110     const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    111     return lattice_scale * Fq; 
     93    q = sqrt(qa*qa + qb*qb + qc*qc); 
     94    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     95    const double Sq = _sq_fcc(qa, qb, qc, dnn, d_factor); 
     96    return _fcc_volume_fraction(radius, dnn) * Pq * Sq; 
    11297} 
Note: See TracChangeset for help on using the changeset viewer.