// one-demo.C : example for using one dimensional grafix displays
// testing of one-dim wavelets : function wt1 of numerical recipes
// (c) wolf 11/96
// 

#define NRANSI
extern "C" {
#include "nr.h" 
#include "nrutil.h"
}

#include "toolbox.h"

float *u,*v,*w,*cv; // the used vectors 

// two floats which values are interactively modified by scrollbars
// the floats must be initialized !! 
float frac = 20; // number of coeffs to use in back transformation 
float fwid = 33; // width of step function (initial value = 33)
int k; // order of wavelets

void wavelet(int n) {  
  int ncen = 333; // centre
  // printf("wavelet %d %d\n",n,k);
  for (int i = 1; i <= n; i++) // construct function to approximate
    // use a span from [-0.5..1.0] to fit into coord system [-1.0..1.0] 
    w[i] = -0.5 + 1.5*(i > ncen-fwid && i < ncen+fwid ?
	       ((float)(i-ncen+fwid)*(float)(ncen+fwid-i))/(fwid*fwid) : 0.0);
  
  // compute coeffs of w -> v, u = |v|
  pwtset(k); // k = order of wavelets = 4,12,20
  for (int i = 1; i <= n; i++) v[i] = w[i];
  wt1(v, n, 1, pwt); // forward trafo
  for (int i=1; i<= n; i++) u[i] = fabs(v[i]);
  int nfrac = irint(frac);   
  float thresh = f_select(n-nfrac,n,u); // threshold
  int i;
  for (i = 1; i <= n; i++) {
    float av = fabs(v[i]);
    if (av <= thresh) cv[i] = 0.0; 
    else cv[i] = v[i]; 
  }
  wt1(cv,n,-1,pwt); // back transform
  float tmp, dmax = 0.0; // get max difference
  for (i = 1; i <= n; i++) if ((tmp=fabs(cv[i]-w[i])) > dmax) dmax=tmp;

  // info_text is a string which is displayed at bottom of the window
  sprintf(info_text,"wavelet of order %d nused %d thresh %f dmax %f",
	  k,nfrac,thresh,dmax);
}

// these solvers that are interactively selectable 
void wave4(int n) { k = 4; wavelet(n); }
void wave12(int n) { k = 12; wavelet(n); }
void wave20(int n) { k = 20; wavelet(n); }

void wave_org(int n) { // compute wavelets itself 
  int i,nfrac = irint(frac); 
  for (i = 1; i <= n; i++) cv[i] = (i==nfrac) ? 1.0 : 0.0;
  wt1(cv, n, -1, pwt); // back transform
  sprintf(info_text,"display wavelet of order %d index %d",k,nfrac);
}

#define NMAX 1024

int main(void) {
    
  u = vector(1,NMAX); // used for determine threshold
  v = vector(1,NMAX); // transformed coeffs
  w = vector(1,NMAX); // original series
  cv = vector(1,NMAX);
  
  // list of all displayed vectors
  disp_vector *dv[] = { new disp_vector(w, "w  : original vector"), 
			new disp_vector(cv,"cv : backtransform  "), 
			0 }; // never forget this terminating NULL !!!!

  // list of interactive controlling variables and functions
  parameter *scr[] = { 
    new parameter("used number of wavelet coefficients", &frac, 0,100), 
    new parameter("width of the peak function to approx", &fwid, 10, 100),
    0 }; // never forget this terminating NULL !!!!

  one_solver *sv[] = { new one_solver(wave4,"order 4"),
		       new one_solver(wave12,"order 12"),
		       new one_solver(wave20,"order 20"),
		       new one_solver(wave_org,"wavelet"),
		       0 }; // !!!

  one_show("one-demo",NMAX,dv,sv,scr);   // this call does the job !

  free_vector(cv,1,NMAX); // what for is this ?
  free_vector(w,1,NMAX);
  free_vector(v,1,NMAX);
  free_vector(u,1,NMAX);

};