//=============================================================
// 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>

const Double_t phi[4] = {0., (1 + TMath::Sqrt(3))/2, 
                         (1 - TMath::Sqrt(3))/2, 0.};

const Double_t h[4] = {(1 + TMath::Sqrt(3))/4, (3 + TMath::Sqrt(3))/4, 
                     (3 - TMath::Sqrt(3))/4, (1 - TMath::Sqrt(3))/4};

// Calculate building block Phi(r) value at r.
Double_t Phi(Double_t 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_t Psi(Double_t 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 daub_ROOT()
{
  Double_t 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;
}