source: sasmodels/sasmodels/models/mass_fractal.c @ 6d96b66

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 6d96b66 was 6d96b66, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

fractal models: handle limits of q=0 and dim=1

  • Property mode set to 100644
File size: 1.0 KB
RevLine 
[6d96b66]1#define INVALID(p) (p.fractal_dim_mass < 1.0)
[5c1d341]2
[6d96b66]3static double
4Iq(double q, double radius, double fractal_dim_mass, double cutoff_length)
[5c1d341]5{
[9c461c7]6    //calculate P(q)
[6d96b66]7    const double pq = square(sph_j1c(q*radius));
[5c1d341]8
9    //calculate S(q)
[6d96b66]10    // S(q) = gamma(D-1) sin((D-1)atan(q c))/q c^(D-1) (1+(q c)^2)^(-(D-1)/2)
11    // lim D->1 [gamma(D-1) sin((D-1) a)] = a
12    // lim q->0 [sin((D-1) atan(q c))/q] = (D-1) c
13    // lim q,D->0,1 [gamma(D-1) sin((D-1)atan(q c))/q] = c
14    double sq;
15    if (q > 0. && fractal_dim_mass > 1.) {
16        const double Dm1 = fractal_dim_mass - 1.0;
17        const double t1 = sas_gamma(Dm1)*sin(Dm1*atan(q*cutoff_length));
18        const double t2 = pow(cutoff_length, Dm1);
19        const double t3 = pow(1.0 + square(q*cutoff_length), -0.5*Dm1);
20        sq = t1 * t2 * t3 / q;
21    } else if (q > 0.) {
22        sq = atan(q*cutoff_length)/q;
23    } else if (fractal_dim_mass > 1.) {
24        const double D = fractal_dim_mass;
25        sq = pow(cutoff_length, D) * sas_gamma(D);
26    } else {
27        sq = cutoff_length;
28    }
[5c1d341]29
[6d96b66]30    return pq * sq;
[5c1d341]31}
Note: See TracBrowser for help on using the repository browser.