"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/F_Math.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

[ To the main getdp source changes report ]

F_Math.cpp  (getdp-3.4.0-source.tgz):F_Math.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
// g++ -std=c++11 on mingw does not define bessel functions // g++ -std=c++11 on mingw does not define bessel functions
#if defined (WIN32) && !defined(__CYGWIN__) #if defined(WIN32) && !defined(__CYGWIN__)
#undef __STRICT_ANSI__ #undef __STRICT_ANSI__
#endif #endif
#include <math.h> #include <math.h>
#include <complex> #include <complex>
#include "ProData.h" #include "ProData.h"
#include "F.h" #include "F.h"
#include "MallocUtils.h" #include "MallocUtils.h"
#include "Message.h" #include "Message.h"
#include "Bessel.h" #include "Bessel.h"
extern struct CurrentData Current ; extern struct CurrentData Current;
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C math functions (scalar, 1 argument, imaginary part set to zero) */ /* C math functions (scalar, 1 argument, imaginary part set to zero) */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
#define scalar_real_1_arg(func, string) \ #define scalar_real_1_arg(func, string) \
int k; \ int k; \
\ \
if(A->Type != SCALAR) \ if(A->Type != SCALAR) \
Message::Error("Non scalar argument for function '" string "'"); \ Message::Error("Non scalar argument for function '" string "'"); \
if(Current.NbrHar > 1 && A->Val[MAX_DIM] != 0.) \ if(Current.NbrHar > 1 && A->Val[MAX_DIM] != 0.) \
Message::Error("Function '" string "' only accepts real arguments");\ Message::Error("Function '" string "' only accepts real arguments"); \
V->Val[0] = func(A->Val[0]) ; \ V->Val[0] = func(A->Val[0]); \
if (Current.NbrHar > 1){ \ if(Current.NbrHar > 1) { \
V->Val[MAX_DIM] = 0. ; \ V->Val[MAX_DIM] = 0.; \
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) \ for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) \
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; \ V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.; \
} \ } \
V->Type = SCALAR; V->Type = SCALAR;
void F_Floor (F_ARG) { scalar_real_1_arg (floor,"Floor") } void F_Floor(F_ARG) { scalar_real_1_arg(floor, "Floor") }
void F_Ceil (F_ARG) { scalar_real_1_arg (ceil, "Ceil") } void F_Ceil(F_ARG) { scalar_real_1_arg(ceil, "Ceil") }
void F_Fabs (F_ARG) { scalar_real_1_arg (fabs, "Fabs") } void F_Fabs(F_ARG) { scalar_real_1_arg(fabs, "Fabs") }
void F_Asin (F_ARG) { scalar_real_1_arg (asin, "Asin") } void F_Asin(F_ARG) { scalar_real_1_arg(asin, "Asin") }
void F_Acos (F_ARG) { scalar_real_1_arg (acos, "Acos") } void F_Acos(F_ARG) { scalar_real_1_arg(acos, "Acos") }
void F_Atan (F_ARG) { scalar_real_1_arg (atan, "Atan") } void F_Atan(F_ARG) { scalar_real_1_arg(atan, "Atan") }
void F_Atanh (F_ARG) { scalar_real_1_arg (atanh, "Atanh")} void F_Atanh(F_ARG) { scalar_real_1_arg(atanh, "Atanh") }
void F_Erf(F_ARG) { scalar_real_1_arg(erf, "Erf") }
#undef scalar_real_1_arg #undef scalar_real_1_arg
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C math functions (complex scalar, 1 argument) */ /* C math functions (complex scalar, 1 argument) */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
#define scalar_cmplx_1_arg(func, string) \ #define scalar_cmplx_1_arg(func, string) \
int k; \ int k; \
std::complex<double> tmp; \ std::complex<double> tmp; \
\ \
if(A->Type != SCALAR) \ if(A->Type != SCALAR) \
Message::Error("Non scalar argument for function '" string "'"); \ Message::Error("Non scalar argument for function '" string "'"); \
\ \
switch(Current.NbrHar){ \ switch(Current.NbrHar) { \
case 1: \ case 1: \
V->Val[0] = func(A->Val[0]) ; \ V->Val[0] = func(A->Val[0]); \
V->Val[MAX_DIM] = 0. ; \ V->Val[MAX_DIM] = 0.; \
break ; \ break; \
case 2: \ case 2: \
tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \ tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \
tmp = func(tmp); \ tmp = func(tmp); \
V->Val[0] = std::real(tmp); \ V->Val[0] = std::real(tmp); \
V->Val[MAX_DIM] = std::imag(tmp); \ V->Val[MAX_DIM] = std::imag(tmp); \
break; \ break; \
default: \ default: \
tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \ tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \
tmp = func(tmp); \ tmp = func(tmp); \
V->Val[0] = std::real(tmp); \ V->Val[0] = std::real(tmp); \
V->Val[MAX_DIM] = std::imag(tmp); \ V->Val[MAX_DIM] = std::imag(tmp); \
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) { \ for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) { \
V->Val[MAX_DIM*k] = std::real(tmp); \ V->Val[MAX_DIM * k] = std::real(tmp); \
V->Val[MAX_DIM*(k+1)] = std::imag(tmp); \ V->Val[MAX_DIM * (k + 1)] = std::imag(tmp); \
} \ } \
} \ } \
V->Type = SCALAR; \ V->Type = SCALAR;
void F_Exp (F_ARG) { scalar_cmplx_1_arg (exp, "Exp") } void F_Exp(F_ARG) { scalar_cmplx_1_arg(exp, "Exp") }
void F_Log (F_ARG) { scalar_cmplx_1_arg (log, "Log") } void F_Log(F_ARG) { scalar_cmplx_1_arg(log, "Log") }
void F_Log10 (F_ARG) { scalar_cmplx_1_arg (log10,"Log10") } void F_Log10(F_ARG) { scalar_cmplx_1_arg(log10, "Log10") }
void F_Sqrt (F_ARG) { scalar_cmplx_1_arg (sqrt, "Sqrt") } void F_Sqrt(F_ARG) { scalar_cmplx_1_arg(sqrt, "Sqrt") }
void F_Sin (F_ARG) { scalar_cmplx_1_arg (sin, "Sin") } void F_Sin(F_ARG) { scalar_cmplx_1_arg(sin, "Sin") }
void F_Cos (F_ARG) { scalar_cmplx_1_arg (cos, "Cos" ) } void F_Cos(F_ARG) { scalar_cmplx_1_arg(cos, "Cos") }
void F_Tan (F_ARG) { scalar_cmplx_1_arg (tan, "Tan") } void F_Tan(F_ARG) { scalar_cmplx_1_arg(tan, "Tan") }
void F_Sinh (F_ARG) { scalar_cmplx_1_arg (sinh, "Sinh") } void F_Sinh(F_ARG) { scalar_cmplx_1_arg(sinh, "Sinh") }
void F_Cosh (F_ARG) { scalar_cmplx_1_arg (cosh, "Cosh") } void F_Cosh(F_ARG) { scalar_cmplx_1_arg(cosh, "Cosh") }
void F_Tanh (F_ARG) { scalar_cmplx_1_arg (tanh, "Tanh") } void F_Tanh(F_ARG) { scalar_cmplx_1_arg(tanh, "Tanh") }
void F_Abs (F_ARG) { scalar_cmplx_1_arg (std::abs, "Abs") } void F_Abs(F_ARG) { scalar_cmplx_1_arg(std::abs, "Abs") }
#undef scalar_cmplx_1_arg #undef scalar_cmplx_1_arg
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* C math functions (scalar, 2 arguments, imaginary part set to zero) */ /* C math functions (scalar, 2 arguments, imaginary part set to zero) */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
#define scalar_real_2_arg(func, string) \ #define scalar_real_2_arg(func, string) \
int k; \ int k; \
\ \
if(A->Type != SCALAR || (A+1)->Type != SCALAR) \ if(A->Type != SCALAR || (A + 1)->Type != SCALAR) \
Message::Error("Non scalar argument(s) for function '" string "'"); \ Message::Error("Non scalar argument(s) for function '" string "'"); \
if(Current.NbrHar > 1 && (A->Val[MAX_DIM] != 0. || \ if(Current.NbrHar > 1 && \
(A+1)->Val[MAX_DIM] != 0.)) \ (A->Val[MAX_DIM] != 0. || (A + 1)->Val[MAX_DIM] != 0.)) \
Message::Error("Function '" string "' only accepts real arguments");\ Message::Error("Function '" string "' only accepts real arguments"); \
\ \
V->Val[0] = func(A->Val[0], (A+1)->Val[0]) ; \ V->Val[0] = func(A->Val[0], (A + 1)->Val[0]); \
if (Current.NbrHar > 1){ \ if(Current.NbrHar > 1) { \
V->Val[MAX_DIM] = 0. ; \ V->Val[MAX_DIM] = 0.; \
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) { \ for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) { \
V->Val[MAX_DIM*k] = func(A->Val[0], (A+1)->Val[0]) ; \ V->Val[MAX_DIM * k] = func(A->Val[0], (A + 1)->Val[0]); \
V->Val[MAX_DIM*(k+1)] = 0. ; \ V->Val[MAX_DIM * (k + 1)] = 0.; \
} \ } \
} \ } \
V->Type = SCALAR; V->Type = SCALAR;
void F_Atan2 (F_ARG) { scalar_real_2_arg (atan2, "Atan2") } void F_Atan2(F_ARG) { scalar_real_2_arg(atan2, "Atan2") }
void F_Fmod (F_ARG) { scalar_real_2_arg (fmod, "Fmod") } void F_Fmod(F_ARG) { scalar_real_2_arg(fmod, "Fmod") }
#undef scalar_real_2_arg #undef scalar_real_2_arg
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* Sign function */ /* Sign function */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void F_Sign(F_ARG) void F_Sign(F_ARG)
{ {
int k; int k;
double x; double x;
if(A->Type != SCALAR) if(A->Type != SCALAR)
Message::Error("Non scalar argument for function 'Sign'"); Message::Error("Non scalar argument for function 'Sign'");
x = A->Val[0]; x = A->Val[0];
if(x >= 0.) if(x >= 0.)
V->Val[0] = 1.; V->Val[0] = 1.;
else if(x < 0.) else if(x < 0.)
V->Val[0] = -1.; V->Val[0] = -1.;
else else
V->Val[0] = 0.; V->Val[0] = 0.;
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* Min, Max */ /* Min, Max */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void F_Min(F_ARG) void F_Min(F_ARG)
{ {
int k; int k;
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function Min"); Message::Error("Non scalar argument(s) for function Min");
V->Val[0] = (A->Val[0] < (A+1)->Val[0]) ? A->Val[0] : (A+1)->Val[0]; V->Val[0] = (A->Val[0] < (A + 1)->Val[0]) ? A->Val[0] : (A + 1)->Val[0];
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_Max(F_ARG) void F_Max(F_ARG)
{ {
int k; int k;
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function Max"); Message::Error("Non scalar argument(s) for function Max");
V->Val[0] = (A->Val[0] > (A+1)->Val[0]) ? A->Val[0] : (A+1)->Val[0]; V->Val[0] = (A->Val[0] > (A + 1)->Val[0]) ? A->Val[0] : (A + 1)->Val[0];
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* Bessel functions jn, yn and their derivatives */ /* Bessel functions jn, yn and their derivatives */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void F_Jn(F_ARG) void F_Jn(F_ARG)
{ {
int k, n; int k, n;
double x; double x;
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for Bessel function of the first kind Message::Error(
'Jn'"); "Non scalar argument(s) for Bessel function of the first kind 'Jn'");
n = (int)A->Val[0]; n = (int)A->Val[0];
x = (A+1)->Val[0]; x = (A + 1)->Val[0];
V->Val[0] = jn(n, x); V->Val[0] = jn(n, x);
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_JnComplex(F_ARG) void F_JnComplex(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for Bessel function of the first kind Message::Error("Non scalar argument(s) for Bessel function of the first "
'JnComplex'"); "kind 'JnComplex'");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double xr = (A+1)->Val[0]; double xr = (A + 1)->Val[0];
double xi = (A+1)->Val[MAX_DIM]; double xi = (A + 1)->Val[MAX_DIM];
double valr, vali; double valr, vali;
BesselJnComplex(n, 1, xr, xi, &valr, &vali); BesselJnComplex(n, 1, xr, xi, &valr, &vali);
V->Val[0] = valr; V->Val[0] = valr;
V->Val[MAX_DIM] = vali; V->Val[MAX_DIM] = vali;
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(int k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_KnComplex(F_ARG) void F_KnComplex(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for modified Bessel function of the s Message::Error("Non scalar argument(s) for modified Bessel function of the "
econd kind 'KnComplex'"); "second kind 'KnComplex'");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double xr = (A+1)->Val[0]; double xr = (A + 1)->Val[0];
double xi = (A+1)->Val[MAX_DIM]; double xi = (A + 1)->Val[MAX_DIM];
double valr, vali; double valr, vali;
BesselKnComplex(n, 1, xr, xi, &valr, &vali); BesselKnComplex(n, 1, xr, xi, &valr, &vali);
V->Val[0] = valr; V->Val[0] = valr;
V->Val[MAX_DIM] = vali; V->Val[MAX_DIM] = vali;
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(int k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_Yn(F_ARG) void F_Yn(F_ARG)
{ {
int k, n; int k, n;
double x; double x;
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for Bessel function of the second 'Yn Message::Error(
'"); "Non scalar argument(s) for Bessel function of the second 'Yn'");
n = (int)A->Val[0]; n = (int)A->Val[0];
x = (A+1)->Val[0]; x = (A + 1)->Val[0];
V->Val[0] = yn(n, x); V->Val[0] = yn(n, x);
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
double dBessel(double *tab, int n, double x) double dBessel(double *tab, int n, double x)
{ {
if(n == 0){ if(n == 0) { return -tab[1]; }
return - tab[1]; else {
} return tab[n - 1] - (double)n / x * tab[n];
else{
return tab[n-1] - (double)n/x * tab[n];
} }
} }
void F_dJn(F_ARG) void F_dJn(F_ARG)
{ {
int k, n;
double x, *jntab;
int k, n; if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
double x, *jntab;
if(A->Type != SCALAR || (A+1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for the derivative of the Bessel " Message::Error("Non scalar argument(s) for the derivative of the Bessel "
"function of the first kind 'dJn'"); "function of the first kind 'dJn'");
n = (int)A->Val[0]; n = (int)A->Val[0];
x = (A+1)->Val[0]; x = (A + 1)->Val[0];
jntab = (double*)Malloc((n + 2) * sizeof(double)); jntab = (double *)Malloc((n + 2) * sizeof(double));
for(k = 0; k < n + 2; k++){ for(k = 0; k < n + 2; k++) { jntab[k] = jn(k, x); }
jntab[k] = jn(k, x);
}
V->Val[0] = dBessel(jntab, n, x); V->Val[0] = dBessel(jntab, n, x);
Free(jntab); Free(jntab);
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
void F_dYn(F_ARG) void F_dYn(F_ARG)
{ {
int k, n; int k, n;
double x, *yntab; double x, *yntab;
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for the derivative of the Bessel " Message::Error("Non scalar argument(s) for the derivative of the Bessel "
"function of the second kind 'dYn'"); "function of the second kind 'dYn'");
n = (int)A->Val[0]; n = (int)A->Val[0];
x = (A+1)->Val[0]; x = (A + 1)->Val[0];
yntab = (double*)Malloc((n + 2) * sizeof(double)); yntab = (double *)Malloc((n + 2) * sizeof(double));
for(k = 0; k < n + 2; k++){ for(k = 0; k < n + 2; k++) { yntab[k] = yn(k, x); }
yntab[k] = yn(k, x);
}
V->Val[0] = dBessel(yntab, n, x); V->Val[0] = dBessel(yntab, n, x);
Free(yntab); Free(yntab);
if (Current.NbrHar > 1){ if(Current.NbrHar > 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR; V->Type = SCALAR;
} }
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
/* Spherical Bessel functions jn, yn and their derivatives */ /* Spherical Bessel functions jn, yn and their derivatives */
/* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */
void F_JnSph(F_ARG) void F_JnSph(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function 'JnSph' (spherical Besse Message::Error("Non scalar argument(s) for function 'JnSph' (spherical "
l function)"); "Bessel function)");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double x = (A+1)->Val[0]; double x = (A + 1)->Val[0];
V->Type = SCALAR; V->Type = SCALAR;
V->Val[0] = Spherical_j_n(n, x); V->Val[0] = Spherical_j_n(n, x);
} }
void F_YnSph(F_ARG) void F_YnSph(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function 'YnSph' (spherical Besse Message::Error("Non scalar argument(s) for function 'YnSph' (spherical "
l function)"); "Bessel function)");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double x = (A+1)->Val[0]; double x = (A + 1)->Val[0];
V->Type = SCALAR; V->Type = SCALAR;
V->Val[0] = Spherical_y_n(n, x); V->Val[0] = Spherical_y_n(n, x);
} }
void F_dJnSph(F_ARG) void F_dJnSph(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function 'dJnSph' (derivative of Message::Error("Non scalar argument(s) for function 'dJnSph' (derivative "
spherical Bessel function)"); "of spherical Bessel function)");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double x = (A+1)->Val[0]; double x = (A + 1)->Val[0];
V->Type = SCALAR; V->Type = SCALAR;
V->Val[0] = (n/x) * Spherical_j_n(n, x) - Spherical_j_n(n+1, x); V->Val[0] = (n / x) * Spherical_j_n(n, x) - Spherical_j_n(n + 1, x);
} }
void F_dYnSph(F_ARG) void F_dYnSph(F_ARG)
{ {
if(A->Type != SCALAR || (A+1)->Type != SCALAR) if(A->Type != SCALAR || (A + 1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for function 'dYnSph' (derivative of Message::Error("Non scalar argument(s) for function 'dYnSph' (derivative "
spherical Bessel function)"); "of spherical Bessel function)");
int n = (int)A->Val[0]; int n = (int)A->Val[0];
double x = (A+1)->Val[0]; double x = (A + 1)->Val[0];
V->Type = SCALAR; V->Type = SCALAR;
V->Val[0] = (n/x) * Spherical_y_n(n, x) - Spherical_y_n(n+1, x); V->Val[0] = (n / x) * Spherical_y_n(n, x) - Spherical_y_n(n + 1, x);
} }
 End of changes. 50 change blocks. 
189 lines changed or deleted 183 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)