#include <iostream>
#include <complex>
#include "besselh.h"
extern "C" {
#include "nr.h"
}
using namespace std;
/* 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-07-30
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----
*/

std::complex<double> besselh( double nu, int kind, double z ){
  float rj,ry,rjp,ryp;

  bessjy( (float) z, (float) nu, &rj, &ry, &rjp, &ryp ); 

  switch( kind ){
  case 1:
    return complex<double> ( rj, +ry );
  case 2:
    return complex<double> ( rj, -ry );
  default:
    std::cout << "besselh: Unknown hankel function kind" << std::endl; 
    exit(1); 
  }
}

#ifdef _BUILD_TESTS_
int main(){
  double nu, z;
  int kind; 

  nu = 1; kind = 1; z = 1; 
  std::cout << "besselh( " << nu << ", " << kind << ", " << z << ") = " << 
    besselh( nu, kind, z ) << std::endl; 

  nu = 1.5; kind = 1; z = 2;
  std::cout << "besselh( " << nu << ", " << kind << ", " << z << ") = " << 
    besselh( nu, kind, z ) << std::endl; 

  nu = 1; kind = 2; z = 1; 
  std::cout << "besselh( " << nu << ", " << kind << ", " << z << ") = " << 
    besselh( nu, kind, z ) << std::endl; 

  nu = 1.5; kind = 2; z = 2; 
  std::cout << "besselh( " << nu << ", " << kind << ", " << z << ") = " << 
    besselh( nu, kind, z ) << std::endl; 
}
#endif