#!/usr/bin/perl
#
# High-order Mandelbrot program.
# Written by Christopher Thomas.
#
# This program draws mandelbrot and julia set analogues using relations of
# the form "z <- z^n + c", for arbitrary n.
#
# Usage: mandel-high <mode> <exp> <imax> <region> <resolution> <output>
#
#  <mode> is "julia" for julia sets, "mandel" for mandelbrot sets.
#  <exp> is the exponent n to use. Non-zero real number.
#  <imax> is the maximum number of iterations to perform per point.
#  <region> is six coordinate pairs of the form "a,b" (for a + bi).
#    The first pair defines the upper-left corner of the scanning region.
#    The second pair defines the lower-right corner of the region.
#    The third pair defines the fixed parameter.
#    For Julia sets, the fixed parameter is C, and the scanned parameter
#    is the starting value of Z.
#    For Mandelbrot sets, the fixed parameter is the starting value of Z
#    (set to 0 for the classic set), and the scanned parameter is C.
#  <resolution> is the resolution of the image to render, of the form
#    "axb".
#  <output> is the filename to write the image to. This is a P3 PPM
#    file (human-readable colour bitmap).
#


#
# Libraries
#

use strict;



#
# Constants
#

# Picture palette info.

my ($palsize);
my (@palette);

if(0)
{
# Light/dark colour banded palette.
# NOTE: This looks ugly, probably because the dark colours look muddy.
$palsize = 16;
@palette =
  ( "  255   0   0", "    0 112 112", "  255 128   0", "    0   0 128",
    "  224 224   0", "   64   0  96", "    0 255   0", "   96   0  64",
    "    0 224 224", "  128   0   0", "    0   0 255", "  128  64   0",
    "  128   0 192", "  112 112   0", "  192   0 128", "    0 128   0" );
}
else
{
# 8-colour rainbow palette.
$palsize = 8;
@palette =
  ( "  255   0   0", "  255 128   0",
    "  224 224   0", "    0 255   0",
    "    0 224 224", "    0   0 255",
    "  128   0 192", "  192   0 128" );
}


#
# Mandelbrot Functions
#

# Generates an iteration count for a given point in the Mandelbrot/Julia
# super-set.
# Arg 1 is the exponent to use.
# Arg 2 is the maximum number of iterations.
# Arg 3 is the starting real Z value.
# Arg 4 is the starting imaginary Z value.
# Arg 5 is the real C value.
# Arg 6 is the imaginary C value.

sub do_mandel
{
  my ($exponent, $imax);
  my ($zr, $zi, $cr, $ci);
  my ($rad, $radmax, $angle);
  my ($icount);

  # Get args.
  $exponent = $_[0];
  $imax = $_[1];
  $zr = $_[2];
  $zi = $_[3];
  $cr = $_[4];
  $ci = $_[5];

  # Proceed if args are valid.
  if ((defined $exponent) && (defined $imax)
    && (defined $zr) && (defined $zi) && (defined $cr) && (defined $ci))
  {
    # Compute limiting radius.
    # Anything outside this grows more than C can correct for, for
    # C within the same limit.
    $radmax = exp(log(2) / ($exponent - 1.0));

    # Initialize iteration count.
    $icount = 0;

    # Calculate radius, for while check.
    $rad = $zr * $zr + $zi * $zi;
    $rad = sqrt($rad);

    # Move this point.
    while (($rad <= $radmax) && ($icount < $imax))
    {
      # Handle Z <- Z^n.
      # We have real N, so do this in polar coordinates.

      # Calculate the angle. We already have radius.
      if (($zr < 1.0e-20) && ($zr > 1.0e-20))
      {
        # Handle the limiting cases of tangent.
        if ($zi < 0.0)
        {
          $angle = -3.1415926535 * 0.5;
        }
        else
        {
          $angle = 3.1415926535 * 0.5;
        }
      }
      else
      {
        # This gives the full -pi..pi, supposedly.
        $angle = atan2($zi, $zr);
      }

      # Do the exponentiation.
      $angle *= $exponent;
      if (1.0e-20 < $rad)
      {
        $rad = exp($exponent * log($rad));
      }


      # Handle Z += C.
      # Convert back to cartesian coordinates for this.

      # Convert to cartesian.
      $zr = $rad * cos($angle);
      $zi = $rad * sin($angle);

      # Add.
      $zr += $cr;
      $zi += $ci;


      # Increment the iteration count.
      $icount++;

      # Calculate radius, as we need it for the while check.
      $rad = $zr * $zr + $zi * $zi;
      $rad = sqrt($rad);
    }
  }

  # Return the iteration count.
  return $icount;
}


# Generates a picture of the requested set.
# Arg 1 points to the argument hash.
# Returns a pointer to the plot data.

sub gen_pic
{
  my ($args_p);
  my ($pic_p);
  my ($h, $v, $width, $height);
  my ($x, $xstart, $xinc, $y, $ystart, $yinc, $cx, $cy);
  my ($exponent, $imax);
  my ($is_mandel);

  # Default to failure.
  $pic_p = undef;

  # Get args.
  $args_p = $_[0];
  if (defined $args_p)
  {
    # Process arguments.
    $width = $$args_p{width};
    $height = $$args_p{height};
    $exponent = $$args_p{exponent};
    $imax = $$args_p{imax};
    $xstart = $$args_p{scan1r};
    $ystart = $$args_p{scan1i};
    $xinc = $$args_p{scan2r} - $xstart;
    $xinc /= ($width - 1);
    $yinc = $$args_p{scan2i} - $ystart;
    $yinc /= ($height - 1);
    $cx = $$args_p{fixedr};
    $cy = $$args_p{fixedi};
    $is_mandel = ("mandel" eq $$args_p{mode} ? 1 : 0);

    # Produce the image.
    $pic_p = [];
    $y = $ystart;
    for ($v = 0; $v < $height; $v++)
    {
      $x = $xstart;
      for ($h = 0; $h < $width; $h++)
      {
        if ($is_mandel)
        {
          # Mandelbrot set.
          # Fixed Z0, scan C.
          $$pic_p[$v][$h] =
            do_mandel($exponent, $imax, $cx, $cy, $x, $y);
        }
        else
        {
          # Julia set.
          # Scan Z0, fix C.
          $$pic_p[$v][$h] =
            do_mandel($exponent, $imax, $x, $y, $cx, $cy);
        }

        $x += $xinc;
      }

      $y += $yinc;
    }
  }

  # Return the plot data.
  return $pic_p;
}



#
# IO Functions
#

# Displays a help screen.
# No arguments.
# No return value.

sub show_help
{
  print << "Endofblock";

 High-order Mandelbrot program.
 Written by Christopher Thomas.

 This program draws mandelbrot and julia set analogues using relations of
 the form "z <- z^n + c", for arbitrary n.

 Usage: mandel-high <mode> <exp> <imax> <region> <resolution> <output>

  <mode> is "julia" for julia sets, "mandel" for mandelbrot sets.
  <exp> is the exponent n to use. Non-zero real number.
  <imax> is the maximum number of iterations to perform per point.
  <region> is six coordinate pairs of the form "a,b" (for a + bi).
    The first pair defines the upper-left corner of the scanning region.
    The second pair defines the lower-right corner of the region.
    The third pair defines the fixed parameter.
    For Julia sets, the fixed parameter is C, and the scanned parameter
    is the starting value of Z.
    For Mandelbrot sets, the fixed parameter is the starting value of Z
    (set to 0 for the classic set), and the scanned parameter is C.
  <resolution> is the resolution of the image to render, of the form
    "axb".
  <output> is the filename to write the image to. This is a P3 PPM
    file (human-readable colour bitmap).

Endofblock
}


# Reads command-line arguments.
# Arg 1 points to the argument array.
# Returns a hash containing argument values, or undefined on error.

sub get_args
{
  my ($argv_p);
  my ($hash_p);
  my ($mode, $exponent, $imax, $loc1, $loc2, $loc3, $res, $oname);

  # Default to failure.
  $hash_p = undef;

  # Get argv pointer.
  $argv_p = $_[0];

  if (defined $argv_p)
  {
    # Read arguments.
    $mode = $$argv_p[0];
    $exponent = $$argv_p[1];
    $imax = $$argv_p[2];
    $loc1 = $$argv_p[3];
    $loc2 = $$argv_p[4];
    $loc3 = $$argv_p[5];
    $res = $$argv_p[6];
    $oname = $$argv_p[7];

    # Check args.
    if ( (("julia" eq $mode) || ("mandel" eq $mode))
      && ((1e-4 < $exponent) || (-1e-4 > $exponent))
      && (0 < $imax)
      && ($loc1 =~ m/^\-?\d+(\.\d+)?\,\-?\d+(\.\d+)?$/)
      && ($loc2 =~ m/^\-?\d+(\.\d+)?\,\-?\d+(\.\d+)?$/)
      && ($loc3 =~ m/^\-?\d+(\.\d+)?\,\-?\d+(\.\d+)?$/)
      && ($res =~ m/^[123456789]\d*x[123456789]\d*$/)
      && ($oname =~ m/\S/) )
    {
      # Args seem ok. Store.

      $hash_p = {};

      $$hash_p{mode} = $mode;
      $$hash_p{exponent} = $exponent;
      $$hash_p{imax} = $imax;

      $loc1 =~ m/(\S+)\,(\S+)/;
      $$hash_p{scan1r} = $1;
      $$hash_p{scan1i} = $2;
      $loc2 =~ m/(\S+)\,(\S+)/;
      $$hash_p{scan2r} = $1;
      $$hash_p{scan2i} = $2;
      $loc3 =~ m/(\S+)\,(\S+)/;
      $$hash_p{fixedr} = $1;
      $$hash_p{fixedi} = $2;

      $res =~ m/(\S+)x(\S+)/;
      $$hash_p{width} = $1;
      $$hash_p{height} = $2;

      $$hash_p{oname} = $oname;
    }
  }

  # Return the argument hash.
  return $hash_p;
}


# Saves the generated picture data to an image file.
# Arg 1 points to the argument hash.
# Arg 2 points to the plot data.
# No return value.

sub write_pic
{
  my ($args_p, $pic_p);
  my ($width, $height, $imax);
  my ($x, $y, $icount);

  # Get args.
  $args_p = $_[0];
  $pic_p = $_[1];

  # Draw the picture if we can.
  if ((defined $args_p) && (defined $pic_p))
  {
    if (!open(OUTFILE, ">$$args_p{oname}"))
    {
      print "### Couldn't write to \"$$args_p{oname}\".\n";
    }
    else
    {
      # Get parameters.
      $width = $$args_p{width};
      $height = $$args_p{height};
      $imax = $$args_p{imax};

      # Header.
      print OUTFILE "P3\n$width $height\n255\n";

      # Draw the picture.
      for ($y = 0; $y < $height; $y++)
      {
        for ($x = 0; $x < $width; $x++)
        {
          $icount = $$pic_p[$y][$x];
          if ($imax == $icount)
          {
            print OUTFILE "    0   0   0";
          }
          else
          {
            print OUTFILE $palette[$icount % $palsize];
          }

          # Wrap after every 4th pixel.
          if (3 == ($x & 0x03))
          {
            print OUTFILE "\n";
          }
        }

        if ($width & 0x03)
        {
          # Width not a multiple of 4, so we need another newline.
          print "\n";
        }
      }
    }
  }
}



#
# Main Program
#

my ($args_p);
my ($pic_p);


# Read arguments, and proceed if we have them.
$args_p = get_args(\@ARGV);
if (!(defined $args_p))
{
  show_help();
}
else
{
  $pic_p = gen_pic($args_p);
  if (!(defined $pic_p))
  {
    print "### ERROR: Unable to generate picture!\n";
  }
  else
  {
    write_pic($args_p, $pic_p);
  }
}


#
# This is the end of the file.
#