Changeset 54bcd4a in sasmodels


Ignore:
Timestamp:
Aug 4, 2016 4:48:37 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:
745b7bb
Parents:
bd49c79
Message:

spherical sld: simplify code so that it works on AMD GPUs

Location:
sasmodels
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/conversion_table.py

    rf67f26c r54bcd4a  
    688688        } 
    689689    ], 
    690     "magsphere": [ 
    691         "SphereModel", 
    692         { 
    693             "sld": "sldSph", 
    694             "radius": "radius", 
    695             "sld_solvent": "sldSolv", 
    696             "sld_M0": "M0_sld_sph", 
    697             "sphere_theta": "M_theta_sph", 
    698             "sphere_phi": "M_phi_sph", 
    699             "sld_solvent_M0": "M0_sld_solv", 
    700             "solvent_theta": "M_theta_solv", 
    701             "solvent_phi": "M_phi_solv", 
    702             "in_spin": "Up_frac_i", 
    703             "out_spin": "Up_frac_f", 
    704             "spin_theta": "Up_theta" 
    705         } 
    706     ], 
    707     "_spherepy": [ 
    708         "SphereModel", 
    709         { 
    710             "sld": "sldSph", 
    711             "radius": "radius", 
    712             "sld_solvent": "sldSolv" 
    713         } 
    714     ], 
    715690    "spherical_sld": [ 
    716691        "SphericalSLDModel", 
    717         { 
    718             "radius_core": "rad_core0", 
    719             "sld_core": "sld_core0", 
    720             "sld_solvent": "sld_solv" 
    721         } 
     692        # Be lazy and use a generator expression to define 
     693        #    sld1: sld_flat0, ... 
     694        #    thickness1: thick_flat0, ... 
     695        #    interface1: thick_inter0, ... 
     696        #    shape1: func_inter0, ... 
     697        #    nu1: nu_inter0, ... 
     698        # but override thickness1 => rad_cor0 and sld1 => sld_core0. 
     699        # Note: explicit key,value pairs given by **{...} override the 
     700        # keys from the gnerator expression ((k,v) for k,v in seq) when 
     701        # used as dict((generator), **{...}) 
     702        dict(((field_new+str(index+1),field_old+str(index)) 
     703             for field_new, field_old in [("sld","sld_flat"), 
     704                                          ("thickness","thick_flat"), 
     705                                          ("interface","thick_inter"), 
     706                                          ("shape","func_inter"), 
     707                                          ("nu","nu_inter"),] 
     708             for index in range(11)), 
     709             **{ 
     710            "n_shells": "n_shells", 
     711            "n_steps": "npts_inter", 
     712            "sld_solvent": "sld_solv", 
     713            "thickness1": "rad_core0", 
     714            "sld1": "sld_core0", 
     715        }) 
    722716    ], 
    723717    "squarewell": [ 
  • sasmodels/convert.py

    rbd49c79 r54bcd4a  
    229229                oldpars['ndensity'] *= 1e15 
    230230        elif name == 'spherical_sld': 
    231             for k in range(1, int(pars['n_shells'])+1): 
    232                 _remove_pd(oldpars, 'thick_flat'+str(k), 'thick_flat') 
    233                 _remove_pd(oldpars, 'thick_inter'+str(k), 'thick_inter') 
     231            oldpars["CONTROL"] -= 1 
     232            # remove polydispersity from shells 
     233            for k in range(1, 11): 
     234                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness') 
     235                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface') 
     236            # remove extra shells 
     237            for k in range(int(pars['n_shells']), 11): 
     238                oldpars.pop('sld_flat'+str(k), 0) 
     239                oldpars.pop('thick_flat'+str(k), 0) 
     240                oldpars.pop('thick_inter'+str(k), 0) 
     241                oldpars.pop('func_inter'+str(k), 0) 
     242                oldpars.pop('nu_inter'+str(k), 0) 
    234243        elif name == 'core_multi_shell': 
    235244            # kill extra shells 
     
    265274                    oldpars.pop(k, None) 
    266275 
     276    #print("convert from",list(sorted(pars))) 
     277    #print("convert to",list(sorted(oldpars.items()))) 
    267278    return oldpars 
    268279 
     
    323334        elif name == 'spherical_sld': 
    324335            pars['n_shells'] = math.ceil(pars['n_shells']) 
    325             pars['npts_inter'] = math.ceil(pars['npts_inter']) 
    326             pars['func_inter0'] = math.trunc(pars['func_inter0']+0.5) 
    327             for k in range(1, 11): 
    328                 pars['func_inter%d'%k] = math.trunc(pars['func_inter%d'%k]+0.5) 
    329                 pars['thick_flat%d_pd_n'%k] = 0 
    330                 pars['thick_inter%d_pd_n'%k] = 0 
    331  
     336            pars['n_steps'] = math.ceil(pars['n_steps']) 
     337            for k in range(1, 12): 
     338                pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5) 
     339            for k in range(2, 12): 
     340                pars['thickness%d_pd_n'%k] = 0 
     341                pars['interface%d_pd_n'%k] = 0 
     342 
  • sasmodels/models/spherical_sld.c

    rdd71228 r54bcd4a  
    11static double form_volume( 
    22    int n_shells, 
    3     double radius_core, 
    4     double thick_inter0, 
    5     double thick_flat[], 
    6     double thick_inter[]) 
     3    double thickness[], 
     4    double interface[]) 
    75{ 
    8     int i; 
    9     double r = radius_core; 
    10     r += thick_inter0; 
    11     for (i=0; i < n_shells; i++) { 
    12         r += thick_inter[i]; 
    13         r += thick_flat[i]; 
     6    double r = 0.0; 
     7    for (int i=0; i < n_shells; i++) { 
     8        r += thickness[i] + interface[i]; 
    149    } 
    1510    return M_4PI_3*cube(r); 
    1611} 
    1712 
     13static double blend(int shape, double nu, double z) 
     14{ 
     15    if (shape==0) { 
     16        const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0)); 
     17        const double denom = 2.0 * sas_erf(nu * M_SQRT1_2); 
     18        return num/denom + 0.5; 
     19    } else if (shape==1) { 
     20        return pow(z, nu); 
     21    } else if (shape==2) { 
     22        return 1.0 - pow(1. - z, nu); 
     23    } else if (shape==3) { 
     24        return expm1(-nu*z)/expm1(-nu); 
     25    } else if (shape==4) { 
     26        return expm1(nu*z)/expm1(nu); 
     27    } else { 
     28        return NAN; 
     29    } 
     30} 
    1831 
    19 static double sphere_sld_kernel( 
     32static double f_linear(double q, double r, double contrast, double slope) 
     33{ 
     34    const double qr = q * r; 
     35    const double qrsq = qr * qr; 
     36    const double bes = sph_j1c(qr); 
     37    double sinqr, cosqr; 
     38    SINCOS(qr, sinqr, cosqr); 
     39    const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); 
     40    const double vol = M_4PI_3 * cube(r); 
     41    return vol*(bes*contrast + fun*slope); 
     42} 
     43 
     44static double Iq( 
    2045    double q, 
    2146    int n_shells, 
    22     int npts_inter, 
    23     double radius_core, 
    24     double sld_core, 
    2547    double sld_solvent, 
    26     double func_inter_core, 
    27     double thick_inter_core, 
    28     double nu_inter_core, 
    29     double sld_flat[], 
    30     double thick_flat[], 
    31     double func_inter[], 
    32     double thick_inter[], 
    33     double nu_inter[] ) { 
     48    double sld[], 
     49    double thickness[], 
     50    double interface[], 
     51    double shape[], 
     52    double nu[], 
     53    int n_steps) 
     54{ 
     55    // iteration for # of shells + core + solvent 
     56    double f=0.0; 
     57    double r=0.0; 
     58    for (int shell=0; shell<n_shells; shell++){ 
     59        const double sld_l = sld[shell]; 
    3460 
    35     int i,j,k; 
    36     int n_s; 
     61        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well. 
     62        f -= M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r); 
     63        r += thickness[shell]; 
     64        f += M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r); 
    3765 
    38     double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2; 
    39     double sign,slope=0.0; 
    40     double pi; 
     66        // iterate over sub_shells in the interface 
     67        const double dr = interface[shell]/n_steps; 
     68        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l; 
     69        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14); 
     70        const int shape_shell = (int)(shape[shell]); 
    4171 
    42     double total_thick=0.0; 
     72        // if there is no interface the equations don't work 
     73        if (dr == 0.) continue; 
    4374 
    44     //TODO: This part can be further simplified 
    45     int fun_type[12]; 
    46     double sld[12]; 
    47     double thick_internal[12]; 
    48     double thick[12]; 
    49     double fun_coef[12]; 
     75        double sld_in = sld_l; 
     76        for (int step=1; step <= n_steps; step++) { 
     77            // find sld_i at the outer boundary of sub-shell step 
     78            //const double z = (double)step/(double)n_steps; 
     79            const double z = (double)step/(double)n_steps; 
     80            const double fraction = blend(shape_shell, nu_shell, z); 
     81            const double sld_out = fraction*delta + sld_l; 
     82            // calculate slope 
     83            const double slope = (sld_out - sld_in)/dr; 
     84            const double contrast = sld_in - slope*r; 
    5085 
    51     fun_type[0] = func_inter_core; 
    52     fun_coef[0] = fabs(nu_inter_core); 
    53     sld[0] = sld_core; 
    54     thick[0] = radius_core; 
    55     thick_internal[0] = thick_inter_core; 
    56  
    57     for (i =1; i<=n_shells; i++){ 
    58         sld[i] = sld_flat[i-1]; 
    59         thick_internal[i]= thick_inter[i-1]; 
    60         thick[i] = thick_flat[i-1]; 
    61         fun_type[i] = func_inter[i-1]; 
    62         fun_coef[i] = fabs(nu_inter[i-1]); 
    63         total_thick += thick[i]; 
    64         total_thick += thick_internal[i]; //doesn't account for core layer 
    65     } 
    66  
    67     sld[n_shells+1] = sld_solvent; 
    68     thick[n_shells+1] = total_thick/5.0; 
    69     thick_internal[n_shells+1] = 0.0; 
    70     fun_coef[n_shells+1] = 0.0; 
    71     fun_type[n_shells+1] = 0; 
    72  
    73     pi = 4.0*atan(1.0); 
    74     f = 0.0; 
    75     r = 0.0; 
    76     vol = 0.0; 
    77     sld_f = sld_core; 
    78  
    79     dz = 0.0; 
    80     // iteration for # of shells + core + solvent 
    81     for (i=0;i<=n_shells+1; i++){ 
    82         // iteration for flat and interface 
    83         for (j=0;j<2;j++){ 
    84             // iteration for sub_shells in the interface 
    85             // starts from #1 sub-layer 
    86             for (n_s=1;n_s<=npts_inter; n_s++){ 
    87                 // for solvent, it doesn't have an interface 
    88                 if (i==n_shells+1 && j==1) 
    89                     break; 
    90                 // for flat layers 
    91                 if (j==0){ 
    92                     dz = thick[i]; 
    93                     sld_i = sld[i]; 
    94                     slope = 0.0; 
    95                 } 
    96                 // for interfacial sub_shells 
    97                 else{ 
    98                     dz = thick_internal[i]/npts_inter; 
    99                     // find sld_i at the outer boundary of sub-layer #n_s 
    100                     sld_i = intersldfunc(fun_type[i], npts_inter, n_s, 
    101                             fun_coef[i], sld[i], sld[i+1]); 
    102                     // calculate slope 
    103                     slope= (sld_i -sld_f)/dz; 
    104                 } 
    105                 contr = sld_f-slope*r; 
    106                 // iteration for the left and right boundary of the shells 
    107                 for (k=0; k<2; k++){ 
    108                     // At r=0, the contribution to I is zero, so skip it. 
    109                     if ( i == 0 && j == 0 && k == 0){ 
    110                         continue; 
    111                     } 
    112                     // On the top of slovent there is no interface; skip it. 
    113                     if (i == n_shells+1 && k == 1){ 
    114                         continue; 
    115                     } 
    116                     // At the right side (outer) boundary 
    117                     if ( k == 1){ 
    118                         sign = 1.0; 
    119                         r += dz; 
    120                     } 
    121                     // At the left side (inner) boundary 
    122                     else{ 
    123                         sign = -1.0; 
    124                     } 
    125  
    126                     qr = q * r; 
    127                     fun = 0.0; 
    128  
    129                     if(qr == 0.0){ 
    130                         bes = sign * 1.0; 
    131                     } 
    132                     else{ 
    133                         // for flat sub-layer 
    134                         //TODO: Single precision calculation fails here 
    135                         bes = sign *  sph_j1c(qr); 
    136                         if (fabs(slope) > 0.0 ){ 
    137                             const double qrsq = qr*qr; 
    138                             double sinqr, cosqr; 
    139                             SINCOS(qr, sinqr, cosqr); 
    140                             fun = sign * 3.0 * r * 
    141                             (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq); 
    142                             // In the onion model Jae-He's formula is rederived 
    143                             // and gives following: 
    144                             //fun = 3.0 * sign * r * 
    145                             //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq); 
    146                             //But this seems not to be working in this case... 
    147                         } 
    148                     } 
    149  
    150                     // update total volume 
    151                     vol = M_4PI_3 * cube(r); 
    152                     f += vol * (bes * contr + fun * slope); 
    153                 } 
    154                 sld_f = sld_i; 
    155                 // no sub-layer iteration (n_s loop) for the flat layer 
    156                 if (j==0) 
    157                     break; 
    158             } 
     86            // iteration for the left and right boundary of the shells 
     87            f -= f_linear(q, r, contrast, slope); 
     88            r += dr; 
     89            f += f_linear(q, r, contrast, slope); 
     90            sld_in = sld_out; 
    15991        } 
    16092    } 
     93    // add in solvent effect 
     94    f -= M_4PI_3 * cube(r) * sld_solvent * sph_j1c(q*r); 
    16195 
    162     f2 = f * f * 1.0e-4; 
    163     return (f2); 
     96    const double f2 = f * f * 1.0e-4; 
     97    return f2; 
    16498} 
    16599 
    166  
    167 /** 
    168  * Function to evaluate 1D SphereSLD function 
    169  * @param q: q-value 
    170  * @return: function value 
    171  */ 
    172 static double Iq(double q, 
    173     int n_shells, 
    174     int npts_inter, 
    175     double radius_core, 
    176     double sld_core, 
    177     double sld_solvent, 
    178     double func_inter0, 
    179     double thick_inter0, 
    180     double nu_inter0, 
    181     double sld_flat[], 
    182     double thick_flat[], 
    183     double func_inter[], 
    184     double thick_inter[], 
    185     double nu_inter[] ) { 
    186  
    187     double intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core, 
    188                 sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0, 
    189                 sld_flat, thick_flat, func_inter, thick_inter, nu_inter); 
    190  
    191     return intensity; 
    192 } 
    193  
    194 /** 
    195  * Function to evaluate 2D SphereSLD function 
    196  * @param q_x: value of Q along x 
    197  * @param q_y: value of Q along y 
    198  * @return: function value 
    199  */ 
    200  
    201 /*static double Iqxy(double qx, double qy, 
    202     int n_shells, 
    203     int npts_inter, 
    204     double radius_core 
    205     double sld_core, 
    206     double sld_solvent, 
    207     double sld_flat[], 
    208     double thick_flat[], 
    209     double func_inter[], 
    210     double thick_inter[], 
    211     double nu_inter[], 
    212     ) { 
    213  
    214     double q = sqrt(qx*qx + qy*qy); 
    215     return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent, 
    216     sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[]) 
    217 }*/ 
    218  
  • sasmodels/models/spherical_sld.py

    r4fd2c63 r54bcd4a  
    190190 
    191191import numpy as np 
    192 from numpy import inf 
     192from numpy import inf, expm1, sqrt 
     193from scipy.special import erf 
    193194 
    194195name = "spherical_sld" 
     
    200201category = "shape:sphere" 
    201202 
     203SHAPES = ["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)", 
     204          "Rexp(-|nu|z)", "Lexp(-|nu|z)"], 
     205 
    202206# pylint: disable=bad-whitespace, line-too-long 
    203207#            ["name", "units", default, [lower, upper], "type", "description"], 
    204 parameters = [["n_shells",          "",               1,      [0, 10],        "volume", "number of shells"], 
    205               ["npts_inter",        "",               35,     [0, inf],       "", "number of points in each sublayer Must be odd number"], 
    206               ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
    207               ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "sld", "sld function flat"], 
    208               ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "sld", "sld function solvent"], 
    209               ["func_inter0",       "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"], 
    210               ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"], 
    211               ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"], 
    212               ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "sld", "sld function flat"], 
    213               ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"], 
    214               ["func_inter[n_shells]",    "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"], 
    215               ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"], 
    216               ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"], 
     208parameters = [["n_shells",             "",           1,      [1, 11],        "volume", "number of shells"], 
     209              ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"], 
     210              ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"], 
     211              ["thickness[n_shells]",  "Ang",        100.0,  [0, inf],       "volume", "thickness shell"], 
     212              ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"], 
     213              ["shape[n_shells]",      "",           0,      SHAPES,         "", "interface shape"], 
     214              ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"], 
     215              ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"], 
    217216              ] 
    218217# pylint: enable=bad-whitespace, line-too-long 
    219 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     218source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sph_j1c.c", "spherical_sld.c"] 
    220219single = False  # TODO: fix low q behaviour 
    221220 
    222221profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    223 def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat, 
    224             thick_flat, func_inter, thick_inter, nu_inter, npts_inter): 
     222 
     223SHAPE_FUNCTIONS = [ 
     224    lambda z, nu: erf(nu/sqrt(2)*(2*z-1))/(2*erf(nu/sqrt(2))) + 0.5,  # erf 
     225    lambda z, nu: z**nu,                    # Rpow 
     226    lambda z, nu: 1 - (1-z)**nu,            # Lpow 
     227    lambda z, nu: expm1(-nu*z)/expm1(-nu),  # Rexp 
     228    lambda z, nu: expm1(nu*z)/expm1(nu),    # Lexp 
     229] 
     230 
     231def profile(n_shells, sld_solvent, sld, thickness, 
     232            interface, shape, nu, n_steps): 
    225233    """ 
    226234    Returns shape profile with x=radius, y=SLD. 
     
    228236 
    229237    z = [] 
    230     beta = [] 
     238    rho = [] 
    231239    z0 = 0 
    232240    # two sld points for core 
    233241    z.append(0) 
    234     beta.append(sld_core) 
    235     z.append(radius_core) 
    236     beta.append(sld_core) 
    237     z0 += radius_core 
    238  
    239     for i in range(1, n_shells+2): 
    240         dz = thick_inter[i-1]/npts_inter 
    241         # j=0 for interface, j=1 for flat layer 
    242         for j in range(0, 2): 
    243             # interation for sub-layers 
    244             for n_s in range(0, npts_inter+1): 
    245                 if j == 1: 
    246                     if i == n_shells+1: 
    247                         break 
    248                     # shift half sub thickness for the first point 
    249                     z0 -= dz#/2.0 
    250                     z.append(z0) 
    251                     #z0 -= dz/2.0 
    252                     z0 += thick_flat[i] 
    253                     sld_i = sld_flat[i] 
    254                     beta.append(sld_flat[i]) 
    255                     dz = 0 
    256                 else: 
    257                     nu = nu_inter[i-1] 
    258                     # decide which sld is which, sld_r or sld_l 
    259                     if i == 1: 
    260                         sld_l = sld_core 
    261                     else: 
    262                         sld_l = sld_flat[i-1] 
    263                     if i == n_shells+1: 
    264                         sld_r = sld_solvent 
    265                     else: 
    266                         sld_r = sld_flat[i] 
    267                     # get function type 
    268                     func_idx = func_inter[i-1] 
    269                     # calculate the sld 
    270                     sld_i = intersldfunc(func_idx, npts_inter, n_s, nu, 
    271                                             sld_l, sld_r) 
    272                 # append to the list 
    273                 z.append(z0) 
    274                 beta.append(sld_i) 
    275                 z0 += dz 
    276                 if j == 1: 
    277                     break 
    278     z.append(z0) 
    279     beta.append(sld_solvent) 
    280     z_ext = z0/5.0 
    281     z.append(z0+z_ext) 
    282     beta.append(sld_solvent) 
     242    rho.append(sld[0]) 
     243 
     244    for i in range(0, n_shells): 
     245        z0 += thickness[i] 
     246        z.append(z0) 
     247        rho.append(sld[i]) 
     248        dz = interface[i]/n_steps 
     249        sld_l = sld[i] 
     250        sld_r = sld[i+1] if i < n_shells-1 else sld_solvent 
     251        interface = SHAPE_FUNCTIONS[int(np.clip(shape[i], 0, len(SHAPES)-1))] 
     252        for step in range(1, n_steps+1): 
     253            portion = interface(float(step)/n_steps, max(abs(nu[i]), 1e-14)) 
     254            z0 += dz 
     255            z.append(z0) 
     256            rho.append((sld_r - sld_l)*portion + sld_l) 
     257    z.append(z0*1.2) 
     258    rho.append(sld_solvent) 
    283259    # return sld profile (r, beta) 
    284     return np.asarray(z), np.asarray(beta)*1e-6 
    285  
    286 def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat): 
     260    return np.asarray(z), np.asarray(rho)*1e-6 
     261 
     262 
     263def ER(n_shells, thickness, interface): 
    287264    n_shells = int(n_shells) 
    288     total_thickness = thick_inter0 
    289     total_thickness += np.sum(thick_inter[:n_shells], axis=0) 
    290     total_thickness += np.sum(thick_flat[:n_shells], axis=0) 
    291     return total_thickness + radius_core 
     265    total = (np.sum(thickness[:n_shells], axis=1) 
     266             + np.sum(interface[:n_shells], axis=1)) 
     267    return total 
    292268 
    293269 
    294270demo = { 
    295     "n_shells": 4, 
    296     "npts_inter": 35.0, 
    297     "radius_core": 50.0, 
    298     "sld_core": 2.07, 
     271    "n_shells": 5, 
     272    "n_steps": 35.0, 
    299273    "sld_solvent": 1.0, 
    300     "thick_inter0": 50.0, 
    301     "func_inter0": 0, 
    302     "nu_inter0": 2.5, 
    303     "sld_flat":[4.0,3.5,4.0,3.5], 
    304     "thick_flat":[100.0,100.0,100.0,100.0], 
    305     "func_inter":[0,0,0,0], 
    306     "thick_inter":[50.0,50.0,50.0,50.0], 
    307     "nu_inter":[2.5,2.5,2.5,2.5], 
     274    "sld":[2.07,4.0,3.5,4.0,3.5], 
     275    "thickness":[50.0,100.0,100.0,100.0,100.0], 
     276    "interface":[50.0,50.0,50.0,50.0], 
     277    "shape": [0,0,0,0,0], 
     278    "nu":[2.5,2.5,2.5,2.5,2.5], 
    308279    } 
    309280 
     
    312283tests = [ 
    313284    # Accuracy tests based on content in test/utest_extra_models.py 
    314     [{"n_shells":4, 
    315         'npts_inter':35, 
    316         "radius_core":50.0, 
    317         "sld_core":2.07, 
     285    [{"n_shells": 5, 
     286        "n_steps": 35, 
    318287        "sld_solvent": 1.0, 
    319         "sld_flat":[4.0,3.5,4.0,3.5], 
    320         "thick_flat":[100.0,100.0,100.0,100.0], 
    321         "func_inter":[0,0,0,0], 
    322         "thick_inter":[50.0,50.0,50.0,50.0], 
    323         "nu_inter":[2.5,2.5,2.5,2.5] 
     288        "sld": [2.07, 4.0, 3.5, 4.0, 3.5], 
     289        "thickness": [50.0, 100.0, 100.0, 100.0, 100.0], 
     290        "interface": [50]*5, 
     291        "shape": [0]*5, 
     292        "nu": [2.5]*5, 
    324293    }, 0.001, 0.001], 
    325294] 
Note: See TracChangeset for help on using the changeset viewer.