source: sasmodels/sasmodels/models/core_shell_cylinder.c @ d42dd4a

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d42dd4a was d42dd4a, checked in by pkienzle, 5 years ago

fix compiler warnings for CUDA

  • Property mode set to 100644
File size: 3.5 KB
RevLine 
[9aa4881]1// vd = volume * delta_rho
[2a0b2b1]2// besarg = q * R * sin(theta)
3// siarg = q * L/2 * cos(theta)
4static double _cyl(double vd, double besarg, double siarg)
[5d4777d]5{
[592343f]6    return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg);
[5d4777d]7}
8
[2a0b2b1]9static double
10form_volume(double radius, double thickness, double length)
[5d4777d]11{
[2a0b2b1]12    return M_PI*square(radius+thickness)*(length+2.0*thickness);
[5d4777d]13}
14
[d277229]15static double
16radius_from_volume(double radius, double thickness, double length)
17{
18    const double volume_outer_cyl = form_volume(radius,thickness,length);
[6d5601c]19    return cbrt(volume_outer_cyl/M_4PI_3);
[d277229]20}
21
22static double
23radius_from_diagonal(double radius, double thickness, double length)
24{
25    const double radius_outer = radius + thickness;
[a94046f]26    const double length_outer = length + 2.0*thickness;
[d277229]27    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer);
28}
29
30static double
31effective_radius(int mode, double radius, double thickness, double length)
32{
[ee60aa7]33    switch (mode) {
[d42dd4a]34    default:
[ee60aa7]35    case 1: // equivalent sphere
[d277229]36        return radius_from_volume(radius, thickness, length);
[ee60aa7]37    case 2: // outer radius
[d277229]38        return radius + thickness;
[ee60aa7]39    case 3: // half outer length
[d277229]40        return 0.5*length + thickness;
[ee60aa7]41    case 4: // half min outer length
[d277229]42        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness);
[ee60aa7]43    case 5: // half max outer length
[d277229]44        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness);
[ee60aa7]45    case 6: // half outer diagonal
[d277229]46        return radius_from_diagonal(radius,thickness,length);
47    }
48}
49
[71b751d]50static void
51Fq(double q,
52    double *F1,
53    double *F2,
[994d77f]54    double core_sld,
55    double shell_sld,
56    double solvent_sld,
57    double radius,
58    double thickness,
59    double length)
[5d4777d]60{
61    // precalculate constants
[2a0b2b1]62    const double core_r = radius;
63    const double core_h = 0.5*length;
[9aa4881]64    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
[2a0b2b1]65    const double shell_r = (radius + thickness);
66    const double shell_h = (0.5*length + thickness);
[9aa4881]67    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
[71b751d]68    double total_F1 = 0.0;
69    double total_F2 = 0.0;
[74768cb]70    for (int i=0; i<GAUSS_N ;i++) {
[2a0b2b1]71        // translate a point in [-1,1] to a point in [0, pi/2]
[74768cb]72        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0;
[2a0b2b1]73        double sin_theta, cos_theta;
[74768cb]74        const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4;
[2a0b2b1]75        SINCOS(theta, sin_theta,  cos_theta);
76        const double qab = q*sin_theta;
77        const double qc = q*cos_theta;
78        const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
79            + _cyl(shell_vd, shell_r*qab, shell_h*qc);
[71b751d]80        total_F1 += GAUSS_W[i] * fq * sin_theta;
81        total_F2 += GAUSS_W[i] * fq * fq * sin_theta;
[5d4777d]82    }
83    // translate dx in [-1,1] to dx in [lower,upper]
[994d77f]84    //const double form = (upper-lower)/2.0*total;
[71b751d]85    *F1 = 1.0e-2 * total_F1 * M_PI_4;
86    *F2 = 1.0e-4 * total_F2 * M_PI_4;
[5d4777d]87}
88
[d86f0fc]89static double
90Iqac(double qab, double qc,
[994d77f]91    double core_sld,
92    double shell_sld,
93    double solvent_sld,
94    double radius,
95    double thickness,
[becded3]96    double length)
[5d4777d]97{
[2a0b2b1]98    const double core_r = radius;
99    const double core_h = 0.5*length;
[9aa4881]100    const double core_vd = form_volume(radius,0,length) * (core_sld-shell_sld);
[2a0b2b1]101    const double shell_r = (radius + thickness);
102    const double shell_h = (0.5*length + thickness);
[9aa4881]103    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld);
[5d4777d]104
[2a0b2b1]105    const double fq = _cyl(core_vd, core_r*qab, core_h*qc)
106        + _cyl(shell_vd, shell_r*qab, shell_h*qc);
[994d77f]107    return 1.0e-4 * fq * fq;
[5d4777d]108}
Note: See TracBrowser for help on using the repository browser.