#!/usr/bin/perl
#
# Test my ability to compute fourier coefficients 
#    for a simple function.
#
# MWR 2/2/2019

$debug = 0;

$mypi = 3.1415926;
# period of the function to fit
$P = 10.0;
# maximum order of Fourier series to compute
$max_order = 5;


for ($order = 0; $order <= $max_order; $order++) {

  $b_coeff[$order] = calc_b_coeff($order);
  if ($debug > -1) {
    printf "#  order %5d  b_coeff %12.5f \n", $order, $b_coeff[$order];
  }

}

exit 0;


$dt = 0.05;
for ($t = 0; $t <= $P; $t += $dt) {


  $f = func_to_fit($t);

  # now compute the Fourier sum to approximate this value
  for ($n = 0; $n <= $max_order; $n++) {
    $sum_val[$n] = 0;
  }
  for ($n = 0; $n <= $max_order; $n++) {
    $omega_n = (2.0*$mypi*$n)/$P;
    $val[$n] = $b_coeff[$n]*sin($omega_n*$t);
    if ($n == 0) {
      $sum_val[$n] = $val[$n];
    }
    else {
      $sum_val[$n] = $sum_val[$n - 1] + $val[$n];
    }
  }


  printf " t %12.5f  f %12.5f ", $t, $f;
  for ($n = 0; $n <= $max_order; $n++) {
     printf " %2d v %9.4f  s %9.4f ", $n, $val[$n], $sum_val[$n];
  }
  printf "\n";

}





exit 0;


###########################################################
# PROCEDURE: func_to_fit
#
# DESCRIPTION: Given some value of time t, return the value
#              of this function at that time.
#
#              This is a simple function: linear rise from
#              0 to 1 at a time P (the period of the function).
#
# USAGE:       func_to_fit   t
#
#              where
#                       t         is the time, in range 0 to P
#
# RETURNS:     the value of the function 
#
sub func_to_fit {

  my($t, $frac, $max, $min, $val);
  $t = $_[0];

  $min = -1.0;
  $max = 1.0;

  # determine what fraction of a period is represented by time t
  $frac = ( ($t/$P) - int($t/$P) );

  # the value of the function at that fraction of a period
  $val = $min + ($max - $min)*$frac;
  if ($debug > 0) {
    printf " func_to_fit:  t %12.5f  frac %12.5f  val %12.5f \n",
             $t, $frac, $val;
  }
  
  return($val);

}


###########################################################
# PROCEDURE: calc_b_coeff
#
# DESCRIPTION: Compute the value of the N'th coefficient
#              of the sine series for our function.
#
# USAGE:       calc_b_coeff   n
#
#              where
#                       n         is the order of coeff to compute
#
# RETURNS:     the value of the coefficient 
#
sub calc_b_coeff {

  my($n, $x);
  my($n_piece, $dx, $omega_n, $sum, $this_piece);

  $n = $_[0];

  $n_piece = 100;
  $dx = $P/$n_piece;
  $omega_n = ((2.0*$mypi*$n)/$P);
  if ($debug > 0) {
    printf "  dx %9.4f  omega_n %9.4f \n", $dx, $omega_n;
  }
  
  $sum = 0;

  # do a numerical integration over one full period
  for ($x = 0; $x < $P; $x += $dx) {
    
    $this_1 = func_to_fit($x);
    $this_2 = sin($omega_n*$x);
    $this_piece = $this_1*$this_2*$dx;
    $sum += $this_piece;
    if ($debug > 1) {
      printf " calc_b_coeff: x %9.4f  func %9.4f  sin %9.4f  dx %9.4f  piece %9.4f  sum %9.4f \n",
           $x, $this_1, $this_2, $dx, $piece, $sum;
    }
  }

  $sum *= (2.0/$P);

  return($sum);

}

