Changes in sasmodels/models/core_shell_ellipsoid.c [71b751d:3c60146] in sasmodels
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_shell_ellipsoid.c
r71b751d r3c60146 36 36 double vol = M_4PI_3*equat_shell*equat_shell*polar_shell; 37 37 return vol; 38 } 39 40 static double 41 radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 42 { 43 const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 44 return cbrt(0.75*volume_ellipsoid/M_PI); 45 } 46 47 static double 48 radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 49 { 50 // Trivial cases 51 if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell; 52 if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.) return 0.; 53 54 // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 55 const double radius_equat_tot = radius_equat_core + thick_shell; 56 const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 57 const double ratio = (radius_polar_tot < radius_equat_tot 58 ? radius_polar_tot / radius_equat_tot 59 : radius_equat_tot / radius_polar_tot); 60 const double e1 = sqrt(1.0 - ratio*ratio); 61 const double b1 = 1.0 + asin(e1) / (e1 * ratio); 62 const double bL = (1.0 + e1) / (1.0 - e1); 63 const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 64 const double delta = 0.75 * b1 * b2; 65 const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot; 66 return 0.5 * cbrt(ddd); 67 } 68 69 static double 70 effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 71 { 72 const double radius_equat_tot = radius_equat_core + thick_shell; 73 const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 74 75 if (mode == 1) { 76 return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 77 } else if (mode == 2) { 78 return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 79 } else if (mode == 3) { 80 return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 81 } else { 82 return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 83 } 38 84 } 39 85
Note: See TracChangeset
for help on using the changeset viewer.