specfun_stdcplx.cpp

Go to the documentation of this file.

// $Id: specfun_stdcplx.cpp,v 1.1.2.10 2019/05/28 11:13:02 garloff Exp $

//#define NO_NS

#include "tbci/specfun_stdcplx.h"
#include "tbci/constants.h"
#include <stdio.h>
#include "tbci/specfun/prototypes2.h"

NAMESPACE_TBCI

CPLX__ complex<double> HypergeometricU(const CPLX__ complex<double>& a,
                                       const CPLX__ complex<double>& b,
                                       const CPLX__ complex<double>& z)
{
  printf ("????\n");
  doublecomplex sin_erg;
  doublecomplex pow_erg;
  doublecomplex temp;
  CPLX__ complex<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__ complex<double>(sin_erg.r, sin_erg.i)
        *(HypergeometricM(a,b,z)/(gamma(a-b+1.0)*gamma(b))
          -CPLX__ complex<double>(pow_erg.r, pow_erg.i)
           *HypergeometricM(a-b+1.0,-b+2.0,z)/(gamma(a)*gamma(-b+2.0))
         );

  return res;
}


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

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


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

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

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

CPLX__ complex<double> besselh1(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  static const CPLX__ complex<double> I(0,1);
  return CPLX__ complex<double>(cyr[0], cyi[0])* MATH__ exp(pi*fnu*I);
}


CPLX__ complex<double> besselh2(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  static const CPLX__ complex<double> I(0,1);
  return CPLX__ complex<double> (cyr[0], cyi[0]) * MATH__ exp(-pi*fnu*I);
}


CPLX__ complex<double> besselj(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  //J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU)
  return CPLX__ complex<double>(cyr[0], cyi[0])* MATH__ cos(pi*fnu)
                                        - bessely(fnu, z)* MATH__ sin(pi*fnu);
}



CPLX__ complex<double> besseli(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  //I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
  return CPLX__ complex<double> (cyr[0], cyi[0])
                                + 2.0/pi*MATH__ sin(pi*fnu)*besselk(fnu,z);
}



CPLX__ complex<double> besselk(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  //o.k.
  return CPLX__ complex<double> (cyr[0], cyi[0]);
}


CPLX__ complex<double> bessely(double order, const CPLX__ complex<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__ complex<double> (cyr[0], cyi[0]);
  //else
  //Y(-FNU,Z) = Y(FNU,Z)*COS(PI*FNU) + J(FNU,Z)*SIN(PI*FNU)
  return CPLX__ complex<double> (cyr[0], cyi[0])*MATH__ cos(pi*fnu)
                                         + besselj(fnu,z)*MATH__ sin(pi*fnu);
}

NAMESPACE_END