Changeset 6d96b66 in sasmodels


Ignore:
Timestamp:
Oct 18, 2016 8:32:18 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
f0afad2
Parents:
3e8ea5d
Message:

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

Location:
sasmodels/models
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/fractal.c

    r2c74c11 r6d96b66  
    1 double Iq(double q, 
    2           double volfraction, 
    3           double radius, 
    4           double fractal_dim, 
    5           double cor_length, 
    6           double sld_block, 
    7           double sld_solvent); 
     1#define INVALID(p) (p.fractal_dim <= 0.0) 
    82 
    9 double Iq(double q, 
    10           double volfraction, 
    11           double radius, 
    12           double fractal_dim, 
    13           double cor_length, 
    14           double sld_block, 
    15           double sld_solvent) 
     3static double 
     4Iq(double q, 
     5   double volfraction, 
     6   double radius, 
     7   double fractal_dim, 
     8   double cor_length, 
     9   double sld_block, 
     10   double sld_solvent) 
    1611{ 
    17     double qr,r0,Df,corr,phi,sldp,sldm; 
    18     double pq,sq,inten; 
     12    //calculate P(q) for the spherical subunits 
     13    const double pq = M_4PI_3*cube(radius) * square((sld_block-sld_solvent)*sph_j1c(q*radius)); 
    1914     
    20      // Actively check the argument - needed for mass fractal - is it needie 
    21      //here? 
    22     if (fractal_dim <= 0.0){ 
    23        return 0.0; 
     15    //calculate S(q),  using Teixeira, Eq(15) 
     16    double sq; 
     17    if (q > 0. && fractal_dim > 1.) { 
     18        // q>0, D>0 
     19        const double D = fractal_dim; 
     20        const double Dm1 = fractal_dim - 1.0; 
     21        // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative 
     22        const double t1 = D*sas_gamma(Dm1)*sin(Dm1*atan(q*cor_length)); 
     23        const double t2 = pow(q*radius, -D); 
     24        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1); 
     25        sq = 1.0 + t1 * t2 * t3; 
     26    } else if (q > 0.) { 
     27        // q>0, D=1 
     28        sq = 1.0 + atan(q*cor_length) / (q*radius); 
     29    } else if (fractal_dim > 1.) { 
     30        // q=0, D>1 
     31        const double D = fractal_dim; 
     32        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0); 
     33    } else { 
     34        // q=0, D=1 
     35        sq = 1.0 + cor_length/radius; 
    2436    } 
    25     
    26     phi = volfraction;        // volume fraction of building block spheres... 
    27     r0 = radius;     //  radius of building block 
    28     Df = fractal_dim;     //  fractal dimension 
    29     corr = cor_length;       //  correlation length of fractal-like aggregates 
    30     sldp = sld_block;       // SLD of building block 
    31     sldm = sld_solvent;       // SLD of matrix or solution 
    32   
    33      qr=q*r0; 
    34      
    35     //calculate P(q) for the spherical subunits 
    36     pq = phi*M_4PI_3*r0*r0*r0*(sldp-sldm)*(sldp-sldm)*sph_j1c(qr)*sph_j1c(qr); 
    37      
    38     //calculate S(q) 
    39     sq = Df*sas_gamma(Df-1.0)*sin((Df-1.0)*atan(q*corr)); 
    40     sq /= pow(qr,Df) * pow((1.0 + 1.0/(q*corr)/(q*corr)),((Df-1.0)/2.0)); 
    41     sq += 1.0; 
    42      
    43     //combine, scale to units cm-1 sr-1 (assuming data on absolute scale) 
    44     //and return 
    45     inten = pq*sq; 
    46     // convert I(1/A) to (1/cm) 
    47     inten *= 1.0e8; 
    48     //convert rho^2 in 10^-6 1/A to 1/A 
    49     inten *= 1.0e-12;     
    50      
    51      
    52     return(inten); 
     37 
     38    // scale to units cm-1 sr-1 (assuming data on absolute scale) 
     39    //    convert I(1/A) to (1/cm)  => 1e8 * I(q) 
     40    //    convert rho^2 in 10^-6 1/A to 1/A  => 1e-12 * I(q) 
     41    //    combined: 1e-4 * I(q) 
     42 
     43    return 1.e-4 * volfraction * pq * sq; 
    5344} 
    5445 

  • sasmodels/models/fractal.py

    r785cbec r6d96b66  
    8484              ["radius",    "Ang",  5.0, [0.0, inf], "", 
    8585               "radius of particles"], 
    86               ["fractal_dim",      "",  2.0, [0.0, 6.0], "", 
     86              ["fractal_dim",      "",  2.0, [1.0, 6.0], "", 
    8787               "fractal dimension"], 
    8888              ["cor_length", "Ang", 100.0, [0.0, inf], "", 

  • sasmodels/models/fractal_core_shell.c

    r3a48772 r6d96b66  
    2727 
    2828 
    29    double intensity = core_shell_kernel(q, 
    30                               radius, 
    31                               thickness, 
    32                               core_sld, 
    33                               shell_sld, 
    34                               solvent_sld); 
     29    const double pq = core_shell_kernel(q, radius, thickness, 
     30                                        core_sld, shell_sld, solvent_sld); 
     31 
     32 
    3533    //calculate S(q) 
    36     double frac_1 = fractal_dim-1.0; 
    37     double qr = q*radius; 
     34    double sq; 
     35    if (q > 0. && fractal_dim > 1.) { 
     36        // q>0, D>0 
     37        const double D = fractal_dim; 
     38        const double Dm1 = fractal_dim - 1.0; 
     39        const double t1 = D*sas_gamma(Dm1)*sin((Dm1)*atan(q*cor_length)); 
     40        const double t2 = pow(q*radius, -D); 
     41        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1); 
     42        sq = 1.0 + t1 * t2 * t3; 
     43    } else if (q > 0.) { 
     44        // q>0, D=1 
     45        sq = 1.0 + atan(q*cor_length) / (q*radius); 
     46    } else if (fractal_dim > 1.) { 
     47        // q=0, D>1 
     48        const double D = fractal_dim; 
     49        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0); 
     50    } else { 
     51        // q=0, D=1 
     52        sq = 1.0 + cor_length/radius; 
     53    } 
    3854 
    39     double t1 = fractal_dim*sas_gamma(frac_1)*sin(frac_1*atan(q*cor_length)); 
    40     double t2 = (1.0 + 1.0/(q*cor_length)/(q*cor_length)); 
    41     double t3 = pow(qr, fractal_dim) * pow(t2, (frac_1/2.0)); 
    42     double sq = t1/t3; 
    43     sq += 1.0; 
    44  
    45     return sq*intensity*volfraction; 
     55    // Note: core_shell_kernel already performs the 1e-4 unit conversion 
     56    return volfraction * sq * pq; 
    4657} 
    4758 

  • sasmodels/models/fractal_core_shell.py

    r4962519 r6d96b66  
    6363#   ["name", "units", default, [lower, upper], "type","description"], 
    6464parameters = [ 
    65     ["radius",      "Ang",        60.0, [0, inf],    "volume", "Sphere core radius"], 
    66     ["thickness",   "Ang",        10.0, [0, inf],    "volume", "Sphere shell thickness"], 
     65    ["radius",      "Ang",        60.0, [0.0, inf],  "volume", "Sphere core radius"], 
     66    ["thickness",   "Ang",        10.0, [0.0, inf],  "volume", "Sphere shell thickness"], 
    6767    ["sld_core",    "1e-6/Ang^2", 1.0,  [-inf, inf], "sld",    "Sphere core scattering length density"], 
    6868    ["sld_shell",   "1e-6/Ang^2", 2.0,  [-inf, inf], "sld",    "Sphere shell scattering length density"], 
    6969    ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"], 
    70     ["volfraction", "",           1.0,  [0, inf],    "",       "Volume fraction of building block spheres"], 
    71     ["fractal_dim",    "",           2.0,  [-inf, inf], "",       "Fractal dimension"], 
    72     ["cor_length",  "Ang",      100.0,  [0, inf],    "",       "Correlation length of fractal-like aggregates"]] 
     70    ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"], 
     71    ["fractal_dim",    "",        2.0,  [1.0, 6.0],  "",       "Fractal dimension"], 
     72    ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"], 
     73] 
    7374# pylint: enable=bad-whitespace, line-too-long 
    7475 

  • sasmodels/models/mass_fractal.c

    r3a48772 r6d96b66  
    1 double form_volume(double radius); 
     1#define INVALID(p) (p.fractal_dim_mass < 1.0) 
    22 
    3 double Iq(double q, 
    4           double radius, 
    5           double fractal_dim_mass, 
    6           double cutoff_length); 
     3static double 
     4Iq(double q, double radius, double fractal_dim_mass, double cutoff_length) 
     5{ 
     6    //calculate P(q) 
     7    const double pq = square(sph_j1c(q*radius)); 
    78 
    8 static double _mass_fractal_kernel(double q, 
    9           double radius, 
    10           double fractal_dim_mass, 
    11           double cutoff_length) 
    12 { 
    13     // Actively check the argument. 
    14     if (fractal_dim_mass <= 1.0){ 
    15        return 0.0; 
     9    //calculate S(q) 
     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; 
    1628    } 
    1729 
    18     //calculate P(q) 
    19     double pq = sph_j1c(q*radius); 
    20     pq = pq*pq; 
    21  
    22     //calculate S(q) 
    23     double mmo = fractal_dim_mass-1.0; 
    24     double sq = sas_gamma(mmo)*sin((mmo)*atan(q*cutoff_length)); 
    25     sq *= pow(cutoff_length, mmo); 
    26     sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0)); 
    27     sq /= q; 
    28  
    29     //combine and return 
    30     double result = pq * sq; 
    31  
    32     return result; 
     30    return pq * sq; 
    3331} 
    34 double form_volume(double radius){ 
    35  
    36     return M_4PI_3*radius*radius*radius; 
    37 } 
    38  
    39 double Iq(double q, 
    40           double radius, 
    41           double fractal_dim_mass, 
    42           double cutoff_length) 
    43 { 
    44     return _mass_fractal_kernel(q, 
    45            radius, 
    46            fractal_dim_mass, 
    47            cutoff_length); 
    48 } 

  • sasmodels/models/mass_fractal.py

    ra807206 r6d96b66  
    7878 
    7979# pylint: disable=bad-whitespace, line-too-long 
    80 #             ["name", "units", default, [lower, upper], "type","description"], 
    81 parameters = [["radius",        "Ang",  10.0, [0.0, inf], "", "Particle radius"], 
    82               ["fractal_dim_mass",      "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"], 
    83               ["cutoff_length", "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 
    84              ] 
     80#   ["name", "units", default, [lower, upper], "type","description"], 
     81parameters = [ 
     82    ["radius",           "Ang",  10.0, [0.0, inf], "", "Particle radius"], 
     83    ["fractal_dim_mass", "",      1.9, [1.0, 6.0], "", "Mass fractal dimension"], 
     84    ["cutoff_length",    "Ang", 100.0, [0.0, inf], "", "Cut-off length"], 
     85] 
    8586# pylint: enable=bad-whitespace, line-too-long 
    8687 

  • sasmodels/models/mass_surface_fractal.c

    r3a48772 r6d96b66  
    1 double form_volume(double radius); 
    2  
    3 double Iq(double q, 
    4           double fractal_dim_mass, 
    5           double fractal_dim_surf, 
    6           double rg_cluster, 
    7           double rg_primary); 
    8  
    9 static double _mass_surface_fractal_kernel(double q, 
     1#define INVALID(p) (p.fractal_dim_mass + p.fractal_dim_surf > 6.) 
     2static double 
     3Iq(double q, 
    104          double fractal_dim_mass, 
    115          double fractal_dim_surf, 
    … …  
    148{ 
    159     //computation 
    16     double tot_dim = 6.0 - fractal_dim_surf - fractal_dim_mass; 
    17     fractal_dim_mass /= 2.0; 
    18     tot_dim /= 2.0; 
     10    const double Dm = 0.5*fractal_dim_mass; 
     11    const double Dt = 0.5*(6.0 - (fractal_dim_mass + fractal_dim_surf)); 
    1912 
    20     double rc_norm = rg_cluster * rg_cluster / (3.0 * fractal_dim_mass); 
    21     double rp_norm = rg_primary * rg_primary / (3.0 * tot_dim); 
     13    const double t1 = Dm==0. ? 1.0 : pow(1.0 + square(q*rg_cluster)/(3.0*Dm), -Dm); 
     14    const double t2 = Dt==0. ? 1.0 : pow(1.0 + square(q*rg_primary)/(3.0*Dt), -Dt); 
     15    const double form = t1*t2; 
    2216 
    23     //x for P 
    24     double x_val1 = 1.0 +  q * q * rc_norm; 
    25     double x_val2 = 1.0 +  q * q * rp_norm; 
    26  
    27     double inv_form = pow(x_val1, fractal_dim_mass) * pow(x_val2, tot_dim); 
    28  
    29     //another singular 
    30     if (inv_form == 0.0) return 0.0; 
    31  
    32     double form_factor = 1.0; 
    33     form_factor /= inv_form; 
    34  
    35     return (form_factor); 
     17    return form; 
    3618} 
    37 double form_volume(double radius){ 
    38  
    39     return M_4PI_3*radius*radius*radius; 
    40 } 
    41  
    42 double Iq(double q, 
    43           double fractal_dim_mass, 
    44           double fractal_dim_surf, 
    45           double rg_cluster, 
    46           double rg_primary) 
    47 { 
    48     return _mass_surface_fractal_kernel(q, 
    49             fractal_dim_mass, 
    50             fractal_dim_surf, 
    51             rg_cluster, 
    52             rg_primary); 
    53 } 

  • sasmodels/models/mass_surface_fractal.py

    r30fbe2e r6d96b66  
    7878 
    7979# pylint: disable=bad-whitespace, line-too-long 
    80 #             ["name", "units", default, [lower, upper], "type","description"], 
    81 parameters = [["fractal_dim_mass",      "",    1.8, [1e-16, 6.0], "", 
    82                "Mass fractal dimension"], 
    83               ["fractal_dim_surf",   "",    2.3, [1e-16, 6.0], "", 
    84                "Surface fractal dimension"], 
    85               ["rg_cluster", "Ang",   86.7, [0.0, inf], "", 
    86                "Cluster radius of gyration"], 
    87               ["rg_primary", "Ang", 4000.,  [0.0, inf], "", 
    88                "Primary particle radius of gyration"], 
    89              ] 
     80#   ["name", "units", default, [lower, upper], "type","description"], 
     81parameters = [ 
     82    ["fractal_dim_mass", "",      1.8, [0.0, 6.0], "", "Mass fractal dimension"], 
     83    ["fractal_dim_surf", "",      2.3, [0.0, 6.0], "", "Surface fractal dimension"], 
     84    ["rg_cluster",       "Ang",  86.7, [0.0, inf], "", "Cluster radius of gyration"], 
     85    ["rg_primary",       "Ang", 4000., [0.0, inf], "", "Primary particle radius of gyration"], 
     86] 
    9087# pylint: enable=bad-whitespace, line-too-long 
    9188 
Note: See TracChangeset for help on using the changeset viewer.