/* dht/test_dht.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_dht.h>


/* Test exact small transform.
 */
int
test_dht_exact(void)
{
  int stat = 0;
  double f_in[3] = { 1.0, 2.0, 3.0 };
  double f_out[3];
  gsl_dht * t = gsl_dht_new(3, 1.0, 1.0);
  gsl_dht_apply(t, f_in, f_out);

  /* Check values. */
  if(fabs( f_out[0]-( 0.375254649407520))/0.375254649407520 > 1.0e-14) stat++;
  if(fabs( f_out[1]-(-0.133507872695560))/0.133507872695560 > 1.0e-14) stat++;
  if(fabs( f_out[2]-( 0.044679925143840))/0.044679925143840 > 1.0e-14) stat++;


  /* Check inverse.
   * We have to adjust the normalization
   * so we can use the same precalculated transform.
   */
  gsl_dht_apply(t, f_out, f_in);
  f_in[0] *= 13.323691936314223*13.323691936314223;  /* jzero[1,4]^2 */
  f_in[1] *= 13.323691936314223*13.323691936314223;
  f_in[2] *= 13.323691936314223*13.323691936314223;

  /* The loss of precision on the inverse
   * is a little surprising. However, this
   * thing is quite tricky since the band-limited
   * function represented by the samples {1,2,3}
   * need not be very nice. Like in any spectral
   * application, you really have to have some
   * a-priori knowledge of the underlying function.
   */
  if(fabs( f_in[0]-1.0)/1.0 > 2.0e-05) stat++;
  if(fabs( f_in[1]-2.0)/2.0 > 2.0e-05) stat++;
  if(fabs( f_in[2]-3.0)/3.0 > 2.0e-05) stat++;

  gsl_dht_free(t);

  return stat;
}



/* Test the transform
 * Integrate[x J_0(a x) / (x^2 + 1), {x,0,Inf}] = K_0(a)
 */
int
test_dht_simple(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 0.0, 100.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = 1.0/(1.0+x*x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* This is a difficult transform to calculate this way,
   * since it does not satisfy the boundary condition and
   * it dies quite slowly. So it is not meaningful to
   * compare this to high accuracy. We only check
   * that it seems to be working.
   */
  if(fabs( f_out[0]-4.00)/4.00 > 0.02) stat++;
  if(fabs( f_out[5]-1.84)/1.84 > 0.02) stat++;
  if(fabs(f_out[10]-1.27)/1.27 > 0.02) stat++;
  if(fabs(f_out[35]-0.352)/0.352 > 0.02) stat++;
  if(fabs(f_out[100]-0.0237)/0.0237 > 0.02) stat++;

  gsl_dht_free(t);

  return stat;
}


/* Test the transform
 * Integrate[ x exp(-x) J_1(a x), {x,0,Inf}] = a F(3/2, 2; 2; -a^2)
 */
int
test_dht_exp1(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 1.0, 20.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = exp(-x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* Spot check.
   * Note that the systematic errors in the calculation
   * are quite large, so it is meaningless to compare
   * to a high accuracy.
   */
  if(fabs( f_out[0]-0.181)/0.181 > 0.02) stat++;
  if(fabs( f_out[5]-0.357)/0.357 > 0.02) stat++;
  if(fabs(f_out[10]-0.211)/0.211 > 0.02) stat++;
  if(fabs(f_out[35]-0.0289)/0.0289 > 0.02) stat++;
  if(fabs(f_out[100]-0.00221)/0.00211 > 0.02) stat++;

  gsl_dht_free(t);

  return stat;
}


/* Test the transform
 * Integrate[ x^2 (1-x^2) J_1(a x), {x,0,1}] = 2/a^2 J_3(a)
 */
int
test_dht_poly1(void)
{
  int stat = 0;
  int n;
  double f_in[128];
  double f_out[128];
  gsl_dht * t = gsl_dht_new(128, 1.0, 1.0);

  for(n=0; n<128; n++) {
    const double x = gsl_dht_x_sample(t, n);
    f_in[n] = x * (1.0 - x*x);
  }

  gsl_dht_apply(t, f_in, f_out);

  /* Spot check. This function satisfies the boundary condition,
   * so the accuracy should be ok.
   */
  if(fabs( f_out[0]-0.057274214)/0.057274214    > 1.0e-07) stat++;
  if(fabs( f_out[5]-(-0.000190850))/0.000190850 > 1.0e-05) stat++;
  if(fabs(f_out[10]-0.000024342)/0.000024342    > 1.0e-04) stat++;
  if(fabs(f_out[35]-(-4.04e-07))/4.04e-07       > 1.0e-03) stat++;
  if(fabs(f_out[100]-1.0e-08)/1.0e-08           > 0.25)    stat++;

  gsl_dht_free(t);

  return stat;
}


int main()
{
  gsl_ieee_env_setup ();

  gsl_test( test_dht_exact(),   "Small Exact DHT");
  gsl_test( test_dht_simple(),  "Simple  DHT");
  gsl_test( test_dht_exp1(),    "Exp  J1 DHT");
  gsl_test( test_dht_poly1(),   "Poly J1 DHT");

  exit (gsl_test_summary());
}