source: sasmodels/sasmodels/models/triaxial_ellipsoid.c @ d507c3a

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d507c3a was 9c461c7, checked in by piotr, 8 years ago
  1. Code review from PK: renamed J1c to sph_j1c
  2. Fixed mass_fractal
  • Property mode set to 100644
File size: 2.9 KB
Line 
1double form_volume(double req_minor, double req_major, double rpolar);
2double Iq(double q, double sld, double solvent_sld,
3    double req_minor, double req_major, double rpolar);
4double Iqxy(double qx, double qy, double sld, double solvent_sld,
5    double req_minor, double req_major, double rpolar, double theta, double phi, double psi);
6
7double form_volume(double req_minor, double req_major, double rpolar)
8{
9    return 1.333333333333333*M_PI*req_minor*req_major*rpolar;
10}
11
12double Iq(double q,
13    double sld,
14    double solvent_sld,
15    double req_minor,
16    double req_major,
17    double rpolar)
18{
19    // if (req_minor > req_major || req_major > rpolar) return NAN;  // Exclude invalid region
20
21    double sn, cn;
22    //double st, ct;
23    //const double lower = 0.0;
24    //const double upper = 1.0;
25    double outer = 0.0;
26    for (int i=0;i<76;i++) {
27        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;
28        const double x = 0.5*(Gauss76Z[i] + 1.0);
29        SINCOS(M_PI_2*x, sn, cn);
30        const double acosx2 = req_minor*req_minor*cn*cn;
31        const double bsinx2 = req_major*req_major*sn*sn;
32        const double c2 = rpolar*rpolar;
33
34        double inner = 0.0;
35        for (int j=0;j<76;j++) {
36            const double y = 0.5*(Gauss76Z[j] + 1.0);
37            const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y);
38            const double fq = sph_j1c(t);
39            inner += Gauss76Wt[j] * fq * fq ;
40        }
41        outer += Gauss76Wt[i] * 0.5 * inner;
42    }
43    //const double fq2 = (upper-lower)/2*outer;
44    const double fq2 = 0.5*outer;
45    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
46    return 1.0e-4 * fq2 * s * s;
47}
48
49double Iqxy(double qx, double qy,
50    double sld,
51    double solvent_sld,
52    double req_minor,
53    double req_major,
54    double rpolar,
55    double theta,
56    double phi,
57    double psi)
58{
59    // if (req_minor > req_major || req_major > rpolar) return NAN;  // Exclude invalid region
60
61    double stheta, ctheta;
62    double sphi, cphi;
63    double spsi, cpsi;
64    double st, ct;
65
66    const double q = sqrt(qx*qx + qy*qy);
67    const double qxhat = qx/q;
68    const double qyhat = qy/q;
69    SINCOS(theta*M_PI_180, stheta, ctheta);
70    SINCOS(phi*M_PI_180, sphi, cphi);
71    SINCOS(psi*M_PI_180, spsi, cpsi);
72    const double calpha = ctheta*cphi*qxhat + stheta*qyhat;
73    const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;
74    const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;
75    const double t = q*sqrt(req_minor*req_minor*cnu*cnu
76                          + req_major*req_major*cmu*cmu
77                          + rpolar*rpolar*calpha*calpha);
78    SINCOS(t, st, ct);
79    const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) );
80    const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);
81
82    return 1.0e-4 * fq * fq * s * s;
83}
84
Note: See TracBrowser for help on using the repository browser.