Changes in / [4f0c993:c935936] in sasmodels


Ignore:
Files:
6 deleted
28 edited

Legend:

Unmodified
Added
Removed

  • explore/jitter.py

    rd4c33d6 r85190c2  
    33dispersity and possible replacement algorithms. 
    44""" 
    5 import sys 
    6  
    75import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    86import matplotlib.pyplot as plt 
    97from matplotlib.widgets import Slider, CheckButtons 
    108from matplotlib import cm 
     9 
    1110import numpy as np 
    1211from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1312 
    14 def draw_beam(ax, view=(0, 0)): 
     13def draw_beam(ax): 
    1514    #ax.plot([0,0],[0,0],[1,-1]) 
    1615    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    … …  
    2322    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2423    y = r*np.outer(np.sin(u), np.ones_like(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] 
     24    z = np.outer(np.ones_like(u), v) 
    3225 
    3326    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    3427 
    35 def draw_jitter(ax, view, jitter): 
    36     size = [0.1, 0.4, 1.0] 
    37     draw_shape = draw_parallelepiped 
    38     #draw_shape = draw_ellipsoid 
     28def 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 
    3934 
    4035    #np.random.seed(10) 
    … …  
    6964        [ 1,  1,  1], 
    7065    ] 
    71     dtheta, dphi, dpsi = jitter 
    7266    if dtheta == 0: 
    7367        cloud = [v for v in cloud if v[0] == 0] 
    … …  
    7670    if dpsi == 0: 
    7771        cloud = [v for v in cloud if v[2] == 0] 
    78     draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     72    draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
    7973    for point in cloud: 
    80         delta = [dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    81         draw_shape(ax, size, view, delta, alpha=0.8) 
     74        shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
     75        draw_shape(ax, size, view, shimmy, alpha=0.8) 
    8276    for v in 'xyz': 
    8377        a, b, c = size 
    … …  
    8680        getattr(ax, v+'axis').label.set_text(v) 
    8781 
    88 def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     82def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
    8983    a,b,c = size 
     84    theta, phi, psi = view 
     85    dtheta, dphi, dpsi = shimmy 
     86 
    9087    u = np.linspace(0, 2 * np.pi, steps) 
    9188    v = np.linspace(0, np.pi, steps) 
    … …  
    9390    y = b*np.outer(np.sin(u), np.sin(v)) 
    9491    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    95     x, y, z = transform_xyz(view, jitter, x, y, z) 
     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] 
    9698 
    9799    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    98100 
    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  
    108 def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     101def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
    109102    a,b,c = size 
     103    theta, phi, psi = view 
     104    dtheta, dphi, dpsi = shimmy 
     105 
    110106    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    111107    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    122118    ]) 
    123119 
    124     x, y, z = transform_xyz(view, jitter, x, y, z) 
     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] 
    125125    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126126 
    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  
    136 def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gauss'): 
    137     theta, phi, psi = view 
    138     dtheta, dphi, dpsi = jitter 
     127def 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 
     136def 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) 
    139145    if dist == 'gauss': 
    140         t = np.linspace(-3, 3, n) 
    141146        weights = exp(-0.5*t**2) 
    142147    elif dist == 'rect': 
    143         t = np.linspace(0, 1, n) 
    144148        weights = np.ones_like(t) 
    145149    else: 
    146150        raise ValueError("expected dist to be 'gauss' or 'rect'") 
    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.) 
     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 
     165def 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 
    160169 
    161170def Rx(angle): 
    … …  
    179188         [0., 0., 1.]] 
    180189    return np.matrix(R) 
    181  
    182 def 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  
    191 def apply_jitter(jitter, points): 
    192     dtheta, dphi, dpsi = jitter 
    193     points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
    194     return points 
    195  
    196 def orient_relative_to_beam(view, points): 
    197     theta, phi, psi = view 
    198     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    199     return points 
    200  
    201 def 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  
    213 def 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') 
    221190 
    222191def main(): 
    … …  
    237206 
    238207    axcolor = 'lightgoldenrodyellow' 
    239  
    240208    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    241209    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    … …  
    244212    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    245213    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
    246  
    247214    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    248215    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    … …  
    250217    sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    251218    sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    252     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dpsi) 
     219    sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
    253220 
    254221    def update(val, axis=None): 
    255         view = stheta.val, sphi.val, spsi.val 
    256         jitter = sdtheta.val, sdphi.val, sdpsi.val 
     222        theta, phi, psi = stheta.val, sphi.val, spsi.val 
     223        dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
    257224        ax.cla() 
    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) 
     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) 
    262229        plt.gcf().canvas.draw() 
    263230 

  • sasmodels/details.py

    rf39759c rccd5f01  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin, radians 
     17from numpy import pi, cos, sin 
    1818 
    1919try: 
    … …  
    2929 
    3030try: 
    31     from typing import List, Tuple, Sequence 
     31    from typing import List 
    3232except ImportError: 
    3333    pass 
    3434else: 
    3535    from .modelinfo import ModelInfo 
    36     from .modelinfo import ParameterTable 
    3736 
    3837 
    … …  
    5453    coordinates, the total circumference decreases as latitude varies from 
    5554    pi r^2 at the equator to 0 at the pole, and the weight associated 
    56     with a range of latitude values needs to be scaled by this circumference. 
     55    with a range of phi values needs to be scaled by this circumference. 
    5756    This scale factor needs to be updated each time the theta value 
    5857    changes.  *theta_par* indicates which of the values in the parameter 
    … …  
    232231    nvalues = kernel.info.parameters.nvalues 
    233232    scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    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) 
     233    values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
    237234    length = np.array([len(w) for w in weights]) 
    238235    offset = np.cumsum(np.hstack((0, length))) 
    … …  
    247244    return call_details, data, is_magnetic 
    248245 
    249 def 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  
    274246 
    275247def convert_magnetism(parameters, values): 
    276     # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    277248    """ 
    278249    Convert magnetism values from polar to rectangular coordinates. 
    … …  
    284255    scale = mag[:,0] 
    285256    if np.any(scale): 
    286         theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
     257        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
    287258        cos_theta = cos(theta) 
    288259        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    … …  
    298269    """ 
    299270    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. 
    308271 
    309272    Returns [p1,p2,...],w where pj is a vector of values for parameter j 

  • sasmodels/kernel_header.c

    rbb4b509 rbb4b509  
    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

    rd4c33d6 rbde38b5  
    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 (not used) 
     27    int32_t theta_par;          // id of spherical correction variable 
    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 
    175186  int step = pd_start; 
    176187 
    … …  
    209220#endif 
    210221#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  } 
    211225  while(i0 < n0) { 
    212226    local_values.vector[p0] = v0[i0]; 
    213227    double weight0 = w0[i0] * weight1; 
    214228//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    } 
    215232#else 
    216233    const double weight0 = 1.0; 
    … …  
    227244      // Note: weight==0 must always be excluded 
    228245      if (weight0 > cutoff) { 
    229         pd_norm += weight0 * CALL_VOLUME(local_values.table); 
     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); 
    230250 
    231251        #ifdef USE_OPENMP 
    … …  
    284304#endif // !MAGNETIC 
    285305//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 
    286           result[q_index] += weight0 * scattering; 
     306          result[q_index] += weight * scattering; 
    287307        } 
    288308      } 

  • sasmodels/kernel_iq.cl

    rd4c33d6 rbde38b5  
    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 (not used) 
     27    int32_t theta_par;          // id of spherical correction variable 
    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 
    171182  int step = pd_start; 
    172183 
    … …  
    206217#endif 
    207218#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  } 
    208222  while(i0 < n0) { 
    209223    local_values.vector[p0] = v0[i0]; 
    210224    double weight0 = w0[i0] * weight1; 
    211225//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    } 
    212229#else 
    213230    const double weight0 = 1.0; 
    … …  
    215232 
    216233//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); 
    217235 
    218236    #ifdef INVALID 
    … …  
    223241      // Note: weight==0 must always be excluded 
    224242      if (weight0 > cutoff) { 
    225         pd_norm += weight0 * CALL_VOLUME(local_values.table); 
     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); 
    226247 
    227248#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    … …  
    275296        const double scattering = CALL_IQ(q, q_index, local_values.table); 
    276297#endif // !MAGNETIC 
    277         this_result += weight0 * scattering; 
     298        this_result += weight * scattering; 
    278299      } 
    279300    } 

  • sasmodels/kernelpy.py

    rd4c33d6 rdbfd471  
    197197 
    198198    pd_norm = 0.0 
     199    spherical_correction = 1.0 
    199200    partial_weight = np.NaN 
    200201    weight = np.NaN 
    201202 
    202203    p0_par = call_details.pd_par[0] 
     204    p0_is_theta = (p0_par == call_details.theta_par) 
    203205    p0_length = call_details.pd_length[0] 
    204206    p0_index = p0_length 
    … …  
    217219            parameters[pd_par] = pd_value[pd_offset+pd_index] 
    218220            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) 
    219224            p0_index = loop_index%p0_length 
    220225 
    221226        weight = partial_weight * pd_weight[p0_offset + p0_index] 
    222227        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) 
    223231        p0_index += 1 
    224232        if weight > cutoff: 
    … …  
    231239 
    232240            # update value and norm 
     241            weight *= spherical_correction 
    233242            total += weight * Iq 
    234243            pd_norm += weight * form_volume() 

  • sasmodels/model_test.py

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

  • sasmodels/models/barbell.c

    r2a0b2b1 r592343f  
    1010//barbell kernel - same as dumbell 
    1111static double 
    12 _bell_kernel(double qab, double qc, double h, double radius_bell, 
    13              double half_length) 
     12_bell_kernel(double q, double h, double radius_bell, 
     13             double half_length, double sin_alpha, double cos_alpha) 
    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 = 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) 
     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) 
    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(qab_r*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    3737        total += Gauss76Wt[i] * Fq; 
    … …  
    4444 
    4545static double 
    46 _fq(double qab, double qc, double h, 
    47     double radius_bell, double radius, double half_length) 
     46_fq(double q, double h, 
     47    double radius_bell, double radius, double half_length, 
     48    double sin_alpha, double cos_alpha) 
    4849{ 
    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); 
     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); 
    5253    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5354    const double Aq = bell_fq + cyl_fq; 
    … …  
    8384        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8485        SINCOS(alpha, sin_alpha, cos_alpha); 
    85         const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
     86        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
    8687        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8788    } 
    … …  
    102103    double q, sin_alpha, cos_alpha; 
    103104    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; 
    106105 
    107106    const double h = -sqrt(square(radius_bell) - square(radius)); 
    108     const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
     107    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
    109108 
    110109    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r7e0b281 r50beefe  
    1 static double 
    2 bcc_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 
     1double form_volume(double radius); 
     2double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
     3double 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); 
    136 
    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 
     7double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
     8double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
     9double _sphereform(double q, double radius, double sld, double solvent_sld); 
    3610 
    37     return Zq; 
     11 
     12double _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); 
    3820} 
    3921 
     22double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    4023 
    41 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    42 static double 
    43 bcc_volume_fraction(double radius, double dnn) 
    44 { 
    45     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
     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); 
    4641} 
    4742 
    48 static double 
    49 form_volume(double radius) 
    50 { 
     43double form_volume(double radius){ 
    5144    return sphere_volume(radius); 
    5245} 
    5346 
    5447 
    55 static double Iq(double q, double dnn, 
     48double Iq(double q, double dnn, 
    5649  double d_factor, double radius, 
    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; 
     50  double sld, double solvent_sld){ 
    6551 
    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; 
     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; 
    9084} 
    9185 
    9286 
    93 static double Iqxy(double qx, double qy, 
     87double Iqxy(double qx, double qy, 
    9488    double dnn, double d_factor, double radius, 
    9589    double sld, double solvent_sld, 
    … …  
    9892    double q, zhat, yhat, xhat; 
    9993    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; 
    10394 
    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; 
     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; 
    108111} 

  • sasmodels/models/bcc_paracrystal.py

    r2a0b2b1 r69e1afc  
    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

    r2a0b2b1 r592343f  
    1414//   radius_cap is the radius of the lens 
    1515//   length is the cylinder length, or the separation between the lens halves 
    16 //   theta is the angle of the cylinder wrt q. 
     16//   alpha is the angle of the cylinder wrt q. 
    1717static double 
    18 _cap_kernel(double qab, double qc, double h, double radius_cap, 
    19             double half_length) 
     18_cap_kernel(double q, double h, double radius_cap, 
     19                      double half_length, double sin_alpha, double cos_alpha) 
    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(theta)) 
     28    //    cos (q (R t - h + L/2) cos(alpha)) 
    2929    // so turn it into: 
    3030    //    cos (m t + b) 
    3131    // where: 
    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) 
     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) 
    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(qab_r*sqrt(radical)); 
     41        const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
    4242        const double Fq = cos(m*t + b) * radical * bj; 
    4343        total += Gauss76Wt[i] * Fq; 
    … …  
    5050 
    5151static double 
    52 _fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
     52_fq(double q, double h, double radius_cap, double radius, double half_length, 
     53    double sin_alpha, double cos_alpha) 
    5354{ 
    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); 
     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); 
    5758    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5859    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    100101    double total = 0.0; 
    101102    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    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; 
    129127 
    130128    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    131     const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
     129    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
    132130 
    133131    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

    r2a0b2b1 rb260926  
    1 static double 
    2 form_volume(double radius, double thick_rim, double thick_face, double length) 
     1double form_volume(double radius, double thick_rim, double thick_face, double length); 
     2double 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 
     13double 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 
     26double form_volume(double radius, double thick_rim, double thick_face, double length) 
    327{ 
    4     return M_PI*square(radius+thick_rim)*(length+2.0*thick_face); 
     28    return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2.0*thick_face); 
    529} 
    630 
    731static double 
    8 bicelle_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) 
     32bicelle_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) 
    1843{ 
    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); 
     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); 
    2550 
    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); 
     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); 
    3055 
    3156    const double t = vol1*dr1*si1*be1 + 
    … …  
    3358                     vol3*dr3*si2*be1; 
    3459 
    35     return t; 
     60    const double retval = t*t; 
     61 
     62    return retval; 
     63 
    3664} 
    3765 
    3866static double 
    39 Iq(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) 
     67bicelle_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) 
    4876{ 
    4977    // set up the integration end points 
    … …  
    5179    const double halflength = 0.5*length; 
    5280 
    53     double total = 0.0; 
     81    double summ = 0.0; 
    5482    for(int i=0;i<N_POINTS_76;i++) { 
    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; 
     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; 
    6190    } 
    6291 
    6392    // calculate value of integral to return 
    64     double answer = total*uplim; 
     93    double answer = uplim*summ; 
     94    return answer; 
     95} 
     96 
     97static double 
     98bicelle_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) 
     109{ 
     110    double q, sin_alpha, cos_alpha; 
     111    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
     112 
     113    double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 
     114                           0.5*length, core_sld, face_sld, rim_sld, 
     115                           solvent_sld, sin_alpha, cos_alpha); 
    65116    return 1.0e-4*answer; 
    66117} 
    67118 
    68 static double 
    69 Iqxy(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) 
     119double 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) 
    80128{ 
    81     double q, sin_alpha, cos_alpha; 
    82     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; 
     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} 
    85133 
    86     double fq = bicelle_kernel(qab, qc, radius, thick_rim, thick_face, 
    87                            0.5*length, core_sld, face_sld, rim_sld, 
    88                            solvent_sld); 
    89     return 1.0e-4*fq*fq; 
     134 
     135double 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; 
    90160} 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r2a0b2b1 rdedcf34  
    1717        double thick_face, 
    1818        double length, 
    19         double sld_core, 
    20         double sld_face, 
    21         double sld_rim, 
    22         double sld_solvent) 
     19        double rhoc, 
     20        double rhoh, 
     21        double rhor, 
     22        double rhosolv) 
    2323{ 
     24    double si1,si2,be1,be2; 
    2425     // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 
    2526     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle 
     27     //    const double uplim = M_PI_4; 
    2628    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 
    2735    const double r_major = r_minor * x_core; 
    2836    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    2937    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
    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); 
     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); 
    3644 
    3745    //initialize integral 
    … …  
    3947    for(int i=0;i<76;i++) { 
    4048        //setup inner integral over the ellipsoidal cross-section 
    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; 
     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); 
    5158        for(int j=0;j<76;j++) { 
    5259            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    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; 
     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); 
    6470        } 
    6571        //now calculate outer integral 
    … …  
    7783          double thick_face, 
    7884          double length, 
    79           double sld_core, 
    80           double sld_face, 
    81           double sld_rim, 
    82           double sld_solvent, 
     85          double rhoc, 
     86          double rhoh, 
     87          double rhor, 
     88          double rhosolv, 
    8389          double theta, 
    8490          double phi, 
    8591          double psi) 
    8692{ 
     93       // THIS NEEDS TESTING 
    8794    double q, xhat, yhat, zhat; 
    8895    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    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; 
     96    const double dr1 = rhoc-rhoh; 
     97    const double dr2 = rhor-rhosolv; 
     98    const double dr3 = rhoh-rhor; 
    9699    const double r_major = r_minor*x_core; 
    97100    const double halfheight = 0.5*length; 
    … …  
    101104 
    102105    // Compute effective radius in rotated coordinates 
    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; 
     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; 
    112115} 
     116 

  • sasmodels/models/core_shell_cylinder.c

    r2a0b2b1 r592343f  
     1double form_volume(double radius, double thickness, double length); 
     2double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
     3    double radius, double thickness, double length); 
     4double 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 
    17// vd = volume * delta_rho 
    2 // besarg = q * R * sin(theta) 
    3 // siarg = q * L/2 * cos(theta) 
    4 static double _cyl(double vd, double besarg, double siarg) 
     8// besarg = q * R * sin(alpha) 
     9// siarg = q * L/2 * cos(alpha) 
     10double _cyl(double vd, double besarg, double siarg); 
     11double _cyl(double vd, double besarg, double siarg) 
    512{ 
    613    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    714} 
    815 
    9 static double 
    10 form_volume(double radius, double thickness, double length) 
     16double form_volume(double radius, double thickness, double length) 
    1117{ 
    12     return M_PI*square(radius+thickness)*(length+2.0*thickness); 
     18    return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
    1319} 
    1420 
    15 static double 
    16 Iq(double q, 
     21double Iq(double q, 
    1722    double core_sld, 
    1823    double shell_sld, 
    … …  
    2328{ 
    2429    // precalculate constants 
    25     const double core_r = radius; 
    26     const double core_h = 0.5*length; 
     30    const double core_qr = q*radius; 
     31    const double core_qh = q*0.5*length; 
    2732    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    28     const double shell_r = (radius + thickness); 
    29     const double shell_h = (0.5*length + thickness); 
     33    const double shell_qr = q*(radius + thickness); 
     34    const double shell_qh = q*(0.5*length + thickness); 
    3035    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3136    double total = 0.0; 
     37    // double lower=0, upper=M_PI_2; 
    3238    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    4347    } 
    4448    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6064    double q, sin_alpha, cos_alpha; 
    6165    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; 
    6466 
    65     const double core_r = radius; 
    66     const double core_h = 0.5*length; 
     67    const double core_qr = q*radius; 
     68    const double core_qh = q*0.5*length; 
    6769    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    68     const double shell_r = (radius + thickness); 
    69     const double shell_h = (0.5*length + thickness); 
     70    const double shell_qr = q*(radius + thickness); 
     71    const double shell_qh = q*(0.5*length + thickness); 
    7072    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7173 
    72     const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    73         + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     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); 
    7476    return 1.0e-4 * fq * fq; 
    7577} 

  • sasmodels/models/core_shell_ellipsoid.c

    r2a0b2b1 r0a3d9b2  
     1double form_volume(double radius_equat_core, 
     2                   double polar_core, 
     3                   double equat_shell, 
     4                   double polar_shell); 
     5double 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); 
    113 
    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> 
    9 static 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; 
    1914 
    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; 
     15double 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); 
    2425 
    25     return fq_core + fq_shell; 
    26 } 
    2726 
    28 static double 
    29 form_volume(double radius_equat_core, 
    30     double x_core, 
    31     double thick_shell, 
    32     double x_polar_shell) 
     27double form_volume(double radius_equat_core, 
     28                   double x_core, 
     29                   double thick_shell, 
     30                   double x_polar_shell) 
    3331{ 
    3432    const double equat_shell = radius_equat_core + thick_shell; 
    … …  
    3937 
    4038static double 
    41 Iq(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) 
     39core_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) 
    4947{ 
    50     const double sld_core_shell = core_sld - shell_sld; 
    51     const double sld_shell_solvent = shell_sld - solvent_sld; 
     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 
    5255 
    5356    const double polar_core = radius_equat_core*x_core; 
    … …  
    5558    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    5659 
    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 
     60    double summ = 0.0;   //initialize intergral 
    6161    for(int i=0;i<76;i++) { 
    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; 
     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; 
    6966    } 
    70     total *= m; 
     67    summ *= 0.5*(uplim-lolim); 
    7168 
    7269    // convert to [cm-1] 
    73     return 1.0e-4 * total; 
     70    return 1.0e-4 * summ; 
    7471} 
    7572 
    7673static double 
    77 Iqxy(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) 
     74core_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) 
    8784{ 
    8885    double q, sin_alpha, cos_alpha; 
    8986    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; 
    9287 
    93     const double sld_core_shell = core_sld - shell_sld; 
    94     const double sld_shell_solvent = shell_sld - solvent_sld; 
     88    const double sldcs = core_sld - shell_sld; 
     89    const double sldss = shell_sld- solvent_sld; 
    9590 
    9691    const double polar_core = radius_equat_core*x_core; 
    … …  
    9893    const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 
    9994 
    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); 
     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); 
    104105 
    105106    //convert to [cm-1] 
    106     return 1.0e-4 * fq * fq; 
     107    answer *= 1.0e-4; 
     108 
     109    return answer; 
    107110} 
     111 
     112double 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 
     134double 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

    r2a0b2b1 r9802ab3  
    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/gauss76.c", "core_shell_ellipsoid.c"] 
     141source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
     142          "core_shell_ellipsoid.c"] 
    142143 
    143144def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    r2a0b2b1 r92dfe0c  
    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; 
    141138 
    142139    // cspkernel in csparallelepiped recoded here 
    … …  
    163160    double tc = length_a + 2.0*thick_rim_c; 
    164161    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    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  
     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     
    172169 
    173170    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    177174               + drB*siA*(siBt-siB)*siC*V2 
    178175               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    179  
     176    
    180177    return 1.0e-4 * f * f; 
    181178} 

  • sasmodels/models/cylinder.c

    rb34fc77 r592343f  
     1double form_volume(double radius, double length); 
     2double fq(double q, double sn, double cn,double radius, double length); 
     3double orient_avg_1D(double q, double radius, double length); 
     4double Iq(double q, double sld, double solvent_sld, double radius, double length); 
     5double Iqxy(double qx, double qy, double sld, double solvent_sld, 
     6    double radius, double length, double theta, double phi); 
     7 
    18#define INVALID(v) (v.radius<0 || v.length<0) 
    29 
    3 static double 
    4 form_volume(double radius, double length) 
     10double form_volume(double radius, double length) 
    511{ 
    612    return M_PI*radius*radius*length; 
    713} 
    814 
    9 static double 
    10 fq(double qab, double qc, double radius, double length) 
     15double fq(double q, double sn, double cn, double radius, double length) 
    1116{ 
    12     return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     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); 
    1321} 
    1422 
    15 static double 
    16 orient_avg_1D(double q, double radius, double length) 
     23double orient_avg_1D(double q, double radius, double length) 
    1724{ 
    1825    // translate a point in [-1,1] to a point in [0, pi/2] 
    1926    const double zm = M_PI_4; 
    20     const double zb = M_PI_4; 
     27    const double zb = M_PI_4;  
    2128 
    2229    double total = 0.0; 
    2330    for (int i=0; i<76 ;i++) { 
    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; 
     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; 
    3036    } 
    3137    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3339} 
    3440 
    35 static double 
    36 Iq(double q, 
     41double Iq(double q, 
    3742    double sld, 
    3843    double solvent_sld, 
    … …  
    4449} 
    4550 
    46 static double 
    47 Iqxy(double qx, double qy, 
     51 
     52double Iqxy(double qx, double qy, 
    4853    double sld, 
    4954    double solvent_sld, 
    … …  
    5560    double q, sin_alpha, cos_alpha; 
    5661    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    57     const double qab = q*sin_alpha; 
    58     const double qc = q*cos_alpha; 
    59  
     62    //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha); 
    6063    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    61     const double form = fq(qab, qc, radius, length); 
     64    const double form = fq(q, sin_alpha, cos_alpha, radius, length); 
    6265    return 1.0e-4 * square(s * form); 
    6366} 

  • sasmodels/models/ellipsoid.c

    r2a0b2b1 r3b571ae  
    1 static double 
    2 form_volume(double radius_polar, double radius_equatorial) 
     1double form_volume(double radius_polar, double radius_equatorial); 
     2double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
     3double Iqxy(double qx, double qy, double sld, double sld_solvent, 
     4    double radius_polar, double radius_equatorial, double theta, double phi); 
     5 
     6double form_volume(double radius_polar, double radius_equatorial) 
    37{ 
    48    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    59} 
    610 
    7 static  double 
    8 Iq(double q, 
     11double Iq(double q, 
    912    double sld, 
    1013    double sld_solvent, 
    … …  
    3841} 
    3942 
    40 static double 
    41 Iqxy(double qx, double qy, 
     43double Iqxy(double qx, double qy, 
    4244    double sld, 
    4345    double sld_solvent, 
    … …  
    4951    double q, sin_alpha, cos_alpha; 
    5052    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    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); 
     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); 
    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

    r2a0b2b1 r61104c8  
    1 static double 
     1double form_volume(double radius_minor, double r_ratio, double length); 
     2double Iq(double q, double radius_minor, double r_ratio, double length, 
     3          double sld, double solvent_sld); 
     4double 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 
     8double 
    29form_volume(double radius_minor, double r_ratio, double length) 
    310{ 
    … …  
    512} 
    613 
    7 static double 
     14double 
    815Iq(double q, double radius_minor, double r_ratio, double length, 
    916   double sld, double solvent_sld) 
    … …  
    5461 
    5562 
    56 static double 
     63double 
    5764Iqxy(double qx, double qy, 
    5865     double radius_minor, double r_ratio, double length, 
    … …  
    6269    double q, xhat, yhat, zhat; 
    6370    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; 
    6771 
    6872    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 
    6973    // Given:    radius_major = r_ratio * radius_minor 
    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); 
     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); 
    7377    const double Aq = be * si; 
    7478    const double delrho = sld - solvent_sld; 

  • sasmodels/models/fcc_paracrystal.c

    r7e0b281 r50beefe  
    1 static double 
    2 fcc_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 
     1double form_volume(double radius); 
     2double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
     3double 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); 
    136 
    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)); 
     7double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
     8double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
    269 
    27     return Zq; 
     10 
     11double _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); 
    2819} 
    2920 
     21double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    3022 
    31 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    32 static double 
    33 fcc_volume_fraction(double radius, double dnn) 
    34 { 
    35     return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
     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); 
    3640} 
    3741 
    38 static double 
    39 form_volume(double radius) 
    40 { 
     42double form_volume(double radius){ 
    4143    return sphere_volume(radius); 
    4244} 
    4345 
    4446 
    45 static double Iq(double q, double dnn, 
     47double Iq(double q, double dnn, 
    4648  double d_factor, double radius, 
    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; 
     49  double sld, double solvent_sld){ 
    5550 
    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); 
     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); 
    7954 
    80     return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     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; 
     83 
     84 
    8185} 
    8286 
    83  
    84 static double Iqxy(double qx, double qy, 
     87double Iqxy(double qx, double qy, 
    8588    double dnn, double d_factor, double radius, 
    8689    double sld, double solvent_sld, 
    … …  
    8992    double q, zhat, yhat, xhat; 
    9093    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     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; 
     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; 
    99112} 

  • sasmodels/models/hollow_cylinder.c

    r2a0b2b1 r592343f  
    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 
    1617static double 
    17 _fq(double qab, double qc, 
    18     double radius, double thickness, double length) 
     18_hollow_cylinder_kernel(double q, 
     19    double radius, double thickness, double length, double sin_val, double cos_val) 
    1920{ 
    20     const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
    21     const double lam2 = sas_2J1x_x(radius*qab); 
     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); 
    2224    const double gamma_sq = square(radius/(radius+thickness)); 
    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); 
     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); 
    2729    return psi*t2; 
    2830} 
    … …  
    4143{ 
    4244    const double lower = 0.0; 
    43     const double upper = 1.0;        //limits of numerical integral 
     45    const double upper = 1.0;           //limits of numerical integral 
    4446 
    45     double summ = 0.0;            //initialize intergral 
     47    double summ = 0.0;                  //initialize intergral 
    4648    for (int i=0;i<76;i++) { 
    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; 
     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; 
    5254    } 
    5355 
    … …  
    6466    double q, sin_alpha, cos_alpha; 
    6567    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, 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); 
     68    const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
     69        sin_alpha, cos_alpha); 
    7070 
    7171    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
     72    return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
    7373} 
     74 

  • sasmodels/models/parallelepiped.c

    r2a0b2b1 rd605080  
    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; 
    7471 
    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); 
     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); 
    7875    const double V = form_volume(length_a, length_b, length_c); 
    7976    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r7e0b281 r50beefe  
    1 static double 
    2 sc_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; 
     1double form_volume(double radius); 
    72 
    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)); 
     3double Iq(double q, 
     4          double dnn, 
     5          double d_factor, 
     6          double radius, 
     7          double sphere_sld, 
     8          double solvent_sld); 
    209 
    21     return Zq; 
    22 } 
     10double 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); 
    2319 
    24 // occupied volume fraction calculated from lattice symmetry and sphere radius 
    25 static double 
    26 sc_volume_fraction(double radius, double dnn) 
    27 { 
    28     return sphere_volume(radius/dnn); 
    29 } 
    30  
    31 static double 
    32 form_volume(double radius) 
     20double form_volume(double radius) 
    3321{ 
    3422    return sphere_volume(radius); 
    3523} 
    3624 
     25static double 
     26sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     27{ 
     28    double cnt, snt; 
     29    SINCOS(theta, cnt, snt); 
    3730 
    38 static double Iq(double q, double dnn, 
    39   double d_factor, double radius, 
    40   double sld, double solvent_sld) 
    41 { 
    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; 
     31    double cnp, snp; 
     32    SINCOS(phi, cnp, snp); 
    4833 
    49  
    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; 
     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); 
    7542} 
    7643 
     44static double 
     45sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
     46{ 
     47    //Function to calculate integrand values for simple cubic structure 
    7748 
    78 static 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) 
     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 
     61double 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 
     103double 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) 
    82112{ 
    83113    double q, zhat, yhat, xhat; 
    84114    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; 
    88115 
    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; 
     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; 
    93128} 

  • sasmodels/models/sc_paracrystal.py

    r2a0b2b1 r69e1afc  
    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

    rb34fc77 r19f996b  
    11static double stacked_disks_kernel( 
    2     double qab, 
    3     double qc, 
     2    double q, 
    43    double halfheight, 
    54    double thick_layer, 
    … …  
    109    double layer_sld, 
    1110    double solvent_sld, 
     11    double sin_alpha, 
     12    double cos_alpha, 
    1213    double d) 
    1314 
    … …  
    1920    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2021 
    21     const double besarg1 = radius*qab; 
    22     //const double besarg2 = radius*qab; 
     22    const double besarg1 = q*radius*sin_alpha; 
     23    //const double besarg2 = q*radius*sin_alpha; 
    2324 
    24     const double sinarg1 = halfheight*qc; 
    25     const double sinarg2 = (halfheight+thick_layer)*qc; 
     25    const double sinarg1 = q*halfheight*cos_alpha; 
     26    const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
    2627 
    2728    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4243 
    4344    // loop for the structure factor S(q) 
    44     double qd_cos_alpha = d*qc; 
     45    double qd_cos_alpha = q*d*cos_alpha; 
    4546    //d*cos_alpha is the projection of d onto q (in other words the component 
    4647    //of d that is parallel to q. 
    … …  
    8384        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8485        SINCOS(zi, sin_alpha, cos_alpha); 
    85         double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
     86        double yyy = stacked_disks_kernel(q, 
    8687                           halfheight, 
    8788                           thick_layer, 
    … …  
    9293                           layer_sld, 
    9394                           solvent_sld, 
     95                           sin_alpha, 
     96                           cos_alpha, 
    9497                           d); 
    9598        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    149152    double phi) 
    150153{ 
     154    int n_stacking = (int)(fp_n_stacking + 0.5); 
    151155    double q, sin_alpha, cos_alpha; 
    152156    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; 
    155157 
    156     int n_stacking = (int)(fp_n_stacking + 0.5); 
    157158    double d = 2.0 * thick_layer + thick_core; 
    158159    double halfheight = 0.5*thick_core; 
    159     double answer = stacked_disks_kernel(qab, qc, 
     160    double answer = stacked_disks_kernel(q, 
    160161                     halfheight, 
    161162                     thick_layer, 
    … …  
    166167                     layer_sld, 
    167168                     solvent_sld, 
     169                     sin_alpha, 
     170                     cos_alpha, 
    168171                     d); 
    169172 

  • sasmodels/models/triaxial_ellipsoid.c

    r2a0b2b1 r68dd6a9  
     1double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
     2double Iq(double q, double sld, double sld_solvent, 
     3    double radius_equat_minor, double radius_equat_major, double radius_polar); 
     4double 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 
    17//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    28 
    3 static double 
    4 form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     9 
     10double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    511{ 
    612    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    713} 
    814 
    9 static double 
    10 Iq(double q, 
     15double Iq(double q, 
    1116    double sld, 
    1217    double sld_solvent, 
    … …  
    4045    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4146    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    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; 
     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; 
    4549} 
    4650 
    47 static double 
    48 Iqxy(double qx, double qy, 
     51double Iqxy(double qx, double qy, 
    4952    double sld, 
    5053    double sld_solvent, 
    … …  
    5861    double q, xhat, yhat, zhat; 
    5962    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 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); 
     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); 
    7069 
    71     return 1.0e-4 * square(vol * drho * fq); 
     70    return 1.0e-4 * square(s * fq); 
    7271} 
     72 

  • sasmodels/weights.py

    rd4c33d6 rf1a8811  
    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 
    6157        if sigma == 0 or self.npts < 2: 
    6258            if lb <= center <= ub: 
Note: See TracChangeset for help on using the changeset viewer.