Changeset a1c5758 in sasmodels for sasmodels/models


Ignore:
Timestamp:
Oct 23, 2017 10:25:27 AM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
34823bc
Parents:
ea60e08 (diff), 55d88b4 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket-776-orientation

Location:
sasmodels/models
Files:
1 added
21 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/flexible_cylinder.py

    r31df0c9 r2573fa1  
    9090    length = 10**np.random.uniform(2, 6) 
    9191    radius = 10**np.random.uniform(1, 3) 
    92     kuhn_length = 10**np.random.uniform(-2, -0.7)*length  # at least 10 segments 
     92    kuhn_length = 10**np.random.uniform(-2, 0)*length 
    9393    pars = dict( 
    9494        length=length, 
    … …  
    100100tests = [ 
    101101    # Accuracy tests based on content in test/utest_other_models.py 
    102     # Currently fails in OCL 
    103     # [{'length':     1000.0, 
    104     #  'kuhn_length': 100.0, 
    105     #  'radius':       20.0, 
    106     #  'sld':           1.0, 
    107     #  'sld_solvent':   6.3, 
    108     #  'background':    0.0001, 
    109     #  }, 0.001, 3509.2187], 
     102    [{'length':     1000.0,  # test T1 
     103      'kuhn_length': 100.0, 
     104      'radius':       20.0, 
     105      'sld':           1.0, 
     106      'sld_solvent':   6.3, 
     107      'background':    0.0001, 
     108     }, 0.001, 3509.2187], 
    110109 
    111110    # Additional tests with larger range of parameters 
    112     [{'length':    1000.0, 
     111    [{'length':    1000.0,  # test T2 
    113112      'kuhn_length': 100.0, 
    114113      'radius':       20.0, 
    … …  
    117116      'background':    0.0001, 
    118117     }, 1.0, 0.000595345], 
    119     [{'length':        10.0, 
     118    [{'length':        10.0,  # test T3 
    120119      'kuhn_length': 800.0, 
    121120      'radius':        2.0, 
    … …  
    124123      'background':    0.001, 
    125124     }, 0.1, 1.55228], 
    126     [{'length':        100.0, 
     125    [{'length':        100.0,  # test T4 
    127126      'kuhn_length': 800.0, 
    128127      'radius':       50.0, 
    … …  
    133132    ] 
    134133 
     134# There are a few branches in the code that ought to have test values: 
     135# 
     136# For length > 4 * kuhn_length 
     137#        if length > 10 * kuhn_length then C is scaled by 3.06 (L/b)^(-0.44) 
     138#        q*kuhn_length <= 3.1  => Sexv_new 
     139#           dS/dQ < 0 has different behaviour from dS/dQ >= 0 
     140#  T2    q*kuhn_length > 3.1   => a_long 
     141# 
     142# For length <= 4 * kuhn_length 
     143#        q*kuhn_length <= max(1.9/Rg_short, 3.0)  => Sdebye((q*Rg)^2) 
     144#           q*Rg < 0.5 uses Pade approx, q*Rg > 1.0 uses math lib 
     145#  T3,T4 q*kuhn_length > max(1.9/Rg_short, 3.0)   => a_short 
     146# 
     147# Note that the transitions between branches may be abrupt.  You can see a 
     148# several percent change around length=10*kuhn_length and length=4*kuhn_length 
     149# using the following: 
     150# 
     151#    sascomp flexible_cylinder -calc=double -sets=10 length=10*kuhn_length,10.000001*kuhn_length 
     152#    sascomp flexible_cylinder -calc=double -sets=10 length=4*kuhn_length,4.000001*kuhn_length 
     153# 
     154# The transition between low q and high q around q*kuhn_length = 3 seems 
     155# to be good to 4 digits or better.  This was tested by computing the value 
     156# on each branches near the transition point and reporting the relative error 
     157# for kuhn lengths of 10, 100 and 1000 and a variety of length:kuhn_length 
     158# ratios. 

  • sasmodels/models/lib/wrc_cyl.c

    r92ce163 r18a2bfc  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
    4 double Sk_WR(double q, double L, double b); 
    5  
    6  
    7 static double 
    8 AlphaSquare(double x) 
    9 { 
    10     // Potentially faster. Needs proper benchmarking. 
    11     // add native_powr to kernel_template 
    12     //double t = native_powr( (1.0 + (x/3.12)*(x/3.12) + 
    13     //     (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
    14     //return t; 
    15  
    16     return pow(1.0+square(x/3.12)+cube(x/8.67), 0.176/3.0); 
    17 } 
    18  
    19 // 
    20 static double 
    21 Rgsquarezero(double q, double L, double b) 
    22 { 
    23     const double r = b/L; 
    24     return (L*b/6.0) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
    25  
    26 } 
    27  
    28 // 
    29 static double 
    30 Rgsquareshort(double q, double L, double b) 
    31 { 
    32     return AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
    33 } 
    34  
    35 // 
    36 static double 
    37 Rgsquare(double q, double L, double b) 
    38 { 
    39     return AlphaSquare(L/b)*L*b/6.0; 
    40 } 
    41  
    42 static double 
    43 sech_WR(double x) 
    44 { 
    45     return(1/cosh(x)); 
    46 } 
    47  
    48 static double 
    49 a1long(double q, double L, double b, double p1, double p2, double q0) 
    50 { 
    51     double C; 
    52  
    53     if( L/b > 10.0) { 
    54         C = 3.06/pow((L/b),0.44); 
     4 
     5static double 
     6Rgsquare(double L, double b) 
     7{ 
     8    const double x = L/b; 
     9    // Use Horner's method to evaluate Pedersen eq 15: 
     10    //     alpha^2 = [1.0 + (x/3.12)^2 + (x/8.67)^3] ^ (0.176/3) 
     11    const double alphasq = 
     12        pow(1.0 + x*x*(1.027284681130835e-01 + 1.534414548417740e-03*x), 
     13            5.866666666666667e-02); 
     14    return alphasq*L*b/6.0; 
     15} 
     16 
     17static double 
     18Rgsquareshort(double L, double b) 
     19{ 
     20    const double r = b/L;  // = 1/n_b in Pedersen ref. 
     21    return Rgsquare(L, b) * (1.0 + r*(-1.5 + r*(1.5 + r*0.75*expm1(-2.0/r)))); 
     22} 
     23 
     24static double 
     25w_WR(double x) 
     26{ 
     27    // Pedersen eq. 16: 
     28    //    w = [1 + tanh((x-C4)/C5)]/2 
     29    const double C4 = 1.523; 
     30    const double C5 = 0.1477; 
     31    return 0.5 + 0.5*tanh((x - C4)/C5); 
     32} 
     33 
     34static double 
     35Sdebye(double qsq) 
     36{ 
     37#if FLOAT_SIZE>4 
     38#  define DEBYE_CUTOFF 0.25  // 1e-15 error 
     39#else 
     40#  define DEBYE_CUTOFF 1.0  // 4e-7 error 
     41#endif 
     42 
     43    if (qsq < DEBYE_CUTOFF) { 
     44        const double x = qsq; 
     45        // mathematica: PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     46        const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     47        const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 
     48        return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 
     49             / ((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 
    5550    } else { 
    56         C = 1.0; 
     51        return 2.*(expm1(-qsq) + qsq)/(qsq*qsq); 
    5752    } 
    58  
     53} 
     54 
     55static double 
     56a_long(double q, double L, double b/*, double p1, double p2, double q0*/) 
     57{ 
     58    const double p1 = 4.12; 
     59    const double p2 = 4.42; 
     60    const double q0 = 3.1; 
     61 
     62    // Constants C1, ..., C5 derived from least squares fit (Pedersen, eq 13,16) 
    5963    const double C1 = 1.22; 
    6064    const double C2 = 0.4288; 
    … …  
    6468    const double miu = 0.585; 
    6569 
    66     const double Rg2 = Rgsquare(q,L,b); 
    67     const double Rg22 = Rg2*Rg2; 
    68     const double Rg = sqrt(Rg2); 
    69     const double Rgb = Rg*q0/b; 
    70  
    71     const double b2 = b*b; 
     70    const double C = (L/b>10.0 ? 3.06*pow(L/b, -0.44) : 1.0); 
     71    //printf("branch B-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     72    const double r2 = Rgsquare(L,b); 
     73    const double r = sqrt(r2); 
     74    const double qr_b = q0*r/b; 
     75    const double qr_b_sq = qr_b*qr_b; 
     76    const double qr_b_4 = qr_b_sq*qr_b_sq; 
     77    const double qr_b_miu = pow(qr_b, -1.0/miu); 
     78    const double em1_qr_b_sq = expm1(-qr_b_sq); 
     79    const double sech2 = 1.0/square(cosh((qr_b-C4)/C5)); 
     80    const double w = w_WR(qr_b); 
     81 
     82    const double t1 = pow(q0, 1.0 + p1 + p2)/(b*(p1-p2)); 
     83    const double t2 = C/(15.0*L) * ( 
     84        + 14.0*b*b*em1_qr_b_sq/(q0*qr_b_sq) 
     85        + 2.0*q0*r2*exp(-qr_b_sq)*(11.0 + 7.0/qr_b_sq)); 
     86    const double t11 = ((C3*qr_b_miu + C2)*qr_b_miu + C1)*qr_b_miu; 
     87    const double t3 = r*sech2/(2.*C5)*t11; 
     88    const double t4 = r*(em1_qr_b_sq + qr_b_sq)*sech2 / (C5*qr_b_4); 
     89    const double t5 = -4.0*r*qr_b*em1_qr_b_sq/qr_b_4 * (1.0 - w); 
     90    const double t10 = 2.0*(em1_qr_b_sq + qr_b_sq)/qr_b_4 * (1.0 - w); //=Sdebye*(1-w) 
     91    const double t6 = 4.0*b/q0 * t10; 
     92    const double t7 = r*((-3.0*C3*qr_b_miu -2.0*C2)*qr_b_miu -1.0*C1)*qr_b_miu/(miu*qr_b); 
     93    const double t9 = C*b/L * (4.0 - exp(-qr_b_sq) * (11.0 + 7.0/qr_b_sq) + 7.0/qr_b_sq)/15.0; 
     94    const double t12 = b*b*M_PI/(L*q0*q0) + t2 + t3 - t4 + t5 - t6 + t7*w; 
     95    const double t13 = -b*M_PI/(L*q0) + t9 + t10 + t11*w; 
     96 
     97    const double a1 = pow(q0,p1)*t13 - t1*pow(q0,-p2)*(t12 + b*p1/q0*t13); 
     98    const double a2 = t1*pow(q0,-p1)*(t12 + b*p1/q0*t13); 
     99 
     100    const double ans = a1*pow(q*b, -p1) + a2*pow(q*b, -p2) + M_PI/(q*L); 
     101    return ans; 
     102} 
     103 
     104static double 
     105_short(double r2, double exp_qr_b, double L, double b, double p1short, 
     106        double p2short, double q0) 
     107{ 
     108    const double qr2 = q0*q0 * r2; 
    72109    const double b3 = b*b*b; 
    73     const double b4 = b3*b; 
    74     const double q02 = q0*q0; 
    75     const double q03 = q0*q0*q0; 
    76     const double q04 = q03*q0; 
    77     const double q05 = q04*q0; 
    78  
    79     const double Rg02 = Rg2*q02; 
    80  
    81     const double t1 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2))) * 
    82          ((11.0/15.0 + (7.0*b2)/(15.0*Rg02))) + 
    83          (7.0*b2)/(15.0*Rg02)))); 
    84  
    85     const double t2 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    86          Rg02/b2))*((1.0 + 0.5*(((-1.0) - 
    87          tanh((-C4 + Rgb/C5))))))); 
    88  
    89     const double t3 = ((C3*pow(Rgb,((-3.0)/miu)) + 
    90          C2*pow(Rgb,((-2.0)/miu)) + 
    91          C1*pow(Rgb,((-1.0)/miu)))); 
    92  
    93     const double t4 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    94  
    95     const double t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - 
    96          b*p2*pow(q0,((-1.0) - p1 - p2)))); 
    97  
    98     const double t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + 
    99          (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 
    100          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 
    101          (7.0*b2)/(15.0*Rg02)))*Rg2)/b))); 
    102  
    103     const double t7 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    104          C2*pow(((Rgb)),((-2.0)/miu)) + 
    105          C1*pow(((Rgb)),((-1.0)/miu))))*pow(sech_WR(((-C4) + 
    106          Rgb)/C5),2)); 
    107  
    108     const double t8 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    109          Rg02/b2))*pow(sech_WR(((-C4) + Rgb)/C5),2)); 
    110  
    111     const double t9 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    112          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))*((1.0 + 0.5*(((-1.0) - 
    113          tanh(((-C4) + Rgb)/C5)))))); 
    114  
    115     const double t10 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    116          Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    117          Rgb)/C5)))))); 
    118  
    119     const double t11 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 
    120           3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)),((-1.0) - 
    121           2.0/miu)))/miu - (C1*Rg*pow(((Rgb)),((-1.0) - 
    122           1.0/miu)))/miu)); 
    123  
    124     const double t12 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    125  
    126     const double t13 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 
    127           (7.0*b2)/(15.0*q02* Rg2))) + 
    128           (7.0*b2)/(15.0*Rg02)))); 
    129  
    130     const double t14 = (2.0*b4*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    131           Rg02/b2))*((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    132           Rgb)/C5)))))); 
    133  
    134     const double t15 = ((C3*pow(((Rgb)),((-3.0)/miu)) + 
    135         C2*pow(((Rgb)),((-2.0)/miu)) + 
    136         C1*pow(((Rgb)),((-1.0)/miu)))); 
    137  
    138  
    139     double yy = (pow(q0,p1)*(((-((b*M_PI)/(L*q0))) +t1/L +t2/(q04*Rg22) + 
    140         0.5*t3*t4)) + (t5*((pow(q0,(p1 - p2))* 
    141         (((-pow(q0,(-p1)))*(((b2*M_PI)/(L*q02) +t6/L +t7/(2.0*C5) - 
    142         t8/(C5*q04*Rg22) + t9/(q04*Rg22) -t10/(q05*Rg22) + 0.5*t11*t12)) - 
    143         b*p1*pow(q0,((-1.0) - p1))*(((-((b*M_PI)/(L*q0))) + t13/L + 
    144         t14/(q04*Rg22) + 0.5*t15*((1.0 + tanh(((-C4) + 
    145         Rgb)/C5))))))))))); 
    146  
    147     return (yy); 
    148 } 
    149  
    150 static double 
    151 a2long(double q, double L, double b, double p1, double p2, double q0) 
    152 { 
    153     double C; 
    154  
    155     if( L/b > 10.0) { 
    156         C = 3.06/pow((L/b),0.44); 
    157     } else { 
    158         C = 1.0; 
    159     } 
    160  
    161     const double C1 = 1.22; 
    162     const double C2 = 0.4288; 
    163     const double C3 = -1.651; 
    164     const double C4 = 1.523; 
    165     const double C5 = 0.1477; 
    166     const double miu = 0.585; 
    167  
    168     const double Rg2 = Rgsquare(q,L,b); 
    169     const double Rg22 = Rg2*Rg2; 
    170     const double b2 = b*b; 
    171     const double b3 = b*b*b; 
    172     const double b4 = b3*b; 
    173     const double q02 = q0*q0; 
    174     const double q03 = q0*q0*q0; 
    175     const double q04 = q03*q0; 
    176     const double q05 = q04*q0; 
    177     const double Rg = sqrt(Rg2); 
    178     const double Rgb = Rg*q0/b; 
    179     const double Rg02 = Rg2*q02; 
    180  
    181     const double t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - 
    182          b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
    183  
    184     const double t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + 
    185          (14.0*b3*pow((double)M_E,(-(Rg02/b2))))/(15.0*q03*Rg2) + 
    186          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*((11.0/15.0 + 
    187          (7*b2)/(15.0*Rg02)))*Rg2)/b)))/L; 
    188  
    189     const double t3 = (Rg*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    190          C2*pow(((Rgb)),((-2.0)/miu)) + 
    191          C1*pow(((Rgb)),((-1.0)/miu))))* 
    192          pow(sech_WR(((-C4) +Rgb)/C5),2.0))/(2.0*C5); 
    193  
    194     const double t4 = (b4*Rg*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    195          Rg02/b2))*pow(sech_WR(((-C4) + 
    196          Rgb)/C5),2))/(C5*q04*Rg22); 
    197  
    198     const double t5 = (2.0*b4*(((2.0*q0*Rg2)/b - 
    199          (2.0*pow((double)M_E,(-(Rg02/b2)))*q0*Rg2)/b))* 
    200          ((1.0 + 0.5*(((-1.0) - tanh(((-C4) + 
    201          Rgb)/C5))))))/(q04*Rg22); 
    202  
    203     const double t6 = (8.0*b4*b*(((-1.0) + pow((double)M_E,(-(Rg02/b2))) + 
    204          Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    205          Rgb)/C5))))))/(q05*Rg22); 
    206  
    207     const double t7 = (((-((3.0*C3*Rg*pow(((Rgb)),((-1.0) - 
    208          3.0/miu)))/miu)) - (2.0*C2*Rg*pow(((Rgb)), 
    209          ((-1.0) - 2.0/miu)))/miu - (C1*Rg*pow(((Rgb)), 
    210          ((-1.0) - 1.0/miu)))/miu)); 
    211  
    212     const double t8 = ((1.0 + tanh(((-C4) + Rgb)/C5))); 
    213  
    214     const double t9 = (b*C*((4.0/15.0 - pow((double)M_E,(-(Rg02/b2)))*((11.0/15.0 + 
    215          (7.0*b2)/(15*Rg02))) + (7.0*b2)/(15.0*Rg02))))/L; 
    216  
    217     const double t10 = (2.0*b4*(((-1) + pow((double)M_E,(-(Rg02/b2))) + 
    218           Rg02/b2))*((1.0 + 0.5*(((-1) - tanh(((-C4) + 
    219           Rgb)/C5))))))/(q04*Rg22); 
    220  
    221     const double yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*M_PI)/(L*q02) + 
    222          t2 + t3 - t4 + t5 - t6 + 0.5*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))* 
    223          (((-((b*M_PI)/(L*q0))) + t9 + t10 + 
    224          0.5*((C3*pow(((Rgb)),((-3.0)/miu)) + 
    225          C2*pow(((Rgb)),((-2.0)/miu)) + 
    226          C1*pow(((Rgb)),((-1.0)/miu))))* 
    227          ((1.0 + tanh(((-C4) + Rgb)/C5)))))))))); 
    228  
    229     return (yy); 
    230 } 
    231  
    232 // 
    233 static double 
    234 a_short(double q, double L, double b, double p1short, 
    235         double p2short, double factor, double pdiff, double q0) 
    236 { 
    237     const double Rg2_sh = Rgsquareshort(q,L,b); 
    238     const double Rg2_sh2 = Rg2_sh*Rg2_sh; 
    239     const double b3 = b*b*b; 
    240     const double t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    241     const double Et1 = exp(t1); 
    242     const double Emt1 = 1.0/Et1; 
    243     const double q02 = q0*q0; 
    244     const double q0p = pow(q0,(-4.0 + p1short) ); 
    245  
    246     double yy = ((factor/(L*(pdiff)*Rg2_sh2)* 
    247         ((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 
    248         2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 
    249         2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*M_PI*q02*q0*Rg2_sh2 + 
    250         Et1*p2short*M_PI*q02*q0*Rg2_sh2)))))); 
    251  
    252     return(yy); 
    253 } 
    254 static double 
    255 a1short(double q, double L, double b, double p1short, double p2short, double q0) 
    256 { 
    257  
    258     double factor = 1.0; 
    259     return a_short(q, L, b, p1short, p2short, factor, p1short - p2short, q0); 
    260 } 
    261  
    262 static double 
    263 a2short(double q, double L, double b, double p1short, double p2short, double q0) 
    264 { 
    265     double factor = -1.0; 
    266     return a_short(q, L, b, p2short, p1short, factor, p1short-p2short, q0); 
    267 } 
    268  
    269 //WR named this w (too generic) 
    270 static double 
    271 w_WR(double x) 
    272 { 
    273     return 0.5*(1 + tanh((x - 1.523)/0.1477)); 
    274 } 
    275  
    276 // 
    277 static double 
    278 u1(double q, double L, double b) 
    279 { 
    280     return Rgsquareshort(q,L,b)*q*q; 
    281 } 
    282  
    283 static double 
    284 u_WR(double q, double L, double b) 
    285 { 
    286     return Rgsquare(q,L,b)*q*q; 
    287 } 
    288  
    289 static double 
    290 Sdebye_kernel(double arg) 
    291 { 
    292     // ORIGINAL 
    293     double result = 2.0*(exp(-arg) + arg -1.0)/(arg*arg); 
    294  
    295     // CONVERSION 1 from http://herbie.uwplse.org/ 
    296     // 
    297     // exhibits discontinuity - needs more investigation 
    298     //double a1 = 1.0/6.0; 
    299     //double a2 = 1.0/72.0; 
    300     //double a3 = 1.0/24.0; 
    301     //double result = pow((1.0 - a1*arg - (a2+a3)*arg*arg), 2); 
    302  
    303     return result; 
    304 } 
    305 static double 
    306 Sdebye(double q, double L, double b) 
    307 { 
    308     double arg = u_WR(q,L,b); 
    309     return Sdebye_kernel(arg); 
    310 } 
    311  
    312 // 
    313 static double 
    314 Sdebye1(double q, double L, double b) 
    315 { 
    316     double arg = u1(q,L,b); 
    317     return Sdebye_kernel(arg); 
    318  
    319 } 
    320  
    321 // 
     110    const double q0p = pow(q0, -4.0 + p1short); 
     111 
     112    double yy = 1.0/(L*r2*r2) * b/exp_qr_b*q0p 
     113        * (8.0*b3*L 
     114           - 8.0*b3*exp_qr_b*L 
     115           + 2.0*b3*exp_qr_b*L*p2short 
     116           - 2.0*b*exp_qr_b*L*p2short*qr2 
     117           + 4.0*b*exp_qr_b*L*qr2 
     118           - 2.0*b3*L*p2short 
     119           + 4.0*b*L*qr2 
     120           - M_PI*exp_qr_b*qr2*q0*r2 
     121           + M_PI*exp_qr_b*p2short*qr2*q0*r2); 
     122 
     123    return yy; 
     124} 
     125static double 
     126a_short(double qp, double L, double b 
     127        /*double p1short, double p2short*/, double q0) 
     128{ 
     129    const double p1short = 5.36; 
     130    const double p2short = 5.62; 
     131 
     132    const double r2 = Rgsquareshort(L,b); 
     133    const double exp_qr_b = exp(r2*square(q0/b)); 
     134    const double pdiff = p1short - p2short; 
     135    const double a1 = _short(r2,exp_qr_b,L,b,p1short,p2short,q0)/pdiff; 
     136    const double a2= -_short(r2,exp_qr_b,L,b,p2short,p1short,q0)/pdiff; 
     137    const double ans = a1*pow(qp*b, -p1short) + a2*pow(qp*b, -p2short) + M_PI/(qp*L); 
     138    return ans; 
     139} 
     140 
     141 
    322142static double 
    323143Sexv(double q, double L, double b) 
    324144{ 
    325  
     145    // Pedersen eq 13, corrected by Chen eq A.5, swapping w and 1-w 
    326146    const double C1=1.22; 
    327147    const double C2=0.4288; 
    328148    const double C3=-1.651; 
    329149    const double miu = 0.585; 
    330     const double qRg = q*sqrt(Rgsquare(q,L,b)); 
    331     const double x = pow(qRg, -1.0/miu); 
    332  
    333     double yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 
    334             w_WR(qRg)*x*(C1 + x*(C2 + x*C3)); 
    335     return (yy); 
    336 } 
    337  
    338  
    339 static double 
    340 Sexvnew(double q, double L, double b) 
    341 { 
    342     double yy; 
    343  
    344     const double C1 =1.22; 
    345     const double C2 =0.4288; 
    346     const double C3 =-1.651; 
    347     const double miu = 0.585; 
     150    const double qr = q*sqrt(Rgsquare(L,b)); 
     151    const double qr_miu = pow(qr, -1.0/miu); 
     152    const double w = w_WR(qr); 
     153    const double t10 = Sdebye(qr*qr)*(1.0 - w); 
     154    const double t11 = ((C3*qr_miu + C2)*qr_miu + C1)*qr_miu; 
     155 
     156    return t10 + w*t11; 
     157} 
     158 
     159// Modified by Yun on Oct. 15, 
     160static double 
     161Sexv_new(double q, double L, double b) 
     162{ 
     163    const double qr = q*sqrt(Rgsquare(L,b)); 
     164    const double qr2 = qr*qr; 
     165    const double C = (L/b > 10.0) ? 3.06*pow(L/b, -0.44) : 1.0; 
     166    const double t9 = C*b/L * (4.0 - exp(-qr2) * (11.0 + 7.0/qr2) + 7.0/qr2)/15.0; 
     167 
     168    const double Sexv_orig = Sexv(q, L, b); 
     169 
     170    // calculating the derivative to decide on the correction (cutoff) term? 
     171    // Note: this is modified from WRs original code 
    348172    const double del=1.05; 
    349     const double qRg = q*sqrt(Rgsquare(q,L,b)); 
    350     const double x = pow(qRg, -1.0/miu); 
    351  
    352  
    353     //calculating the derivative to decide on the corection (cutoff) term? 
    354     // I have modified this from WRs original code 
    355     const double qdel = (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q); 
    356     const double C_star2 =(qdel >= 0.0) ? 0.0 : 1.0; 
    357  
    358     yy = (1.0 - w_WR(qRg))*Sdebye(q,L,b) + 
    359          C_star2*w_WR(qRg)* 
    360          x*(C1 + x*(C2 + x*C3)); 
    361  
    362     return (yy); 
    363 } 
    364  
    365 double Sk_WR(double q, double L, double b) 
    366 { 
    367     // 
    368     const double p1 = 4.12; 
    369     const double p2 = 4.42; 
    370     const double p1short = 5.36; 
    371     const double p2short = 5.62; 
    372     const double q0 = 3.1; 
    373  
    374     double q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
    375     double Sexvmodify, ans; 
    376  
    377     const double C =(L/b > 10.0) ? 3.06/pow((L/b),0.44) : 1.0; 
    378  
    379     if( L > 4*b ) { // Longer Chains 
    380        if (q*b <= 3.1) {   //Modified by Yun on Oct. 15, 
    381          Sexvmodify = Sexvnew(q, L, b); 
    382          ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - 
    383             (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
    384  
    385        } else { //q(i)*b > 3.1 
    386          ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + 
    387                a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + M_PI/(q*L); 
    388        } 
    389     } else { //L <= 4*b Shorter Chains 
    390        if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
    391          if (q*b<=0.01) { 
    392             ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
    393          } else { 
    394             ans = Sdebye1(q,L,b); 
    395          } 
    396        } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
    397          ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + 
    398                a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + 
    399                M_PI/(q*L); 
    400        } 
     173    const double qdel = (Sexv(q*del,L,b) - Sexv_orig)/(q*(del - 1.0)); 
     174 
     175    if (qdel < 0) { 
     176        //printf("branch A1-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     177        return t9 + Sexv_orig; 
     178    } else { 
     179        //printf("branch A2-%d q=%g L=%g b=%g\n", C==1.0, q, L, b); 
     180        const double w = w_WR(qr); 
     181        const double t10 = Sdebye(qr*qr)*(1.0 - w); 
     182        return t9 + t10; 
    401183    } 
    402  
    403   return(ans); 
    404 } 
     184} 
     185 
     186 
     187static double 
     188Sk_WR(double q, double L, double b) 
     189{ 
     190    const double Rg_short = sqrt(Rgsquareshort(L, b)); 
     191    double q0short = fmax(1.9/Rg_short, 3.0); 
     192    double ans; 
     193 
     194 
     195    if( L > 4*b ) { // L > 4*b : Longer Chains 
     196        if (q*b <= 3.1) { 
     197            ans = Sexv_new(q, L, b); 
     198        } else { //q(i)*b > 3.1 
     199            ans = a_long(q, L, b /*, p1, p2, q0*/); 
     200        } 
     201    } else { // L <= 4*b : Shorter Chains 
     202        if (q*b <= q0short) { // q*b <= fmax(1.9/Rg_short, 3) 
     203            //printf("branch C-%d q=%g L=%g b=%g\n", square(q*Rg_short)<DEBYE_CUTOFF, q, L, b); 
     204            // Note that q0short is usually 3, but it will be greater than 3 
     205            // small enough b, depending on the L/b ratio: 
     206            //     L/b == 1 => b < 2.37 
     207            //     L/b == 2 => b < 1.36 
     208            //     L/b == 3 => b < 1.00 
     209            //     L/b == 4 => b < 0.816 
     210            // 2017-10-01 pkienzle: moved low q approximation into Sdebye() 
     211            ans = Sdebye(square(q*Rg_short)); 
     212        } else {  // q*b > max(1.9/Rg_short, 3) 
     213            //printf("branch D q=%g L=%g b=%g\n", q, L, b); 
     214            ans = a_short(q, L, b /*, p1short, p2short*/, q0short); 
     215        } 
     216    } 
     217 
     218    return ans; 
     219} 

  • sasmodels/models/barbell.c

    r592343f rbecded3  
    1 double form_volume(double radius_bell, double radius, double length); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3         double radius_bell, double radius, double length); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5         double radius_bell, double radius, double length, 
    6         double theta, double phi); 
    7  
    81#define INVALID(v) (v.radius_bell < v.radius) 
    92 
    103//barbell kernel - same as dumbell 
    114static double 
    12 _bell_kernel(double q, double h, double radius_bell, 
    13              double half_length, double sin_alpha, double cos_alpha) 
     5_bell_kernel(double qab, double qc, double h, double radius_bell, 
     6             double half_length) 
    147{ 
    158    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2619    //    m = q R cos(alpha) 
    2720    //    b = q(L/2-h) cos(alpha) 
    28     const double m = q*radius_bell*cos_alpha; // cos argument slope 
    29     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    30     const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 
     21    const double m = radius_bell*qc; // cos argument slope 
     22    const double b = (half_length-h)*qc; // cos argument intercept 
     23    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    3124    double total = 0.0; 
    3225    for (int i = 0; i < 76; i++){ 
    3326        const double t = Gauss76Z[i]*zm + zb; 
    3427        const double radical = 1.0 - t*t; 
    35         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     28        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3629        const double Fq = cos(m*t + b) * radical * bj; 
    3730        total += Gauss76Wt[i] * Fq; 
    … …  
    4437 
    4538static double 
    46 _fq(double q, double h, 
    47     double radius_bell, double radius, double half_length, 
    48     double sin_alpha, double cos_alpha) 
     39_fq(double qab, double qc, double h, 
     40    double radius_bell, double radius, double half_length) 
    4941{ 
    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); 
     42    const double bell_fq = _bell_kernel(qab, qc, h, radius_bell, half_length); 
     43    const double bj = sas_2J1x_x(radius*qab); 
     44    const double si = sas_sinx_x(half_length*qc); 
    5345    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5446    const double Aq = bell_fq + cyl_fq; 
    … …  
    5648} 
    5749 
    58  
    59 double form_volume(double radius_bell, 
    60         double radius, 
    61         double length) 
     50static double 
     51form_volume(double radius_bell, 
     52    double radius, 
     53    double length) 
    6254{ 
    6355    // bell radius should never be less than radius when this is called 
    … …  
    7062} 
    7163 
    72 double Iq(double q, double sld, double solvent_sld, 
    73           double radius_bell, double radius, double length) 
     64static double 
     65Iq(double q, double sld, double solvent_sld, 
     66    double radius_bell, double radius, double length) 
    7467{ 
    7568    const double h = -sqrt(radius_bell*radius_bell - radius*radius); 
    … …  
    8477        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8578        SINCOS(alpha, sin_alpha, cos_alpha); 
    86         const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha); 
     79        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    8780        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8881    } 
    … …  
    9689 
    9790 
    98 double Iqxy(double qx, double qy, 
    99         double sld, double solvent_sld, 
    100         double radius_bell, double radius, double length, 
    101         double theta, double phi) 
     91static double 
     92Iqxy(double qab, double qc, 
     93    double sld, double solvent_sld, 
     94    double radius_bell, double radius, double length) 
    10295{ 
    103     double q, sin_alpha, cos_alpha; 
    104     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    105  
    10696    const double h = -sqrt(square(radius_bell) - square(radius)); 
    107     const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha); 
     97    const double Aq = _fq(qab, qc, h, radius_bell, radius, 0.5*length); 
    10898 
    10999    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/bcc_paracrystal.c

    r50beefe rea60e08  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 26-27-28, multiplied by |q| 
     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; 
    68 
    7 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _BCCeval(double Theta, double Phi, double temp1, double temp3); 
    9 double _sphereform(double q, double radius, double sld, double solvent_sld); 
     9#if 1 
     10    // Matsuoka 29-30-31 
     11    //     Z_k numerator: 1 - exp(a)^2 
     12    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) 
     13    // Rewriting numerator 
     14    //         => -(exp(2a) - 1) 
     15    //         => -expm1(2a) 
     16    // Rewriting denominator 
     17    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     18    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
     19    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     20    const double exp_arg = exp(arg); 
     21    const double Zq = -cube(expm1(2.0*arg)) 
     22        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     24          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     25 
     26#elif 0 
     27    // ** Alternate form, which perhaps is more approachable 
     28    //     Z_k numerator   => -[(exp(2a) - 1) / 2.exp(a)] 2.exp(a) 
     29    //                     => -[sinh(a)] exp(a) 
     30    //     Z_k denominator => [(exp(2a) + 1) / 2.exp(a) - cos(d a_k)] 2.exp(a) 
     31    //                     => [cosh(a) - cos(d a_k)] 2.exp(a) 
     32    //     => Z_k = -sinh(a) / [cosh(a) - cos(d a_k)] 
     33    //            = sinh(-a) / [cosh(-a) - cos(d a_k)] 
     34    // 
     35    // One more step leads to the form in sasview 3.x for 2d models 
     36    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)] 
     37    // 
     38    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     39    const double sinh_qd = sinh(arg); 
     40    const double cosh_qd = cosh(arg); 
     41    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1)) 
     42                    * sinh_qd/(cosh_qd - cos(dnn*a2)) 
     43                    * sinh_qd/(cosh_qd - cos(dnn*a3)); 
     44#else 
     45    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     46    const double tanh_qd = tanh(arg); 
     47    const double cosh_qd = cosh(arg); 
     48    const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd) 
     49                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd) 
     50                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd); 
     51#endif 
     52 
     53    return Zq; 
     54} 
    1055 
    1156 
    12 double _BCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    13  
    14         const double Da = d_factor*dnn; 
    15         const double temp1 = q*q*Da*Da; 
    16         const double temp3 = q*dnn; 
    17  
    18         double retVal = _BCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    19         return(retVal); 
     57// occupied volume fraction calculated from lattice symmetry and sphere radius 
     58static double 
     59bcc_volume_fraction(double radius, double dnn) 
     60{ 
     61    return 2.0*sphere_volume(sqrt(0.75)*radius/dnn); 
    2062} 
    2163 
    22 double _BCCeval(double Theta, double Phi, double temp1, double temp3) { 
    23  
    24         double result; 
    25         double sin_theta,cos_theta,sin_phi,cos_phi; 
    26         SINCOS(Theta, sin_theta, cos_theta); 
    27         SINCOS(Phi, sin_phi, cos_phi); 
    28  
    29         const double temp6 =  sin_theta; 
    30         const double temp7 =  sin_theta*cos_phi + sin_theta*sin_phi + cos_theta; 
    31         const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta; 
    32         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta; 
    33  
    34         const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9)); 
    35         result = cube(1.0 - (temp10*temp10))*temp6 
    36             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    38               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    39  
    40         return (result); 
    41 } 
    42  
    43 double form_volume(double radius){ 
     64static double 
     65form_volume(double radius) 
     66{ 
    4467    return sphere_volume(radius); 
    4568} 
    4669 
    4770 
    48 double Iq(double q, double dnn, 
    49   double d_factor, double radius, 
    50   double sld, double solvent_sld){ 
     71static double Iq(double q, double dnn, 
     72    double d_factor, double radius, 
     73    double sld, double solvent_sld) 
     74{ 
     75    // translate a point in [-1,1] to a point in [0, 2 pi] 
     76    const double phi_m = M_PI; 
     77    const double phi_b = M_PI; 
     78    // translate a point in [-1,1] to a point in [0, pi] 
     79    const double theta_m = M_PI_2; 
     80    const double theta_b = M_PI_2; 
    5181 
    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; 
     82    double outer_sum = 0.0; 
     83    for(int i=0; i<150; i++) { 
     84        double inner_sum = 0.0; 
     85        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     86        double sin_theta, cos_theta; 
     87        SINCOS(theta, sin_theta, cos_theta); 
     88        const double qc = q*cos_theta; 
     89        const double qab = q*sin_theta; 
     90        for(int j=0;j<150;j++) { 
     91            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     92            double sin_phi, cos_phi; 
     93            SINCOS(phi, sin_phi, cos_phi); 
     94            const double qa = qab*cos_phi; 
     95            const double qb = qab*sin_phi; 
     96            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     97            inner_sum += Gauss150Wt[j] * form; 
     98        } 
     99        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     100        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     101    } 
     102    outer_sum *= theta_m; 
     103    const double Zq = outer_sum/(4.0*M_PI); 
     104    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     105    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    84106} 
    85107 
    86108 
    87 double Iqxy(double qx, double qy, 
     109static double Iqxy(double qa, double qb, double qc, 
    88110    double dnn, double d_factor, double radius, 
    89     double sld, double solvent_sld, 
    90     double theta, double phi, double psi) 
     111    double sld, double solvent_sld) 
    91112{ 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    94  
    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; 
     113    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     114    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor); 
     115    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     116    return bcc_volume_fraction(radius, dnn) * Pq * Zq; 
    111117} 

  • sasmodels/models/capped_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double radius_cap, double length); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3     double radius, double radius_cap, double length); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5     double radius, double radius_cap, double length, double theta, double phi); 
    6  
    71#define INVALID(v) (v.radius_cap < v.radius) 
    82 
    … …  
    148//   radius_cap is the radius of the lens 
    159//   length is the cylinder length, or the separation between the lens halves 
    16 //   alpha is the angle of the cylinder wrt q. 
     10//   theta is the angle of the cylinder wrt q. 
    1711static double 
    18 _cap_kernel(double q, double h, double radius_cap, 
    19                       double half_length, double sin_alpha, double cos_alpha) 
     12_cap_kernel(double qab, double qc, double h, double radius_cap, 
     13    double half_length) 
    2014{ 
    2115    // translate a point in [-1,1] to a point in [lower,upper] 
    … …  
    2620 
    2721    // cos term in integral is: 
    28     //    cos (q (R t - h + L/2) cos(alpha)) 
     22    //    cos (q (R t - h + L/2) cos(theta)) 
    2923    // so turn it into: 
    3024    //    cos (m t + b) 
    3125    // where: 
    32     //    m = q R cos(alpha) 
    33     //    b = q(L/2-h) cos(alpha) 
    34     const double m = q*radius_cap*cos_alpha; // cos argument slope 
    35     const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 
    36     const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 
     26    //    m = q R cos(theta) 
     27    //    b = q(L/2-h) cos(theta) 
     28    const double m = radius_cap*qc; // cos argument slope 
     29    const double b = (half_length-h)*qc; // cos argument intercept 
     30    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3731    double total = 0.0; 
    3832    for (int i=0; i<76 ;i++) { 
    3933        const double t = Gauss76Z[i]*zm + zb; 
    4034        const double radical = 1.0 - t*t; 
    41         const double bj = sas_2J1x_x(qrst*sqrt(radical)); 
     35        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    4236        const double Fq = cos(m*t + b) * radical * bj; 
    4337        total += Gauss76Wt[i] * Fq; 
    … …  
    5044 
    5145static double 
    52 _fq(double q, double h, double radius_cap, double radius, double half_length, 
    53     double sin_alpha, double cos_alpha) 
     46_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length) 
    5447{ 
    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); 
     48    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length); 
     49    const double bj = sas_2J1x_x(radius*qab); 
     50    const double si = sas_sinx_x(half_length*qc); 
    5851    const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 
    5952    const double Aq = cap_Fq + cyl_Fq; 
    … …  
    6154} 
    6255 
    63 double form_volume(double radius, double radius_cap, double length) 
     56static double 
     57form_volume(double radius, double radius_cap, double length) 
    6458{ 
    6559    // cap radius should never be less than radius when this is called 
    … …  
    9084} 
    9185 
    92 double Iq(double q, double sld, double solvent_sld, 
    93           double radius, double radius_cap, double length) 
     86static double 
     87Iq(double q, double sld, double solvent_sld, 
     88    double radius, double radius_cap, double length) 
    9489{ 
    9590    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    … …  
    10196    double total = 0.0; 
    10297    for (int i=0; i<76 ;i++) { 
    103         const double alpha= Gauss76Z[i]*zm + zb; 
    104         double sin_alpha, cos_alpha; // slots to hold sincos function output 
    105         SINCOS(alpha, sin_alpha, cos_alpha); 
    106  
    107         const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha); 
    108         // sin_alpha for spherical coord integration 
    109         total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     98        const double theta = Gauss76Z[i]*zm + zb; 
     99        double sin_theta, cos_theta; // slots to hold sincos function output 
     100        SINCOS(theta, sin_theta, cos_theta); 
     101        const double qab = q*sin_theta; 
     102        const double qc = q*cos_theta; 
     103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
     104        // scale by sin_theta for spherical coord integration 
     105        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    110106    } 
    111107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    118114 
    119115 
    120 double Iqxy(double qx, double qy, 
     116static double 
     117Iqxy(double qab, double qc, 
    121118    double sld, double solvent_sld, double radius, 
    122     double radius_cap, double length, 
    123     double theta, double phi) 
     119    double radius_cap, double length) 
    124120{ 
    125     double q, sin_alpha, cos_alpha; 
    126     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    127  
    128121    const double h = sqrt(radius_cap*radius_cap - radius*radius); 
    129     const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha); 
     122    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length); 
    130123 
    131124    // Multiply by contrast^2 and convert to cm-1 

  • sasmodels/models/core_shell_bicelle.c

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

  • sasmodels/models/core_shell_bicelle_elliptical.c

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

  • sasmodels/models/core_shell_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double core_sld, double shell_sld, double solvent_sld, 
    3     double radius, double thickness, double length); 
    4 double Iqxy(double qx, double qy, double core_sld, double shell_sld, double solvent_sld, 
    5     double radius, double thickness, double length, double theta, double phi); 
    6  
    71// vd = volume * delta_rho 
    8 // besarg = q * R * sin(alpha) 
    9 // siarg = q * L/2 * cos(alpha) 
    10 double _cyl(double vd, double besarg, double siarg); 
    11 double _cyl(double vd, double besarg, double siarg) 
     2// besarg = q * R * sin(theta) 
     3// siarg = q * L/2 * cos(theta) 
     4static double _cyl(double vd, double besarg, double siarg) 
    125{ 
    136    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 
    147} 
    158 
    16 double form_volume(double radius, double thickness, double length) 
     9static double 
     10form_volume(double radius, double thickness, double length) 
    1711{ 
    18     return M_PI*(radius+thickness)*(radius+thickness)*(length+2.0*thickness); 
     12    return M_PI*square(radius+thickness)*(length+2.0*thickness); 
    1913} 
    2014 
    21 double Iq(double q, 
     15static double 
     16Iq(double q, 
    2217    double core_sld, 
    2318    double shell_sld, 
    … …  
    2823{ 
    2924    // precalculate constants 
    30     const double core_qr = q*radius; 
    31     const double core_qh = q*0.5*length; 
     25    const double core_r = radius; 
     26    const double core_h = 0.5*length; 
    3227    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    33     const double shell_qr = q*(radius + thickness); 
    34     const double shell_qh = q*(0.5*length + thickness); 
     28    const double shell_r = (radius + thickness); 
     29    const double shell_h = (0.5*length + thickness); 
    3530    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3631    double total = 0.0; 
    37     // double lower=0, upper=M_PI_2; 
    3832    for (int i=0; i<76 ;i++) { 
    39         // translate a point in [-1,1] to a point in [lower,upper] 
    40         //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    41         double sn, cn; 
    42         const double alpha = 0.5*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    43         SINCOS(alpha, sn, cn); 
    44         const double fq = _cyl(core_vd, core_qr*sn, core_qh*cn) 
    45             + _cyl(shell_vd, shell_qr*sn, shell_qh*cn); 
    46         total += Gauss76Wt[i] * fq * fq * sn; 
     33        // translate a point in [-1,1] to a point in [0, pi/2] 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     35        double sin_theta, cos_theta; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     37        SINCOS(theta, sin_theta,  cos_theta); 
     38        const double qab = q*sin_theta; 
     39        const double qc = q*cos_theta; 
     40        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     41            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4743    } 
    4844    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5248 
    5349 
    54 double Iqxy(double qx, double qy, 
     50double Iqxy(double qab, double qc, 
    5551    double core_sld, 
    5652    double shell_sld, 
    … …  
    5854    double radius, 
    5955    double thickness, 
    60     double length, 
    61     double theta, 
    62     double phi) 
     56    double length) 
    6357{ 
    64     double q, sin_alpha, cos_alpha; 
    65     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    66  
    67     const double core_qr = q*radius; 
    68     const double core_qh = q*0.5*length; 
     58    const double core_r = radius; 
     59    const double core_h = 0.5*length; 
    6960    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld); 
    70     const double shell_qr = q*(radius + thickness); 
    71     const double shell_qh = q*(0.5*length + thickness); 
     61    const double shell_r = (radius + thickness); 
     62    const double shell_h = (0.5*length + thickness); 
    7263    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    7364 
    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); 
     65    const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
     66        + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    7667    return 1.0e-4 * fq * fq; 
    7768} 

  • sasmodels/models/core_shell_ellipsoid.c

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

  • sasmodels/models/core_shell_ellipsoid.py

    r30b60d2 r8db25bf  
    141141# pylint: enable=bad-whitespace, line-too-long 
    142142 
    143 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 
    144           "core_shell_ellipsoid.c"] 
     143source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    145144 
    146145def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c rbecded3  
    1 double form_volume(double length_a, double length_b, double length_c,  
    2                    double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    3 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    4           double solvent_sld, double length_a, double length_b, double length_c, 
    5           double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    6 double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, 
    7             double crim_sld, double solvent_sld, double length_a, double length_b, 
    8             double length_c, double thick_rim_a, double thick_rim_b, 
    9             double thick_rim_c, double theta, double phi, double psi); 
    10  
    11 double form_volume(double length_a, double length_b, double length_c,  
    12                    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     1static double 
     2form_volume(double length_a, double length_b, double length_c, 
     3    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    134{ 
    145    //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 +  
     6    return length_a * length_b * length_c + 
     7           2.0 * thick_rim_a * length_b * length_c + 
    178           2.0 * thick_rim_b * length_a * length_c + 
    189           2.0 * thick_rim_c * length_a * length_b; 
    1910} 
    2011 
    21 double Iq(double q, 
     12static double 
     13Iq(double q, 
    2214    double core_sld, 
    2315    double arim_sld, 
    … …  
    3426    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3527    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     28 
    3729    const double mu = 0.5 * q * length_b; 
    38      
     30 
    3931    //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 +  
     32    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     33    //double vol = length_a * length_b * length_c + 
     34    //       2.0 * thick_rim_a * length_b * length_c + 
    4335    //       2.0 * thick_rim_b * length_a * length_c + 
    4436    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     37 
    4638    // Scale sides by B 
    4739    const double a_scaled = length_a / length_b; 
    … …  
    10193            //   ( 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 
    10294            // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     95            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c 
     96 
    10597            //  correct FF : sum of square of phase factors 
    10698            inner_total += Gauss76Wt[j] * form * form; 
    … …  
    118110} 
    119111 
    120 double Iqxy(double qx, double qy, 
     112static double 
     113Iqxy(double qa, double qb, double qc, 
    121114    double core_sld, 
    122115    double arim_sld, 
    … …  
    129122    double thick_rim_a, 
    130123    double thick_rim_b, 
    131     double thick_rim_c, 
    132     double theta, 
    133     double phi, 
    134     double psi) 
     124    double thick_rim_c) 
    135125{ 
    136     double q, zhat, yhat, xhat; 
    137     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    138  
    139126    // cspkernel in csparallelepiped recoded here 
    140127    const double dr0 = core_sld-solvent_sld; 
    … …  
    160147    double tc = length_a + 2.0*thick_rim_c; 
    161148    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162     double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    163     double siB = sas_sinx_x(0.5*q*length_b*yhat); 
    164     double siC = sas_sinx_x(0.5*q*length_c*zhat); 
    165     double siAt = sas_sinx_x(0.5*q*ta*xhat); 
    166     double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167     double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     149    double siA = sas_sinx_x(0.5*length_a*qa); 
     150    double siB = sas_sinx_x(0.5*length_b*qb); 
     151    double siC = sas_sinx_x(0.5*length_c*qc); 
     152    double siAt = sas_sinx_x(0.5*ta*qa); 
     153    double siBt = sas_sinx_x(0.5*tb*qb); 
     154    double siCt = sas_sinx_x(0.5*tc*qc); 
     155 
    169156 
    170157    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    174161               + drB*siA*(siBt-siB)*siC*V2 
    175162               + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     163 
    177164    return 1.0e-4 * f * f; 
    178165} 

  • sasmodels/models/cylinder.c

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

  • sasmodels/models/ellipsoid.c

    r3b571ae rbecded3  
    1 double form_volume(double radius_polar, double radius_equatorial); 
    2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 
    3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    4     double radius_polar, double radius_equatorial, double theta, double phi); 
    5  
    6 double form_volume(double radius_polar, double radius_equatorial) 
     1static double 
     2form_volume(double radius_polar, double radius_equatorial) 
    73{ 
    84    return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 
    95} 
    106 
    11 double Iq(double q, 
     7static  double 
     8Iq(double q, 
    129    double sld, 
    1310    double sld_solvent, 
    … …  
    4138} 
    4239 
    43 double Iqxy(double qx, double qy, 
     40static double 
     41Iqxy(double qab, double qc, 
    4442    double sld, 
    4543    double sld_solvent, 
    4644    double radius_polar, 
    47     double radius_equatorial, 
    48     double theta, 
    49     double phi) 
     45    double radius_equatorial) 
    5046{ 
    51     double q, sin_alpha, cos_alpha; 
    52     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    53     const double r = sqrt(square(radius_equatorial*sin_alpha) 
    54                           + square(radius_polar*cos_alpha)); 
    55     const double f = sas_3j1x_x(q*r); 
     47    const double qr = sqrt(square(radius_equatorial*qab) + square(radius_polar*qc)); 
     48    const double f = sas_3j1x_x(qr); 
    5649    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    5750 
    5851    return 1.0e-4 * square(f * s); 
    5952} 
    60  

  • sasmodels/models/elliptical_cylinder.c

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

  • sasmodels/models/fcc_paracrystal.c

    r50beefe rf728001  
    1 double form_volume(double radius); 
    2 double Iq(double q,double dnn,double d_factor, double radius,double sld, double solvent_sld); 
    3 double Iqxy(double qx, double qy, double dnn, 
    4     double d_factor, double radius,double sld, double solvent_sld, 
    5     double theta, double phi, double psi); 
     1static double 
     2fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 17-18-19, multiplied by |q| 
     5    const double a1 = ( qa + qb)/2.0; 
     6    const double a2 = ( qa + qc)/2.0; 
     7    const double a3 = ( qb + qc)/2.0; 
    68 
    7 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi); 
    8 double _FCCeval(double Theta, double Phi, double temp1, double temp3); 
     9    // Matsuoka 23-24-25 
     10    //     Z_k numerator: 1 - exp(a)^2 
     11    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) 
     12    // Rewriting numerator 
     13    //         => -(exp(2a) - 1) 
     14    //         => -expm1(2a) 
     15    // Rewriting denominator 
     16    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     17    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
     18    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     19    const double exp_arg = exp(arg); 
     20    const double Zq = -cube(expm1(2.0*arg)) 
     21        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
     24 
     25    return Zq; 
     26} 
    927 
    1028 
    11 double _FCC_Integrand(double q, double dnn, double d_factor, double theta, double phi) { 
    12  
    13         const double Da = d_factor*dnn; 
    14         const double temp1 = q*q*Da*Da; 
    15         const double temp3 = q*dnn; 
    16  
    17         double retVal = _FCCeval(theta,phi,temp1,temp3)/(4.0*M_PI); 
    18         return(retVal); 
     29// occupied volume fraction calculated from lattice symmetry and sphere radius 
     30static double 
     31fcc_volume_fraction(double radius, double dnn) 
     32{ 
     33    return 4.0*sphere_volume(M_SQRT1_2*radius/dnn); 
    1934} 
    2035 
    21 double _FCCeval(double Theta, double Phi, double temp1, double temp3) { 
    22  
    23         double result; 
    24         double sin_theta,cos_theta,sin_phi,cos_phi; 
    25         SINCOS(Theta, sin_theta, cos_theta); 
    26         SINCOS(Phi, sin_phi, cos_phi); 
    27  
    28         const double temp6 =  sin_theta; 
    29         const double temp7 =  sin_theta*sin_phi + cos_theta; 
    30         const double temp8 = -sin_theta*cos_phi + cos_theta; 
    31         const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi; 
    32  
    33         const double temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
    34         result = cube(1.0-(temp10*temp10))*temp6 
    35             / ( (1.0 - 2.0*temp10*cos(0.5*temp3*temp7) + temp10*temp10) 
    36               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp8) + temp10*temp10) 
    37               * (1.0 - 2.0*temp10*cos(0.5*temp3*temp9) + temp10*temp10)); 
    38  
    39         return (result); 
    40 } 
    41  
    42 double form_volume(double radius){ 
     36static double 
     37form_volume(double radius) 
     38{ 
    4339    return sphere_volume(radius); 
    4440} 
    4541 
    4642 
    47 double Iq(double q, double dnn, 
     43static double Iq(double q, double dnn, 
    4844  double d_factor, double radius, 
    49   double sld, double solvent_sld){ 
     45  double sld, double solvent_sld) 
     46{ 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI; 
     49    const double phi_b = M_PI; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_2; 
     52    const double theta_b = M_PI_2; 
    5053 
    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); 
     54    double outer_sum = 0.0; 
     55    for(int i=0; i<150; i++) { 
     56        double inner_sum = 0.0; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        double sin_theta, cos_theta; 
     59        SINCOS(theta, sin_theta, cos_theta); 
     60        const double qc = q*cos_theta; 
     61        const double qab = q*sin_theta; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     64            double sin_phi, cos_phi; 
     65            SINCOS(phi, sin_phi, cos_phi); 
     66            const double qa = qab*cos_phi; 
     67            const double qb = qab*sin_phi; 
     68            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     69            inner_sum += Gauss150Wt[j] * form; 
     70        } 
     71        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73    } 
     74    outer_sum *= theta_m; 
     75    const double Zq = outer_sum/(4.0*M_PI); 
     76    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
    5477 
    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; 
     78    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
     79} 
    8380 
    8481 
     82static double Iqxy(double qa, double qb, double qc, 
     83    double dnn, double d_factor, double radius, 
     84    double sld, double solvent_sld) 
     85{ 
     86    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     87    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     88    const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor); 
     89    return fcc_volume_fraction(radius, dnn) * Pq * Zq; 
    8590} 
    86  
    87 double Iqxy(double qx, double qy, 
    88     double dnn, double d_factor, double radius, 
    89     double sld, double solvent_sld, 
    90     double theta, double phi, double psi) 
    91 { 
    92     double q, zhat, yhat, xhat; 
    93     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    94  
    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; 
    112 } 

  • sasmodels/models/hollow_cylinder.c

    r592343f rbecded3  
    1 double form_volume(double radius, double thickness, double length); 
    2 double Iq(double q, double radius, double thickness, double length, double sld, 
    3         double solvent_sld); 
    4 double Iqxy(double qx, double qy, double radius, double thickness, double length, double sld, 
    5         double solvent_sld, double theta, double phi); 
    6  
    71//#define INVALID(v) (v.radius_core >= v.radius) 
    82 
    … …  
    148} 
    159 
    16  
    1710static double 
    18 _hollow_cylinder_kernel(double q, 
    19     double radius, double thickness, double length, double sin_val, double cos_val) 
     11_fq(double qab, double qc, 
     12    double radius, double thickness, double length) 
    2013{ 
    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); 
     14    const double lam1 = sas_2J1x_x((radius+thickness)*qab); 
     15    const double lam2 = sas_2J1x_x(radius*qab); 
    2416    const double gamma_sq = square(radius/(radius+thickness)); 
    25     //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 
    26     //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 
    27     const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
    28     const double t2 = sas_sinx_x(0.5*q*length*cos_val); 
     17    //Note: lim_{thickness -> 0} psi = sas_J0(radius*qab) 
     18    //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qab) 
     19    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq);    //SRK 10/19/00 
     20    const double t2 = sas_sinx_x(0.5*length*qc); 
    2921    return psi*t2; 
    3022} 
    3123 
    32 double 
     24static double 
    3325form_volume(double radius, double thickness, double length) 
    3426{ 
    … …  
    3830 
    3931 
    40 double 
     32static double 
    4133Iq(double q, double radius, double thickness, double length, 
    4234    double sld, double solvent_sld) 
    4335{ 
    4436    const double lower = 0.0; 
    45     const double upper = 1.0;           //limits of numerical integral 
     37    const double upper = 1.0;        //limits of numerical integral 
    4638 
    47     double summ = 0.0;                  //initialize intergral 
     39    double summ = 0.0;            //initialize intergral 
    4840    for (int i=0;i<76;i++) { 
    49         const double cos_val = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    50         const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    51         const double inter = _hollow_cylinder_kernel(q, radius, thickness, length, 
    52                                                      sin_val, cos_val); 
    53         summ += Gauss76Wt[i] * inter * inter; 
     41        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
     42        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     43        const double form = _fq(q*sin_theta, q*cos_theta, 
     44                                radius, thickness, length); 
     45        summ += Gauss76Wt[i] * form * form; 
    5446    } 
    5547 
    … …  
    5951} 
    6052 
    61 double 
    62 Iqxy(double qx, double qy, 
     53static double 
     54Iqxy(double qab, double qc, 
    6355    double radius, double thickness, double length, 
    64     double sld, double solvent_sld, double theta, double phi) 
     56    double sld, double solvent_sld) 
    6557{ 
    66     double q, sin_alpha, cos_alpha; 
    67     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    68     const double Aq = _hollow_cylinder_kernel(q, radius, thickness, length, 
    69         sin_alpha, cos_alpha); 
     58    const double form = _fq(qab, qc, radius, thickness, length); 
    7059 
    7160    const double vol = form_volume(radius, thickness, length); 
    72     return _hollow_cylinder_scaling(Aq*Aq, solvent_sld-sld, vol); 
     61    return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
    7362} 
    74  

  • sasmodels/models/parallelepiped.c

    rd605080 r9b7b23f  
    1 double form_volume(double length_a, double length_b, double length_c); 
    2 double Iq(double q, double sld, double solvent_sld, 
    3     double length_a, double length_b, double length_c); 
    4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    5     double length_a, double length_b, double length_c, 
    6     double theta, double phi, double psi); 
    7  
    8 double form_volume(double length_a, double length_b, double length_c) 
     1static double 
     2form_volume(double length_a, double length_b, double length_c) 
    93{ 
    104    return length_a * length_b * length_c; 
    … …  
    126 
    137 
    14 double Iq(double q, 
     8static double 
     9Iq(double q, 
    1510    double sld, 
    1611    double solvent_sld, 
    … …  
    2015{ 
    2116    const double mu = 0.5 * q * length_b; 
    22      
     17 
    2318    // Scale sides by B 
    2419    const double a_scaled = length_a / length_b; 
    2520    const double c_scaled = length_c / length_b; 
    26          
     21 
    2722    // outer integral (with gauss points), integration limits = 0, 1 
    2823    double outer_total = 0; //initialize integral 
    … …  
    5752 
    5853 
    59 double Iqxy(double qx, double qy, 
     54static double 
     55Iqxy(double qa, double qb, double qc, 
    6056    double sld, 
    6157    double solvent_sld, 
    6258    double length_a, 
    6359    double length_b, 
    64     double length_c, 
    65     double theta, 
    66     double phi, 
    67     double psi) 
     60    double length_c) 
    6861{ 
    69     double q, xhat, yhat, zhat; 
    70     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    71  
    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); 
     62    const double siA = sas_sinx_x(0.5*length_a*qa); 
     63    const double siB = sas_sinx_x(0.5*length_b*qb); 
     64    const double siC = sas_sinx_x(0.5*length_c*qc); 
    7565    const double V = form_volume(length_a, length_b, length_c); 
    7666    const double drho = (sld - solvent_sld); 

  • sasmodels/models/sc_paracrystal.c

    r50beefe rf728001  
    1 double form_volume(double radius); 
     1static double 
     2sc_Zq(double qa, double qb, double qc, double dnn, double d_factor) 
     3{ 
     4    // Equations from Matsuoka 9-10-11, multiplied by |q| 
     5    const double a1 = qa; 
     6    const double a2 = qb; 
     7    const double a3 = qc; 
    28 
    3 double Iq(double q, 
    4           double dnn, 
    5           double d_factor, 
    6           double radius, 
    7           double sphere_sld, 
    8           double solvent_sld); 
     9    // Matsuoka 13-14-15 
     10    //     Z_k numerator: 1 - exp(a)^2 
     11    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a) 
     12    // Rewriting numerator 
     13    //         => -(exp(2a) - 1) 
     14    //         => -expm1(2a) 
     15    // Rewriting denominator 
     16    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 
     17    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1 
     18    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3); 
     19    const double exp_arg = exp(arg); 
     20    const double Zq = -cube(expm1(2.0*arg)) 
     21        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0) 
     22          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0) 
     23          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0)); 
    924 
    10 double Iqxy(double qx, double qy, 
    11             double dnn, 
    12             double d_factor, 
    13             double radius, 
    14             double sphere_sld, 
    15             double solvent_sld, 
    16             double theta, 
    17             double phi, 
    18             double psi); 
     25    return Zq; 
     26} 
    1927 
    20 double form_volume(double radius) 
     28// occupied volume fraction calculated from lattice symmetry and sphere radius 
     29static double 
     30sc_volume_fraction(double radius, double dnn) 
     31{ 
     32    return sphere_volume(radius/dnn); 
     33} 
     34 
     35static double 
     36form_volume(double radius) 
    2137{ 
    2238    return sphere_volume(radius); 
    2339} 
    2440 
     41 
    2542static double 
    26 sc_eval(double theta, double phi, double temp3, double temp4, double temp5) 
     43Iq(double q, double dnn, 
     44    double d_factor, double radius, 
     45    double sld, double solvent_sld) 
    2746{ 
    28     double cnt, snt; 
    29     SINCOS(theta, cnt, snt); 
     47    // translate a point in [-1,1] to a point in [0, 2 pi] 
     48    const double phi_m = M_PI_4; 
     49    const double phi_b = M_PI_4; 
     50    // translate a point in [-1,1] to a point in [0, pi] 
     51    const double theta_m = M_PI_4; 
     52    const double theta_b = M_PI_4; 
    3053 
    31     double cnp, snp; 
    32     SINCOS(phi, cnp, snp); 
    3354 
    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); 
     55    double outer_sum = 0.0; 
     56    for(int i=0; i<150; i++) { 
     57        double inner_sum = 0.0; 
     58        const double theta = Gauss150Z[i]*theta_m + theta_b; 
     59        double sin_theta, cos_theta; 
     60        SINCOS(theta, sin_theta, cos_theta); 
     61        const double qc = q*cos_theta; 
     62        const double qab = q*sin_theta; 
     63        for(int j=0;j<150;j++) { 
     64            const double phi = Gauss150Z[j]*phi_m + phi_b; 
     65            double sin_phi, cos_phi; 
     66            SINCOS(phi, sin_phi, cos_phi); 
     67            const double qa = qab*cos_phi; 
     68            const double qb = qab*sin_phi; 
     69            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
     70            inner_sum += Gauss150Wt[j] * form; 
     71        } 
     72        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
     73        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     74    } 
     75    outer_sum *= theta_m; 
     76    const double Zq = outer_sum/M_PI_2; 
     77    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     78 
     79    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    4280} 
    4381 
     82 
    4483static double 
    45 sc_integrand(double dnn, double d_factor, double qq, double xx, double yy) 
     84Iqxy(double qa, double qb, double qc, 
     85    double dnn, double d_factor, double radius, 
     86    double sld, double solvent_sld) 
    4687{ 
    47     //Function to calculate integrand values for simple cubic structure 
    48  
    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); 
     88    const double q = sqrt(qa*qa + qb*qb + qc*qc); 
     89    const double Pq = sphere_form(q, radius, sld, solvent_sld); 
     90    const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor); 
     91    return sc_volume_fraction(radius, dnn) * Pq * Zq; 
    5992} 
    60  
    61 double Iq(double q, 
    62           double dnn, 
    63           double d_factor, 
    64           double radius, 
    65           double sphere_sld, 
    66           double solvent_sld) 
    67 { 
    68         const double va = 0.0; 
    69         const double vb = M_PI_2; //orientation average, outer integral 
    70  
    71     double summ=0.0; 
    72     double answer=0.0; 
    73  
    74         for(int i=0;i<150;i++) { 
    75                 //setup inner integral over the ellipsoidal cross-section 
    76                 double summj=0.0; 
    77                 double zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0; 
    78                 for(int j=0;j<150;j++) { 
    79                         //150 gauss points for the inner integral 
    80                         double zij = ( Gauss150Z[j]*(vb-va) + va + vb )/2.0; 
    81                         double tmp = Gauss150Wt[j] * sc_integrand(dnn,d_factor,q,zi,zij); 
    82                         summj += tmp; 
    83                 } 
    84                 //now calculate the value of the inner integral 
    85                 answer = (vb-va)/2.0*summj; 
    86  
    87                 //now calculate outer integral 
    88                 double tmp = Gauss150Wt[i] * answer; 
    89                 summ += tmp; 
    90         }               //final scaling is done at the end of the function, after the NT_FP64 case 
    91  
    92         answer = (vb-va)/2.0*summ; 
    93  
    94         //Volume fraction calculated from lattice symmetry and sphere radius 
    95         // NB: 4/3 pi r^3 / dnn^3 = 4/3 pi(r/dnn)^3 
    96         const double latticeScale = sphere_volume(radius/dnn); 
    97  
    98         answer *= sphere_form(q, radius, sphere_sld, solvent_sld)*latticeScale; 
    99  
    100         return answer; 
    101 } 
    102  
    103 double Iqxy(double qx, double qy, 
    104           double dnn, 
    105           double d_factor, 
    106           double radius, 
    107           double sphere_sld, 
    108           double solvent_sld, 
    109           double theta, 
    110           double phi, 
    111           double psi) 
    112 { 
    113     double q, zhat, yhat, xhat; 
    114     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
    115  
    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; 
    128 } 

  • sasmodels/models/sc_paracrystal.py

    r8f04da4 r9bc4882  
    161161    [{'theta': 10.0, 'phi': 20, 'psi': 30.0}, (0.023, 0.045), 0.0177333171285], 
    162162    ] 
    163  
    164  

  • sasmodels/models/stacked_disks.c

    r19f996b rbecded3  
    1 static double stacked_disks_kernel( 
    2     double q, 
     1static double 
     2stacked_disks_kernel( 
     3    double qab, 
     4    double qc, 
    35    double halfheight, 
    46    double thick_layer, 
    … …  
    911    double layer_sld, 
    1012    double solvent_sld, 
    11     double sin_alpha, 
    12     double cos_alpha, 
    1313    double d) 
    1414 
    … …  
    2020    // zi is the dummy variable for the integration (x in Feigin's notation) 
    2121 
    22     const double besarg1 = q*radius*sin_alpha; 
    23     //const double besarg2 = q*radius*sin_alpha; 
     22    const double besarg1 = radius*qab; 
     23    //const double besarg2 = radius*qab; 
    2424 
    25     const double sinarg1 = q*halfheight*cos_alpha; 
    26     const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 
     25    const double sinarg1 = halfheight*qc; 
     26    const double sinarg2 = (halfheight+thick_layer)*qc; 
    2727 
    2828    const double be1 = sas_2J1x_x(besarg1); 
    … …  
    4343 
    4444    // loop for the structure factor S(q) 
    45     double qd_cos_alpha = q*d*cos_alpha; 
     45    double qd_cos_alpha = d*qc; 
    4646    //d*cos_alpha is the projection of d onto q (in other words the component 
    4747    //of d that is parallel to q. 
    … …  
    6161 
    6262 
    63 static double stacked_disks_1d( 
     63static double 
     64stacked_disks_1d( 
    6465    double q, 
    6566    double thick_core, 
    … …  
    8485        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8586        SINCOS(zi, sin_alpha, cos_alpha); 
    86         double yyy = stacked_disks_kernel(q, 
     87        double yyy = stacked_disks_kernel(q*sin_alpha, q*cos_alpha, 
    8788                           halfheight, 
    8889                           thick_layer, 
    … …  
    9394                           layer_sld, 
    9495                           solvent_sld, 
    95                            sin_alpha, 
    96                            cos_alpha, 
    9796                           d); 
    9897        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    … …  
    105104} 
    106105 
    107 static double form_volume( 
     106static double 
     107form_volume( 
    108108    double thick_core, 
    109109    double thick_layer, 
    … …  
    116116} 
    117117 
    118 static double Iq( 
     118static double 
     119Iq( 
    119120    double q, 
    120121    double thick_core, 
    … …  
    140141 
    141142 
    142 static double Iqxy(double qx, double qy, 
     143static double 
     144Iqxy(double qab, double qc, 
    143145    double thick_core, 
    144146    double thick_layer, 
    … …  
    148150    double core_sld, 
    149151    double layer_sld, 
    150     double solvent_sld, 
    151     double theta, 
    152     double phi) 
     152    double solvent_sld) 
    153153{ 
    154154    int n_stacking = (int)(fp_n_stacking + 0.5); 
    155     double q, sin_alpha, cos_alpha; 
    156     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    157  
    158155    double d = 2.0 * thick_layer + thick_core; 
    159156    double halfheight = 0.5*thick_core; 
    160     double answer = stacked_disks_kernel(q, 
     157    double answer = stacked_disks_kernel(qab, qc, 
    161158                     halfheight, 
    162159                     thick_layer, 
    … …  
    167164                     layer_sld, 
    168165                     solvent_sld, 
    169                      sin_alpha, 
    170                      cos_alpha, 
    171166                     d); 
    172167 

  • sasmodels/models/triaxial_ellipsoid.c

    r68dd6a9 rbecded3  
    1 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar); 
    2 double Iq(double q, double sld, double sld_solvent, 
    3     double radius_equat_minor, double radius_equat_major, double radius_polar); 
    4 double Iqxy(double qx, double qy, double sld, double sld_solvent, 
    5     double radius_equat_minor, double radius_equat_major, double radius_polar, double theta, double phi, double psi); 
    6  
    71//#define INVALID(v) (v.radius_equat_minor > v.radius_equat_major || v.radius_equat_major > v.radius_polar) 
    82 
    9  
    10 double form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     3static double 
     4form_volume(double radius_equat_minor, double radius_equat_major, double radius_polar) 
    115{ 
    126    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
    137} 
    148 
    15 double Iq(double q, 
     9static double 
     10Iq(double q, 
    1611    double sld, 
    1712    double sld_solvent, 
    … …  
    4540    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    4641    const double fqsq = outer/4.0;  // = outer*um*zm*8.0/(4.0*M_PI); 
    47     const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
    48     return 1.0e-4 * s * s * fqsq; 
     42    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     43    const double drho = (sld - sld_solvent); 
     44    return 1.0e-4 * square(vol*drho) * fqsq; 
    4945} 
    5046 
    51 double Iqxy(double qx, double qy, 
     47static double 
     48Iqxy(double qa, double qb, double qc, 
    5249    double sld, 
    5350    double sld_solvent, 
    5451    double radius_equat_minor, 
    5552    double radius_equat_major, 
    56     double radius_polar, 
    57     double theta, 
    58     double phi, 
    59     double psi) 
     53    double radius_polar) 
    6054{ 
    61     double q, xhat, yhat, zhat; 
    62     ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 
     55    const double qr = sqrt(square(radius_equat_minor*qa) 
     56                           + square(radius_equat_major*qb) 
     57                           + square(radius_polar*qc)); 
     58    const double fq = sas_3j1x_x(qr); 
     59    const double vol = form_volume(radius_equat_minor, radius_equat_major, radius_polar); 
     60    const double drho = (sld - sld_solvent); 
    6361 
    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); 
    69  
    70     return 1.0e-4 * square(s * fq); 
     62    return 1.0e-4 * square(vol * drho * fq); 
    7163} 
    72  
Note: See TracChangeset for help on using the changeset viewer.