specfun.cpp

Go to the documentation of this file.

// $Id: specfun.cpp,v 1.3.2.9 2008/07/17 17:07:08 garloff Exp $

//#define NO_NS

#include "specfun.h"
#include <stdio.h>
#include "prototypes2.h"

NAMESPACE_TBCI

cplx<double> HypergeometricU(const cplx<double>& a, const cplx<double>& b,
                             const cplx<double>& z)
{
  printf ("????\n");
  doublecomplex sin_erg;
  doublecomplex pow_erg;
  doublecomplex temp;
  cplx<double> res;

  temp.r = pi*b.real(); temp.i = pi*b.imag();
  z_sin(&sin_erg, &temp);

  temp.r = -b.real() + 1.0; temp.i = -b.imag();
  pow_zz(&pow_erg, (doublecomplex*)(&z), &temp);
 
  res = pi/cplx<double>(sin_erg.r, sin_erg.i)
          *(HypergeometricM(a,b,z)/(gamma(a-b+1)*gamma(b))
            -cplx<double>(pow_erg.r, pow_erg.i)
             *HypergeometricM(a-b+1,-b+2,z)/(gamma(a)*gamma(-b+2))
           );
  return res;
}


cplx<double> gamma(const cplx<double>& z)
{
  floatcomplex arg = {(float)z.real(), (float)z.imag()};
  floatcomplex erg;
  cgamma_(&erg, &arg);

  return cplx<double>(erg.r, erg.i);
}


cplx<double> HypergeometricM(const cplx<double>& a, const cplx<double>& b,
                             const cplx<double>& z)
{
  doublecomplex erg = {0, 0};
  static const LONG_ int lnchf = 0; //Log-Ausgabe (1)
  static const LONG_ int ip = 10;   //Anzahl der genauen Stellen

  conhyp_(&erg, (doublecomplex*)&a, (doublecomplex*)&b, 
          (doublecomplex*)&z, &lnchf, &ip);

  return cplx<double>(erg.r, erg.i);
}

cplx<double> besselh1(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, m, n; 
  double cyr[1], cyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  m = 1;    // Art H1
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesh_(&zr, &zi, &fnu, &kode, &m, &n, cyr, cyi, &nz, &ierr);
  if (ierr || nz) 
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  static const cplx<double> I(0,1);
  return cplx<double>(cyr[0], cyi[0])* MATH__ exp(pi*fnu*I);
}


cplx<double> besselh2(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, m, n; 
  double cyr[1], cyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  m = 2;    // Art H2
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesh_(&zr, &zi, &fnu, &kode, &m, &n, cyr, cyi, &nz, &ierr);
  if (ierr || nz)
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  static const cplx<double> I(0,1);
  return cplx<double> (cyr[0], cyi[0]) * MATH__ exp(-pi*fnu*I);
}


cplx<double> besselj(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, n; 
  double cyr[1], cyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesj_(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, &ierr);
  if (ierr || nz)
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  //J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU)
  return cplx<double>(cyr[0], cyi[0]) * MATH__ cos(pi*fnu) 
                              - bessely(fnu, z)* MATH__ sin(pi*fnu);
}



cplx<double> besseli(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, n; 
  double cyr[1], cyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesi_(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, &ierr);
  if (ierr || nz)
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  //I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
  return cplx<double> (cyr[0], cyi[0]) 
                        + 2.0/pi*MATH__ sin(pi*fnu)*besselk(fnu,z);
}



cplx<double> besselk(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, n; 
  double cyr[1], cyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesk_(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, &ierr);
  if (ierr || nz)
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  //o.k.
  return cplx<double> (cyr[0], cyi[0]);
}


cplx<double> bessely(double order, const cplx<double>& z)
{
  double zr, zi, fnu;
  LONG_ int   kode, n; 
  double cyr[1], cyi[1];
  double wcyr[1], wcyi[1];
  LONG_ int   nz, ierr;

  zr = z.real(); zi = z.imag();

  kode = 1; // Unskaliert 
  n = 1;    // Laenge der Sequenz
  fnu = MATH__ fabs(order);
  zbesy_(&zr, &zi, &fnu, &kode, &n, cyr, cyi, &nz, wcyr, wcyi, &ierr);
  if (ierr || nz)
    fprintf (stderr, "Error computing Hankel function\n");
  if (order >= 0.0)
    return cplx<double> (cyr[0], cyi[0]);
  //else
  //Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
  return cplx<double> (cyr[0], cyi[0])*MATH__ cos(pi*fnu) 
                               + besselj(fnu,z)*MATH__ sin(pi*fnu);
}

NAMESPACE_END

---