Changeset d119f34 in sasmodels


Ignore:
Timestamp:
Aug 2, 2016 12:13:29 PM (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:
0dd4199, 3a45c2c
Parents:
f95556f
Message:

onion: now matches sasview

Location:
sasmodels
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/convert.py

    r2c74c11 rd119f34  
    316316        elif name == 'core_multi_shell': 
    317317            pars['n'] = min(math.ceil(pars['n']), 4) 
     318        elif name == 'onion': 
     319            pars['n_shells'] = math.ceil(pars['n_shells']) 
    318320        elif name == 'spherical_sld': 
    319321            for k in range(1, 11): 
  • 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 
Note: See TracChangeset for help on using the changeset viewer.