#ifndef CYLFUNC_H
#define CYLFUNC_H

/*
 * METRIC --- Mode expansion modeling in integrated optics / photonics
 * http://metric.computational-photonics.eu/ 
 */

/*
 * cylfunc.h
 * Bessel-  and Hankel-functions 
 */

#include"complex.h"

/* ---------------------------------------------------------------------- */

/*
 * Bessel and Hankel functions of integer order, real argument
 */

/* Bessel function of the first kind J_0, real argument, modified version of
   float bessj0(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_J0(const double& x);

/* Bessel function of the first kind J_1, real argument, modified version of
   float bessj1(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_J1(const double& x);

/* Bessel function of the first kind J_n, integer order n>=2, real argument, modified version of
   float bessj(int n, float x)
   (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_J(int n, const double& x);

/* Bessel function of the second kind Y_0, real argument, modified version of
   float bessy0(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_Y0(const double& x);

/* Bessel function of the second kind Y_1, real argument, modified version of
   float bessy1(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_Y1(const double& x);

/* Bessel function of the second kind Y_n, integer order n>=2, real argument, modified version of
   float bessy(int n, float x)
   (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_Y(int n, const double& x);

/* modified Bessel function of the first kind I_0, real argument, modified version of
   float bessi0(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_I0(const double& x);

/* modified Bessel function of the first kind I_1, real argument, modified version of
   float bessi1(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_I1(const double& x);

/* modified Bessel function of the first kind I_n, integer order n>=2, real argument, modified version of
   float bessi(int n, float x)
   (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_I(int n, const double& x);

/* modified Bessel function of the second kind K_0, real argument, modified version of
   float bessk0(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_K0(const double& x);

/* modified Bessel function of the second kind K_1, real argument, modified version of
   float bessk1(float x), (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_K1(const double& x);

/* modified Bessel function of the second kind K_n, integer order n>=2, real argument, modified version of
   float bessk(int n, float x)
   (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */
double bessel_K(int n, const double& x);

/* Hankel function H_n^(2), integer order n>=0, real argument */
Complex hankel_H2(int n, const double& x);

/* Hankel function H_n^(2), integer order n>=0, real argument,
   first derivative  */
Complex d_hankel_H2(int n, const double& x);

/* ---------------------------------------------------------------------- */

/*
 * Bessel- and Hankel-functions of complex order, for complex arguments
 * based on: File mod_zbes.f90
 *
 * Author: Masao Kodama
 * Address: 21-20, Gakuen 1 chome, Mastue-shi, Shimane-ken, 690-0825 Japan
 * Email: mkodama@mable.ne.jp
 *
 * Masao Kodama, "Algorithm 912: A module for calculating cylindrical
 * functions of complex order and complex argument," ACM Transactions on
 * Mathematical Software 37 (4), Article No. 47 
 *
 * & references cited therein
 */

/* status code info is set upon return:
    0: normal completion 
   10: unreliable output: potential underflow/overflow, insufficient precision, 
       theoretically indefinite result (e.g. at origin)
   20: output out of range */

/* Bessel function J_n, complex order znu, complex argument zz */
Complex bessel_J(const Complex& znu, const Complex& zz, int& info);

/* Bessel function Y_n, complex order znu, complex argument zz */
Complex bessel_Y(const Complex& znu, const Complex& zz, int& info);

/* Hankel function H_n^(1), complex order znu, complex argument zz */
Complex hankel_H1(const Complex& znu, const Complex& zz, int& info);

/* Hankel function H_n^(2), complex order znu, complex argument zz */
Complex hankel_H2(const Complex& znu, const Complex& zz, int& info);

/* & derivatives of these with respect to the argument */
Complex  d_bessel_J(const Complex& znu, const Complex& zz, int& info);
Complex  d_bessel_Y(const Complex& znu, const Complex& zz, int& info);
Complex d_hankel_H1(const Complex& znu, const Complex& zz, int& info);
Complex d_hankel_H2(const Complex& znu, const Complex& zz, int& info);

/* ---------------------------------------------------------------------- */

/*
 * Bessel and Hankel functions of complex order and complex arguments
 * 
 * Circurs - Circular optical microresonators simulation tools
 * Kirankumar R. Hiremath (k.r.hiremath@ieee.org)
 * University of Twente, Department of Applied Mathematics 
 * P.O. Box 217, 7500AE Enschede, The Netherlands
 * (2005)
 * 
 * - modified -
 */

/*
  common computational parameters & indicators:
	
  KODE = 1 => computation without scaling
  KODE = 2 => computation with scaling

  Error codes:
  IERR =  0 => :-) Success
  IERR =  1 => inout error
  IERR =  2 => overflow. no computation. Re(\zeta) too large for KODE = 1
  IERR =  3 => |z| too large. ans may wrong.
  IERR =  4 => |z| too large. no computation.
  IERR =  5 => no computation, algorithm termination condition not met
  IERR = 11 => situation beyond scope of routines
  IERR = 12 => division by zero

  NZ: Underflow indicator. 0: Normal return, 1: underflow
*/

/* Bessel function of the first kind */
Complex BesselJ(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* Bessel function of the second kind */
Complex BesselY(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* derivative of the Bessel function of the first kind */
Complex DBesselJ(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* derivative of the Bessel function of the second kind */
Complex DBesselY(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* Hankel function of the first kind */
Complex HankelH1(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* Hankel function of the second kind */
Complex HankelH2(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* derivative of Hankel function of the first kind */
Complex DHankelH1(Complex order, Complex argument, int KODE, int *NZ, int *IERR);
/* derivative of Hankel function of the second kind */ 
Complex DHankelH2(Complex order, Complex argument, int KODE, int *NZ, int *IERR);

#endif   // CYLFUNC_H