Changeset 99658f6 in sasmodels for sasmodels/models/triaxial_ellipsoid.c


Ignore:
Timestamp:
Nov 6, 2018 4:10:43 PM (21 months ago)
Author:
grethevj
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
cf3d0ce
Parents:
5024a56
Message:

updated ER functions including cylinder excluded volume, to match 4.x

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/triaxial_ellipsoid.c

    rd42dd4a r99658f6  
    55{ 
    66    return M_4PI_3*radius_equat_minor*radius_equat_major*radius_polar; 
     7} 
     8 
     9static double 
     10radius_from_curvature(double radius_equat_minor, double radius_equat_major, double radius_polar) 
     11{ 
     12    // Trivial cases 
     13    if (radius_equat_minor == radius_equat_major == radius_polar) return radius_polar; 
     14    if (radius_equat_minor * radius_equat_major * radius_polar == 0.)  return 0.; 
     15 
     16 
     17    double r_equat_equiv, r_polar_equiv; 
     18    double radii[3] = {radius_equat_minor, radius_equat_major, radius_polar}; 
     19    double radmax = fmax(radii[0],fmax(radii[1],radii[2])); 
     20 
     21    double radius_1 = radmax; 
     22    double radius_2 = radmax; 
     23    double radius_3 = radmax; 
     24 
     25    for(int irad=0; irad<3; irad++) { 
     26        if (radii[irad] < radius_1) { 
     27            radius_3 = radius_2; 
     28            radius_2 = radius_1; 
     29            radius_1 = radii[irad]; 
     30            } else { 
     31                if (radii[irad] < radius_2) { 
     32                        radius_2 = radii[irad]; 
     33                } 
     34            } 
     35    } 
     36    if(radius_2-radius_1 > radius_3-radius_2) { 
     37        r_equat_equiv = sqrt(radius_2*radius_3); 
     38        r_polar_equiv = radius_1; 
     39    } else  { 
     40        r_equat_equiv = sqrt(radius_1*radius_2); 
     41        r_polar_equiv = radius_3; 
     42    } 
     43 
     44    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     45    const double ratio = (r_polar_equiv < r_equat_equiv 
     46                          ? r_polar_equiv / r_equat_equiv 
     47                          : r_equat_equiv / r_polar_equiv); 
     48    const double e1 = sqrt(1.0 - ratio*ratio); 
     49    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     50    const double bL = (1.0 + e1) / (1.0 - e1); 
     51    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     52    const double delta = 0.75 * b1 * b2; 
     53    const double ddd = 2.0 * (delta + 1.0) * r_polar_equiv * r_equat_equiv * r_equat_equiv; 
     54    return 0.5 * cbrt(ddd); 
    755} 
    856 
     
    3280    switch (mode) { 
    3381    default: 
    34     case 1: // equivalent sphere 
     82    case 1: // equivalent biaxial ellipsoid average curvature 
     83        return radius_from_curvature(radius_equat_minor,radius_equat_major, radius_polar); 
     84    case 2: // equivalent volume sphere 
    3585        return radius_from_volume(radius_equat_minor,radius_equat_major, radius_polar); 
    36     case 2: // min radius 
     86    case 3: // min radius 
    3787        return radius_from_min_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    38     case 3: // max radius 
     88    case 4: // max radius 
    3989        return radius_from_max_dimension(radius_equat_minor,radius_equat_major, radius_polar); 
    4090    } 
Note: See TracChangeset for help on using the changeset viewer.