Changeset 3ac4e1b in sasmodels


Ignore:
Timestamp:
Oct 18, 2016 8:32:55 AM (7 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:
0a3d9b2
Parents:
158cee4
git-author:
Paul Kienzle <pkienzle@…> (10/18/16 08:32:20)
git-committer:
Paul Kienzle <pkienzle@…> (10/18/16 08:32:55)
Message:

stacked_disks: code tidying

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/stacked_disks.c

    r6831fa0 r3ac4e1b  
    2727 
    2828static 
    29 double _kernel(double qq, 
     29double _kernel(double q, 
    3030               double radius, 
    3131               double core_sld, 
     
    4141 
    4242{ 
    43         // qq is the q-value for the calculation (1/A) 
    44         // radius is the core radius of the cylinder (A) 
    45         // *_sld are the respective SLD's 
    46         // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A) 
    47         // zi is the dummy variable for the integration (x in Feigin's notation) 
     43    // q is the q-value for the calculation (1/A) 
     44    // radius is the core radius of the cylinder (A) 
     45    // *_sld are the respective SLD's 
     46    // halfheight is the *Half* CORE-LENGTH of the cylinder = L (A) 
     47    // zi is the dummy variable for the integration (x in Feigin's notation) 
    4848 
    49         const double besarg1 = qq*radius*sin_alpha; 
    50         //const double besarg2 = qq*radius*sin_alpha; 
     49    const double besarg1 = q*radius*sin_alpha; 
     50    //const double besarg2 = q*radius*sin_alpha; 
    5151 
    52         const double sinarg1 = qq*halfheight*cos_alpha; 
    53         const double sinarg2 = qq*(halfheight+thick_layer)*cos_alpha; 
     52    const double sinarg1 = q*halfheight*cos_alpha; 
     53    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
    5454 
    55         const double be1 = sas_J1c(besarg1); 
    56         //const double be2 = sas_J1c(besarg2); 
    57         const double be2 = be1; 
    58         const double si1 = sinc(sinarg1); 
    59         const double si2 = sinc(sinarg2); 
     55    const double be1 = sas_J1c(besarg1); 
     56    //const double be2 = sas_J1c(besarg2); 
     57    const double be2 = be1; 
     58    const double si1 = sinc(sinarg1); 
     59    const double si2 = sinc(sinarg2); 
    6060 
    61         const double dr1 = (core_sld-solvent_sld); 
    62         const double dr2 = (layer_sld-solvent_sld); 
    63         const double area = M_PI*radius*radius; 
    64         const double totald = 2.0*(thick_layer+halfheight); 
     61    const double dr1 = core_sld - solvent_sld; 
     62    const double dr2 = layer_sld - solvent_sld; 
     63    const double area = M_PI*radius*radius; 
     64    const double totald = 2.0*(thick_layer + halfheight); 
    6565 
    66         const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1); 
    67         const double t2 = area*dr2*(totald*si2-2.0*halfheight*si1)*(be2); 
     66    const double t1 = area * (2.0*halfheight) * dr1 * si1 * be1; 
     67    const double t2 = area * dr2 * (totald*si2 - 2.0*halfheight*si1) * be2; 
    6868 
     69    double pq = square(t1 + t2); 
    6970 
    70         double retval =((t1+t2)*(t1+t2)); 
     71    // loop for the structure factor S(q) 
     72    double qd_cos_alpha = q*d*cos_alpha; 
     73    double debye_arg = -0.5*square(qd_cos_alpha*sigma_dnn); 
     74    double sq=0.0; 
     75    for (int kk=1; kk<n_stacking; kk++) { 
     76        sq += (n_stacking-kk) * cos(qd_cos_alpha*kk) * exp(debye_arg*kk); 
     77    } 
     78    // end of loop for S(q) 
     79    sq = 1.0 + 2.0*sq/n_stacking; 
    7180 
    72         // loop for the structure facture S(q) 
    73         double sqq=0.0; 
    74         for(int kk=1;kk<n_stacking;kk+=1) { 
    75                 double qd_cos_alpha = qq*d*cos_alpha; 
    76                 double dexpt=square(qd_cos_alpha*sigma_dnn)*kk/2.0; 
    77                 sqq += (n_stacking-kk)*cos(qd_cos_alpha*kk)*exp(-1.*dexpt); 
    78         } 
    79         // end of loop for S(q) 
    80         sqq=1.0+2.0*sqq/n_stacking; 
    81  
    82         return retval * sqq; 
     81    return pq * sq; 
    8382} 
    8483 
     
    9594                            double solvent_sld) 
    9695{ 
    97 /*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks 
     96/*    StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks 
    9897like clay platelets that are not exfoliated 
    9998*/ 
    100         double summ = 0.0;      //initialize integral 
     99    double summ = 0.0;    //initialize integral 
    101100 
    102         double d = 2.0*thick_layer+thick_core; 
    103         double halfheight = 0.5*thick_core; 
     101    double d = 2.0*thick_layer+thick_core; 
     102    double halfheight = 0.5*thick_core; 
    104103 
    105         for(int i=0; i<N_POINTS_76; i++) { 
    106                 double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
    107                 double sin_alpha, cos_alpha; // slots to hold sincos function output 
    108                 SINCOS(zi, sin_alpha, cos_alpha); 
    109                 double yyy = _kernel(q, 
    110                                    radius, 
    111                                    core_sld, 
    112                                    layer_sld, 
    113                                    solvent_sld, 
    114                                    halfheight, 
    115                                    thick_layer, 
    116                                    sin_alpha, 
    117                                    cos_alpha, 
    118                                    sigma_dnn, 
    119                                    d, 
    120                                    n_stacking); 
    121                 summ += Gauss76Wt[i] * yyy * sin_alpha; 
    122         } 
     104    for(int i=0; i<N_POINTS_76; i++) { 
     105        double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     106        double sin_alpha, cos_alpha; // slots to hold sincos function output 
     107        SINCOS(zi, sin_alpha, cos_alpha); 
     108        double yyy = _kernel(q, 
     109                           radius, 
     110                           core_sld, 
     111                           layer_sld, 
     112                           solvent_sld, 
     113                           halfheight, 
     114                           thick_layer, 
     115                           sin_alpha, 
     116                           cos_alpha, 
     117                           sigma_dnn, 
     118                           d, 
     119                           n_stacking); 
     120        summ += Gauss76Wt[i] * yyy * sin_alpha; 
     121    } 
    123122 
    124         double answer = M_PI_4*summ; 
     123    double answer = M_PI_4*summ; 
    125124 
    126         //Convert to [cm-1] 
    127         answer *= 1.0e-4; 
    128  
    129         return answer; 
     125    //Convert to [cm-1] 
     126    return 1.0e-4*answer; 
    130127} 
    131128 
Note: See TracChangeset for help on using the changeset viewer.