Plugins
Loading...
Searching...
No Matches
entire_incomplete_gamma_function.h
Go to the documentation of this file.
1/*
2 * SPDX-FileCopyrightText: 2025 Gesellschaft fuer Anlagen- und Reaktorsicherheit gGmbH
3 * SPDX-License-Identifier: EUPL-1.2
4 * SPDX-FileContributor: Sabine Stichel
5 * SPDX-FileContributor: Goethe Universität Frankfurt
6 * SPDX-FileType: SOURCE
7 *
8 * This file is part of d3f++.
9 * d3f++ is an extension for UG4. Licensing information and citation requirements of UG4 are provided in LICENSES/UG4-LGPL_2.1
10 */
11
12
13/*
14 * Original by J.Musuuza
15 */
16
17// extern headers
18#include <iostream>
19#include <sstream>
20#include <fstream>
21#include <string>
22#include <cmath>
23#include <vector>
24
25
26#define _USE_MATH_DEFINES
27#include <math.h>
28#include <float.h> //for DBL_MAX
29#include <iomanip>
30#include <cstdlib>
31#include <time.h>
32#include "common/common.h"
33
34#include "lib_disc/common/groups_util.h"
35#include "lib_disc/common/function_group.h"
36#include "lib_grid/tools/subset_group.h"
37#include "bridge/bridge.h"
38#include "lib_disc/domain.h"
39#include "lib_disc/domain_traits.h"
40#include "lib_algebra/cpu_algebra_types.h"
41#include "lib_disc/quadrature/quadrature.h"
42#include "lib_disc/local_finite_element/local_finite_element_provider.h"
43#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
44#include "lib_disc/spatial_disc/user_data/user_data.h"
45#include "common/util/provider.h"
46#include "lib_disc/spatial_disc/disc_util/fv1_geom.h"
47#include "lib_disc/spatial_disc/disc_util/hfv1_geom.h"
48#include "common/parser/rapidxml/rapidxml_print.hpp"
49#include "common/parser/rapidxml/rapidxml.hpp"
50#include "../plugins/experimental/biogas/save_load.h"
51
52
53//using namespace std;
54//using namespace rapidxml;
55#ifdef UG_FOR_LUA
56#include "bindings/lua/lua_user_data.h"
57#endif
58
59namespace ug{
60namespace d3f{
61
62
63/*functions copied from: entire_incomplete_gamma_function.c*/
64
65
66// double Entire_Incomplete_Gamma_Function(double x, double nu)
67//
68// Description:
69// The entire incomplete gamma function, also called the regularized
70// incomplete gamma function, is defined as the integral from 0 to x of
71// the integrand t^(nu-1) exp(-t) / gamma(nu) dt. The parameter nu is
72// sometimes referred to as the shape parameter.
73//
74// Arguments:
75// double x Upper limit of the integral with integrand given above.
76// double nu The shape parameter of the entire incomplete gamma function.
77//
78// Return Values:
79//
80// Example:
81// double x, g, nu;
82//
83// g = Entire_Incomplete_Gamma_Function( x, nu );
84template <typename TData, int dim, typename TRet>
85double FractalField<TData,dim,TRet>::Entire_Incomplete_Gamma_Function(double x, double nu)
86{
87 return (double) xEntire_Incomplete_Gamma_Function((long double)x,(long double)nu);
88}
89
90
92// long double xEntire_Incomplete_Gamma_Function(long double x, //
93// long double nu) //
94// //
95// Description: //
96// The entire incomplete gamma function, also called the regularized //
97// incomplete gamma function, is defined as the integral from 0 to x of //
98// the integrand t^(nu-1) exp(-t) / gamma(nu) dt. The parameter nu is //
99// sometimes referred to as the shape parameter. //
100// //
101// Arguments: //
102// long double x Upper limit of the integral with integrand given above.//
103// long double nu The shape parameter of the entire incomplete gamma //
104// function. //
105// //
106// Return Values: //
107// //
108// Example: //
109// long double x, g, nu; //
110// //
111// g = xEntire_Incomplete_Gamma_Function( x, nu ); //
113template <typename TData, int dim, typename TRet>
114long double FractalField<TData,dim,TRet>::xEntire_Incomplete_Gamma_Function(long double x, long double nu)
115{
116
117 if (x == 0.0L) return 0.0L;
118 if (fabsl(x) <= 1.0L) return xSmall_x(x, nu);
119 if (fabsl(x) < (nu + 1.0L) ) return xMedium_x(x, nu);
120 return xLarge_x(x, nu);
121}
122
123
125// static long double xSmall_x(long double x, long double nu) //
126// //
127// Description: //
128// This function approximates the entire incomplete gamma function for //
129// x, where -1 <= x <= 1. //
130// //
131// Arguments: //
132// long double x Upper limit of the integral with integrand described //
133// in the section under Entire_Incomplete_Gamma_Function. //
134// long double nu The shape parameter of the entire incomplete gamma //
135// function. //
136// //
137// Return Values: //
138// The entire incomplete gamma function: //
139// I(0,x) t^(nu-1) Exp(-t) dt / Gamma(nu). //
140// //
141// Example: //
142// long double x, g, nu; //
143// //
144// g = xSmall_x( x, nu); //
146#define Nterms 20
147template <typename TData, int dim, typename TRet>
148long double FractalField<TData,dim,TRet>::xSmall_x(long double x, long double nu)
149{
150 long double terms[Nterms];
151 long double x_term = -x;
152 long double x_power = 1.0L;
153 long double sum;
154 int i;
155 long double buffer;
156
157 for (i = 0; i < Nterms; i++)
158 {
159 buffer= xFactorial(i);
160 terms[i] = (x_power / xFactorial(i)) / (i + nu);
161 x_power *= x_term;
162 }
163 sum = terms[Nterms-1];
164 for (i = Nterms-2; i >= 0; i--)
165 sum += terms[i];
166 if ( nu <= Gamma_Function_Max_Arg() )
167 return powl(x,nu) * sum / xGamma_Function(nu);
168 else
169 return expl(nu * logl(x) + logl(sum) - xLn_Gamma_Function(nu));
170}
171
172
174// static long double xMedium_x(long double x, long double nu) //
175// //
176// Description: //
177// This function approximates the entire incomplete gamma function for //
178// x, where 1 < x < nu + 1. //
179// //
180// If nu + 1 < x, then one should use xLarge_x(x,nu). //
181// //
182// Arguments: //
183// long double x Upper limit of the integral with integrand described //
184// in the section under Entire_Incomplete_Gamma_Function. //
185// long double nu The shape parameter of the entire incomplete gamma //
186// function. //
187// //
188// Return Values: //
189// The entire incomplete gamma function: //
190// I(0,x) t^(nu-1) exp(-t) dt / gamma(nu). //
191// //
192// Example: //
193// long double x, g, nu; //
194// //
195// g = xMedium_x( x, nu); //
197template <typename TData, int dim, typename TRet>
198long double FractalField<TData,dim,TRet>::xMedium_x(long double x, long double nu)
199{
200 long double coef;
201 long double term = 1.0L / nu;
202 long double corrected_term = term;
203 long double temp_sum = term;
204 long double correction = -temp_sum + corrected_term;
205 long double sum1 = temp_sum;
206 long double sum2;
207 long double epsilon = 0.0L;
208 int i;
209
210 if (nu > Gamma_Function_Max_Arg())
211 {
212 coef = expl( nu * logl(x) - x - xLn_Gamma_Function(nu) );
213 if (coef > 0.0L) epsilon = DBL_EPSILON/coef;
214 }
215 else
216 {
217 coef = powl(x, nu) * expl(-x) / xGamma_Function(nu);
218 epsilon = DBL_EPSILON/coef;
219 }
220 if (epsilon <= 0.0L) epsilon = (long double) DBL_EPSILON;
221
222 for (i = 1; term > epsilon * sum1; i++)
223 {
224 term *= x / (nu + i);
225 corrected_term = term + correction;
226 temp_sum = sum1 + corrected_term;
227 correction = (sum1 - temp_sum) + corrected_term;
228 sum1 = temp_sum;
229 }
230 sum2 = sum1;
231 sum1 *= coef;
232 correction += sum2 - sum1 / coef;
233 term *= x / (nu + i);
234 sum2 = term + correction;
235 for (i++; (term + correction) > epsilon * sum2; i++)
236 {
237 term *= x / (nu + i);
238 corrected_term = term + correction;
239 temp_sum = sum2 + corrected_term;
240 correction = (sum2 - temp_sum) + corrected_term;
241 sum2 = temp_sum;
242 }
243
244 sum2 += correction;
245 sum2 *= coef;
246 return sum1 + sum2;
247}
248
249
251// static long double xLarge_x(long double x, long double nu) //
252// //
253// Description: //
254// This function approximates the entire incomplete gamma function for //
255// x, where nu + 1 <= x. //
256// //
257// If 0 <= x < nu + 1, then one should use xSmall_x(x,nu). //
258// //
259// Arguments: //
260// long double x Upper limit of the integral with integrand described //
261// in the section under Entire_Incomplete_Gamma_Function. //
262// long double nu The shape parameter of the entire incomplete gamma //
263// function. //
264// //
265// Return Values: //
266// If x is positive and is less than 171 then Gamma(x) is returned and //
267// if x > 171 then DBL_MAX is returned. //
268// //
269// Example: //
270// long double x, g, nu; //
271// //
272// g = xLarge_x( x, nu); //
274template <typename TData, int dim, typename TRet>
275long double FractalField<TData,dim,TRet>::xLarge_x(long double x, long double nu)
276{
277 long double temp = 1.0L / nu;
278 long double sum = temp;
279 long double coef;
280 int i = 0;
281 int n;
282
283 n = (int)(x - nu - 1.0L) + 1;
284 for (i = 1; i < n; i++)
285 {
286 temp *= x / (nu + i);
287 sum += temp;
288 }
289 if ( nu <= Gamma_Function_Max_Arg() )
290 {
291 coef = powl(x, nu) * expl(-x) / xGamma_Function(nu);
292 return xMedium_x(x, nu + n) + coef * sum;
293 }
294 else
295 {
296 return expl(logl(sum) + nu * logl(x) - x - xLn_Gamma_Function(nu)) + xMedium_x(x, nu + n);
297 }
298}
299
300
301
302}
303}
ug::d3f::FractalField::xLarge_x
long double xLarge_x(long double x, long double nu)
Definition entire_incomplete_gamma_function.h:275
ug::d3f::FractalField::Entire_Incomplete_Gamma_Function
double Entire_Incomplete_Gamma_Function(double x, double nu)
Definition entire_incomplete_gamma_function.h:85
ug::d3f::FractalField::xEntire_Incomplete_Gamma_Function
long double xEntire_Incomplete_Gamma_Function(long double x, long double nu)
Definition entire_incomplete_gamma_function.h:114
ug::d3f::FractalField::xMedium_x
long double xMedium_x(long double x, long double nu)
Definition entire_incomplete_gamma_function.h:198
ug::d3f::FractalField::xSmall_x
long double xSmall_x(long double x, long double nu)
Definition entire_incomplete_gamma_function.h:148
cpu_algebra_types.h
domain_traits.h
Nterms
#define Nterms
Definition entire_incomplete_gamma_function.h:146
function_group.h
fv1_geom.h
groups_util.h
hfv1_geom.h
user_data.h
local_finite_element_provider.h
lua_user_data.h
subset_group.h