Changeset 11ca2ab in sasmodels


Ignore:
Timestamp:
Oct 14, 2016 1:13:55 AM (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:
5abbad7
Parents:
a80e64c
git-author:
Paul Kienzle <pkienzle@…> (10/13/16 23:39:19)
git-committer:
Paul Kienzle <pkienzle@…> (10/14/16 01:13:55)
Message:

add macros for theta-phi-psi orientation

Location:
sasmodels
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/kernel_header.c

    re1d6983 r11ca2ab  
    148148inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
     150 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
     151     SINCOS(phi*M_PI_180, sn, cn); \ 
     152     q = sqrt(qx*qx + qy*qy); \ 
     153     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \ 
     154     sn = sqrt(1 - cn*cn); \ 
     155 } while (0) 
     156 
     157 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 
     158    q = sqrt(qx*qx + qy*qy); \ 
     159    const double qxhat = qx/q; \ 
     160    const double qyhat = qy/q; \ 
     161    double sin_theta, cos_theta; \ 
     162    double sin_phi, cos_phi; \ 
     163    double sin_psi, cos_psi; \ 
     164    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 
     165    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 
     166    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 
     167    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 
     168    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 
     169    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 
     170 } while (0) 

  • sasmodels/models/barbell.c

    r0d6e865 r11ca2ab  
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    4040    const double integral = total*zm; 
    41     const double bell_Fq = 2*M_PI*cube(radius_bell)*integral; 
    42     return bell_Fq; 
     41    const double bell_fq = 2*M_PI*cube(radius_bell)*integral; 
     42    return bell_fq; 
    4343} 
     44 
     45static double 
     46_fq(double q, double h, 
     47    double radius_bell, double radius, double half_length, 
     48    double sin_alpha, double cos_alpha) 
     49{ 
     50    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
     51    const double bj = sas_J1c(q*radius*sin_alpha); 
     52    const double si = sinc(q*half_length*cos_alpha); 
     53    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
     54    const double Aq = bell_fq + cyl_fq; 
     55    return Aq; 
     56} 
     57 
    4458 
    4559double form_volume(double radius_bell, 
    … …  
    4761        double length) 
    4862{ 
    49  
    5063    // bell radius should never be less than radius when this is called 
    5164    const double hdist = sqrt(square(radius_bell) - square(radius)); 
    … …  
    7184        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7285        SINCOS(alpha, sin_alpha, cos_alpha); 
    73  
    74         const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    75         const double bj = sas_J1c(q*radius*sin_alpha); 
    76         const double si = sinc(q*half_length*cos_alpha); 
    77         const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    78         const double Aq = bell_Fq + cyl_Fq; 
     86        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
    7987        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8088    } 
    … …  
    93101        double theta, double phi) 
    94102{ 
    95     // Compute angle alpha between q and the cylinder axis 
    96     double sn, cn; // slots to hold sincos function output 
    97     SINCOS(phi*M_PI_180, sn, cn); 
    98     const double q = sqrt(qx*qx + qy*qy); 
    99     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    100     const double alpha = acos(cos_val); // rod angle relative to q 
     103    double q, sin_alpha, cos_alpha; 
     104    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    101105 
    102106    const double h = -sqrt(square(radius_bell) - square(radius)); 
    103     const double half_length = 0.5*length; 
    104  
    105     double sin_alpha, cos_alpha; // slots to hold sincos function output 
    106     SINCOS(alpha, sin_alpha, cos_alpha); 
    107     const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    108     const double bj = sas_J1c(q*radius*sin_alpha); 
    109     const double si = sinc(q*half_length*cos_alpha); 
    110     const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    111     const double Aq = cyl_Fq + bell_Fq; 
     107    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
    112108 
    113109    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r0bef47b r11ca2ab  
    3131        const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    3232        const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
     33 
    3334        const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    3435        result = pow(1.0-(temp10*temp10),3)*temp6 
    … …  
    8182 
    8283    return answer; 
    83  
    84  
    8584} 
    8685 
    8786 
    88 double Iqxy(double qx, double qy, double dnn, 
    89     double d_factor, double radius,double sld, double solvent_sld, 
    90     double theta, double phi, double psi){ 
     87double Iqxy(double qx, double qy, 
     88    double dnn, double d_factor, double radius, 
     89    double sld, double solvent_sld, 
     90    double theta, double phi, double psi) 
     91{ 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    9194 
    92   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 
    93   //double q_z; 
    94   double cos_val_b3, cos_val_b2, cos_val_b1; 
    95   double a1_dot_q, a2_dot_q,a3_dot_q; 
    96   double answer; 
    97   double Zq, Fkq, Fkq_2; 
     95    const double a1 = +cos_a3 - cos_a1 + cos_a2; 
     96    const double a2 = +cos_a3 + cos_a1 - cos_a2; 
     97    const double a3 = -cos_a3 + cos_a1 + cos_a2; 
    9898 
    99   //convert to q and make scaled values 
    100   double q = sqrt(qx*qx+qy*qy); 
    101   double q_x = qx/q; 
    102   double q_y = qy/q; 
     99    const double qd = 0.5*q*dnn; 
     100    const double exp_qd = exp(0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3)); 
     101    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd; 
     102    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd; 
    103103 
    104   //convert angle degree to radian 
    105   theta = theta * M_PI_180; 
    106   phi = phi * M_PI_180; 
    107   psi = psi * M_PI_180; 
     104    const double Zq = sinh_qd/(cosh_qd - cos(qd*a1)) 
     105                    * sinh_qd/(cosh_qd - cos(qd*a2)) 
     106                    * sinh_qd/(cosh_qd - cos(qd*a3)); 
    108107 
    109   const double Da = d_factor*dnn; 
    110   const double s1 = dnn/sqrt(0.75); 
    111  
    112  
    113   //the occupied volume of the lattice 
    114   const double latticescale = 2.0*sphere_volume(radius/s1); 
    115   // q vector 
    116   //q_z = 0.0; // for SANS; assuming qz is negligible 
    117   /// Angles here are respect to detector coordinate 
    118   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 
    119     // b3 axis orientation 
    120     b3_x = cos(theta) * cos(phi); 
    121     b3_y = sin(theta); 
    122     //b3_z = -cos(theta) * sin(phi); 
    123     cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 
    124  
    125     //alpha = acos(cos_val_b3); 
    126     // b1 axis orientation 
    127     b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    128     b1_y = sin(psi)*cos(theta); 
    129     cos_val_b1 = b1_x*q_x + b1_y*q_y; 
    130     // b2 axis orientation 
    131     b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    132         b2_y = cos(theta)*cos(psi); 
    133     cos_val_b2 = b2_x*q_x + b2_y*q_y; 
    134  
    135     // The following test should always pass 
    136     if (fabs(cos_val_b3)>1.0) { 
    137       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    138       cos_val_b3 = 1.0; 
    139     } 
    140     if (fabs(cos_val_b2)>1.0) { 
    141       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    142       cos_val_b2 = 1.0; 
    143     } 
    144     if (fabs(cos_val_b1)>1.0) { 
    145       //printf("bcc_ana_2D: Unexpected error: cos()>1\n"); 
    146       cos_val_b1 = 1.0; 
    147     } 
    148     // Compute the angle btw vector q and the a3 axis 
    149     a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 
    150  
    151     // a1 axis 
    152     a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 
    153  
    154     // a2 axis 
    155     a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 
    156  
    157  
    158     // Get Fkq and Fkq_2 
    159     Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
    160     Fkq_2 = Fkq*Fkq; 
    161     // Call Zq=Z1*Z2*Z3 
    162     Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
    163     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
    164     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
    165  
    166   // Use SphereForm directly from libigor 
    167   answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 
    168  
    169   return answer; 
    170  } 
     108    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
     109    //the occupied volume of the lattice 
     110    const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
     111    return lattice_scale * Fq; 
     112} 

  • sasmodels/models/cylinder.c

    r0d6e865 r11ca2ab  
    3333        // alpha(theta,phi) the projection of the cylinder on the detector plane 
    3434        SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn; 
     35        total += Gauss76Wt[i] * square(fq(q, sn, cn, radius, length)) * sn; 
    3636    } 
    3737    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5858    double phi) 
    5959{ 
    60     double sn, cn; // slots to hold sincos function output 
    61  
    62     // Compute angle alpha between q and the cylinder axis 
    63     SINCOS(phi*M_PI_180, sn, cn); 
    64     const double q = sqrt(qx*qx + qy*qy); 
    65     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    66  
    67     const double alpha = acos(cos_val); 
    68  
    69     SINCOS(alpha, sn, cn); 
     60    double q, sin_alpha, cos_alpha; 
     61    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    7062    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    71     return 1.0e-4 * square(s * fq(q, sn, cn, radius, length)); 
     63    return 1.0e-4 * square(s * fq(q, sin_alpha, cos_alpha, radius, length)); 
    7264} 

  • sasmodels/models/fcc_paracrystal.c

    r0bef47b r11ca2ab  
    8585} 
    8686 
     87double Iqxy(double qx, double qy, 
     88    double dnn, double d_factor, double radius, 
     89    double sld, double solvent_sld, 
     90    double theta, double phi, double psi) 
     91{ 
     92    double q, cos_a1, cos_a2, cos_a3; 
     93    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    8794 
    88 double Iqxy(double qx, double qy, double dnn, 
    89     double d_factor, double radius,double sld, double solvent_sld, 
    90     double theta, double phi, double psi){ 
     95    const double a1 = +cos_a3 + cos_a2; 
     96    const double a2 = +cos_a3 + cos_a1; 
     97    const double a3 = -cos_a3 + cos_a2; 
     98    const double qd = 0.5*q*dnn; 
     99    const double exp_qd = exp(0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3)); 
     100    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd; 
     101    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd; 
    91102 
    92   double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z, 
    93   // double q_z; 
    94   double cos_val_b3, cos_val_b2, cos_val_b1; 
    95   double a1_dot_q, a2_dot_q,a3_dot_q; 
    96   double answer; 
    97   double Zq, Fkq, Fkq_2; 
     103    const double Zq = sinh_qd/(cosh_qd - cos(qd*a1)) 
     104                    * sinh_qd/(cosh_qd - cos(qd*a2)) 
     105                    * sinh_qd/(cosh_qd - cos(qd*a3)); 
    98106 
    99   //convert to q and make scaled values 
    100   double q = sqrt(qx*qx+qy*qy); 
    101   double q_x = qx/q; 
    102   double q_y = qy/q; 
    103  
    104   //convert angle degree to radian 
    105   theta = theta * M_PI_180; 
    106   phi = phi * M_PI_180; 
    107   psi = psi * M_PI_180; 
    108  
    109   const double Da = d_factor*dnn; 
    110   const double s1 = dnn/sqrt(0.75); 
    111  
    112  
    113   //the occupied volume of the lattice 
    114   const double latticescale = 2.0*sphere_volume(radius/s1); 
    115   // q vector 
    116   // q_z = 0.0; // for SANS; assuming qz is negligible 
    117   /// Angles here are respect to detector coordinate 
    118   ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854)) 
    119     // b3 axis orientation 
    120     b3_x = cos(theta) * cos(phi); 
    121     b3_y = sin(theta); 
    122     //b3_z = -cos(theta) * sin(phi); 
    123     cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z; 
    124  
    125     //alpha = acos(cos_val_b3); 
    126     // b1 axis orientation 
    127     b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi); 
    128     b1_y = sin(psi)*cos(theta); 
    129     cos_val_b1 = b1_x*q_x + b1_y*q_y; 
    130     // b2 axis orientation 
    131     b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi); 
    132         b2_y = cos(theta)*cos(psi); 
    133     cos_val_b2 = b2_x*q_x + b2_y*q_y; 
    134  
    135     // The following test should always pass 
    136     if (fabs(cos_val_b3)>1.0) { 
    137       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    138       cos_val_b3 = 1.0; 
    139     } 
    140     if (fabs(cos_val_b2)>1.0) { 
    141       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    142       cos_val_b2 = 1.0; 
    143     } 
    144     if (fabs(cos_val_b1)>1.0) { 
    145       //printf("FCC_ana_2D: Unexpected error: cos()>1\n"); 
    146       cos_val_b1 = 1.0; 
    147     } 
    148     // Compute the angle btw vector q and the a3 axis 
    149     a3_dot_q = 0.5*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3); 
    150  
    151     // a1 axis 
    152     a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1); 
    153  
    154     // a2 axis 
    155     a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2); 
    156  
    157  
    158     // Get Fkq and Fkq_2 
    159     Fkq = exp(-0.5*pow(Da/dnn,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q)); 
    160     Fkq_2 = Fkq*Fkq; 
    161     // Call Zq=Z1*Z2*Z3 
    162     Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2); 
    163     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2); 
    164     Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2); 
    165  
    166   // Use SphereForm directly from libigor 
    167   answer = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale; 
    168  
    169   return answer; 
    170  } 
     107    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
     108    //the occupied volume of the lattice 
     109    const double lattice_scale = 2.0*sphere_volume(M_SQRT1_2*radius/dnn); 
     110    return lattice_scale * Fq; 
     111} 

  • sasmodels/models/sc_paracrystal.c

    r0bef47b r11ca2ab  
    6060} 
    6161 
    62 static 
    63 double sc_crystal_kernel(double q, 
     62double Iq(double q, 
    6463          double dnn, 
    6564          double d_factor, 
    … …  
    103102} 
    104103 
    105 static 
    106 double sc_crystal_kernel_2d(double q, double q_x, double q_y, 
     104double Iqxy(double qx, double qy, 
    107105          double dnn, 
    108106          double d_factor, 
    … …  
    114112          double psi) 
    115113{ 
    116     //convert angle degree to radian 
    117     theta = theta * M_PI_180; 
    118     phi = phi * M_PI_180; 
    119     psi = psi * M_PI_180; 
     114    double q, cos_a1, cos_a2, cos_a3; 
     115    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 
    120116 
    121     const double qda_2 = pow(q*d_factor*dnn,2.0); 
     117    const double qd = q*dnn; 
     118    const double exp_qd = exp(0.5*square(qd*d_factor)); 
     119    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd; 
     120    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd; 
    122121 
    123     double snt, cnt; 
    124     SINCOS(theta, snt, cnt); 
     122    const double Zq = sinh_qd/(cosh_qd - cos(qd*cos_a1)) 
     123                    * sinh_qd/(cosh_qd - cos(qd*cos_a2)) 
     124                    * sinh_qd/(cosh_qd - cos(qd*cos_a3)); 
    125125 
    126     double snp, cnp; 
    127     SINCOS(phi, snp, cnp); 
    128  
    129     double sns, cns; 
    130     SINCOS(psi, sns, cns); 
    131  
    132     /// Angles here are respect to detector coordinate instead of against 
    133     //  q coordinate(PRB 36, 3, 1754) 
    134     // a3 axis orientation 
    135  
    136     const double a3_x = cnt * cnp; 
    137     const double a3_y = snt; 
    138  
    139     // Compute the angle btw vector q and the a3 axis 
    140     double cos_val_a3 = a3_x*q_x + a3_y*q_y; 
    141  
    142     // a1 axis orientation 
    143     const double a1_x = -cnp*sns * snt+snp*cns; 
    144     const double a1_y = sns*cnt; 
    145  
    146     double cos_val_a1 = a1_x*q_x + a1_y*q_y; 
    147  
    148     // a2 axis orientation 
    149     const double a2_x = -snt*cns*cnp-sns*snp; 
    150     const double a2_y = cnt*cns; 
    151  
    152     // a2 axis 
    153     const double cos_val_a2 =  a2_x*q_x + a2_y*q_y; 
    154  
    155     // The following test should always pass 
    156     if (fabs(cos_val_a3)>1.0) { 
    157         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    158         cos_val_a3 = 1.0; 
    159     } 
    160     if (fabs(cos_val_a1)>1.0) { 
    161         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    162         cos_val_a1 = 1.0; 
    163     } 
    164     if (fabs(cos_val_a2)>1.0) { 
    165         //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n"); 
    166         cos_val_a3 = 1.0; 
    167     } 
    168  
    169     const double a3_dot_q = dnn*q*cos_val_a3; 
    170     const double a1_dot_q = dnn*q*cos_val_a1; 
    171     const double a2_dot_q = dnn*q*cos_val_a2; 
    172  
    173     // Call Zq=Z1*Z2*Z3 
    174     double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_dot_q)+exp(-qda_2)); 
    175     Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a2_dot_q)+exp(-qda_2)); 
    176     Zq *= (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a3_dot_q)+exp(-qda_2)); 
    177  
    178     // Use SphereForm directly from libigor 
    179     double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    180  
    181     //consider scales 
    182     const double latticeScale = sphere_volume(radius/dnn); 
    183     answer *= latticeScale; 
    184  
    185     return answer; 
     126    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
     127    //the occupied volume of the lattice 
     128    const double lattice_scale = sphere_volume(radius/dnn); 
     129    return lattice_scale * Fq; 
    186130} 
    187  
    188 double Iq(double q, 
    189           double dnn, 
    190           double d_factor, 
    191           double radius, 
    192           double sphere_sld, 
    193           double solvent_sld) 
    194 { 
    195     return sc_crystal_kernel(q, 
    196               dnn, 
    197               d_factor, 
    198               radius, 
    199               sphere_sld, 
    200               solvent_sld); 
    201 } 
    202  
    203 // Iqxy is never called since no orientation or magnetic parameters. 
    204 double Iqxy(double qx, double qy, 
    205             double dnn, 
    206             double d_factor, 
    207             double radius, 
    208             double sphere_sld, 
    209             double solvent_sld, 
    210             double theta, 
    211             double phi, 
    212             double psi) 
    213 { 
    214     double q = sqrt(qx*qx + qy*qy); 
    215  
    216  
    217     return sc_crystal_kernel_2d(q, qx/q, qy/q, 
    218                   dnn, 
    219                   d_factor, 
    220                   radius, 
    221                   sphere_sld, 
    222                   solvent_sld, 
    223                   theta, 
    224                   phi, 
    225                   psi); 
    226  
    227 } 
    228  

  • sasmodels/models/triaxial_ellipsoid.c

    ra807206 r11ca2ab  
    5858    double psi) 
    5959{ 
    60     double stheta, ctheta; 
    61     double sphi, cphi; 
    62     double spsi, cpsi; 
     60    double q, calpha, cmu, cnu; 
     61    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 
    6362 
    64     const double q = sqrt(qx*qx + qy*qy); 
    65     const double qxhat = qx/q; 
    66     const double qyhat = qy/q; 
    67     SINCOS(theta*M_PI_180, stheta, ctheta); 
    68     SINCOS(phi*M_PI_180, sphi, cphi); 
    69     SINCOS(psi*M_PI_180, spsi, cpsi); 
    70     const double calpha = ctheta*cphi*qxhat + stheta*qyhat; 
    71     const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 
    72     const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 
    7363    const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 
    7464                          + radius_equat_major*radius_equat_major*cmu*cmu 
Note: See TracChangeset for help on using the changeset viewer.