/**********************************************************
* This program tests the subroutine BESSK to calculate the*
* modified Bessel function of the third kind of order N   *
* for any real positive argument X.                       *
* ------------------------------------------------------- *
* SAMPLE RUN:                                             *
*                                                         *
*  X = 1.209700                                           *
*                                                         *
*  N = 0                                                  *
*  Y = 0.314324                                           *
*                                                         *
*  N = 1                                                  *
*  Y = 0.428051                                           *
*                                                         *
*  N = 2                                                  *
*  Y = 1.022022                                           *
*                                                         *
*  N = 3                                                  *
*  Y = 3.807473                                           *
*                                                         *
*  N = 4                                                  *
*  Y = 19.906737                                          *
*                                                         *
* ------------------------------------------------------- *
* Reference:   From Numath Library By Tuan Dang Trong     *
* in Fortran 77 [BIBLI 18].                               *
*                                                         *
*                      C++ Version By J-P Moreau, Paris.  *
*                             (www.jpmoreau.fr)           *
**********************************************************/
#include <stdio.h>
#include <math.h>

#define BIG  1e30

    double X, Y;
    int i, N;

    double BESSK0(double X);
    double BESSK1(double X);
    double BESSI0(double X);
    double BESSI1(double X);


    double BESSK(int N, double X) {
/* -----------------------------------------------------------------------
!     CE SOUS-PROGRAMME CALCULE LA FONCTION BESSEL MODIFIFIEE 3E ESPECE
!     D'ORDRE N ENTIER POUR TOUT X REEL POSITIF > 0.  ON UTILISE ICI LA
!     FORMULE DE RECURRENCE CLASSIQUE EN PARTANT DE BESSK0 ET BESSK1.
!
!     THIS ROUTINE CALCULATES THE MODIFIED BESSEL FUNCTION OF THE THIRD
!     KIND OF INTEGER ORDER, N FOR ANY POSITIVE REAL ARGUMENT, X. THE
!     CLASSICAL RECURSION FORMULA IS USED, STARTING FROM BESSK0 AND BESSK1.
! ------------------------------------------------------------------------ 
!     REFERENCE:
!     C.W.CLENSHAW, CHEBYSHEV SERIES FOR MATHEMATICAL FUNCTIONS,
!     MATHEMATICAL TABLES, VOL.5, 1962.
! ----------------------------------------------------------------------*/
      double TOX,BK,BKM,BKP;
      int J;
      if (N == 0) return BESSK0(X);
      if (N == 1) return BESSK1(X);
      if (X == 0.0) return BIG;  //arbitrary big value
      TOX = 2.0/X;
      BK  = BESSK1(X);
      BKM = BESSK0(X);
      for (J=1; J<N; J++) {
        BKP = BKM + J*TOX*BK;
        BKM = BK;
        BK  = BKP;
      }
      return BK;
    }

    double BESSK0(double X)  {
/* ---------------------------------------------------------------------
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DU 3EME ESPECE D'ORDRE 0
!     POUR TOUT X REEL NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ZERO FOR ANY POSITIVE REAL ARGUMENT, X.
! --------------------------------------------------------------------*/
      double Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,tmp;
      P1=-0.57721566; P2= 0.42278420; P3=0.23069756; P4=0.3488590e-1;
      P5= 0.262698e-2; P6=0.10750e-3; P7=0.74e-5;
      Q1= 1.25331414; Q2=-0.7832358e-1; Q3=0.2189568e-1; Q4=-0.1062446e-1;
      Q5= 0.587872e-2; Q6=-0.251540e-2; Q7=0.53208e-3;
      if (X == 0.0) return BIG;  //arbitrary big value
      if (X <= 2.0) {
        Y=X*X/4.0;
        AX=-log(X/2.0)*BESSI0(X);
        tmp = AX+(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))));
	return tmp;
      }
      else {
        Y=2.0/X;
        AX=exp(-X)/sqrt(X);
        tmp = AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))));
	return tmp;
      }
    }

    double BESSK1(double X)  {
/* ------------------------------------------------------------------------
!     CALCUL DE LA FONCTION BESSEL MODIFIEE DE 3EME ESPECE D'ORDRE 1
!     POUR TOUT X REEL POSITF NON NUL.
!
!     CALCULATES THE THE MODIFIED BESSEL FUNCTION OF THE THIRD KIND OF 
!     ORDER ONE FOR ANY POSITIVE REAL ARGUMENT, X.
! -----------------------------------------------------------------------*/
      double Y,AX,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7;
      P1=1.0; P2=0.15443144; P3=-0.67278579; P4=-0.18156897;
      P5=-0.1919402e-1; P6=-0.110404e-2; P7=-0.4686e-4;
      Q1=1.25331414; Q2=0.23498619; Q3=-0.3655620e-1; Q4=0.1504268e-1;
      Q5=-0.780353e-2; Q6=0.325614e-2; Q7=-0.68245e-3;
      if (X == 0.0) return BIG;
      if (X <= 2.0) {
        Y=X*X/4.0;
        AX=log(X/2.0)*BESSI1(X);
        return (AX+(1.0/X)*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))));
      }
      else {
        Y=2.0/X;
        AX=exp(-X)/sqrt(X);
        return (AX*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7)))))));
      }
    }

//  Bessel Function of the 1st kind of order zero

    double BESSI0(double X)  {
      double Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX;
      P1=1.0; P2=3.5156229; P3=3.0899424; P4=1.2067429;
      P5=0.2659732; P6=0.360768e-1; P7=0.45813e-2;
      Q1=0.39894228; Q2=0.1328592e-1; Q3=0.225319e-2; Q4=-0.157565e-2;
      Q5=0.916281e-2; Q6=-0.2057706e-1; Q7=0.2635537e-1;
      Q8=-0.1647633e-1; Q9=0.392377e-2;
      if (fabs(X) < 3.75) {
        Y=pow((X/3.75),2);
        return (P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))));
      }
      else {
        AX=fabs(X);
        Y=3.75/AX;
        BX=exp(AX)/sqrt(AX);
        AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))));
        return AX*BX;
      }
    }

//  Bessel Function of the 1st kind of order one

    double BESSI1(double X)  {
    double Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,AX,BX;
      P1=0.5; P2=0.87890594; P3=0.51498869; P4=0.15084934;
      P5=0.2658733e-1; P6=0.301532e-2; P7=0.32411e-3;
      Q1=0.39894228; Q2=-0.3988024e-1; Q3=-0.362018e-2;
      Q4=0.163801e-2;Q5=-0.1031555e-1; Q6= 0.2282967e-1;
      Q7=-0.2895312e-1; Q8=0.1787654e-1; Q9=-0.420059e-2;
      if (fabs(X) < 3.75) {
        Y=pow((X/3.75),2);
        return (X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))));
      }
      else {
        AX=fabs(X);
        Y=3.75/AX;
        BX=exp(AX)/sqrt(AX);
        AX=Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))));
        return AX*BX;
      }
    }


void main()  {

   N=5;           // number of integer orders
   X=1.2097;      // value of real argument

   printf("\n X = %f\n\n", X);
   for (i=0; i<N; i++) {
     Y=BESSK(i,X);
     printf(" N = %d\n", i);
     printf(" Y = %f\n\n", Y);
   }
   printf("\n");
}

// End of file Tbessk.cpp