swf.h (4477B)
1 /* Copyright (C) 2024 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #ifndef SWF_H 17 #define SWF_H 18 19 #include <rsys/rsys.h> 20 21 /* Library symbol management */ 22 #if defined(SWF_SHARED_BUILD) 23 #define SWF_API extern EXPORT_SYM 24 #else 25 #define SWF_API extern 26 #endif 27 28 #ifndef NDEBUG 29 #define SWF(Func) ASSERT(swf_ ## Func == RES_OK) 30 #else 31 #define SWF(Func) swf_ ## Func 32 #endif 33 34 enum swf_prediction { 35 SWF_LINEAR, 36 SWF_QUADRATIC 37 }; 38 39 /* Forward declarations of external data types */ 40 struct mem_allocator; 41 42 struct swf_H_tabulate_args { 43 double x_min; 44 double x_max; 45 46 /* The step member variable is used to calculate the x arguments to be 47 * tabulated. Each x to be tabulated is calculated by adding its product with 48 * the step value to the previous one. The delta x is therefore relative to 49 * each x value: the smaller the x, the smaller the tabulation step. And this 50 * is what we're looking for, since the function is more difficult to estimate 51 * near 0 */ 52 double step; /* x_next = x + x*step: */ 53 54 /* Force renormalization of the function. Note that H is already normalized to 55 * the correct interval, so normalization may be unnecessary. Worse still! It 56 * could introduce numerical inaccuracy. So don't normalize if you don't 57 * have to. */ 58 int normalize; 59 60 struct mem_allocator* allocator; /* NULL <=> use default allocator */ 61 }; 62 63 /* The H function's default tab values guarantee both numerical accuracy and a 64 * small memory footprint */ 65 #define SWF_H2D_TABULATE_ARGS_DEFAULT__ {8.0e-3, 6.05, 1e-3, 0, NULL} 66 static const struct swf_H_tabulate_args SWF_H2D_TABULATE_ARGS_DEFAULT = 67 SWF_H2D_TABULATE_ARGS_DEFAULT__; 68 #define SWF_H3D_TABULATE_ARGS_DEFAULT__ {7.5e-3, 3.337, 1e-3, 0, NULL} 69 static const struct swf_H_tabulate_args SWF_H3D_TABULATE_ARGS_DEFAULT = 70 SWF_H3D_TABULATE_ARGS_DEFAULT__; 71 72 /* Forward declarations of opaque data types */ 73 struct swf_tabulation; 74 75 BEGIN_DECLS 76 77 /* These functions evaluate the H function, as specified in: 78 * 79 * "The floating random walk and its application to Monte Carlo solutions of 80 * heat equations" - Haji-Sheikh A. and Sparrow E.M., SIAM Journal on Applied 81 * Mathematics, 14(2): 370-389, 1966 82 * 83 * These function are expressed as an infinite series, which is numerically 84 * evaluated by summing all terms, until convergence is reached. 85 * Two functions are defined here, H2d and H3d, which are used in 2D and 3D 86 * respectively: 87 * +INF 88 * __ 89 * H2d(x) = 1 - 2 >_ exp(-x*ak^2) / (ak*J1(ak)); ak is the k^th root of J0 90 * k=1 91 * 92 * +INF 93 * __ 94 * H3d(x) = 1 + 2 >_ (-1)^k * exp(-(PI*k)^2*x) 95 * k=1 96 * 97 * The H function is the cumulated function of the probability density used to 98 * sample passage times for a given sphere radius. Its parameter x is positive 99 * and is equal to alpha*tau/r^2, with: 100 * - alpha the diffusivity of the solid material: lambda/(rho*C) [m^2/s] 101 * - tau > 0 the time interval [s] 102 * - r > 0 is the radius of the sphere [m] */ 103 SWF_API double 104 swf_H2d_eval 105 (const double x); 106 107 SWF_API double 108 swf_H3d_eval 109 (const double x); 110 111 SWF_API double 112 swf_H3d_inverse 113 (const double y); /* in [0, 1[ */ 114 115 SWF_API double 116 swf_H2d_inverse 117 (const double y); /* in [0, 1[ */ 118 119 SWF_API res_T 120 swf_H2d_tabulate 121 (const struct swf_H_tabulate_args* args, 122 struct swf_tabulation** tab); 123 124 SWF_API res_T 125 swf_H3d_tabulate 126 (const struct swf_H_tabulate_args* args, 127 struct swf_tabulation** tab); 128 129 SWF_API res_T 130 swf_tabulation_ref_get 131 (struct swf_tabulation* tab); 132 133 SWF_API res_T 134 swf_tabulation_ref_put 135 (struct swf_tabulation* tab); 136 137 SWF_API double /* x */ 138 swf_tabulation_inverse 139 (const struct swf_tabulation* tab, 140 const enum swf_prediction prediction, 141 const double y); /* in [0, 1[ */ 142 143 END_DECLS 144 145 #endif /* SWF_H */