Changeset ec87470 in sasmodels for sasmodels/models


Ignore:
Timestamp:
Aug 3, 2016 10:26:20 AM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c, 3d9001f
Parents:
edf06e1 (diff), d119f34 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'origin/master' into polydisp

Location:
sasmodels/models
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_ellipsoid.py

    rec45c4f r7f1ee79  
    9999category = "shape:ellipsoid" 
    100100 
    101 single = False  # TODO: maybe using sph_j1c inside gfn would help? 
    102101# pylint: disable=bad-whitespace, line-too-long 
    103102#             ["name", "units", default, [lower, upper], "type", "description"], 

  • sasmodels/models/gaussian_peak.py

    r2c74c11 r7f1ee79  
    3838category = "shape-independent" 
    3939 
    40 single = False 
    4140#             ["name", "units", default, [lower, upper], "type","description"], 
    4241parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"], 

  • sasmodels/models/hardsphere.py

    r9eb3632 r7f1ee79  
    5858category = "structure-factor" 
    5959structure_factor = True 
    60 single = False 
     60single = False # TODO: check 
    6161 
    6262#             ["name", "units", default, [lower, upper], "type","description"], 

  • sasmodels/models/hollow_cylinder.c

    r2f5c6d4 rf95556f  
    66        double solvent_sld, double theta, double phi); 
    77 
    8 #define INVALID(v) (v.core_radius >= v.radius || v.radius >= v.length) 
     8#define INVALID(v) (v.core_radius >= v.radius) 
    99 
    1010// From Igor library 
    11 static double hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     11static double hollow_cylinder_scaling( 
     12    double integrand, double delrho, double volume) 
    1213{ 
    13         double answer; 
    14         // Multiply by contrast^2 
    15         answer = integrand*delrho*delrho; 
     14    double answer; 
     15    // Multiply by contrast^2 
     16    answer = integrand*delrho*delrho; 
    1617 
    17         //normalize by cylinder volume 
    18         answer *= volume*volume; 
     18    //normalize by cylinder volume 
     19    answer *= volume*volume; 
    1920 
    20         //convert to [cm-1] 
    21         answer *= 1.0e-4; 
     21    //convert to [cm-1] 
     22    answer *= 1.0e-4; 
    2223 
    23         return answer; 
     24    return answer; 
    2425} 
    2526 
    26 static double _hollow_cylinder_kernel(double q, double core_radius, double radius, 
    27         double length, double dum) 
     27 
     28static double _hollow_cylinder_kernel( 
     29    double q, double core_radius, double radius, double length, double dum) 
    2830{ 
    29     double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables 
    30      
    31     gamma = core_radius/radius; 
    32     arg1 = q*radius*sqrt(1.0-dum*dum);          //1=shell (outer radius) 
    33     arg2 = q*core_radius*sqrt(1.0-dum*dum);                     //2=core (inner radius) 
    34     if (arg1 == 0.0){ 
    35         lam1 = 1.0; 
    36     }else{ 
    37         lam1 = sas_J1c(arg1); 
    38     } 
    39     if (arg2 == 0.0){ 
    40         lam2 = 1.0; 
    41     }else{ 
    42         lam2 = sas_J1c(arg2); 
    43     } 
    44     //Todo: Need to check psi behavior as gamma goes to 1. 
    45     psi = (lam1 -  gamma*gamma*lam2)/(1.0-gamma*gamma);         //SRK 10/19/00 
    46     sinarg = q*length*dum/2.0; 
    47     if (sinarg == 0.0){ 
    48         t2 = 1.0; 
    49     }else{ 
    50         t2 = sin(sinarg)/sinarg; 
    51     } 
     31    const double qs = q*sqrt(1.0-dum*dum); 
     32    const double lam1 = sas_J1c(radius*qs); 
     33    const double lam2 = sas_J1c(core_radius*qs); 
     34    const double gamma_sq = square(core_radius/radius); 
     35    //Note: lim_{r -> r_c} psi = J0(core_radius*qs) 
     36    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); 
     39} 
    5240 
    53     retval = psi*psi*t2*t2; 
    54      
    55     return(retval); 
    56 } 
    57 static double hollow_cylinder_analytical_2D_scaled(double q, double q_x, double q_y, double radius, double core_radius, double length, double sld, 
    58         double solvent_sld, double theta, double phi) { 
    59         double cyl_x, cyl_y; //, cyl_z 
    60         //double q_z; 
    61         double vol, cos_val, delrho; 
    62         double answer; 
    63         //convert angle degree to radian 
    64         double pi = 4.0*atan(1.0); 
    65         theta = theta * pi/180.0; 
    66         phi = phi * pi/180.0; 
    67         delrho = solvent_sld - sld; 
    6841 
    69         // Cylinder orientation 
    70         cyl_x = cos(theta) * cos(phi); 
    71         cyl_y = sin(theta); 
    72         //cyl_z = -cos(theta) * sin(phi); 
     42static double hollow_cylinder_analytical_2D_scaled( 
     43    double q, double q_x, double q_y, double radius, double core_radius, 
     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; 
    7354 
    74         // q vector 
    75         //q_z = 0; 
     55    // Cylinder orientation 
     56    cyl_x = cos(theta) * cos(phi); 
     57    cyl_y = sin(theta); 
     58    //cyl_z = -cos(theta) * sin(phi); 
    7659 
    77         // Compute the angle btw vector q and the 
    78         // axis of the cylinder 
    79         cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
     60    // q vector 
     61    //q_z = 0; 
    8062 
    81         answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val); 
     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; 
    8266 
    83         vol = form_volume(radius, core_radius, length); 
    84         answer = hollow_cylinder_scaling(answer, delrho, vol); 
     67    answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val); 
    8568 
    86         return answer; 
     69    vol = form_volume(radius, core_radius, length); 
     70    answer = hollow_cylinder_scaling(answer, delrho, vol); 
     71 
     72    return answer; 
    8773} 
    8874 
    … …  
    9076double form_volume(double radius, double core_radius, double length) 
    9177{ 
    92         double pi = 4.0*atan(1.0); 
    93         double v_shell = pi*length*(radius*radius-core_radius*core_radius); 
    94         return(v_shell); 
     78    double v_shell = M_PI*length*(radius*radius-core_radius*core_radius); 
     79    return(v_shell); 
    9580} 
    9681 
    9782 
    98 double Iq(double q, double radius, double core_radius, double length, double sld, 
    99         double solvent_sld) 
     83double Iq(double q, double radius, double core_radius, double length, 
     84    double sld, double solvent_sld) 
    10085{ 
    10186    int i; 
    102         int nord=76;                    //order of integration 
    103         double lower,upper,zi, inter;           //upper and lower integration limits 
    104         double summ,answer,delrho;                      //running tally of integration 
    105         double norm,volume;     //final calculation variables 
    106          
    107         delrho = solvent_sld - sld; 
    108         lower = 0.0; 
    109         upper = 1.0;            //limits of numerical integral 
     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 
    11090 
    111         summ = 0.0;                     //initialize intergral 
    112         for(i=0;i<nord;i++) { 
    113                 zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 
    114                 inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi); 
    115                 summ += inter; 
    116         } 
    117          
    118         norm = summ*(upper-lower)/2.0; 
    119         volume = form_volume(radius, core_radius, length); 
    120         answer = hollow_cylinder_scaling(norm, delrho, volume); 
    121          
    122         return(answer); 
     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, core_radius, radius, length, zi); 
     98        summ += inter; 
     99    } 
     100 
     101    norm = summ*(upper-lower)/2.0; 
     102    volume = form_volume(radius, core_radius, length); 
     103    delrho = solvent_sld - sld; 
     104    answer = hollow_cylinder_scaling(norm, delrho, volume); 
     105 
     106    return(answer); 
    123107} 
    124108 
    125 double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld, 
    126         double solvent_sld, double theta, double phi) 
     109 
     110double Iqxy(double qx, double qy, double radius, double core_radius, 
     111    double length, double sld, double solvent_sld, double theta, double phi) 
    127112{ 
    128         double q; 
    129         q = sqrt(qx*qx+qy*qy); 
    130         return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi); 
     113    const double q = sqrt(qx*qx+qy*qy); 
     114    return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi); 
    131115} 

  • sasmodels/models/hollow_cylinder.py

    r42356c8 r58210db  
    7979# pylint: enable=bad-whitespace, line-too-long 
    8080 
    81 source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
     81source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    8282 
    8383# pylint: disable=W0613 

  • sasmodels/models/lamellar_hg_stack_caille.py

    r42356c8 r7f1ee79  
    9090category = "shape:lamellae" 
    9191 
    92 single = False 
     92single = False  # TODO: check 
    9393parameters = [ 
    9494    #   [ "name", "units", default, [lower, upper], "type", 

  • sasmodels/models/lamellar_stack_caille.py

    r42356c8 r7f1ee79  
    8383category = "shape:lamellae" 
    8484 
    85 single = False 
     85single = False  # TODO: check 
    8686# pylint: disable=bad-whitespace, line-too-long 
    8787#             ["name", "units", default, [lower, upper], "type","description"], 

  • sasmodels/models/linear_pearls.py

    r42356c8 rd2d6100  
    6161    ] 
    6262# pylint: enable=bad-whitespace, line-too-long 
     63single=False 
    6364 
    6465source = ["linear_pearls.c"] 

  • sasmodels/models/multilayer_vesicle.py

    r42356c8 rf67f26c  
    7474    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    7575    ["radius", "Ang", 60.0, [0.0, inf],  "", "Core radius of the multishell"], 
    76     ["thick_shell", "Ang",        10.0, [0.0, inf],  "sld", "Shell thickness"], 
     76    ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "Shell thickness"], 
    7777    ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "Water thickness"], 
    7878    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "Core scattering length density"], 

  • sasmodels/models/onion.c

    r639c4e3 rd119f34  
    22static double 
    33f_exp(double q, double r, double sld_in, double sld_out, 
    4     double thickness, double A) 
    5 { 
    6   const double vol = M_4PI_3 * cube(r); 
    7   const double qr = q * r; 
    8   const double alpha = A * r/thickness; 
    9   const double bes = sph_j1c(qr); 
    10   const double B = (sld_out - sld_in)/expm1(A); 
    11   const double C = sld_in - B; 
    12   double fun; 
    13   if (qr == 0.0) { 
    14     fun = 1.0; 
    15   } else if (fabs(A) > 0.0) { 
    16     const double qrsq = qr*qr; 
    17     const double alphasq = alpha*alpha; 
    18     const double sumsq = alphasq + qrsq; 
    19     double sinqr, cosqr; 
    20     SINCOS(qr, sinqr, cosqr); 
    21     fun = -3.0*( 
    22             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    23                 - (alpha*sinqr/qr - cosqr) 
    24         ) / sumsq; 
    25   } else { 
    26     fun = bes; 
    27   } 
    28   return vol * (B*fun + C*bes); 
    29 } 
    30  
    31 static double 
    32 f_linear(double q, double r, double sld, double slope) 
     4    double thickness, double A, double side) 
    335{ 
    346  const double vol = M_4PI_3 * cube(r); 
    357  const double qr = q * r; 
    368  const double bes = sph_j1c(qr); 
    37   double fun = 0.0; 
    38   if (qr > 0.0) { 
    39     const double qrsq = qr*qr; 
     9  const double alpha = A * r/thickness; 
     10  double result; 
     11  if (qr == 0.0) { 
     12    result = 1.0; 
     13  } else if (fabs(A) > 0.0) { 
     14    const double qrsq = qr * qr; 
     15    const double alphasq = alpha * alpha; 
     16    const double sumsq = alphasq + qrsq; 
    4017    double sinqr, cosqr; 
    4118    SINCOS(qr, sinqr, cosqr); 
    42     // Jae-He's code seems to simplify to this 
    43     //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); 
    44     // Rederiving the math, we get the following instead: 
    45     fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); 
     19    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr; 
     20    const double t2 = alpha*sinqr/qr - cosqr; 
     21    const double fun = -3.0*(t1/sumsq - t2)/sumsq; 
     22    const double slope = (sld_out - sld_in)/expm1(A); 
     23    const double contrast = slope*exp(A*side); 
     24    result = contrast*fun + (sld_in-slope)*bes; 
     25  } else { 
     26    result = sld_in*bes; 
    4627  } 
    47   return vol * (sld*bes + fun); 
    48 } 
    49  
    50 static double 
    51 f_constant(double q, double r, double sld) 
    52 { 
    53   const double bes = sph_j1c(q * r); 
    54   const double vol = M_4PI_3 * cube(r); 
    55   return sld * vol * bes; 
     28  return vol * result; 
    5629} 
    5730 
    … …  
    6942static double 
    7043Iq(double q, double sld_core, double core_radius, double sld_solvent, 
    71     double n, double sld_in[], double sld_out[], double thickness[], 
     44    double n_shells, double sld_in[], double sld_out[], double thickness[], 
    7245    double A[]) 
    7346{ 
    74   int i; 
    75   double r = core_radius; 
    76   double f = f_constant(q, r, sld_core); 
    77   for (i=0; i<n; i++){ 
    78     const double r0 = r; 
    79     r += thickness[i]; 
    80     if (r == r0) { 
    81       // no thickness, so nothing to add 
    82     } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) { 
    83       f -= f_constant(q, r0, sld_in[i]); 
    84       f += f_constant(q, r, sld_in[i]); 
    85     } else if (fabs(A[i]) < 1.0e-4) { 
    86       const double slope = (sld_out[i] - sld_in[i])/thickness[i]; 
    87       f -= f_linear(q, r0, sld_in[i], slope); 
    88       f += f_linear(q, r, sld_out[i], slope); 
    89     } else { 
    90       f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); 
    91       f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); 
    92     } 
     47  int n = (int)(n_shells+0.5); 
     48  double r_out = core_radius; 
     49  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0); 
     50  for (int i=0; i < n; i++){ 
     51    const double r_in = r_out; 
     52    r_out += thickness[i]; 
     53    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0); 
     54    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0); 
    9355  } 
    94   f -= f_constant(q, r, sld_solvent); 
     56  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0); 
    9557  const double f2 = f * f * 1.0e-4; 
    9658 

  • sasmodels/models/onion.py

    r63c6a08 rd119f34  
    394394    # Could also specify them individually as 
    395395    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1, 
     396    #"core_radius_pd_n": 10, 
     397    #"core_radius_pd": 0.4, 
     398    #"thickness4_pd_n": 10, 
     399    #"thickness4_pd": 0.4, 
    396400    } 
    397401 

  • sasmodels/models/polymer_micelle.py

    r42356c8 rd2d6100  
    5151# pylint: enable=bad-whitespace, line-too-long 
    5252 
     53single = False 
     54 
    5355source = ["lib/sph_j1c.c", "polymer_micelle_kernel.c"] 
    5456 

  • sasmodels/models/polymer_micelle_kernel.c

    r2c74c11 rc915373  
    3636    // Self-correlation term of the chains 
    3737    const double qrg2 = q*radius_gyr*q*radius_gyr; 
    38     const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(exp(-qrg2)-1+qrg2)/(qrg2*qrg2); 
     38    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2); 
    3939    const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain; 
    4040 
    4141    // Interference cross-term between core and chains 
    42     const double chain_ampl = (qrg2 == 0.0) ? 1.0 : (1-exp(-qrg2))/qrg2; 
     42    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; 
    4343    const double bes_corona = sinc(q*(radius_core + d_penetration * radius_gyr)); 
    4444    const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona * 

  • sasmodels/models/rpa.c

    rd680a2b r6351bfa  
    2525  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd; 
    2626  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd; 
    27   double S0da,S0db,S0dc,Pdd,S0dd; 
    28   double Kaa,Kbb,Kcc,Kdd; 
    29   double Kba,Kca,Kda,Kcb,Kdb,Kdc; 
     27  double S0da,S0db,S0dc; 
     28  double Pdd,S0dd; 
     29  double Kaa,Kbb,Kcc; 
     30  double Kba,Kca,Kcb; 
     31  double Kda,Kdb,Kdc,Kdd; 
    3032  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc; 
    3133  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33; 

  • sasmodels/models/surface_fractal.py

    rec45c4f r33875e3  
    8080parameters = [["radius",        "Ang", 10.0, [0, inf],   "", 
    8181               "Particle radius"], 
    82               ["surface_dim",   "",    2.0,  [0, inf],   "", 
     82              ["surface_dim",   "",    2.0,  [1, 3],   "", 
    8383               "Surface fractal dimension"], 
    8484              ["cutoff_length", "Ang", 500., [0.0, inf], "", 

  • sasmodels/models/unified_power_Rg.py

    r2c74c11 rec77322  
    8080 
    8181    with errstate(divide='ignore', invalid='ignore'): 
    82         result = np.empty(q.shape, 'd') 
     82        result = np.zeros(q.shape, 'd') 
    8383        for i in range(ilevel): 
    8484            exp_now = exp(-(q*rg[i])**2/3.) 

  • sasmodels/models/lib/sas_erf.c

    r1557a1e redf06e1  
    280280    else 
    281281        x = a;*/ 
    282  
    283     x = fabsf(a); 
     282    // TODO: tinycc does not support fabsf 
     283    x = fabs(a); 
    284284 
    285285 
    … …  
    302302            return (0.0); 
    303303    } 
    304  
    305304    z = expf(z); 
    306305 
    … …  
    332331    float y, z; 
    333332 
    334     if (fabsf(x) > 1.0) 
     333    // TODO: tinycc does not support fabsf 
     334    if (fabs(x) > 1.0) 
    335335        return (1.0 - erfcf(x)); 
    336336 

  • sasmodels/models/spherical_sld.py

    r63c6a08 redf06e1  
    217217              ] 
    218218# pylint: enable=bad-whitespace, line-too-long 
    219 source = ["lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     219source = ["lib/polevl.c", "lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
    220220single = False  # TODO: fix low q behaviour 
    221221 
Note: See TracChangeset for help on using the changeset viewer.