From fb4312c8c45224dac79e3e24b9b4081f5b64cb8d Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 24 Jul 2025 10:36:14 -0400 Subject: [PATCH 1/2] fixed parametrization of UniformInShell function --- core/PhysiCell_rules.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/core/PhysiCell_rules.cpp b/core/PhysiCell_rules.cpp index 3b73a38a3..470537e85 100644 --- a/core/PhysiCell_rules.cpp +++ b/core/PhysiCell_rules.cpp @@ -2331,13 +2331,13 @@ std::vector UniformInShell( double r1, double r2 ) sqrt_one_minus_T -= T; sqrt_one_minus_T = sqrt( sqrt_one_minus_T ); - double param1 = pow( UniformRandom() , 0.33333333333333333333333333333333333333 ); // xi^(1/3), - // param1 *= (r2-r1); - // param1 += r1; - double param2 = param1; // xi^(1/3) - param2 *= 2.0; // 2 * xi^(1/3) - param2 *= sqrt_T; // 2 * xi(1) * T^(1/2) - param2 *= sqrt_one_minus_T; // 2 * xi(1) * T^(1/2) * (1-T)^(1/2) + double r1_3 = r1*r1*r1; + double r2_3 = r2*r2*r2; + double param1 = pow( r1_3 + (r2_3-r1_3) * UniformRandom() , 0.33333333333333333333333333333333333333 ); // r scaled by xi^(1/3) + double param2 = param1; // r * xi^(1/3) + param2 *= 2.0; // 2 * r * xi^(1/3) + param2 *= sqrt_T; // 2 * r * xi(1) * T^(1/2) + param2 *= sqrt_one_minus_T; // 2 * r * xi(1) * T^(1/2) * (1-T)^(1/2) double theta = UniformRandom(); // U(0,1) theta *= two_pi; // U(0,2*pi) From 6fe2f7ba493baa73ece982b8dc8d4b56f6ce135f Mon Sep 17 00:00:00 2001 From: Heber Lima da Rocha Date: Thu, 24 Jul 2025 13:29:18 -0400 Subject: [PATCH 2/2] UniformInShell function simplified using UniformOnUnitSphere --- core/PhysiCell_rules.cpp | 29 ++++++++++------------------- 1 file changed, 10 insertions(+), 19 deletions(-) diff --git a/core/PhysiCell_rules.cpp b/core/PhysiCell_rules.cpp index 470537e85..a63792d73 100644 --- a/core/PhysiCell_rules.cpp +++ b/core/PhysiCell_rules.cpp @@ -2321,28 +2321,19 @@ std::vector UniformInAnnulus( double r1, double r2 ) return {x,y,0.0}; } -std::vector UniformInShell( double r1, double r2 ) +std::vector UniformInShell(double r1, double r2) { - static double two_pi = 6.283185307179586; + double r1_3 = r1 * r1 * r1; + double r2_3 = r2 * r2 * r2; + + double xi = UniformRandom(); + double r = pow(r1_3 + xi * (r2_3 - r1_3), 1.0 / 3.0); - double T = UniformRandom(); - double sqrt_T = sqrt(T); - double sqrt_one_minus_T = 1.0; - sqrt_one_minus_T -= T; - sqrt_one_minus_T = sqrt( sqrt_one_minus_T ); + std::vector dir = UniformOnUnitSphere(); + for (int i = 0; i < 3; ++i) + dir[i] *= r; - double r1_3 = r1*r1*r1; - double r2_3 = r2*r2*r2; - double param1 = pow( r1_3 + (r2_3-r1_3) * UniformRandom() , 0.33333333333333333333333333333333333333 ); // r scaled by xi^(1/3) - double param2 = param1; // r * xi^(1/3) - param2 *= 2.0; // 2 * r * xi^(1/3) - param2 *= sqrt_T; // 2 * r * xi(1) * T^(1/2) - param2 *= sqrt_one_minus_T; // 2 * r * xi(1) * T^(1/2) * (1-T)^(1/2) - - double theta = UniformRandom(); // U(0,1) - theta *= two_pi; // U(0,2*pi) - - return { param2*sin(theta) , param2*cos(theta), param1*(1-2*T) }; + return dir; } void setup_cell_rules( void )