188 lines
5.0 KiB
C++
188 lines
5.0 KiB
C++
// Copied from Jean-Pierre Moreau's Home Page
|
|
// I would like to thank him for sharing his code.
|
|
/*****************************************************************
|
|
* Complex Bessel Function of the 1st Kind of integer order *
|
|
* -------------------------------------------------------------- *
|
|
* SAMPLE RUN: *
|
|
* *
|
|
* Complex Bessel Function of the 1st Kind of integer order *
|
|
* *
|
|
* Input complex argument (real imaginary): 1 2 *
|
|
* Input integer order: 1 *
|
|
* *
|
|
* Function value: 1.291848 1.010488 *
|
|
* *
|
|
* *
|
|
* C++ Release 1.2 By J-P Moreau, Paris. *
|
|
* (www.jpmoreau.fr) *
|
|
* -------------------------------------------------------------- *
|
|
* Release 1.1: Corrected bug in Function ZLn (atan replaced by *
|
|
* function ATAN2) 11/10/2005. *
|
|
* Release 1.2: ZPower replaced by IZpower (integer exponant). *
|
|
* Limitations increased in CDIV, MAXK=20 *
|
|
*****************************************************************/
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
|
|
#define MAXK 20 //09/21/2009 (15 before)
|
|
#define FALSE 0
|
|
#define TRUE 1
|
|
|
|
double HALF = 0.5, ONE = 1.0, FPF = 5.5;
|
|
double PI=4.0*atan(1.0);
|
|
|
|
int nu; // order of complex Bessel
|
|
double z[2], z1[2]; // Complex numbers
|
|
double x,y;
|
|
|
|
// Z=Z1/Z2
|
|
void CDIV(double *Z1, double *Z2, double *Z) {
|
|
double D;
|
|
D=Z2[0]*Z2[0]+Z2[1]*Z2[1];
|
|
if (D > 1e30) {
|
|
Z[0]=0.0; Z[1]=0.0; return;
|
|
}
|
|
if (D == 0) return;
|
|
Z[0]=(Z1[0]*Z2[0]+Z1[1]*Z2[1])/D;
|
|
Z[1]=(Z1[1]*Z2[0]-Z1[0]*Z2[1])/D;
|
|
}
|
|
|
|
// Z=Z1*Z2
|
|
void CMUL(double *Z1, double *Z2, double *Z) {
|
|
Z[0]=Z1[0]*Z2[0] - Z1[1]*Z2[1];
|
|
Z[1]=Z1[0]*Z2[1] + Z1[1]*Z2[0];
|
|
}
|
|
|
|
// compute Z^N
|
|
void IZPower(double *z, int n, double *z1) {
|
|
double temp[2],temp1[2]; int i;
|
|
if (n==0) {
|
|
z1[0]=1.0;
|
|
z1[1]=0.0;
|
|
}
|
|
else if (n==1) {
|
|
z1[0]=z[0];
|
|
z1[1]=z[1];
|
|
}
|
|
else {
|
|
temp1[0]=z[0]; temp1[1]=z[1];
|
|
for (i=2; i<=n; i++) {
|
|
CMUL(temp1,z,temp);
|
|
temp1[0]=temp[0];
|
|
temp1[1]=temp[1];
|
|
}
|
|
z1[0]=temp[0];
|
|
z1[1]=temp[1];
|
|
}
|
|
}
|
|
|
|
double Fact(int k) {
|
|
int i; double f;
|
|
f=ONE;
|
|
for (i=2; i<=k; i++) f *= ONE*i;
|
|
return f;
|
|
}
|
|
|
|
/******************************************
|
|
* FUNCTION GAMMA(X) *
|
|
* --------------------------------------- *
|
|
* Returns the value of Gamma(x) in double *
|
|
* precision as EXP(LN(GAMMA(X))) for X>0. *
|
|
******************************************/
|
|
double Gamma(double xx) {
|
|
double cof[6];
|
|
double stp,x,tmp,ser;
|
|
int j;
|
|
cof[0]=76.18009173;
|
|
cof[1]=-86.50532033;
|
|
cof[2]=24.01409822;
|
|
cof[3]=-1.231739516;
|
|
cof[4]=0.120858003e-2;
|
|
cof[5]=-0.536382e-5;
|
|
stp=2.50662827465;
|
|
x=xx-ONE;
|
|
tmp=x+FPF;
|
|
tmp=(x+HALF)*log(tmp)-tmp;
|
|
ser=ONE;
|
|
for (j=0; j<6; j++) {
|
|
x += ONE;
|
|
ser += cof[j]/x;
|
|
}
|
|
return (exp(tmp+log(stp*ser)));
|
|
}
|
|
|
|
// main subroutine
|
|
void CBESSJ(double *z, int nu, double *z1) {
|
|
/*--------------------------------------------------
|
|
inf. (-z^2/4)^k
|
|
Jnu(z) = (z/2)^nu x Sum ------------------
|
|
k=0 k! x Gamma(nu+k+1)
|
|
(nu must be >= 0). Here k=15.
|
|
---------------------------------------------------*/
|
|
int k;
|
|
double sum[2],tmp[2],tmp1[2],tmp2[2],tmp3[2];
|
|
// calculate (z/2)^nu in tmp3
|
|
tmp[0]=2.0; tmp[1]=0.0;
|
|
CDIV(z,tmp,tmp1);
|
|
IZPower(tmp1,nu,tmp3);
|
|
sum[0]=0.0; sum[1]=0.0;
|
|
//calculate Sum
|
|
for (k=0; k<=MAXK; k++) {
|
|
// calculate (-z^2/4)^k
|
|
IZPower(z,2,tmp);
|
|
tmp[0]=-tmp[0]; tmp[1]=-tmp[1];
|
|
tmp1[0]=4.0; tmp1[1]=0.0;
|
|
CDIV(tmp,tmp1,tmp2);
|
|
|
|
if (z[1]==0) { //case real number
|
|
tmp[0]=pow(tmp2[0],k);
|
|
tmp[1]=0.0;
|
|
}
|
|
else //case complex number
|
|
IZPower(tmp2,k,tmp);
|
|
|
|
// divide by k!
|
|
tmp1[0]=Fact(k); tmp1[1]=0.0;
|
|
if (z[1]==0) {
|
|
tmp2[0]=tmp[0]/tmp1[0];
|
|
tmp2[1]=0.0;
|
|
}
|
|
else
|
|
CDIV(tmp,tmp1,tmp2);
|
|
// divide by Gamma(nu+k+1)
|
|
tmp1[0]=Gamma(ONE*(nu+k+1)); tmp1[1]=0.0;
|
|
|
|
if (z[1]==0) {
|
|
tmp[0]=tmp2[0]/tmp1[0];
|
|
tmp[1]=0.0;
|
|
}
|
|
else
|
|
CDIV(tmp2,tmp1,tmp);
|
|
// actualize sum
|
|
sum[0] += tmp[0];
|
|
sum[1] += tmp[1];
|
|
|
|
}
|
|
// multiply (z/2)^nu by sum
|
|
CMUL(tmp3,sum,z1);
|
|
}
|
|
|
|
|
|
// int main() {
|
|
|
|
// double z1[2];
|
|
|
|
// printf("\n Complex Bessel Function of the 1st Kind of integer order\n\n");
|
|
// printf(" Input complex argument (real imaginary): ");
|
|
// scanf("%lf %lf", &x, &y);
|
|
// z[0]=x; z[1]=y;
|
|
// printf(" Input integer order: "); scanf("%d", &nu);
|
|
|
|
// CBESSJ(z,nu,z1);
|
|
|
|
// printf("\n Function value: %f %f\n\n", z1[0], z1[1]);
|
|
|
|
// }
|
|
|
|
// end of file cbessj.cpp
|