cbessj_cpp.txt

BESSEL PROGRAMS IN C/C

/*****************************************************************
*    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.1 By J-P Moreau, Paris. *
* -------------------------------------------------------------- *
* Release 1.1: Corrected bug in Function ZLn (atan replaced by   *
*              function ATAN2) 11/10/2005.                       *
*****************************************************************/
#include <stdio.h>
#include <math.h>

#define  MAXK  15
#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;

  void ATAN2(double numerator,double denominator, double *phase) {
  //Returns an argument (phase) between 0 and 2*PI
    if (fabs(denominator)<1.0e-6) {
      if (numerator < 0.0) *phase=-PI/2;
      else *phase=PI/2;
      return;
    }
    else {
      *phase=atan(numerator/denominator);
      if (denominator > 0.0) {
        if (numerator >=0.0) return;
        else {
          *phase = *phase + 2*PI;
          return;
        }
      }
      else
        *phase = *phase + PI;
    } //else
  } //ATAN2()

  // compute z = Ln(z1)
  int ZLn(double *z1, double *z) {
    double r,t;
    if (z1[0] == 0.0 && z1[1] == 0.0)
      return FALSE;
    else {
      r=sqrt(z1[0]*z1[0]+z1[1]*z1[1]);
      if (z1[0] != 0.0)
        //t=atan(z1[1]/z1[0]);
	ATAN2(z1[1],z1[0],&t);
      else
        t=PI/2.0;
      z[0]=log(r);
      z[1]=t;
      return TRUE;
    }
  }

  // Z=Z1/Z2
  void CDIV(double *Z1, double *Z2, double *Z) {
    double D;
    D=Z2[0]*Z2[0]+Z2[1]*Z2[1];
	if (D > 1e12) {
      Z[0]=0.0; Z[1]=0.0; return;
    }
    if (D < 1e-12) 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];
  }

  // z1=exp(z)
  void CEXP(double *Z, double *Z1) {
    double tmp;
    if (exp(Z[0]) > 1e16) tmp=1e16; else tmp=exp(Z[0]);
    Z1[0]=tmp*cos(Z[1]);
    Z1[1]=tmp*sin(Z[1]);
  }

  // compute Z1^Z2
  void ZPower(double *z1, double *z2, double *z) {
    double z3[2];
    if (ZLn(z1,z3)) {
      CMUL(z3,z2,z3);
      CEXP(z3,z);
    }
  }

  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 tmp2
  tmp[0]=2.0; tmp[1]=0.0;
  CDIV(z,tmp,tmp1);
  tmp[0]=(double) ONE*nu;
  ZPower(tmp1,tmp,tmp3);
  sum[0]=0.0; sum[1]=0.0;
  //calculate Sum
  for (k=0; k<=MAXK; k++) {
    // calculate (-z^2/4)^k
    tmp1[0]=2.0; tmp1[1]=0.0;
    ZPower(z,tmp1,tmp);
    tmp[0]=-tmp[0]; tmp[1]=-tmp[1];
    tmp1[0]=4.0; tmp1[1]=0.0;
    CDIV(tmp,tmp1,tmp2);
    tmp1[0]=1.0*k; tmp1[1]=0.0;
    ZPower(tmp2,tmp1,tmp);
    // divide by k!
    tmp1[0]=Fact(k); tmp1[1]=0.0;
    CDIV(tmp,tmp1,tmp2);
    // divide by Gamma(nu+k+1)
    tmp1[0]=Gamma(ONE*(nu+k+1)); tmp1[1]=0.0;
    CDIV(tmp2,tmp1,tmp);
    // actualize sum
    sum[0] += tmp[0];
    sum[1] += tmp[1];
  }
  // multiply (z/2)^nu by sum
  CMUL(tmp3,sum,z1);
}


void main() {

  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