Commit 0acc57db authored by Christophe Geuzaine's avatar Christophe Geuzaine

test Bessel Jn with complex args

parent c5b71447
Pipeline #1521 passed with stage
in 9 minutes 32 seconds
......@@ -948,6 +948,7 @@ struct StringXFunction2Nbr F_Function[] = { /* #Par #Arg */
{"Min" , (CAST)F_Min , 0, 2 },
{"Max" , (CAST)F_Max , 0, 2 },
{"Jn" , (CAST)F_Jn , 0, 2 },
{"JnComplex" , (CAST)F_JnComplex , 0, 2 },
{"Yn" , (CAST)F_Yn , 0, 2 },
{"dJn" , (CAST)F_dJn , 0, 2 },
{"dYn" , (CAST)F_dYn , 0, 2 },
......
......@@ -129,6 +129,7 @@ void F_Sign (F_ARG) ;
void F_Min (F_ARG) ;
void F_Max (F_ARG) ;
void F_Jn (F_ARG) ;
void F_JnComplex (F_ARG) ;
void F_Yn (F_ARG) ;
void F_dJn (F_ARG) ;
void F_dYn (F_ARG) ;
......@@ -341,7 +342,7 @@ void Vector_Find_Jk_K (const double hrk[3],
void Vector_Find_hrk_K(const double Jk[3],
const double Ja, const double ha, const double Jb, const double hb,
double hrk[3]);
void Tensor_dJkdhrk_K(const double hr[3],
void Tensor_dJkdhrk_K(const double hr[3],
const double Ja, const double ha, const double Jb, const double hb,
double mutg[6]);
......
......@@ -46,12 +46,25 @@ void F_Hypot(F_ARG)
void F_TanhC2(F_ARG)
{
// @Jon: check this and the behavior of the overall function for large args
// lim_x\to\inf cosh(x) = +\inf
// lim_x\to\inf sinh(x) = 1
double denom =
SQU(cosh(A->Val[0])*cos(A->Val[MAX_DIM])) +
SQU(sinh(A->Val[0])*sin(A->Val[MAX_DIM]));
printf("arg=%g cosh(arg)=%g\n", A->Val[0], cosh(A->Val[0]));
V->Val[0] = sinh(A->Val[0])*cosh(A->Val[0]) / denom ;
V->Val[MAX_DIM] = sin(A->Val[MAX_DIM])*cos(A->Val[MAX_DIM]) / denom ;
V->Type = SCALAR ;
printf("numer_real=%g, numer_imag=%g, denom = %g\n",
sinh(A->Val[0])*cosh(A->Val[0]) ,
sin(A->Val[MAX_DIM])*cos(A->Val[MAX_DIM]),
denom);
}
/* ------------------------------------------------------------------------ */
......
......@@ -9,6 +9,7 @@
#include "F.h"
#include "MallocUtils.h"
#include "Message.h"
#include "Bessel.h"
extern struct CurrentData Current ;
......@@ -198,6 +199,27 @@ void F_Jn(F_ARG)
V->Type = SCALAR;
}
void F_JnComplex(F_ARG)
{
if(A->Type != SCALAR || (A+1)->Type != SCALAR)
Message::Error("Non scalar argument(s) for Bessel function of the first kind 'JnComplex'");
int n = (int)A->Val[0];
double xr = (A+1)->Val[0];
double xi = (A+1)->Val[MAX_DIM];
double valr, vali;
BesselJnComplex(n, 1, xr, xi, &valr, &vali);
V->Val[0] = valr;
V->Val[MAX_DIM] = vali;
if (Current.NbrHar > 1){
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->Type = SCALAR;
}
void F_Yn(F_ARG)
{
int k, n;
......
......@@ -9,19 +9,19 @@
#if defined(HAVE_NO_FORTRAN)
static void zbesj_(double*, double*, double*, int*, int*, double*,
static void zbesj_(double*, double*, double*, int*, int*, double*,
double*, int*, int*)
{
Message::Fatal("Bessel functions require Fortran compiler");
}
static void zbesy_(double*, double*, double*, int*, int*, double*,
static void zbesy_(double*, double*, double*, int*, int*, double*,
double*, int*, double*, double*, int*)
{
Message::Fatal("Bessel functions require Fortran compiler");
}
static void zbesh_(double*, double*, double*, int*, int*, int*,
static void zbesh_(double*, double*, double*, int*, int*, int*,
double*, double*, int*, int*)
{
Message::Fatal("Bessel functions require Fortran compiler");
......@@ -36,12 +36,12 @@ static void zbesh_(double*, double*, double*, int*, int*, int*,
#endif
extern "C" {
void zbesj_(double*, double*, double*, int*, int*, double*,
void zbesj_(double*, double*, double*, int*, int*, double*,
double*, int*, int*);
void zbesy_(double*, double*, double*, int*, int*, double*,
double*, int*, double*,
double*, int*);
void zbesh_(double*, double*, double*, int*, int*, int*, double*,
void zbesh_(double*, double*, double*, int*, int*, int*, double*,
double*, int*, int*);
}
......@@ -84,14 +84,23 @@ int BesselJn(double n, int num, double x, double *val)
int nz = 0, ierr = 0, kode = 1;
double xi = 0.0;
double* ji = new double[num];
zbesj_(&x, &xi, &n, &kode, &num, val, ji, &nz, &ierr) ;
delete[] ji;
return BesselError(ierr, "BesselJn");
}
int BesselJnComplex(double n, int num, double xr, double xi, double *valr, double *vali)
{
int nz = 0, ierr = 0, kode = 1;
zbesj_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr) ;
return BesselError(ierr, "BesselJnComplex");
}
int BesselSphericalJn(double n, int num, double x, double *val)
{
int ierr = BesselJn(n+0.5, num, x, val);
......@@ -124,7 +133,7 @@ int BesselYn(double n, int num, double x, double *val)
zbesy_(&x, &xi, &n, &kode, &num, val, yi, &nz, auxyr, auxyi, &ierr);
delete[] yi;
delete[] yi;
delete[] auxyr;
delete[] auxyi;
......@@ -158,7 +167,7 @@ int BesselHn(int type, double n, int num, double x, std::complex<double> *val)
int nz = 0, ierr = 0, kode = 1;
double* hr = new double[num];
double* hi = new double[num];
double xi = 0.0;
double xi = 0.0;
zbesh_(&x, &xi, &n, &kode, &type, &num, hr, hi, &nz, &ierr);
......@@ -221,4 +230,3 @@ double AltSpherical_y_n(int n, double x)
BesselAltSphericalYn(n, 1, x, &res);
return res;
}
......@@ -22,6 +22,8 @@ int BesselJn(double n, int num, double x, double *val);
int BesselSphericalJn(double n, int num, double x, double *val);
int BesselAltSphericalJn(double n, int num, double x, double *val);
int BesselJnComplex(double n, int num, double xr, double xi, double *valr, double *vali);
int BesselYn(double n, int num, double x, double *val);
int BesselSphericalYn(double n, int num, double x, double *val);
int BesselAltSphericalYn(double n, int num, double x, double *val);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment