#!/usr/bin/perl
use strict;
use warnings;

use Math::MatrixReal;

use GD;

# Globals.
# Yes, I know, globals considered harmful.
use vars qw(
  $VERSION
  $output_file $image_width $image_height
  $u_x $u_y $x_scale $y_scale
  $foreground_color $background_color
  $x_axis_color $y_axis_color $z_axis_color
);

$VERSION = 1.010;

# I've put some functions into another file.
# Make sure it is in the same directory you invoke the program
# from.
BEGIN {
   require 'functions.res';
}

# Output file
$output_file = 'out.png';

# Image specs:
$image_width  = 1024;
$image_height = 768;

# Where is the 0/0?
$u_x = $image_width  / 2;
$u_y = $image_height / 2;

# Scaling
$x_scale = 10;
$y_scale = 10;

# Colors
my @foreground_color = (0,   0,   0);
my @background_color = (255, 255, 255);
my @x_axis_color     = (255, 0,   0);
my @y_axis_color     = (0,   255, 0);
my @z_axis_color     = (0,   0,   255);

# Prepare new image
my $image = new GD::Image($image_width, $image_height);

# Resolve colors
$foreground_color = $image->colorResolve( @foreground_color );
$background_color = $image->colorResolve( @background_color );
$x_axis_color     = $image->colorResolve( @x_axis_color );
$y_axis_color     = $image->colorResolve( @y_axis_color );
$z_axis_color     = $image->colorResolve( @z_axis_color );

# Background
$image->fill(0, 0, $background_color);

# Parameter range and increment
my $t_lower_range_boundary = -50;
my $t_upper_range_boundary = 50;
my $t_increment            = 0.1;

my $u_lower_range_boundary = 0;
my $u_upper_range_boundary = 0;
my $u_increment            = .1;

# Separate plots for all t's: f(u)
my $separate_plots = 0;

# color range for coloring ||s(t,u)||

my $plot_ranged_color = 1;
my $max_distance      = 15;

my $color_steps       = 150;
my $color_start       = 0;
my $color_end         = 255;

my $color_step_diff   = abs($color_start - $color_end) / $color_steps;
$color_step_diff      *= -1 if $color_start > $color_end;

my $color_form        = '$color, $color, $color';

foreach (1..$color_steps) {
   my $color = $color_start + $color_step_diff * $_;
   $image->colorResolve( eval $color_form );
   die $@ if $@;
}

# Component function for the x coordinate of s
my $x_func = sub {
   my $t = shift;

   return sin($t)*5;
};

# Component function for the y coordinate of s
my $y_func = sub {
   my $t = shift;

   return ($t);
};

# Component function for the z coordinate of s
my $z_func = sub {
   my $t = shift;

   return cos($t)*5;
};

# s_function will become a function that returns a M::MatrixReal
# vector for any t that is passed.
my $s_function = make_vector_function( $x_func, $y_func, $z_func );

# Now we define the plane we want to project onto.

# Define one point on the plane.
my $plane_point = Math::MatrixReal->new_from_cols(
                      [
                        [0, 0, 0 ]
                      ]
                    );

# Define one directional vector.
my $plane_d1 = Math::MatrixReal->new_from_cols(
                      [
                        [ .4, 1, 0 ]
                      ]
                    );

# Define another one.
my $plane_d2 = Math::MatrixReal->new_from_cols(
                      [
                        [ .4, 0, 1 ]
                      ]
                    );

# Get the normal of the plane.
my $normal_vector = $plane_d1->vector_product( $plane_d2 );


# Generate a matrix needed to solve the linear equation system
# (d, e be the directional vectors, n the normal, p the plane point)
# (let i be the solution)
# 
# x(t) + n1*i1 = d1*i2 + e1*i3 + p1
# y(t) + n2*i1 = d2*i2 + e2*i3 + p2
# z(t) + n3*i1 = d3*i2 + e3*i3 + p3
# 
# This is the intersection of the plane and the corresponding orthogonal
# line through s(t). Another form:
# 
# n1*i1 + d1*i2 + e1*i3 = p1 + x(t)
# n2*i1 + d2*i2 + e2*i3 = p2 + y(t)
# n3*i1 + d3*i2 + e3*i3 = p3 + z(t)
# 
# linear_eq_system will be the matrix of n,d,e.
# result_vector will be p+s(t)

my $linear_eq_system = Math::MatrixReal->new_from_cols(
                      [
                        $normal_vector,
                        $plane_d1,
                        $plane_d2,
                      ]
);

# Do the first step of solving the les by decomposing it into
# an LR matrix.
my $LR_matrix = $linear_eq_system->decompose_LR();

draw_coordinate_axis(
  $image,
  $LR_matrix,
  $plane_point,
  [ $x_axis_color, $y_axis_color, $z_axis_color ],
  100   # Length in units
);

# Reset function plotting routine
reset_plot();

# For any point given by $s_function that is in the range
for (my $t = $t_lower_range_boundary; $t <= $t_upper_range_boundary; $t += $t_increment ) {

   for (my $u = $u_lower_range_boundary; $u <= $u_upper_range_boundary; $u += $u_increment ) {
      #($u,$t)=($t,$u);
      # Generate result_vector
      my $result_vector = Math::MatrixReal->new_from_cols(
                            [
                              $plane_point + $s_function->($t, $u)
                            ]
      );
      
      # Solve the l_e_s.
      my ($dimension, $p_vector, undef) = $LR_matrix->solve_LR($result_vector);
      
      # Sanity check
      die "Invalid result of linear equation system: dimension $dimension."
        if $dimension;
      
      # The second and third component of the resulting vector
      # are the linear factors that identify any 2D position on the
      # plane. Hence, if the directional vectors plane_d1 and plane_d2
      # are unit vectors (|plane_d1|=|plane_d2|=1), we may assume
      # those two components to be the screen coordinates of our
      # projection.
      
      # print join( ', ', $p_vector->element(2,1), $p_vector->element(3,1) ), "\n";

      my $color;
      if ( $plot_ranged_color ) {
         my $distance = abs( $p_vector->element(1,1) * $normal_vector->length() );
         if ( $distance > $max_distance ) {
            $color = $color_end;
         } else {
            $color = $color_start + $distance / $max_distance * ($color_end - $color_start);
         }
         $color = $image->colorClosest( eval $color_form );
         die $@ if $@;
      } else {
         $color = $foreground_color;
      }

      plot_step( $p_vector->element(2,1), $p_vector->element(3,1), $image, $color );
   }

   # Finish function plotting routine
   plot_finish( $image, $foreground_color ) if $separate_plots;
}

# Write image to file
output_image( $image, $output_file );





sub draw_coordinate_axis {
   my $image       = shift;
   my $LR_matrix   = shift;
   my $plane_point = shift;
   my $colors      = shift;
   my $length      = shift || 10;

   my @axis_names = qw(x y z);

   $colors = [] if ref $colors ne 'ARRAY';
   while (@{$colors} < 3) {
      push @$colors, $foreground_color;
   }

   # Define unit vectors
   my @unit_vectors = (
     Math::MatrixReal->new_from_cols( [ [ 1, 0, 0 ] ] ),
     Math::MatrixReal->new_from_cols( [ [ 0, 1, 0 ] ] ),
     Math::MatrixReal->new_from_cols( [ [ 0, 0, 1 ] ] ),
   );

   foreach my $unit (@unit_vectors) {
      # Generate result_vector
      my $result_vector = Math::MatrixReal->new_from_cols(
                            [
                              $plane_point + $unit
                            ]
      );

      # Solve the l_e_s.
      my ($dimension, $p_vector, undef) = $LR_matrix->solve_LR($result_vector);
      
      # Sanity check
      die "Invalid result of linear equation system: dimension $dimension."
        if $dimension;

      my $axis_color = shift @$colors;

      # Draw axis
      $image->line(
         l_g( $p_vector->element(2,1) * (-1 * $length / 2),
              $p_vector->element(3,1) * (-1 * $length / 2)
         ),
         l_g( $p_vector->element(2,1) * (     $length / 2),
              $p_vector->element(3,1) * (     $length / 2),
         ),
         $axis_color
      );

      my $desc_x = $p_vector->element(2,1) * $length / 2;
      my $desc_y = $p_vector->element(3,1) * $length / 2;

      my ($min_x, $min_y) = g_l(0,0);
      my ($max_x, $max_y) = g_l( $image->getBounds() );

      my ($diff_high) = (g_l(15,0))[0] - (g_l(0,0))[0];
      my ($diff_low)  = (g_l(7, 0))[0] - (g_l(0,0))[0];

      $desc_x = $max_x - $diff_high if $desc_x > $max_x - $diff_high;
      $desc_x = $min_x + $diff_low  if $desc_x < $min_x + $diff_low;
      $desc_y = $max_y - $diff_low  if $desc_y > $max_y - $diff_low;
      $desc_y = $min_y + $diff_high if $desc_y < $min_y + $diff_high;

      # axis descriptors
      $image->string(GD::Font->MediumBold, l_g($desc_x, $desc_y), shift @axis_names, $axis_color);
   }
}