Changeset 4f0c993 in sasmodels


Ignore:
Timestamp:
Jun 28, 2017 10:14:45 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
4f611f1
Parents:
73cbc5b (diff), c935936 (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 branch 'master' into ticket-776-orientation

Files:
13 added
29 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    r2a55a6f r8a5f021  
    2020/.idea 
    2121/sasmodels.egg-info/ 
     22/example/Fit_*/ 
     23/example/batch_fit.csv 

  • explore/jitter.py

    r85190c2 rd4c33d6  
    33dispersity and possible replacement algorithms. 
    44""" 
     5import sys 
     6 
    57import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    68import matplotlib.pyplot as plt 
    79from matplotlib.widgets import Slider, CheckButtons 
    810from matplotlib import cm 
    9  
    1011import numpy as np 
    1112from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1213 
    13 def draw_beam(ax): 
     14def draw_beam(ax, view=(0, 0)): 
    1415    #ax.plot([0,0],[0,0],[1,-1]) 
    1516    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    2223    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2324    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     25    z = 1.3*np.outer(np.ones_like(u), v) 
     26 
     27    theta, phi = view 
     28    shape = x.shape 
     29    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     30    points = Rz(phi)*Ry(theta)*points 
     31    x, y, z = [v.reshape(shape) for v in points] 
    2532 
    2633    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2734 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     35def draw_jitter(ax, view, jitter): 
     36    size = [0.1, 0.4, 1.0] 
     37    draw_shape = draw_parallelepiped 
     38    #draw_shape = draw_ellipsoid 
    3439 
    3540    #np.random.seed(10) 
    … …  
    6469        [ 1,  1,  1], 
    6570    ] 
     71    dtheta, dphi, dpsi = jitter 
    6672    if dtheta == 0: 
    6773        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    7076    if dpsi == 0: 
    7177        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     78    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    7379    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     80        delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     81        draw_shape(ax, size, view, delta, alpha=0.8) 
    7682    for v in 'xyz': 
    7783        a, b, c = size 
    … …  
    8086        getattr(ax, v+'axis').label.set_text(v) 
    8187 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     88def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    8389    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    8790    u = np.linspace(0, 2 * np.pi, steps) 
    8891    v = np.linspace(0, np.pi, steps) 
    … …  
    9093    y = b*np.outer(np.sin(u), np.sin(v)) 
    9194    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     95    x, y, z = transform_xyz(view, jitter, x, y, z) 
    9896 
    9997    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    10098 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     99    draw_labels(ax, view, jitter, [ 
     100         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     101         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     102         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     103         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     104         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     105         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     106    ]) 
     107 
     108def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    102109    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106110    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107111    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    118122    ]) 
    119123 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     124    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126126 
    127 def draw_sphere(ax, radius=10., steps=100): 
    128     u = np.linspace(0, 2 * np.pi, steps) 
    129     v = np.linspace(0, np.pi, steps) 
    130  
    131     x = radius * np.outer(np.cos(u), np.sin(v)) 
    132     y = radius * np.outer(np.sin(u), np.sin(v)) 
    133     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    134     ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135  
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
     127    draw_labels(ax, view, jitter, [ 
     128         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     129         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     130         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     131         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     132         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     133         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     134    ]) 
     135 
     136def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 
     137    theta, phi, psi = view 
     138    dtheta, dphi, dpsi = jitter 
    145139    if dist == 'gauss': 
     140        t = np.linspace(-3, 3, n) 
    146141        weights = exp(-0.5*t**2) 
    147142    elif dist == 'rect': 
     143        t = np.linspace(0, 1, n) 
    148144        weights = np.ones_like(t) 
    149145    else: 
    150146        raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
    163     ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164  
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
     147 
     148    # mesh in theta, phi formed by rotating z 
     149    z = np.matrix([[0], [0], [radius]]) 
     150    points = np.hstack([Rx(phi_i)*Ry(theta_i)*z 
     151                        for theta_i in dtheta*t 
     152                        for phi_i in dphi*t]) 
     153    # rotate relative to beam 
     154    points = orient_relative_to_beam(view, points) 
     155 
     156    w = np.outer(weights, weights) 
     157 
     158    x, y, z = [np.array(v).flatten() for v in points] 
     159    ax.scatter(x, y, z, c=w.flatten(), marker='o', vmin=0., vmax=1.) 
    169160 
    170161def Rx(angle): 
    … …  
    188179         [0., 0., 1.]] 
    189180    return np.matrix(R) 
     181 
     182def transform_xyz(view, jitter, x, y, z): 
     183    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     184    shape = x.shape 
     185    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     186    points = apply_jitter(jitter, points) 
     187    points = orient_relative_to_beam(view, points) 
     188    x, y, z = [np.array(v).reshape(shape) for v in points] 
     189    return x, y, z 
     190 
     191def apply_jitter(jitter, points): 
     192    dtheta, dphi, dpsi = jitter 
     193    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     194    return points 
     195 
     196def orient_relative_to_beam(view, points): 
     197    theta, phi, psi = view 
     198    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     199    return points 
     200 
     201def draw_labels(ax, view, jitter, text): 
     202    labels, locations, orientations = zip(*text) 
     203    px, py, pz = zip(*locations) 
     204    dx, dy, dz = zip(*orientations) 
     205 
     206    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     207    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     208 
     209    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     210        zdir = np.asarray(zdir).flatten() 
     211        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     212 
     213def draw_sphere(ax, radius=10., steps=100): 
     214    u = np.linspace(0, 2 * np.pi, steps) 
     215    v = np.linspace(0, np.pi, steps) 
     216 
     217    x = radius * np.outer(np.cos(u), np.sin(v)) 
     218    y = radius * np.outer(np.sin(u), np.sin(v)) 
     219    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
     220    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    190221 
    191222def main(): 
    … …  
    206237 
    207238    axcolor = 'lightgoldenrodyellow' 
     239 
    208240    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209241    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    212244    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213245    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     246 
    214247    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215248    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    … …  
    217250    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218251    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
     252    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 
    220253 
    221254    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     255        view = stheta.val, sphi.val, spsi.val 
     256        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    224257        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     258        draw_beam(ax, (0, 0)) 
     259        draw_jitter(ax, view, jitter) 
     260        #draw_jitter(ax, view, (0,0,0)) 
     261        draw_mesh(ax, view, jitter) 
    229262        plt.gcf().canvas.draw() 
    230263 

  • sasmodels/details.py

    rccd5f01 rf39759c  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin 
     17from numpy import pi, cos, sin, radians 
    1818 
    1919try: 
    … …  
    2929 
    3030try: 
    31     from typing import List 
     31    from typing import List, Tuple, Sequence 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
     36    from .modelinfo import ParameterTable 
    3637 
    3738 
    … …  
    5354    coordinates, the total circumference decreases as latitude varies from 
    5455    pi r^2 at the equator to 0 at the pole, and the weight associated 
    55     with a range of phi values needs to be scaled by this circumference. 
     56    with a range of latitude values needs to be scaled by this circumference. 
    5657    This scale factor needs to be updated each time the theta value 
    5758    changes.  *theta_par* indicates which of the values in the parameter 
    … …  
    231232    nvalues = kernel.info.parameters.nvalues 
    232233    scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    233     values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
     234    # skipping scale and background when building values and weights 
     235    values, weights = zip(*pairs[2:npars+2]) if npars else ((), ()) 
     236    weights = correct_theta_weights(kernel.info.parameters, values, weights) 
    234237    length = np.array([len(w) for w in weights]) 
    235238    offset = np.cumsum(np.hstack((0, length))) 
    … …  
    244247    return call_details, data, is_magnetic 
    245248 
     249def correct_theta_weights(parameters, values, weights): 
     250    # type: (ParameterTable, Sequence[np.ndarray], Sequence[np.ndarray]) -> Sequence[np.ndarray] 
     251    """ 
     252    If there is a theta parameter, update the weights of that parameter so that 
     253    the cosine weighting required for polar integration is preserved.  Avoid 
     254    evaluation strictly at the pole, which would otherwise send the weight to 
     255    zero. 
     256 
     257    Note: values and weights do not include scale and background 
     258    """ 
     259    # TODO: document code, explaining why scale and background are skipped 
     260    # given that we don't have scale and background in the list, we 
     261    # should be calling the variables something other than values and weights 
     262    # Apparently the parameters.theta_offset similarly skips scale and 
     263    # and background, so the indexing works out. 
     264    if parameters.theta_offset >= 0: 
     265        index = parameters.theta_offset 
     266        theta = values[index] 
     267        theta_weight = np.minimum(abs(cos(radians(theta))), 1e-6) 
     268        # copy the weights list so we can update it 
     269        weights = list(weights) 
     270        weights[index] = theta_weight*np.asarray(weights[index]) 
     271        weights = tuple(weights) 
     272    return weights 
     273 
    246274 
    247275def convert_magnetism(parameters, values): 
     276    # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    248277    """ 
    249278    Convert magnetism values from polar to rectangular coordinates. 
    … …  
    255284    scale = mag[:,0] 
    256285    if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     286        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    258287        cos_theta = cos(theta) 
    259288        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    … …  
    269298    """ 
    270299    Create a mesh grid of dispersion parameters and weights. 
     300 
     301    pars is a list of pairs (values, weights), where the values are the 
     302    individual parameter values at which to evaluate the polydispersity 
     303    distribution and weights are the weights associated with each value. 
     304 
     305    Only the volume parameters should be included in this list.  Orientation 
     306    parameters do not affect the calculation of effective radius or volume 
     307    ratio. 
    271308 
    272309    Returns [p1,p2,...],w where pj is a vector of values for parameter j 

  • sasmodels/kernel_header.c

    rbb4b509 r73cbc5b  
    164164    SINCOS(phi*M_PI_180, sn, cn); \ 
    165165    q = sqrt(qx*qx + qy*qy); \ 
    166     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     166    cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
    167167    sn = sqrt(1 - cn*cn); \ 
    168168    } while (0) 

  • sasmodels/kernel_iq.c

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
    … …  
    173173 
    174174 
    175 #if MAX_PD>0 
    176   const int theta_par = details->theta_par; 
    177   const int fast_theta = (theta_par == p0); 
    178   const int slow_theta = (theta_par >= 0 && !fast_theta); 
    179   double spherical_correction = 1.0; 
    180 #else 
    181   // Note: if not polydisperse the weights cancel and we don't need the 
    182   // spherical correction. 
    183   const double spherical_correction = 1.0; 
    184 #endif 
    185  
    186175  int step = pd_start; 
    187176 
    … …  
    220209#endif 
    221210#if MAX_PD>0 
    222   if (slow_theta) { // Theta is not in inner loop 
    223     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    224   } 
    225211  while(i0 < n0) { 
    226212    local_values.vector[p0] = v0[i0]; 
    227213    double weight0 = w0[i0] * weight1; 
    228214//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
    229     if (fast_theta) { // Theta is in inner loop 
    230       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    231     } 
    232215#else 
    233216    const double weight0 = 1.0; 
    … …  
    244227      // Note: weight==0 must always be excluded 
    245228      if (weight0 > cutoff) { 
    246         // spherical correction is set at a minimum of 1e-6, otherwise there 
    247         // would be problems looking at models with theta=90. 
    248         const double weight = weight0 * spherical_correction; 
    249         pd_norm += weight * CALL_VOLUME(local_values.table); 
     229        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    250230 
    251231        #ifdef USE_OPENMP 
    … …  
    304284#endif // !MAGNETIC 
    305285//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 
    306           result[q_index] += weight * scattering; 
     286          result[q_index] += weight0 * scattering; 
    307287        } 
    308288      } 

  • sasmodels/kernel_iq.cl

    rbde38b5 rd4c33d6  
    2525    int32_t num_weights;        // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    27     int32_t theta_par;          // id of spherical correction variable 
     27    int32_t theta_par;          // id of spherical correction variable (not used) 
    2828} ProblemDetails; 
    2929 
    … …  
    169169 
    170170 
    171 #if MAX_PD>0 
    172   const int theta_par = details->theta_par; 
    173   const bool fast_theta = (theta_par == p0); 
    174   const bool slow_theta = (theta_par >= 0 && !fast_theta); 
    175   double spherical_correction = 1.0; 
    176 #else 
    177   // Note: if not polydisperse the weights cancel and we don't need the 
    178   // spherical correction. 
    179   const double spherical_correction = 1.0; 
    180 #endif 
    181  
    182171  int step = pd_start; 
    183172 
    … …  
    217206#endif 
    218207#if MAX_PD>0 
    219   if (slow_theta) { // Theta is not in inner loop 
    220     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    221   } 
    222208  while(i0 < n0) { 
    223209    local_values.vector[p0] = v0[i0]; 
    224210    double weight0 = w0[i0] * weight1; 
    225211//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
    226     if (fast_theta) { // Theta is in inner loop 
    227       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    228     } 
    229212#else 
    230213    const double weight0 = 1.0; 
    … …  
    232215 
    233216//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); } 
    234 //if (q_index == 0) printf("sphcor: %g\n", spherical_correction); 
    235217 
    236218    #ifdef INVALID 
    … …  
    241223      // Note: weight==0 must always be excluded 
    242224      if (weight0 > cutoff) { 
    243         // spherical correction is set at a minimum of 1e-6, otherwise there 
    244         // would be problems looking at models with theta=90. 
    245         const double weight = weight0 * spherical_correction; 
    246         pd_norm += weight * CALL_VOLUME(local_values.table); 
     225        pd_norm += weight0 * CALL_VOLUME(local_values.table); 
    247226 
    248227#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    … …  
    296275        const double scattering = CALL_IQ(q, q_index, local_values.table); 
    297276#endif // !MAGNETIC 
    298         this_result += weight * scattering; 
     277        this_result += weight0 * scattering; 
    299278      } 
    300279    } 

  • sasmodels/kernelpy.py

    rdbfd471 rd4c33d6  
    197197 
    198198    pd_norm = 0.0 
    199     spherical_correction = 1.0 
    200199    partial_weight = np.NaN 
    201200    weight = np.NaN 
    202201 
    203202    p0_par = call_details.pd_par[0] 
    204     p0_is_theta = (p0_par == call_details.theta_par) 
    205203    p0_length = call_details.pd_length[0] 
    206204    p0_index = p0_length 
    … …  
    219217            parameters[pd_par] = pd_value[pd_offset+pd_index] 
    220218            partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    221             if call_details.theta_par >= 0: 
    222                 cor = sin(pi / 180 * parameters[call_details.theta_par]) 
    223                 spherical_correction = max(abs(cor), 1e-6) 
    224219            p0_index = loop_index%p0_length 
    225220 
    226221        weight = partial_weight * pd_weight[p0_offset + p0_index] 
    227222        parameters[p0_par] = pd_value[p0_offset + p0_index] 
    228         if p0_is_theta: 
    229             cor = cos(pi/180 * parameters[p0_par]) 
    230             spherical_correction = max(abs(cor), 1e-6) 
    231223        p0_index += 1 
    232224        if weight > cutoff: 
    … …  
    239231 
    240232            # update value and norm 
    241             weight *= spherical_correction 
    242233            total += weight * Iq 
    243234            pd_norm += weight * form_volume() 

  • sasmodels/model_test.py

    rbb4b509 r73cbc5b  
    7979    suite = unittest.TestSuite() 
    8080 
    81     if models[0] == 'all': 
     81    if models[0] in core.KINDS: 
    8282        skip = models[1:] 
    83         models = list_models() 
     83        models = list_models(models[0]) 
    8484    else: 
    8585        skip = [] 

  • sasmodels/models/barbell.c

    r592343f r2a0b2b1  
    1010//barbell kernel - same as dumbell 
    1111static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     12_bell_kernel(double qab, double qc, double h, double radius_bell, 
     13             double half_length) 
    1414{ 
    1515    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626    //    m = q R cos(alpha) 
    2727    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     28    const double m = radius_bell*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    3232    for (int i = 0; i < 76; i++){ 
    3333        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
    … …  
    4444 
    4545static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, 
     47    double radius_bell, double radius, double half_length) 
    4948{ 
    50     const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 
    51     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    52     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     49    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     50    const double bj = sas_2J1x_x(radius*qab); 
     51    const double si = sas_sinx_x(half_length*qc); 
    5352    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5453    const double Aq = bell_fq + cyl_fq; 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     85        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8786        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8887    } 
    … …  
    103102    double q, sin_alpha, cos_alpha; 
    104103    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     104    const double qab = q*sin_alpha; 
     105    const double qc = q*cos_alpha; 
    105106 
    106107    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     108    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    108109 
    109110    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     5    const double a1 = (+qa + qb + qc)/2.0; 
     6    const double a2 = (-qa - qb + qc)/2.0; 
     7    const double a3 = (-qa + qb - qc)/2.0; 
     8#else 
     9    const double a1 = (+qa + qb - qc)/2.0; 
     10    const double a2 = (+qa - qb + qc)/2.0; 
     11    const double a3 = (-qa + qb + qc)/2.0; 
     12#endif 
    613 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     14#if 1 
     15    // Numerator: (1 - exp(a)^2)^3 
     16    //         => (-(exp(2a) - 1))^3 
     17    //         => -expm1(2a)^3 
     18    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     19    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     20    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     21    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     22    const double exp_arg = exp(arg); 
     23    const double Zq = -cube(expm1(2.0*arg)) 
     24        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     26          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     27#else 
     28    // Alternate form, which perhaps is more approachable 
     29    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     30    const double sinh_qd = sinh(arg); 
     31    const double cosh_qd = cosh(arg); 
     32    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
     33                    * sinh_qd/(cosh_qd - cos(dnn*a2)) 
     34                    * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     35#endif 
     36 
     37    return Zq; 
     38} 
    1039 
    1140 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     41// occupied volume fraction calculated from lattice symmetry and sphere radius 
     42static double 
     43bcc_volume_fraction(double radius, double dnn) 
     44{ 
     45    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2046} 
    2147 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     48static double 
     49form_volume(double radius) 
     50{ 
    4451    return sphere_volume(radius); 
    4552} 
    4653 
    4754 
    48 double Iq(double q, double dnn, 
     55static double Iq(double q, double dnn, 
    4956  double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     57  double sld, double solvent_sld) 
     58{ 
     59    // translate a point in [-1,1] to a point in [0, 2 pi] 
     60    const double phi_m = M_PI; 
     61    const double phi_b = M_PI; 
     62    // translate a point in [-1,1] to a point in [0, pi] 
     63    const double theta_m = M_PI_2; 
     64    const double theta_b = M_PI_2; 
    5165 
    52         //Volume fraction calculated from lattice symmetry and sphere radius 
    53         const double s1 = dnn/sqrt(0.75); 
    54         const double latticescale = 2.0*sphere_volume(radius/s1); 
    55  
    56     const double va = 0.0; 
    57     const double vb = 2.0*M_PI; 
    58     const double vaj = 0.0; 
    59     const double vbj = M_PI; 
    60  
    61     double summ = 0.0; 
    62     double answer = 0.0; 
    63         for(int i=0; i<150; i++) { 
    64                 //setup inner integral over the ellipsoidal cross-section 
    65                 double summj=0.0; 
    66                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    67                 for(int j=0;j<150;j++) { 
    68                         //20 gauss points for the inner integral 
    69                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    70                         double yyy = Gauss150Wt[j] * _BCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    71                         summj += yyy; 
    72                 } 
    73                 //now calculate the value of the inner integral 
    74                 double answer = (vbj-vaj)/2.0*summj; 
    75  
    76                 //now calculate outer integral 
    77                 summ = summ+(Gauss150Wt[i] * answer); 
    78         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    79  
    80         answer = (vb-va)/2.0*summ; 
    81         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    82  
    83     return answer; 
     66    double outer_sum = 0.0; 
     67    for(int i=0; i<150; i++) { 
     68        double inner_sum = 0.0; 
     69        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     70        double sin_theta, cos_theta; 
     71        SINCOS(theta, sin_theta, cos_theta); 
     72        const double qc = q*cos_theta; 
     73        const double qab = q*sin_theta; 
     74        for(int j=0;j<150;j++) { 
     75            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     76            double sin_phi, cos_phi; 
     77            SINCOS(phi, sin_phi, cos_phi); 
     78            const double qa = qab*cos_phi; 
     79            const double qb = qab*sin_phi; 
     80            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     81            inner_sum += Gauss150Wt[j] * form; 
     82        } 
     83        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     84        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     85    } 
     86    outer_sum *= theta_m; 
     87    const double Zq = outer_sum/(4.0*M_PI); 
     88    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     89    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    8490} 
    8591 
    8692 
    87 double Iqxy(double qx, double qy, 
     93static double Iqxy(double qx, double qy, 
    8894    double dnn, double d_factor, double radius, 
    8995    double sld, double solvent_sld, 
    … …  
    9298    double q, zhat, yhat, xhat; 
    9399    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     100    const double qa = q*xhat; 
     101    const double qb = q*yhat; 
     102    const double qc = q*zhat; 
    94103 
    95     const double a1 = +xhat - zhat + yhat; 
    96     const double a2 = +xhat + zhat - yhat; 
    97     const double a3 = -xhat + zhat + yhat; 
    98  
    99     const double qd = 0.5*q*dnn; 
    100     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    101     const double tanh_qd = tanh(arg); 
    102     const double cosh_qd = cosh(arg); 
    103     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    105                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    106  
    107     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    108     //the occupied volume of the lattice 
    109     const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    110     return lattice_scale * Fq; 
     104    q = sqrt(qa*qa + qb*qb + qc*qc); 
     105    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     106    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     107    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    111108} 

  • sasmodels/models/bcc_paracrystal.py

    r69e1afc r2a0b2b1  
    8181.. figure:: img/parallelepiped_angle_definition.png 
    8282 
    83     Orientation of the crystal with respect to the scattering plane, when  
     83    Orientation of the crystal with respect to the scattering plane, when 
    8484    $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 
    8585 

  • sasmodels/models/capped_cylinder.c

    r592343f r2a0b2b1  
    1414//   radius_cap is the radius of the lens 
    1515//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     16//   theta is the angle of the cylinder wrt q. 
    1717static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     18_cap_kernel(double qab, double qc, double h, double radius_cap, 
     19            double half_length) 
    2020{ 
    2121    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2626 
    2727    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     28    //    cos (q (R t - h + L/2) cos(theta)) 
    2929    // so turn it into: 
    3030    //    cos (m t + b) 
    3131    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     32    //    m = q R cos(theta) 
     33    //    b = q(L/2-h) cos(theta) 
     34    const double m = radius_cap*qc; // cos argument slope 
     35    const double b = (half_length-h)*qc; // cos argument intercept 
     36    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3737    double total = 0.0; 
    3838    for (int i=0; i<76 ;i++) { 
    3939        const double t = Gauss76Z[i]*zm + zb; 
    4040        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
    … …  
    5050 
    5151static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     52_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5453{ 
    55     const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 
    56     const double bj = sas_2J1x_x(q*radius*sin_alpha); 
    57     const double si = sas_sinx_x(q*half_length*cos_alpha); 
     54    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     55    const double bj = sas_2J1x_x(radius*qab); 
     56    const double si = sas_sinx_x(half_length*qc); 
    5857    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5958    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    101100    double total = 0.0; 
    102101    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    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; 
     102        const double theta = Gauss76Z[i]*zm + zb; 
     103        double sin_theta, cos_theta; // slots to hold sincos function output 
     104        SINCOS(theta, sin_theta, cos_theta); 
     105        const double qab = q*sin_theta; 
     106        const double qc = q*cos_theta; 
     107        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     108        // scale by sin_theta for spherical coord integration 
     109        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110110    } 
    111111    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    125125    double q, sin_alpha, cos_alpha; 
    126126    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     127    const double qab = q*sin_alpha; 
     128    const double qc = q*cos_alpha; 
    127129 
    128130    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     131    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130132 
    131133    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    rb260926 r2a0b2b1  
    1 double form_volume(double radius, double thick_rim, double thick_face, double length); 
    2 double Iq(double q, 
    3           double radius, 
    4           double thick_rim, 
    5           double thick_face, 
    6           double length, 
    7           double core_sld, 
    8           double face_sld, 
    9           double rim_sld, 
    10           double solvent_sld); 
    11  
    12  
    13 double Iqxy(double qx, double qy, 
    14           double radius, 
    15           double thick_rim, 
    16           double thick_face, 
    17           double length, 
    18           double core_sld, 
    19           double face_sld, 
    20           double rim_sld, 
    21           double solvent_sld, 
    22           double theta, 
    23           double phi); 
    24  
    25  
    26 double form_volume(double radius, double thick_rim, double thick_face, double length) 
     1static double 
     2form_volume(double radius, double thick_rim, double thick_face, double length) 
    273{ 
    28     return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
     4    return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
    295} 
    306 
    317static double 
    32 bicelle_kernel(double q, 
    33               double rad, 
    34               double radthick, 
    35               double facthick, 
    36               double halflength, 
    37               double rhoc, 
    38               double rhoh, 
    39               double rhor, 
    40               double rhosolv, 
    41               double sin_alpha, 
    42               double cos_alpha) 
     8bicelle_kernel(double qab, 
     9    double qc, 
     10    double radius, 
     11    double thick_radius, 
     12    double thick_face, 
     13    double halflength, 
     14    double sld_core, 
     15    double sld_face, 
     16    double sld_rim, 
     17    double sld_solvent) 
    4318{ 
    44     const double dr1 = rhoc-rhoh; 
    45     const double dr2 = rhor-rhosolv; 
    46     const double dr3 = rhoh-rhor; 
    47     const double vol1 = M_PI*square(rad)*2.0*(halflength); 
    48     const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 
    49     const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 
     19    const double dr1 = sld_core-sld_face; 
     20    const double dr2 = sld_rim-sld_solvent; 
     21    const double dr3 = sld_face-sld_rim; 
     22    const double vol1 = M_PI*square(radius)*2.0*(halflength); 
     23    const double vol2 = M_PI*square(radius+thick_radius)*2.0*(halflength+thick_face); 
     24    const double vol3 = M_PI*square(radius)*2.0*(halflength+thick_face); 
    5025 
    51     const double be1 = sas_2J1x_x(q*(rad)*sin_alpha); 
    52     const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha); 
    53     const double si1 = sas_sinx_x(q*(halflength)*cos_alpha); 
    54     const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha); 
     26    const double be1 = sas_2J1x_x((radius)*qab); 
     27    const double be2 = sas_2J1x_x((radius+thick_radius)*qab); 
     28    const double si1 = sas_sinx_x((halflength)*qc); 
     29    const double si2 = sas_sinx_x((halflength+thick_face)*qc); 
    5530 
    5631    const double t = vol1*dr1*si1*be1 + 
    … …  
    5833                     vol3*dr3*si2*be1; 
    5934 
    60     const double retval = t*t; 
    61  
    62     return retval; 
    63  
     35    return t; 
    6436} 
    6537 
    6638static double 
    67 bicelle_integration(double q, 
    68                    double rad, 
    69                    double radthick, 
    70                    double facthick, 
    71                    double length, 
    72                    double rhoc, 
    73                    double rhoh, 
    74                    double rhor, 
    75                    double rhosolv) 
     39Iq(double q, 
     40    double radius, 
     41    double thick_radius, 
     42    double thick_face, 
     43    double length, 
     44    double sld_core, 
     45    double sld_face, 
     46    double sld_rim, 
     47    double sld_solvent) 
    7648{ 
    7749    // set up the integration end points 
    … …  
    7951    const double halflength = 0.5*length; 
    8052 
    81     double summ = 0.0; 
     53    double total = 0.0; 
    8254    for(int i=0;i<N_POINTS_76;i++) { 
    83         double alpha = (Gauss76Z[i] + 1.0)*uplim; 
    84         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    85         SINCOS(alpha, sin_alpha, cos_alpha); 
    86         double yyy = Gauss76Wt[i] * bicelle_kernel(q, rad, radthick, facthick, 
    87                              halflength, rhoc, rhoh, rhor, rhosolv, 
    88                              sin_alpha, cos_alpha); 
    89         summ += yyy*sin_alpha; 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
     56        double sin_theta, cos_theta; // slots to hold sincos function output 
     57        SINCOS(theta, sin_theta, cos_theta); 
     58        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     59                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    9061    } 
    9162 
    9263    // calculate value of integral to return 
    93     double answer = uplim*summ; 
    94     return answer; 
     64    double answer = total*uplim; 
     65    return 1.0e-4*answer; 
    9566} 
    9667 
    9768static double 
    98 bicelle_kernel_2d(double qx, double qy, 
    99           double radius, 
    100           double thick_rim, 
    101           double thick_face, 
    102           double length, 
    103           double core_sld, 
    104           double face_sld, 
    105           double rim_sld, 
    106           double solvent_sld, 
    107           double theta, 
    108           double phi) 
     69Iqxy(double qx, double qy, 
     70    double radius, 
     71    double thick_rim, 
     72    double thick_face, 
     73    double length, 
     74    double core_sld, 
     75    double face_sld, 
     76    double rim_sld, 
     77    double solvent_sld, 
     78    double theta, 
     79    double phi) 
    10980{ 
    11081    double q, sin_alpha, cos_alpha; 
    11182    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     83    const double qab = q*sin_alpha; 
     84    const double qc = q*cos_alpha; 
    11285 
    113     double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
     86    double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    11487                           0.5*length, core_sld, face_sld, rim_sld, 
    115                            solvent_sld, sin_alpha, cos_alpha); 
    116     return 1.0e-4*answer; 
     88                           solvent_sld); 
     89    return 1.0e-4*fq*fq; 
    11790} 
    118  
    119 double Iq(double q, 
    120           double radius, 
    121           double thick_rim, 
    122           double thick_face, 
    123           double length, 
    124           double core_sld, 
    125           double face_sld, 
    126           double rim_sld, 
    127           double solvent_sld) 
    128 { 
    129     double intensity = bicelle_integration(q, radius, thick_rim, thick_face, 
    130                        length, core_sld, face_sld, rim_sld, solvent_sld); 
    131     return intensity*1.0e-4; 
    132 } 
    133  
    134  
    135 double Iqxy(double qx, double qy, 
    136           double radius, 
    137           double thick_rim, 
    138           double thick_face, 
    139           double length, 
    140           double core_sld, 
    141           double face_sld, 
    142           double rim_sld, 
    143           double solvent_sld, 
    144           double theta, 
    145           double phi) 
    146 { 
    147     double intensity = bicelle_kernel_2d(qx, qy, 
    148                       radius, 
    149                       thick_rim, 
    150                       thick_face, 
    151                       length, 
    152                       core_sld, 
    153                       face_sld, 
    154                       rim_sld, 
    155                       solvent_sld, 
    156                       theta, 
    157                       phi); 
    158  
    159     return intensity; 
    160 } 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    rdedcf34 r2a0b2b1  
    1717        double thick_face, 
    1818        double length, 
    19         double rhoc, 
    20         double rhoh, 
    21         double rhor, 
    22         double rhosolv) 
     19        double sld_core, 
     20        double sld_face, 
     21        double sld_rim, 
     22        double sld_solvent) 
    2323{ 
    24     double si1,si2,be1,be2; 
    2524     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2625     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
    27      //    const double uplim = M_PI_4; 
    2826    const double halfheight = 0.5*length; 
    29     //const double va = 0.0; 
    30     //const double vb = 1.0; 
    31     // inner integral limits 
    32     //const double vaj=0.0; 
    33     //const double vbj=M_PI; 
    34  
    3527    const double r_major = r_minor * x_core; 
    3628    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3729    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    38     const double dr1 = (rhoc-rhoh)   *M_PI*r_minor*r_major*(2.0*halfheight);; 
    39     const double dr2 = (rhor-rhosolv)*M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    40     const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    41     //const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
    42     //const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
    43     //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     30    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     31    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face); 
     32    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
     33    const double dr1 = vol1*(sld_core-sld_face); 
     34    const double dr2 = vol2*(sld_rim-sld_solvent); 
     35    const double dr3 = vol3*(sld_face-sld_rim); 
    4436 
    4537    //initialize integral 
    … …  
    4739    for(int i=0;i<76;i++) { 
    4840        //setup inner integral over the ellipsoidal cross-section 
    49         // since we generate these lots of times, why not store them somewhere? 
    50         //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    51         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    52         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    53         double inner_sum=0; 
    54         double sinarg1 = q*halfheight*cos_alpha; 
    55         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    56         si1 = sas_sinx_x(sinarg1); 
    57         si2 = sas_sinx_x(sinarg2); 
     41        //const double va = 0.0; 
     42        //const double vb = 1.0; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     45        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     46        const double qab = q*sin_theta; 
     47        const double qc = q*cos_theta; 
     48        const double si1 = sas_sinx_x(halfheight*qc); 
     49        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     50        double inner_sum=0.0; 
    5851        for(int j=0;j<76;j++) { 
    5952            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    60             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    61             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    62             const double rr = sqrt(r2A - r2B*cos(beta)); 
    63             double besarg1 = q*rr*sin_alpha; 
    64             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    65             be1 = sas_2J1x_x(besarg1); 
    66             be2 = sas_2J1x_x(besarg2); 
    67             inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    68                                               dr2*si2*be2 + 
    69                                               dr3*si2*be1); 
     53            // inner integral limits 
     54            //const double vaj=0.0; 
     55            //const double vbj=M_PI; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     58            const double rr = sqrt(r2A - r2B*cos(phi)); 
     59            const double be1 = sas_2J1x_x(rr*qab); 
     60            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     61            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     62 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    7064        } 
    7165        //now calculate outer integral 
    … …  
    8377          double thick_face, 
    8478          double length, 
    85           double rhoc, 
    86           double rhoh, 
    87           double rhor, 
    88           double rhosolv, 
     79          double sld_core, 
     80          double sld_face, 
     81          double sld_rim, 
     82          double sld_solvent, 
    8983          double theta, 
    9084          double phi, 
    9185          double psi) 
    9286{ 
    93        // THIS NEEDS TESTING 
    9487    double q, xhat, yhat, zhat; 
    9588    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    96     const double dr1 = rhoc-rhoh; 
    97     const double dr2 = rhor-rhosolv; 
    98     const double dr3 = rhoh-rhor; 
     89    const double qa = q*xhat; 
     90    const double qb = q*yhat; 
     91    const double qc = q*zhat; 
     92 
     93    const double dr1 = sld_core-sld_face; 
     94    const double dr2 = sld_rim-sld_solvent; 
     95    const double dr3 = sld_face-sld_rim; 
    9996    const double r_major = r_minor*x_core; 
    10097    const double halfheight = 0.5*length; 
    … …  
    104101 
    105102    // Compute effective radius in rotated coordinates 
    106     const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat)); 
    107     const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat) 
    108                                    + square((r_minor+thick_rim)*yhat)); 
    109     const double be1 = sas_2J1x_x( q*r_hat ); 
    110     const double be2 = sas_2J1x_x( q*rshell_hat ); 
    111     const double si1 = sas_sinx_x( q*halfheight*zhat ); 
    112     const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq; 
     103    const double qr_hat = sqrt(square(r_major*qa) + square(r_minor*qb)); 
     104    const double qrshell_hat = sqrt(square((r_major+thick_rim)*qa) 
     105                                   + square((r_minor+thick_rim)*qb)); 
     106    const double be1 = sas_2J1x_x( qr_hat ); 
     107    const double be2 = sas_2J1x_x( qrshell_hat ); 
     108    const double si1 = sas_sinx_x( halfheight*qc ); 
     109    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
     110    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1; 
     111    return 1.0e-4 * fq*fq; 
    115112} 
    116  

  • sasmodels/models/core_shell_cylinder.c

    r592343f r2a0b2b1  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
    5     double radius, double thickness, double length, double theta, double phi); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
    … …  
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // translate a point in [-1,1] to a point in [0, pi/2] 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6460    double q, sin_alpha, cos_alpha; 
    6561    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     62    const double qab = q*sin_alpha; 
     63    const double qc = q*cos_alpha; 
    6664 
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     65    const double core_r = radius; 
     66    const double core_h = 0.5*length; 
    6967    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     68    const double shell_r = (radius + thickness); 
     69    const double shell_h = (0.5*length + thickness); 
    7270    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7371 
    74     const double fq = _cyl(core_vd, core_qr*sin_alpha, core_qh*cos_alpha) 
    75         + _cyl(shell_vd, shell_qr*sin_alpha, shell_qh*cos_alpha); 
     72    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     73        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7674    return 1.0e-4 * fq * fq; 
    7775} 

  • sasmodels/models/core_shell_ellipsoid.c

    r0a3d9b2 r2a0b2b1  
    1 double form_volume(double radius_equat_core, 
    2                    double polar_core, 
    3                    double equat_shell, 
    4                    double polar_shell); 
    5 double Iq(double q, 
    6           double radius_equat_core, 
    7           double x_core, 
    8           double thick_shell, 
    9           double x_polar_shell, 
    10           double core_sld, 
    11           double shell_sld, 
    12           double solvent_sld); 
    131 
     2// Converted from Igor function gfn4, using the same pattern as ellipsoid 
     3// for evaluating the parts of the integral. 
     4//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN 
     5//                       BY (53) & (58-59) IN CHEN AND 
     6//                       KOTLARCHYK REFERENCE 
     7// 
     8//       <OBLATE ELLIPSOID> 
     9static double 
     10_cs_ellipsoid_kernel(double qab, double qc, 
     11    double equat_core, double polar_core, 
     12    double equat_shell, double polar_shell, 
     13    double sld_core_shell, double sld_shell_solvent) 
     14{ 
     15    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc)); 
     16    const double si_core = sas_3j1x_x(qr_core); 
     17    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core; 
     18    const double fq_core = si_core*volume_core*sld_core_shell; 
    1419 
    15 double Iqxy(double qx, double qy, 
    16           double radius_equat_core, 
    17           double x_core, 
    18           double thick_shell, 
    19           double x_polar_shell, 
    20           double core_sld, 
    21           double shell_sld, 
    22           double solvent_sld, 
    23           double theta, 
    24           double phi); 
     20    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc)); 
     21    const double si_shell = sas_3j1x_x(qr_shell); 
     22    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell; 
     23    const double fq_shell = si_shell*volume_shell*sld_shell_solvent; 
    2524 
     25    return fq_core + fq_shell; 
     26} 
    2627 
    27 double form_volume(double radius_equat_core, 
    28                    double x_core, 
    29                    double thick_shell, 
    30                    double x_polar_shell) 
     28static double 
     29form_volume(double radius_equat_core, 
     30    double x_core, 
     31    double thick_shell, 
     32    double x_polar_shell) 
    3133{ 
    3234    const double equat_shell = radius_equat_core + thick_shell; 
    … …  
    3739 
    3840static double 
    39 core_shell_ellipsoid_xt_kernel(double q, 
    40           double radius_equat_core, 
    41           double x_core, 
    42           double thick_shell, 
    43           double x_polar_shell, 
    44           double core_sld, 
    45           double shell_sld, 
    46           double solvent_sld) 
     41Iq(double q, 
     42    double radius_equat_core, 
     43    double x_core, 
     44    double thick_shell, 
     45    double x_polar_shell, 
     46    double core_sld, 
     47    double shell_sld, 
     48    double solvent_sld) 
    4749{ 
    48     const double lolim = 0.0; 
    49     const double uplim = 1.0; 
    50  
    51  
    52     const double delpc = core_sld - shell_sld; //core - shell 
    53     const double delps = shell_sld - solvent_sld; //shell - solvent 
    54  
     50    const double sld_core_shell = core_sld - shell_sld; 
     51    const double sld_shell_solvent = shell_sld - solvent_sld; 
    5552 
    5653    const double polar_core = radius_equat_core*x_core; 
    … …  
    5855    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5956 
    60     double summ = 0.0;   //initialize intergral 
     57    // translate from [-1, 1] => [0, 1] 
     58    const double m = 0.5; 
     59    const double b = 0.5; 
     60    double total = 0.0;     //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    62         double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim ); 
    63         double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell, 
    64                           polar_shell, delpc, delps, q); 
    65         summ += Gauss76Wt[i] * yyy; 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
     63        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     64        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
     65            radius_equat_core, polar_core, 
     66            equat_shell, polar_shell, 
     67            sld_core_shell, sld_shell_solvent); 
     68        total += Gauss76Wt[i] * fq * fq; 
    6669    } 
    67     summ *= 0.5*(uplim-lolim); 
     70    total *= m; 
    6871 
    6972    // convert to [cm-1] 
    70     return 1.0e-4 * summ; 
     73    return 1.0e-4 * total; 
    7174} 
    7275 
    7376static double 
    74 core_shell_ellipsoid_xt_kernel_2d(double qx, double qy, 
    75           double radius_equat_core, 
    76           double x_core, 
    77           double thick_shell, 
    78           double x_polar_shell, 
    79           double core_sld, 
    80           double shell_sld, 
    81           double solvent_sld, 
    82           double theta, 
    83           double phi) 
     77Iqxy(double qx, double qy, 
     78    double radius_equat_core, 
     79    double x_core, 
     80    double thick_shell, 
     81    double x_polar_shell, 
     82    double core_sld, 
     83    double shell_sld, 
     84    double solvent_sld, 
     85    double theta, 
     86    double phi) 
    8487{ 
    8588    double q, sin_alpha, cos_alpha; 
    8689    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     90    const double qab = q*sin_alpha; 
     91    const double qc = q*cos_alpha; 
    8792 
    88     const double sldcs = core_sld - shell_sld; 
    89     const double sldss = shell_sld- solvent_sld; 
     93    const double sld_core_shell = core_sld - shell_sld; 
     94    const double sld_shell_solvent = shell_sld - solvent_sld; 
    9095 
    9196    const double polar_core = radius_equat_core*x_core; 
    … …  
    9398    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9499 
    95     // Call the IGOR library function to get the kernel: 
    96     // MUST use gfn4 not gf2 because of the def of params. 
    97     double answer = gfn4(cos_alpha, 
    98                   radius_equat_core, 
    99                   polar_core, 
    100                   equat_shell, 
    101                   polar_shell, 
    102                   sldcs, 
    103                   sldss, 
    104                   q); 
     100    double fq = _cs_ellipsoid_kernel(qab, qc, 
     101                  radius_equat_core, polar_core, 
     102                  equat_shell, polar_shell, 
     103                  sld_core_shell, sld_shell_solvent); 
    105104 
    106105    //convert to [cm-1] 
    107     answer *= 1.0e-4; 
    108  
    109     return answer; 
     106    return 1.0e-4 * fq * fq; 
    110107} 
    111  
    112 double Iq(double q, 
    113           double radius_equat_core, 
    114           double x_core, 
    115           double thick_shell, 
    116           double x_polar_shell, 
    117           double core_sld, 
    118           double shell_sld, 
    119           double solvent_sld) 
    120 { 
    121     double intensity = core_shell_ellipsoid_xt_kernel(q, 
    122            radius_equat_core, 
    123            x_core, 
    124            thick_shell, 
    125            x_polar_shell, 
    126            core_sld, 
    127            shell_sld, 
    128            solvent_sld); 
    129  
    130     return intensity; 
    131 } 
    132  
    133  
    134 double Iqxy(double qx, double qy, 
    135           double radius_equat_core, 
    136           double x_core, 
    137           double thick_shell, 
    138           double x_polar_shell, 
    139           double core_sld, 
    140           double shell_sld, 
    141           double solvent_sld, 
    142           double theta, 
    143           double phi) 
    144 { 
    145     double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy, 
    146                        radius_equat_core, 
    147                        x_core, 
    148                        thick_shell, 
    149                        x_polar_shell, 
    150                        core_sld, 
    151                        shell_sld, 
    152                        solvent_sld, 
    153                        theta, 
    154                        phi); 
    155  
    156     return intensity; 
    157 } 

  • sasmodels/models/core_shell_ellipsoid.py

    r9802ab3 r2a0b2b1  
    2525ellipsoid. This may have some undesirable effects if the aspect ratio of the 
    2626ellipsoid is large (ie, if $X << 1$ or $X >> 1$ ), when the $S(q)$ 
    27 - which assumes spheres - will not in any case be valid.  Generating a  
    28 custom product model will enable separate effective volume fraction and effective  
     27- which assumes spheres - will not in any case be valid.  Generating a 
     28custom product model will enable separate effective volume fraction and effective 
    2929radius in the $S(q)$. 
    3030 
    … …  
    4444 
    4545.. math:: 
    46     \begin{align}     
     46    \begin{align} 
    4747    F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ 
    4848    &+ f(q,radius\_equat\_core + thick\_shell,radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) 
    49     \end{align}  
     49    \end{align} 
    5050 
    5151where 
    52   
     52 
    5353.. math:: 
    5454 
    … …  
    7777   F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 
    7878 
    79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,  
     79For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 
    8080see the :ref:`elliptical-cylinder` model for further information. 
    8181 
    … …  
    139139# pylint: enable=bad-whitespace, line-too-long 
    140140 
    141 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    142           "core_shell_ellipsoid.c"] 
     141source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    143142 
    144143def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c r2a0b2b1  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    … …  
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
    … …  
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36 
    3737    const double mu = 0.5 * q * length_b; 
    38      
     38 
    3939    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     40    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     41    //double vol = length_a * length_b * length_c + 
     42    //       2.0 * thick_rim_a * length_b * length_c + 
    4343    //       2.0 * thick_rim_b * length_a * length_c + 
    4444    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     45 
    4646    // Scale sides by B 
    4747    const double a_scaled = length_a / length_b; 
    … …  
    101101            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102102            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     103            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     104 
    105105            //  correct FF : sum of square of phase factors 
    106106            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    136136    double q, zhat, yhat, xhat; 
    137137    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     138    const double qa = q*xhat; 
     139    const double qb = q*yhat; 
     140    const double qc = q*zhat; 
    138141 
    139142    // cspkernel in csparallelepiped recoded here 
    … …  
    160163    double tc = length_a + 2.0*thick_rim_c; 
    161164    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     165    double siA = sas_sinx_x(0.5*length_a*qa); 
     166    double siB = sas_sinx_x(0.5*length_b*qb); 
     167    double siC = sas_sinx_x(0.5*length_c*qc); 
     168    double siAt = sas_sinx_x(0.5*ta*qa); 
     169    double siBt = sas_sinx_x(0.5*tb*qb); 
     170    double siCt = sas_sinx_x(0.5*tc*qc); 
     171 
    169172 
    170173    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    174177               + drB*siA*(siBt-siB)*siC*V2 
    175178               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     179 
    177180    return 1.0e-4 * f * f; 
    178181} 

  • sasmodels/models/cylinder.c

    r592343f rb34fc77  
    1 double form_volume(double radius, double length); 
    2 double fq(double q, double sn, double cn,double radius, double length); 
    3 double orient_avg_1D(double q, double radius, double length); 
    4 double Iq(double q, double sld, double solvent_sld, double radius, double length); 
    5 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    6     double radius, double length, double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius<0 || v.length<0) 
    92 
    10 double form_volume(double radius, double length) 
     3static double 
     4form_volume(double radius, double length) 
    115{ 
    126    return M_PI*radius*radius*length; 
    137} 
    148 
    15 double fq(double q, double sn, double cn, double radius, double length) 
     9static double 
     10fq(double qab, double qc, double radius, double length) 
    1611{ 
    17     // precompute qr and qh to save time in the loop 
    18     const double qr = q*radius; 
    19     const double qh = q*0.5*length;  
    20     return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 
     12    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
    2113} 
    2214 
    23 double orient_avg_1D(double q, double radius, double length) 
     15static double 
     16orient_avg_1D(double q, double radius, double length) 
    2417{ 
    2518    // translate a point in [-1,1] to a point in [0, pi/2] 
    2619    const double zm = M_PI_4; 
    27     const double zb = M_PI_4;  
     20    const double zb = M_PI_4; 
    2821 
    2922    double total = 0.0; 
    3023    for (int i=0; i<76 ;i++) { 
    31         const double alpha = Gauss76Z[i]*zm + zb; 
    32         double sn, cn; // slots to hold sincos function output 
    33         // alpha(theta,phi) the projection of the cylinder on the detector plane 
    34         SINCOS(alpha, sn, cn); 
    35         total += Gauss76Wt[i] * square( fq(q, sn, cn, radius, length) ) * sn; 
     24        const double theta = Gauss76Z[i]*zm + zb; 
     25        double sin_theta, cos_theta; // slots to hold sincos function output 
     26        // theta (theta,phi) the projection of the cylinder on the detector plane 
     27        SINCOS(theta , sin_theta, cos_theta); 
     28        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3630    } 
    3731    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3933} 
    4034 
    41 double Iq(double q, 
     35static double 
     36Iq(double q, 
    4237    double sld, 
    4338    double solvent_sld, 
    … …  
    4944} 
    5045 
    51  
    52 double Iqxy(double qx, double qy, 
     46static double 
     47Iqxy(double qx, double qy, 
    5348    double sld, 
    5449    double solvent_sld, 
    … …  
    6055    double q, sin_alpha, cos_alpha; 
    6156    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    62     //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 
     57    const double qab = q*sin_alpha; 
     58    const double qc = q*cos_alpha; 
     59 
    6360    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    64     const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
     61    const double form = fq(qab, qc, radius, length); 
    6562    return 1.0e-4 * square(s * form); 
    6663} 

  • sasmodels/models/ellipsoid.c

    r3b571ae r2a0b2b1  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
    … …  
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qx, double qy, 
    4442    double sld, 
    4543    double sld_solvent, 
    … …  
    5149    double q, sin_alpha, cos_alpha; 
    5250    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     51    const double qab = q*sin_alpha; 
     52    const double qc = q*cos_alpha; 
     53 
     54    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     55    const double f = sas_3j1x_x(qr); 
    5656    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5757 
    5858    return 1.0e-4 * square(f * s); 
    5959} 
    60  

  • sasmodels/models/elliptical_cylinder.c

    r61104c8 r2a0b2b1  
    1 double form_volume(double radius_minor, double r_ratio, double length); 
    2 double Iq(double q, double radius_minor, double r_ratio, double length, 
    3           double sld, double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius_minor, double r_ratio, double length, 
    5             double sld, double solvent_sld, double theta, double phi, double psi); 
    6  
    7  
    8 double 
     1static double 
    92form_volume(double radius_minor, double r_ratio, double length) 
    103{ 
    … …  
    125} 
    136 
    14 double 
     7static double 
    158Iq(double q, double radius_minor, double r_ratio, double length, 
    169   double sld, double solvent_sld) 
    … …  
    6154 
    6255 
    63 double 
     56static double 
    6457Iqxy(double qx, double qy, 
    6558     double radius_minor, double r_ratio, double length, 
    … …  
    6962    double q, xhat, yhat, zhat; 
    7063    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     64    const double qa = q*xhat; 
     65    const double qb = q*yhat; 
     66    const double qc = q*zhat; 
    7167 
    7268    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    7369    // Given:    radius_major = r_ratio * radius_minor 
    74     const double r = radius_minor*sqrt(square(r_ratio*xhat) + square(yhat)); 
    75     const double be = sas_2J1x_x(q*r); 
    76     const double si = sas_sinx_x(q*zhat*0.5*length); 
     70    const double qr = radius_minor*sqrt(square(r_ratio*qa) + square(qb)); 
     71    const double be = sas_2J1x_x(qr); 
     72    const double si = sas_sinx_x(qc*0.5*length); 
    7773    const double Aq = be * si; 
    7874    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4#if 0  // Equations as written in Matsuoka 
     5    const double a1 = ( qa + qb)/2.0; 
     6    const double a2 = (-qa + qc)/2.0; 
     7    const double a3 = (-qa + qb)/2.0; 
     8#else 
     9    const double a1 = ( qa + qb)/2.0; 
     10    const double a2 = ( qa + qc)/2.0; 
     11    const double a3 = ( qb + qc)/2.0; 
     12#endif 
    613 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     14    // Numerator: (1 - exp(a)^2)^3 
     15    //         => (-(exp(2a) - 1))^3 
     16    //         => -expm1(2a)^3 
     17    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     18    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     19    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     20    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     21    const double exp_arg = exp(arg); 
     22    const double Zq = -cube(expm1(2.0*arg)) 
     23        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     25          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     26 
     27    return Zq; 
     28} 
    929 
    1030 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     31// occupied volume fraction calculated from lattice symmetry and sphere radius 
     32static double 
     33fcc_volume_fraction(double radius, double dnn) 
     34{ 
     35    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1936} 
    2037 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     38static double 
     39form_volume(double radius) 
     40{ 
    4341    return sphere_volume(radius); 
    4442} 
    4543 
    4644 
    47 double Iq(double q, double dnn, 
     45static double Iq(double q, double dnn, 
    4846  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     47  double sld, double solvent_sld) 
     48{ 
     49    // translate a point in [-1,1] to a point in [0, 2 pi] 
     50    const double phi_m = M_PI; 
     51    const double phi_b = M_PI; 
     52    // translate a point in [-1,1] to a point in [0, pi] 
     53    const double theta_m = M_PI_2; 
     54    const double theta_b = M_PI_2; 
    5055 
    51         //Volume fraction calculated from lattice symmetry and sphere radius 
    52         const double s1 = dnn*sqrt(2.0); 
    53         const double latticescale = 4.0*sphere_volume(radius/s1); 
     56    double outer_sum = 0.0; 
     57    for(int i=0; i<150; i++) { 
     58        double inner_sum = 0.0; 
     59        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     60        double sin_theta, cos_theta; 
     61        SINCOS(theta, sin_theta, cos_theta); 
     62        const double qc = q*cos_theta; 
     63        const double qab = q*sin_theta; 
     64        for(int j=0;j<150;j++) { 
     65            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     66            double sin_phi, cos_phi; 
     67            SINCOS(phi, sin_phi, cos_phi); 
     68            const double qa = qab*cos_phi; 
     69            const double qb = qab*sin_phi; 
     70            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     71            inner_sum += Gauss150Wt[j] * form; 
     72        } 
     73        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     74        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     75    } 
     76    outer_sum *= theta_m; 
     77    const double Zq = outer_sum/(4.0*M_PI); 
     78    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5479 
    55     const double va = 0.0; 
    56     const double vb = 2.0*M_PI; 
    57     const double vaj = 0.0; 
    58     const double vbj = M_PI; 
    59  
    60     double summ = 0.0; 
    61     double answer = 0.0; 
    62         for(int i=0; i<150; i++) { 
    63                 //setup inner integral over the ellipsoidal cross-section 
    64                 double summj=0.0; 
    65                 const double zphi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;             //the outer dummy is phi 
    66                 for(int j=0;j<150;j++) { 
    67                         //20 gauss points for the inner integral 
    68                         double ztheta = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;             //the inner dummy is theta 
    69                         double yyy = Gauss150Wt[j] * _FCC_Integrand(q,dnn,d_factor,ztheta,zphi); 
    70                         summj += yyy; 
    71                 } 
    72                 //now calculate the value of the inner integral 
    73                 double answer = (vbj-vaj)/2.0*summj; 
    74  
    75                 //now calculate outer integral 
    76                 summ = summ+(Gauss150Wt[i] * answer); 
    77         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    78  
    79         answer = (vb-va)/2.0*summ; 
    80         answer = answer*sphere_form(q,radius,sld,solvent_sld)*latticescale; 
    81  
    82     return answer; 
     80    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     81} 
    8382 
    8483 
    85 } 
    86  
    87 double Iqxy(double qx, double qy, 
     84static double Iqxy(double qx, double qy, 
    8885    double dnn, double d_factor, double radius, 
    8986    double sld, double solvent_sld, 
    … …  
    9289    double q, zhat, yhat, xhat; 
    9390    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     91    const double qa = q*xhat; 
     92    const double qb = q*yhat; 
     93    const double qc = q*zhat; 
    9494 
    95     const double a1 = yhat + xhat; 
    96     const double a2 = xhat + zhat; 
    97     const double a3 = yhat + zhat; 
    98     const double qd = 0.5*q*dnn; 
    99     const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
    100     const double tanh_qd = tanh(arg); 
    101     const double cosh_qd = cosh(arg); 
    102     const double Zq = tanh_qd/(1. - cos(qd*a1)/cosh_qd) 
    103                     * tanh_qd/(1. - cos(qd*a2)/cosh_qd) 
    104                     * tanh_qd/(1. - cos(qd*a3)/cosh_qd); 
    105  
    106     //if (isnan(Zq)) printf("q:(%g,%g) qd: %g a1: %g a2: %g a3: %g arg: %g\n", qx, qy, qd, a1, a2, a3, arg); 
    107  
    108     const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq; 
    109     //the occupied volume of the lattice 
    110     const double lattice_scale = 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    111     return lattice_scale * Fq; 
     95    q = sqrt(qa*qa + qb*qb + qc*qc); 
     96    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     97    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     98    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
    11299} 

  • sasmodels/models/hollow_cylinder.c

    r592343f r2a0b2b1  
    11double form_volume(double radius, double thickness, double length); 
    22double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
     3    double solvent_sld); 
    44double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
     5    double solvent_sld, double theta, double phi); 
    66 
    77//#define INVALID(v) (v.radius_core >= v.radius) 
    … …  
    1414} 
    1515 
    16  
    1716static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     17_fq(double qab, double qc, 
     18    double radius, double thickness, double length) 
    2019{ 
    21     const double qs = q*sin_val; 
    22     const double lam1 = sas_2J1x_x((radius+thickness)*qs); 
    23     const double lam2 = sas_2J1x_x(radius*qs); 
     20    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     21    const double lam2 = sas_2J1x_x(radius*qab); 
    2422    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     23    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     24    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     25    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     26    const double t2 = sas_sinx_x(0.5*length*qc); 
    2927    return psi*t2; 
    3028} 
    … …  
    4341{ 
    4442    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     43    const double upper = 1.0;        //limits of numerical integral 
    4644 
    47     double summ = 0.0;                  //initialize intergral 
     45    double summ = 0.0;            //initialize intergral 
    4846    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     47        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     48        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     49        const double form = _fq(q*sin_theta, q*cos_theta, 
     50                                radius, thickness, length); 
     51        summ += Gauss76Wt[i] * form * form; 
    5452    } 
    5553 
    … …  
    6664    double q, sin_alpha, cos_alpha; 
    6765    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    68     const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
    69         sin_alpha, cos_alpha); 
     66    const double qab = q*sin_alpha; 
     67    const double qc = q*cos_alpha; 
     68 
     69    const double form = _fq(qab, qc, radius, thickness, length); 
    7070 
    7171    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7373} 
    74  

  • sasmodels/models/parallelepiped.c

    rd605080 r2a0b2b1  
    2020{ 
    2121    const double mu = 0.5 * q * length_b; 
    22      
     22 
    2323    // Scale sides by B 
    2424    const double a_scaled = length_a / length_b; 
    2525    const double c_scaled = length_c / length_b; 
    26          
     26 
    2727    // outer integral (with gauss points), integration limits = 0, 1 
    2828    double outer_total = 0; //initialize integral 
    … …  
    6969    double q, xhat, yhat, zhat; 
    7070    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     71    const double qa = q*xhat; 
     72    const double qb = q*yhat; 
     73    const double qc = q*zhat; 
    7174 
    72     const double siA = sas_sinx_x(0.5*length_a*q*xhat); 
    73     const double siB = sas_sinx_x(0.5*length_b*q*yhat); 
    74     const double siC = sas_sinx_x(0.5*length_c*q*zhat); 
     75    const double siA = sas_sinx_x(0.5*length_a*qa); 
     76    const double siB = sas_sinx_x(0.5*length_b*qb); 
     77    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7578    const double V = form_volume(length_a, length_b, length_c); 
    7679    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe r7e0b281  
    1 double form_volume(double radius); 
     1static double 
     2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    const double a1 = qa; 
     5    const double a2 = qb; 
     6    const double a3 = qc; 
    27 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     8    // Numerator: (1 - exp(a)^2)^3 
     9    //         => (-(exp(2a) - 1))^3 
     10    //         => -expm1(2a)^3 
     11    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a)) 
     12    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     13    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1) 
     14    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     15    const double exp_arg = exp(arg); 
     16    const double Zq = -cube(expm1(2.0*arg)) 
     17        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     18          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     19          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
    920 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     21    return Zq; 
     22} 
    1923 
    20 double form_volume(double radius) 
     24// occupied volume fraction calculated from lattice symmetry and sphere radius 
     25static double 
     26sc_volume_fraction(double radius, double dnn) 
     27{ 
     28    return sphere_volume(radius/dnn); 
     29} 
     30 
     31static double 
     32form_volume(double radius) 
    2133{ 
    2234    return sphere_volume(radius); 
    2335} 
    2436 
    25 static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     37 
     38static double Iq(double q, double dnn, 
     39  double d_factor, double radius, 
     40  double sld, double solvent_sld) 
    2741{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     42    // translate a point in [-1,1] to a point in [0, 2 pi] 
     43    const double phi_m = M_PI_4; 
     44    const double phi_b = M_PI_4; 
     45    // translate a point in [-1,1] to a point in [0, pi] 
     46    const double theta_m = M_PI_4; 
     47    const double theta_b = M_PI_4; 
    3048 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3349 
    34         double temp6 = snt; 
    35         double temp7 = -1.0*temp3*snt*cnp; 
    36         double temp8 = temp3*snt*snp; 
    37         double temp9 = temp3*cnt; 
    38         double result = temp6/((1.0-temp4*cos((temp7))+temp5)* 
    39                                (1.0-temp4*cos((temp8))+temp5)* 
    40                                (1.0-temp4*cos((temp9))+temp5)); 
    41         return (result); 
     50    double outer_sum = 0.0; 
     51    for(int i=0; i<150; i++) { 
     52        double inner_sum = 0.0; 
     53        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     54        double sin_theta, cos_theta; 
     55        SINCOS(theta, sin_theta, cos_theta); 
     56        const double qc = q*cos_theta; 
     57        const double qab = q*sin_theta; 
     58        for(int j=0;j<150;j++) { 
     59            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     60            double sin_phi, cos_phi; 
     61            SINCOS(phi, sin_phi, cos_phi); 
     62            const double qa = qab*cos_phi; 
     63            const double qb = qab*sin_phi; 
     64            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     65            inner_sum += Gauss150Wt[j] * form; 
     66        } 
     67        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     68        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     69    } 
     70    outer_sum *= theta_m; 
     71    const double Zq = outer_sum/M_PI_2; 
     72    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     73 
     74    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    4275} 
    4376 
    44 static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
    46 { 
    47     //Function to calculate integrand values for simple cubic structure 
    4877 
    49         double da = d_factor*dnn; 
    50         double temp1 = qq*qq*da*da; 
    51         double temp2 = cube(-expm1(-temp1)); 
    52         double temp3 = qq*dnn; 
    53         double temp4 = 2.0*exp(-0.5*temp1); 
    54         double temp5 = exp(-1.0*temp1); 
    55  
    56         double integrand = temp2*sc_eval(yy,xx,temp3,temp4,temp5)/M_PI_2; 
    57  
    58         return(integrand); 
    59 } 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
     78static double Iqxy(double qx, double qy, 
     79    double dnn, double d_factor, double radius, 
     80    double sld, double solvent_sld, 
     81    double theta, double phi, double psi) 
    11282{ 
    11383    double q, zhat, yhat, xhat; 
    11484    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     85    const double qa = q*xhat; 
     86    const double qb = q*yhat; 
     87    const double qc = q*zhat; 
    11588 
    116     const double qd = q*dnn; 
    117     const double arg = 0.5*square(qd*d_factor); 
    118     const double tanh_qd = tanh(arg); 
    119     const double cosh_qd = cosh(arg); 
    120     const double Zq = tanh_qd/(1. - cos(qd*zhat)/cosh_qd) 
    121                     * tanh_qd/(1. - cos(qd*yhat)/cosh_qd) 
    122                     * tanh_qd/(1. - cos(qd*xhat)/cosh_qd); 
    123  
    124     const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; 
    125     //the occupied volume of the lattice 
    126     const double lattice_scale = sphere_volume(radius/dnn); 
    127     return lattice_scale * Fq; 
     89    q = sqrt(qa*qa + qb*qb + qc*qc); 
     90    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     91    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
     92    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    12893} 

  • sasmodels/models/sc_paracrystal.py

    r69e1afc r2a0b2b1  
    152152    [{}, 0.001, 10.3048], 
    153153    [{}, 0.215268, 0.00814889], 
    154     [{}, (0.414467), 0.001313289], 
     154    [{}, 0.414467, 0.001313289], 
    155155    [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 
    156156    [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 
    157157    ] 
    158  
    159  

  • sasmodels/models/stacked_disks.c

    r19f996b rb34fc77  
    11static double stacked_disks_kernel( 
    2     double q, 
     2    double qab, 
     3    double qc, 
    34    double halfheight, 
    45    double thick_layer, 
    … …  
    910    double layer_sld, 
    1011    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1312    double d) 
    1413 
    … …  
    2019    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2120 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     21    const double besarg1 = radius*qab; 
     22    //const double besarg2 = radius*qab; 
    2423 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     24    const double sinarg1 = halfheight*qc; 
     25    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2726 
    2827    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4342 
    4443    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     44    double qd_cos_alpha = d*qc; 
    4645    //d*cos_alpha is the projection of d onto q (in other words the component 
    4746    //of d that is parallel to q. 
    … …  
    8483        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8584        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     85        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8786                           halfheight, 
    8887                           thick_layer, 
    … …  
    9392                           layer_sld, 
    9493                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9794                           d); 
    9895        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    152149    double phi) 
    153150{ 
    154     int n_stacking = (int)(fp_n_stacking + 0.5); 
    155151    double q, sin_alpha, cos_alpha; 
    156152    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     153    const double qab = q*sin_alpha; 
     154    const double qc = q*cos_alpha; 
    157155 
     156    int n_stacking = (int)(fp_n_stacking + 0.5); 
    158157    double d = 2.0 * thick_layer + thick_core; 
    159158    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     159    double answer = stacked_disks_kernel(qab, qc, 
    161160                     halfheight, 
    162161                     thick_layer, 
    … …  
    167166                     layer_sld, 
    168167                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171168                     d); 
    172169 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 r2a0b2b1  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
    … …  
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qx, double qy, 
    5249    double sld, 
    5350    double sld_solvent, 
    … …  
    6158    double q, xhat, yhat, zhat; 
    6259    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     60    const double qa = q*xhat; 
     61    const double qb = q*yhat; 
     62    const double qc = q*zhat; 
    6363 
    64     const double r = sqrt(square(radius_equat_minor*xhat) 
    65                           + square(radius_equat_major*yhat) 
    66                           + square(radius_polar*zhat)); 
    67     const double fq = sas_3j1x_x(q*r); 
    68     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     64    const double qr = sqrt(square(radius_equat_minor*qa) 
     65                           + square(radius_equat_major*qb) 
     66                           + square(radius_polar*qc)); 
     67    const double fq = sas_3j1x_x(qr); 
     68    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     69    const double drho = (sld - sld_solvent); 
    6970 
    70     return 1.0e-4 * square(s * fq); 
     71    return 1.0e-4 * square(vol * drho * fq); 
    7172} 
    72  

  • sasmodels/weights.py

    r41e7f2e rd4c33d6  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            #center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
Note: See TracChangeset for help on using the changeset viewer.