#!/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 = 4.0;
# maximum order of Fourier series to compute
$max_order = 10;
# this is 1 when the program starts, but can be set to 0 later 
#   by the "func_to_fit()" routine
$first_flag = 1;



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

  $a_coeff[$order] = calc_a_coeff($order);
  $b_coeff[$order] = calc_b_coeff($order);

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

}



$dt = 0.01;
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;

    $a_val[$n] = $a_coeff[$n]*cos($omega_n*$t);
    $b_val[$n] = $b_coeff[$n]*sin($omega_n*$t);

    if ($n == 0) {
      $sum_val[$n] = ($a_val[$n] + $b_val[$n]);
    }
    else {
      $sum_val[$n] = $sum_val[$n - 1] + ($a_val[$n] + $b_val[$n]);
    }
  }


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

  printf "#QQ t %12.5f  f %12.5f  sum %12.5f \n", $t, $f, $sum_val[$max_order];

}





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: 
#
#                    3 + sin(2 pi t/ P )
#
# 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];


  if (0 == 1) {
          # this is a sine curve with offset from zero
			 if ($first_flag == 1) {
				printf "#  function is 3 + sin(2 pi t / %.1f) \n", $P;
				$first_flag = 0;
			 }
          $val = 3.0 + sin(2*$mypi*$t / $P);
  }
  
  if (0 == 1) {
          # this is a square wave
			 if ($first_flag == 1) {
				printf "#  function is a square wave, amplitude 1 \n";
				$first_flag = 0;
			 }
			 if ($t < 1.0) {
				$val = 1.0;
			 }
			 else {
				$val = 0.0;
			 }
  }
  
  if (0 == 1) {
          # this is a sawtooth wave
			 if ($first_flag == 1) {
				printf "#  function is a sawtooth wave, amplitude 1 \n";
				$first_flag = 0;
			 }
          $val = 1.0*($t/$P);
  }

  if (1 == 1) {
          # this is a "tooth" wave, which looks like
          #      _/\_  
          # with a period of 4 and amplitude of 1 
			 if ($first_flag == 1) {
				printf "#  function is a tooth wave, amplitude 1 \n";
				$first_flag = 0;
			 }
          if ($t < ($P/4.0)) {
            $val = 0.0;
          }
          elsif ($t < $P/2.0) {
            $val = ($t - 1.0);
          }
          elsif ($t < (3.0*$P/4.0)) {
            $val = (3.0 - $t);
          }
          else {
            $val = 0.0;
          }
  }
  

  


  if ($debug > 0) {
    printf " func_to_fit:  t %12.5f  val %12.5f \n",
             $t, $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];

  if ($n == 0) {
    return(0.0);
  }

  $n_piece = 3000;
  $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);

}


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

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

  $n = $_[0];

  $n_piece = 3000;
  $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 = cos($omega_n*$x);
    $this_piece = $this_1*$this_2*$dx;
    $sum += $this_piece;
    if ($debug > 1) {
      printf " calc_a_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;
    }
  }

  if ($n > 0) {
    $sum *= (2.0/$P);
  }
  else {
    $sum *= (1.0/$P);
  }

  return($sum);

}

