use strict;
use warnings;


# pi
sub pi () { 4 * atan2(1, 1) }
#      ^^
# For constant folding


# make_vector_function
# Returns a sub that will return a vector for any
# value of t that is passed.
# The function is compiled from three component functions.

sub make_vector_function {
   my $x_of_t = shift;
   my $y_of_t = shift;
   my $z_of_t = shift;

   return sub {
      my $t = shift;

      return( Math::MatrixReal->new_from_cols(
        [
          [ $x_of_t->($t, @_), $y_of_t->($t, @_), $z_of_t->($t, @_) ]
        ]
      ) );
   };
}

# output_image
# Takes a GD::Image object and a filename as argument.
# Returns 1 on success.
# 
# If the filename is equal to '', the png image is printed
# to STDOUT, otherwise, the file is opened and the data
# written into it.

sub output_image {
   # Image object
   my $image = shift;
   # file name
   my $file  = shift;

   # Print to file if the filename's existant
   if ($file ne '') {

      # open for writing
      open my $fh, ">", $file
        or die "Could not open file ($file): $!";

      binmode $fh; # Make Winlame happy

      # print data
      print $fh $image->png;

      # close and check for success
      close $fh
        or die "Could not complete write process to file ($file): $!";

   } else {
      binmode STDOUT; # Wooohooo, Windows!
      # output the data to STDOUT
      print $image->png;
   }

   return 1;
}



# Aux function l_g (logical to graphical)
# Takes an x/y pair of logical coordinates as
# argument and returns the corresponding graphical
# coordinates.

sub l_g {
   my $x = shift;
   my $y = shift;

   # A logical unit is a graphical one displaced by u_?
   # and multiplied with the appropriate scaling factor.

   # Move to the right by $u_x (x-Ursprung) pixels.
   $x = $u_x + $x * $x_scale;
   # Move down by $u_y (y-Ursprung) pixels.
   $y = $u_y - $y * $y_scale;

   return $x, $y;
}


# Aux function g_l (graphical to logical)
# Takes an x/y pair of graphical coordinates as
# argument and returns the corresponding
# (later: nearest) logical coordinates.

sub g_l {
   my $x = shift;
   my $y = shift;

   # A graphical unit is a logical one displaced by u_?
   # and divided by the appropriate scaling factor.

   # Move to the left by $u_x (x-Ursprung) pixels.
   $x = ( $x - $u_x ) / $x_scale;
   # Move up by $u_y (y-Ursprung) pixels.
   $y = ( $y - $u_y ) / $y_scale;

   return $x, $y;
}


# Function plotting routines
# plot_step and reset_plot

# Closures ahead!
# private vars used for caching
{
   # Initial graphical coords
   my $init_l_x;
   my $init_l_y;

   # Previous graphical coords
   my $prev_g_x;
   my $prev_g_y;

   # Reset cache
   sub reset_plot {
      $prev_g_x = $prev_g_y = undef;
      $init_l_x = $init_l_y = undef;
   }

   # Plot one step of the function
   sub plot_step {
      # logical coords and image
      my($l_x, $l_y, $image, $color) = @_;

      # If no cache, then set cache
      if ( not defined $prev_g_x or not defined $prev_g_y ) {

         ( $prev_g_x, $prev_g_y ) = l_g( $l_x, $l_y );
         ( $init_l_x, $init_l_y ) = ( $l_x, $l_y );

      } else {

         # plot

         my ($g_x, $g_y) = l_g( $l_x, $l_y );

         $image->line( $prev_g_x, $prev_g_y,
                       $g_x, $g_y,
                       $color
                     );

         # set cache to new values
         ( $prev_g_x, $prev_g_y ) = ( $g_x, $g_y );
      }

      return 1;
   }

   sub plot_finish {
      my $image = shift;
      my $color = shift;
      return if not defined $init_l_x or not defined $init_l_y;

      plot_step($init_l_x, $init_l_y, $image, $color);

      reset_plot();
   }
}



# e (precision???)
sub e { 2.71828182845905 }

# alias for log
sub ln { log(@_) }

# tan
sub tan { sin($_[0]) / cos($_[0]) }

# sec - secant
sub sec { 1 / cos($_[0]) }

# csc - cosecant
sub csc { 1 / sin($_[0]) }

# cot - cotangent
sub cot { cos($_[0]) / sin($_[0]) }

# acos
# Computes the arc cosine acos(z) = -i log(z + sqrt(z*z-1)).
# Above formula is for complex numbers.
# 
sub acos { atan2( sqrt( 1 - $_[0] * $_[0] ), $_[0] ) }

# asin
#
# Computes the arc sine asin(z) = -i log(iz + sqrt(1-z*z)).
# Above formula is for complex numbers.
# 
sub asin { atan2($_[0], sqrt( 1 - $_[0] * $_[0] ) ) }

# atan
#
# Computes the arc tangent atan(z) = i/2 log((i+z) / (i-z)).
# Above formula is for complex numbers.
sub atan { atan2($_[0], 1) }

# asec - arc secant
sub asec { acos(1 / $_[0]) }

# acsc - arc cosecant
sub acsc { asin(1 / $_[0]) }

# acot - arc cotangent
sub acot { atan(1 / $_[0]) }

# cosh - hyperbolic cosine
sub cosh {
   my $ex = exp($_[0]);
   return $ex ? ($ex + 1 / $ex) / 2 : die;
}

# sinh - hyperbolic sine
sub sinh {
   return 0 if $_[0] == 0;
   my $ex = exp($_[0]);
   return $ex ? ($ex - 1 / $ex) / 2 : die;
}

# tanh - hyperbolic tangent
sub tanh { sinh($_[0]) / cosh($_[0]) }

# sech - hyperbolic secant
sub sech { 1 / cosh($_[0]) }

# csch - hyperbolic cosecant
sub csch { 1 / sinh($_[0]) }

# coth - hyperbolic cotangent
sub coth { cosh($_[0]) / sinh($_[0]) }

# acosh
#
# Computes the arc hyperbolic cosine acosh(z) = log(z + sqrt(z*z-1)).
#
sub acosh { log( $_[0] + sqrt( $_[0] * $_[0] - 1 ) ) }

# asinh
#
# Computes the arc hyperbolic sine asinh(z) = log(z + sqrt(z*z+1))
#
sub asinh { log( $_[0] + sqrt( $_[0] * $_[0] + 1 ) ) }

# atanh
#
# Computes the arc hyperbolic tangent atanh(z) = 1/2 log((1+z) / (1-z)).
#
sub atanh { log( (1+$_[0]) / (1-$_[0]) ) / 2 }

# asech
#
# Computes the hyperbolic arc secant asech(z) = acosh(1 / z).
#
sub asech { acosh( 1 / $_[0] ) }

# acsch
#
# Computes the hyperbolic arc cosecant acsch(z) = asinh(1 / z).
#
sub acsch { asinh( 1 / $_[0] ) }

# acoth
#
# Computes the arc hyperbolic cotangent acoth(z) = 1/2 log((1+z) / (z-1)).
#
sub acoth { log( (1+$_[0]) / ($_[0]-1) ) / 2 }







1;