Changeset 5bddd89 in sasmodels


Ignore:
Timestamp:
Oct 14, 2016 4:05:39 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:
9068f4c
Parents:
0b717c5
Message:

use ORIENT macro for remaining symmetric models

Location:
sasmodels/models
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/capped_cylinder.c

    r0d6e865 r5bddd89  
    4949} 
    5050 
     51static double 
     52_fq(double q, double h, double radius_cap, double radius, double half_length, 
     53    double sin_alpha, double cos_alpha) 
     54{ 
     55    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
     56    const double bj = sas_J1c(q*radius*sin_alpha); 
     57    const double si = sinc(q*half_length*cos_alpha); 
     58    const double cyl_Fq = 2.*M_PI*radius*radius*half_length*bj*si; 
     59    const double Aq = cap_Fq + cyl_Fq; 
     60    return Aq; 
     61} 
     62 
    5163double form_volume(double radius, double radius_cap, double length) 
    5264{ 
    … …  
    93105        SINCOS(alpha, sin_alpha, cos_alpha); 
    94106 
    95         const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    96         const double bj = sas_J1c(q*radius*sin_alpha); 
    97         const double si = sinc(q*half_length*cos_alpha); 
    98         const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    99         const double Aq = cap_Fq + cyl_Fq; 
    100         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; // sin_alpha for spherical coord integration 
     107        const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
     108        // sin_alpha for spherical coord integration 
     109        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    101110    } 
    102111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    114123    double theta, double phi) 
    115124{ 
    116     // Compute angle alpha between q and the cylinder axis 
    117     double sn, cn; 
    118     SINCOS(phi*M_PI_180, sn, cn); 
    119     const double q = sqrt(qx*qx+qy*qy); 
    120     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    121     const double alpha = acos(cos_val); // rod angle relative to q 
     125    double q, sin_alpha, cos_alpha; 
     126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    122127 
    123128    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    124     const double half_length = 0.5*length; 
    125  
    126     double sin_alpha, cos_alpha; // slots to hold sincos function output 
    127     SINCOS(alpha, sin_alpha, cos_alpha); 
    128     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    129     const double bj = sas_J1c(q*radius*sin_alpha); 
    130     const double si = sinc(q*half_length*cos_alpha); 
    131     const double cyl_Fq = M_PI*radius*radius*length*bj*si; 
    132     const double Aq = cap_Fq + cyl_Fq; 
     129    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
    133130 
    134131    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    r0d6e865 r5bddd89  
    2626double form_volume(double radius, double thick_rim, double thick_face, double length) 
    2727{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2*thick_face); 
     28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
    2929} 
    3030 
    … …  
    3939              double rhor, 
    4040              double rhosolv, 
    41               double dum) 
     41              double sin_alpha, 
     42              double cos_alpha) 
    4243{ 
    4344    double si1,si2,be1,be2; 
    … …  
    4950    const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 
    5051    const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 
    51     double sn,cn; 
    52     SINCOS(dum, sn, cn); 
    53     double besarg1 = qq*rad*sn; 
    54     double besarg2 = qq*(rad+radthick)*sn; 
    55     double sinarg1 = qq*length*cn; 
    56     double sinarg2 = qq*(length+facthick)*cn; 
     52    double besarg1 = qq*rad*sin_alpha; 
     53    double besarg2 = qq*(rad+radthick)*sin_alpha; 
     54    double sinarg1 = qq*length*cos_alpha; 
     55    double sinarg2 = qq*(length+facthick)*cos_alpha; 
    5756 
    5857    be1 = sas_J1c(besarg1); 
    … …  
    6564                     vol3*dr3*si2*be1; 
    6665 
    67     const double retval = t*t*sn; 
     66    const double retval = t*t*sin_alpha; 
    6867 
    69     return(retval); 
     68    return retval; 
    7069 
    7170} 
    … …  
    8382{ 
    8483    // set up the integration end points 
    85     const double uplim = M_PI/4; 
    86     const double halfheight = length/2.0; 
     84    const double uplim = M_PI_4; 
     85    const double halfheight = 0.5*length; 
    8786 
    8887    double summ = 0.0; 
    8988    for(int i=0;i<N_POINTS_76;i++) { 
    90         double zi = (Gauss76Z[i] + 1.0)*uplim; 
     89        double alpha = (Gauss76Z[i] + 1.0)*uplim; 
     90        double sin_alpha, cos_alpha; // slots to hold sincos function output 
     91        SINCOS(alpha, sin_alpha, cos_alpha); 
    9192        double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 
    92                              halfheight, rhoc, rhoh, rhor,rhosolv, zi); 
     93                             halfheight, rhoc, rhoh, rhor, rhosolv, 
     94                             sin_alpha, cos_alpha); 
    9395        summ += yyy; 
    9496    } 
    … …  
    9698    // calculate value of integral to return 
    9799    double answer = uplim*summ; 
    98     return(answer); 
     100    return answer; 
    99101} 
    100102 
    101103static double 
    102 bicelle_kernel_2d(double q, double q_x, double q_y, 
     104bicelle_kernel_2d(double qx, double qy, 
    103105          double radius, 
    104106          double thick_rim, 
    … …  
    112114          double phi) 
    113115{ 
    114     //convert angle degree to radian 
    115     theta *= M_PI_180; 
    116     phi *= M_PI_180; 
     116    double q, sin_alpha, cos_alpha; 
     117    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    117118 
    118     // Cylinder orientation 
    119     const double cyl_x = sin(theta) * cos(phi); 
    120     const double cyl_y = sin(theta) * sin(phi); 
    121  
    122     // Compute the angle btw vector q and the axis of the cylinder 
    123     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    124     const double alpha = acos( cos_val ); 
    125  
    126     // Get the kernel 
    127119    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
    128                            length/2.0, core_sld, face_sld, rim_sld, 
    129                            solvent_sld, alpha) / fabs(sin(alpha)); 
     120                           0.5*length, core_sld, face_sld, rim_sld, 
     121                           solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 
    130122 
    131123    answer *= 1.0e-4; 
    … …  
    162154          double phi) 
    163155{ 
    164     double q; 
    165     q = sqrt(qx*qx+qy*qy); 
    166     double intensity = bicelle_kernel_2d(q, qx/q, qy/q, 
     156    double intensity = bicelle_kernel_2d(qx, qy, 
    167157                      radius, 
    168158                      thick_rim, 

  • sasmodels/models/core_shell_cylinder.c

    r0d6e865 r5bddd89  
    1111double _cyl(double twovd, double besarg, double siarg) 
    1212{ 
    13     const double bj = (besarg == 0.0 ? 0.5 : 0.5*sas_J1c(besarg)); 
    14     const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg); 
    15     return twovd*si*bj; 
     13    return twovd * sinc(siarg) * sas_J1c(besarg); 
    1614} 
    1715 
    1816double form_volume(double radius, double thickness, double length) 
    1917{ 
    20     return M_PI*(radius+thickness)*(radius+thickness)*(length+2*thickness); 
     18    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    2119} 
    2220 
    … …  
    6664    double phi) 
    6765{ 
    68     double sn, cn; // slots to hold sincos function output 
    69  
    70     // Compute angle alpha between q and the cylinder axis 
    71     SINCOS(phi*M_PI_180, sn, cn); 
    72     // # The following correction factor exists in sasview, but it can't be 
    73     // # right, so we are leaving it out for now. 
    74     // const double correction = fabs(cn)*M_PI_2; 
    75     const double q = sqrt(qx*qx+qy*qy); 
    76     const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    77     const double alpha = acos(cos_val); 
     66    double q, sin_alpha, cos_alpha; 
     67    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    7868 
    7969    const double core_qr = q*radius; 
    … …  
    8676                             * (shell_sld-solvent_sld); 
    8777 
    88     SINCOS(alpha, sn, cn); 
    89     const double fq = _cyl(core_twovd, core_qr*sn, core_qh*cn) 
    90         + _cyl(shell_twovd, shell_qr*sn, shell_qh*cn); 
     78    const double fq = _cyl(core_twovd, core_qr*sin_alpha, core_qh*cos_alpha) 
     79        + _cyl(shell_twovd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
    9180    return 1.0e-4 * fq * fq; 
    9281} 

  • sasmodels/models/core_shell_ellipsoid.c

    r0d6e865 r5bddd89  
    3232    const double equat_shell = radius_equat_core + thick_shell; 
    3333    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    34     double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell; 
     34    double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 
    3535    return vol; 
    3636} 
    … …  
    6060 
    6161    for(int i=0;i<N_POINTS_76;i++) { 
    62         double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     62        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    6363        double yyy = Gauss76Wt[i] * gfn4(zi, 
    6464                                  radius_equat_core, 
    … …  
    7272    } 
    7373 
    74     double answer = (uplim-lolim)/2.0*summ; 
     74    double answer = 0.5*(uplim-lolim)*summ; 
    7575    //convert to [cm-1] 
    7676    answer *= 1.0e-4; 
    … …  
    8080 
    8181static double 
    82 core_shell_ellipsoid_xt_kernel_2d(double q, double q_x, double q_y, 
     82core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    8383          double radius_equat_core, 
    8484          double x_core, 
    … …  
    9191          double phi) 
    9292{ 
    93     //convert angle degree to radian 
    94     theta = theta * M_PI_180; 
    95     phi = phi * M_PI_180; 
    96  
    97     // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis. 
    98     const double cyl_x = sin(theta) * cos(phi); 
    99     const double cyl_y = sin(theta) * sin(phi); 
     93    double q, sin_alpha, cos_alpha; 
     94    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    10095 
    10196    const double sldcs = core_sld - shell_sld; 
    10297    const double sldss = shell_sld- solvent_sld; 
    103  
    104     // Compute the angle btw vector q and the 
    105     // axis of the cylinder 
    106     const double cos_val = cyl_x*q_x + cyl_y*q_y; 
    10798 
    10899    const double polar_core = radius_equat_core*x_core; 
    … …  
    112103    // Call the IGOR library function to get the kernel: 
    113104    // MUST use gfn4 not gf2 because of the def of params. 
    114     double answer = gfn4(cos_val, 
     105    double answer = gfn4(cos_alpha, 
    115106                  radius_equat_core, 
    116107                  polar_core, 
    … …  
    160151          double phi) 
    161152{ 
    162     double q; 
    163     q = sqrt(qx*qx+qy*qy); 
    164     double intensity = core_shell_ellipsoid_xt_kernel_2d(q, qx/q, qy/q, 
     153    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    165154                       radius_equat_core, 
    166155                       x_core, 

  • sasmodels/models/ellipsoid.c

    r0d6e865 r5bddd89  
    4949    double phi) 
    5050{ 
    51     double sn, cn; 
    52  
    53     const double q = sqrt(qx*qx + qy*qy); 
    54     SINCOS(phi*M_PI_180, sn, cn); 
    55     const double cos_alpha = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q); 
    56     const double alpha = acos(cos_alpha); 
    57     SINCOS(alpha, sn, cn); 
    58     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sn); 
     51    double q, sin_alpha, cos_alpha; 
     52    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     53    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
    5954    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6055 

  • sasmodels/models/hollow_cylinder.c

    r0d6e865 r5bddd89  
    11double form_volume(double radius, double thickness, double length); 
    2  
    32double Iq(double q, double radius, double thickness, double length, double sld, 
    43        double solvent_sld); 
    … …  
    98 
    109// From Igor library 
    11 static double hollow_cylinder_scaling( 
    12     double integrand, double delrho, double volume) 
     10static double 
     11_hollow_cylinder_scaling(double integrand, double delrho, double volume) 
    1312{ 
    14     double answer; 
    15     // Multiply by contrast^2 
    16     answer = integrand*delrho*delrho; 
    17  
    18     //normalize by cylinder volume 
    19     answer *= volume*volume; 
    20  
    21     //convert to [cm-1] 
    22     answer *= 1.0e-4; 
    23  
    24     return answer; 
     13    return 1.0e-4 * square(volume * delrho * integrand); 
    2514} 
    2615 
    2716 
    28 static double _hollow_cylinder_kernel( 
    29     double q, double radius, double thickness, double length, double dum) 
     17static double 
     18_hollow_cylinder_kernel(double q, 
     19    double radius, double thickness, double length, double sin_val, double cos_val) 
    3020{ 
    31     const double qs = q*sqrt(1.0-dum*dum); 
     21    const double qs = q*sin_val; 
    3222    const double lam1 = sas_J1c((radius+thickness)*qs); 
    3323    const double lam2 = sas_J1c(radius*qs); 
    … …  
    3525    //Note: lim_{r -> r_c} psi = J0(radius_core*qs) 
    3626    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    37     const double t2 = sinc(q*length*dum/2.0); 
    38     return square(psi*t2); 
     27    const double t2 = sinc(0.5*q*length*cos_val); 
     28    return psi*t2; 
    3929} 
    4030 
    41  
    42 static double hollow_cylinder_analytical_2D_scaled( 
    43     double q, double q_x, double q_y, double radius, double thickness, 
    44     double length, double sld, double solvent_sld, double theta, double phi) 
    45 { 
    46     double cyl_x, cyl_y; //, cyl_z 
    47     //double q_z; 
    48     double vol, cos_val, delrho; 
    49     double answer; 
    50     //convert angle degree to radian 
    51     theta = theta * M_PI_180; 
    52     phi = phi * M_PI_180; 
    53     delrho = solvent_sld - sld; 
    54  
    55     // Cylinder orientation 
    56     cyl_x = sin(theta) * cos(phi); 
    57     cyl_y = sin(theta) * sin(phi); 
    58     //cyl_z = -cos(theta) * sin(phi); 
    59  
    60     // q vector 
    61     //q_z = 0; 
    62  
    63     // Compute the angle btw vector q and the 
    64     // axis of the cylinder 
    65     cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
    66  
    67     answer = _hollow_cylinder_kernel(q, radius, thickness, length, cos_val); 
    68  
    69     vol = form_volume(radius, thickness, length); 
    70     answer = hollow_cylinder_scaling(answer, delrho, vol); 
    71  
    72     return answer; 
    73 } 
    74  
    75  
    76 double form_volume(double radius, double thickness, double length) 
     31double 
     32form_volume(double radius, double thickness, double length) 
    7733{ 
    7834    double v_shell = M_PI*length*((radius+thickness)*(radius+thickness)-radius*radius); 
    … …  
    8137 
    8238 
    83 double Iq(double q, double radius, double thickness, double length, 
     39double 
     40Iq(double q, double radius, double thickness, double length, 
    8441    double sld, double solvent_sld) 
    8542{ 
    86     int i; 
    87     double lower,upper,zi, inter;               //upper and lower integration limits 
    88     double summ,answer,delrho;                  //running tally of integration 
    89     double norm,volume; //final calculation variables 
     43    const double lower = 0.0; 
     44    const double upper = 1.0;           //limits of numerical integral 
    9045 
    91     lower = 0.0; 
    92     upper = 1.0;                //limits of numerical integral 
    93  
    94     summ = 0.0;                 //initialize intergral 
    95     for (i=0;i<76;i++) { 
    96         zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 
    97         inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, radius, thickness, length, zi); 
    98         summ += inter; 
     46    double summ = 0.0;                  //initialize intergral 
     47    for (int i=0;i<76;i++) { 
     48        const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     49        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
     50        const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
     51                                                     sin_val, cos_val); 
     52        summ += Gauss76Wt[i] * inter; 
    9953    } 
    10054 
    101     norm = summ*(upper-lower)/2.0; 
    102     volume = form_volume(radius, thickness, length); 
    103     delrho = solvent_sld - sld; 
    104     answer = hollow_cylinder_scaling(norm, delrho, volume); 
     55    const double Aq = 0.5*summ*(upper-lower); 
     56    const double volume = form_volume(radius, thickness, length); 
     57    return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 
     58} 
    10559 
    106     return(answer); 
     60double 
     61Iqxy(double qx, double qy, 
     62    double radius, double thickness, double length, 
     63    double sld, double solvent_sld, double theta, double phi) 
     64{ 
     65    double q, sin_alpha, cos_alpha; 
     66    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     67    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
     68        sin_alpha, cos_alpha); 
     69 
     70    const double vol = form_volume(radius, thickness, length); 
     71    return _hollow_cylinder_scaling(Aq, solvent_sld-sld, vol); 
    10772} 
    10873 
    10974 
    110 double Iqxy(double qx, double qy, double radius, double thickness, 
    111     double length, double sld, double solvent_sld, double theta, double phi) 
    112 { 
    113     const double q = sqrt(qx*qx+qy*qy); 
    114     return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, thickness, length, sld, solvent_sld, theta, phi); 
    115 } 
Note: See TracChangeset for help on using the changeset viewer.