//=============================================================
// daubechies.C
//
// A macro that calculates Daubechies building blocks and 
// 	wavelets .
//
// ALC Sanchez
// IITHEP
//
// Update record
//	23 Nov 1999	ALC Sanchez	*started writing 
//============================================================

#include <iostream.h>
#include <math.h>

const double phi[4] = {0., (1 + sqrt(3))/2, (1 - sqrt(3))/2, 0.};

const double h[4] = {(1 + sqrt(3))/4, (3 + sqrt(3))/4, 
                     (3 - sqrt(3))/4, (1 - sqrt(3))/4};

// Calculate building block Phi(r) value at r.
double Phi(double r) {

  if ( (r<=0) || (r>=3) ) return 0;
  else if (r==1) return phi[1];
  else if (r==2) return phi[2];
  else { 
    cout << "\n Calculating ... ";
    return (h[0]*Phi(2*r) + h[1]*Phi(2*r-1) + h[2]*Phi(2*r-2) +
            h[3]*Phi(2*r-3));
  }
}

// Calculate Daubechies wavelet function Psi(r) at r.
double Psi(double r) {
  return (-h[0]*Phi(2*r-1) +h[1]*Phi(2*r) -h[2]*Phi(2*r+1) + h[3]*Phi(2*r+2));
}

// Main part
int main()
{
  double r, phi, psi;

  cout << "\n This calculates the value of the basic "
       << "\n Daubechies' building block Phi(r) and"
       << "\n wavelet Psi(r) at r."
       << "\n\n Enter value of r: ";
  cin >> r;

  phi = Phi(r);
  psi = Psi(r);

  cout << "\n\n\t Phi( " << r <<" ) = " << phi << " .";
  cout << "\n\t Psi( " << r <<" ) = " << psi << " .\n\n";
  return 0;
}