#!/usr/bin/perl
##########################################################################
# Program    : JuliaGen                                                  #
# Author     : Apostolos Syropoulos (apostolo@obelix.ee.duth.gr)         #
# Date       : June 1, 1997                                              #
# Language   : Perl                                                      #
# Version    : 1.0                                                       #
# Description: Program that computes a Julia set and outputs either a    #
#              PostScript file or a PDF file that contains the image of  #
#              the set. The program computes the sets with the Dinstance #
#              Estimator method, as this is described in                 #
#              "The Science of Fractal Images",                          #
#              H.-O. Peitgen and D. Saupe editors,                       #
#              Springer-Verlag, New York,                                #
#              1988, pp 198.                                             #
#                                                                        #
##########################################################################
################# sub max computes the maximum of two numbers ############
##########################################################################
sub max
{
    my ($a, $b)=@_;
    return ($a>$b) ? $a : $b;
}
##########################################################################
### transforms a decimal up to 15 to its hexadecimal equivalent ##########
##########################################################################
sub hexdigit
{
    my ($a)=@_;
    
    $HEX++;
    if($a < 10)
    {
	return"$a";
    }
    else
    {
	return "A" if($a==10);
	return "B" if($a==11);
	return "C" if($a==12);
	return "D" if($a==13);
	return "E" if($a==14);
	return "F" if($a==15);
    }
}
#################################################################
############ Computes the distance of x0, y0 from the J-Set #####
#################################################################
sub JSetDist
{
   my ($x_0, $y_0)=@_;
   $dist=0;
   $iter=0;
   $xorbit[0]=$x_0;
   $yorbit[0]=$y_0;
   $x2=$x_0*$x_0;
   $y2=$y_0*$y_0;
   while ($iter < ($maxiter-1) && ($x2+$y2)<$huge)
   {
      $temp=$x2-$y2+$Cx;
      $y_0=2.0*$x_0*$y_0+$Cy;
      $x_0=$temp;
      $iter++;
      $x2=$x_0*$x_0;
      $y2=$y_0*$y_0;
      $xorbit[$iter]= $x_0;
      $yorbit[$iter]= $y_0;
   }
   $flag=0;
   if (($x2+$y2) > $huge)   
   {
      $xder=0.0;
      $yder=0.0;
      $i=0;
      OVERFLOW: while ($i<$iter) 
      {
         $temp=2.0*($xorbit[$i]*$xder-$yorbit[$i]*$yder)+1.0;
         $yder=2.0*($yorbit[$i]*$xder+$xorbit[$i]*$yder);
         $xder=$temp;
         $flag=1, last OVERFLOW if (&max(abs($xder),abs($yder)) > $overflow);
	 $i++;
      }
     
      $dist=log($x2+$y2)*sqrt($x2+$y2)/sqrt($xder*$xder+$yder*$yder) 
                                                       if (!$flag);
   }
   return $dist;
}
###################################################################
############ Emits the PS file header #############################
###################################################################
sub emitPSheader
{
    my($date);
    $date=localtime;
    open(JULIA, ">$FN")|| die "Can't create file $FN\n";
    if($type==1)
    {
	print JULIA <<PSHeader;
%!PS-Adobe-2.0
%%Title: Julia set for z --> z^2 + ($Cx+i$Cy) 
%%Creator: JuliaGen, Copyright 1997 AS
%%For: Apostolos Syropoulos, apostolo\@obelix.ee.duth.gr
%%CreationDate: $date
%%Pages: 1
%%DocumentPaperSizes: $paper
%%EndComments
PSHeader
    }
    elsif($type==2)
    {
       print JULIA <<EPSHeader;
%!PS-Adobe-2.0 EPSF-3.0
%%Title: Julia set for z --> z^2 + ($Cx+i$Cy)
%%Creator: JuliaGen, Copyright 1997 AS
%%For: Apostolos Syropoulos, apostolo\@obelix.ee.duth.gr
%%CreationDate: $date
%%BoundingBox: 85 170 426 511
%%EndComments
EPSHeader
    }       
    print JULIA <<RestHEADER;
%%BeginProlog
/cm { 28.35 mul } def
/picstr $Xscreen 7 add 8 idiv string def
%%EndProlog
%%Page: 1 1
gsave
3 cm 6 cm translate
12 cm 12 cm scale
$Xscreen $Yscreen 1 [$Xscreen 0 0 $Yscreen 0 0]
{  currentfile 
    picstr readhexstring pop }
image
RestHEADER
    close JULIA;
}
#####################################################################
############## returns date in format suitable for PDF ##############
#####################################################################
sub PDFDate
{
    my($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst)=localtime;
    my ($d, @A);
    $mon++;
    $d=localtime;
    @A=split(/ /,$d);
    $year=$A[4];
    $sec="0$sec" if(length($sec)==1);
    $min="0$min" if(length($min)==1);
    $hour="0$hour" if(length($hour)==1);
    $mday="0$mday" if(length($mday)==1);
    $mon="0$mon" if(length($mon)==1);
    return "(D:$year$mon$mday$hour$min$sec)";
}
#####################################################################
####### Emits the PDF file "header" #################################
#####################################################################
sub emitPDFheader
{
    my ($urx, $ury)=(596, 842);
    my ($bindata)=chr(199)x4;
    ($urx, $ury)=(612, 792) if($paper eq "letter");
    open(JULIA, ">$FN")|| die "Can't create file $file_name\n";
    binmode JULIA;
    print JULIA <<PDFstart;
%PDF-1.1
%$bindata
1 0 obj
<<
/Type/Catalog
/Pages 2 0 R
>> 
endobj
2 0 obj
<<
/Type/Pages
/Kids [3 0 R]
/Count 1
>>
endobj
3 0 obj
<<
/Type/Page
/MediaBox [0 0 $urx $ury]
/Parent 2 0 R
/Resources
<</ProcSet 7 0 R /XObject <</Im0 5 0 R>> >>
/Contents 4 0 R
>>
endobj
4 0 obj
<</Length 47 >>
stream
1 0 0 1 100 170  cm
340 0 0 340 0 0 cm
/Im0 Do
endstream
endobj
5 0 obj
<<
/Type/XObject
/Subtype/Image
/Name/Im0
/Width $Xscreen
/Height $Xscreen
/BitsPerComponent 1
/ColorSpace/DeviceGray
/Filter/ASCIIHexDecode
/Length 6 0 R
>>
stream
PDFstart
    close JULIA;
}
############################################################
############ prints PS file trailer ########################
############################################################
sub emitPStrailer
{
    open(JL,">>$FN")||
    die "Something wrong has happened: I can't open file $FN\n"; 
    print JL <<TRAILER;

grestore
showpage
%%EOF
TRAILER
    close JL;
}
############################################################
############ prints PDF file trailer ########################
############################################################
sub emitPDFtrailer
{
    my ($imageLen, $date);
    open(JL, ">>$FN")||
    die "Something wrong has happened: I can't open file $FN\n";
    binmode JL;
############################
### calculate image size
############################
    $imageLen=$HEX + 2 + $newlines;
    $date=&PDFDate;
    print JL <<MORE;
>
endstream
endobj
6 0 obj
$imageLen
endobj
7 0 obj
[/PDF/ImageB]
endobj
8 0 obj
<<
/Author (Apostolos Syropoulos, apostolo\@obelix.ee.duth.gr)
/Creator (JuliaGen)
/CreationDate $date
/Title (Julia set for z --> z^2 + \\($Cx+i$Cy\\))
>>
endobj
MORE
    close JL;
############################
### calculate obj sizes
############################
    $bytes=0;
    $i=0;
    open(JL,"$FN")||
    die "Something wrong has happened: I can't open file $FN\n"; 
    binmode JL;
    while(<JL>)
    {
        if(/^(\d+) (\d+) obj$/)
        {
           $BYTES[$i]=$bytes;
           $i++;
           $bytes +=length($_);
        }
        else
        {
           $bytes += length($_);
        }
    }
    close JL;
    $BYTES[$i]=$bytes;
    open(JL, ">>$FN")||
    die "Something wrong has happened: I can't open file $FN\n";
    binmode JL;
    print JL "xref\n";
    print JL "0 ",$i+1,"\n";
    print JL "0"x10," 65535 f \n";
    for ($j=0; $j<$i; $j++)
    {
        print JL "0" x (10-length($BYTES[$j])), $BYTES[$j], " 00000 n \n";
    }
    print JL "trailer\n";
    print JL "<<\n/Size $i\n";
    print JL "/Root 1 0 R\n";
    print JL "/Info 8 0 R\n";
    print JL ">>\n";
    print JL "startxref\n";
    print JL $BYTES[$i],"\n";
    print JL"%%EOF";
    close JL;
}
############################################################
#################### Main Program ##########################
############################################################
if(!@ARGV)
{
    print <<USAGE;
This is JuliaGen, version 1.0 <24/07/1997>.
Generating images of Julia Sets.
Usage: juliagen [-tn] [-sm] [-pP] [-r] <output file name>
-tn is used to specify the output file format, possible values for n:
n=1 => PostScript (default)
n=2 => EPSF
n=3 => PDF

-sm is used to specify the image size in pixels, possible values for m:
m=0 => 64x64
m=1 => 128x128
m=2 => 256x256 (default)
m=3 => 512x512
m=4 => 1024x1024
m=5 => 2048x2048

-pP is used to specify the paper size format, possible values for P: 
A4 (default) and letter 

-r is used by the user who wants to enter his/her values for xmin, 
xmax, ymin and ymax (default values: xmin=ymin=-1.6 and xmax=ymax=1.6).

USAGE
    exit(0);
}
####################################################################
## get command line arguments ######################################
####################################################################
$type=1;     #default behavior: emit PostScript
$size=2;     #default image size: 256x256
$paper="A4"; #default paper size
$rectagle=0; #default: do not enter xmin, xmax, ymin, and ymax.
SWITCHES: while($_ = $ARGV[0], /^-/)
{
    shift;
    if(/^-t(\d)/)
    {
	$type= ($1<1 || $1>3) ? 1 : $1;
    }
    elsif(/^-s(\d)/)
    {
	$size= ($1>5 || $1<0) ? 0 : $1;
    }
    elsif(/^-r$/)
    {
        $rectangle=1;
    }
    elsif(/^-p(\w*)/)
    {
	if($1 eq "A4" || $1 eq "letter")
        {
	    $paper=$1;
	}
        else
        {
            die "Not supported paper size format $1\n";
	}
    }
    elsif(/^-[^tspr].*/)
    {
	die "Invalid switch: $_\n";
    }
    elsif (!@ARGV)
    {
	last SWITCHES;
    }
}    
die "No output file name specified!\n" if(!@ARGV);
$FN=$ARGV[0];
warn "\n****File $FN exists and it will be overwriten.\n\n" if(-e $FN);
############################################################
################# Data Input ###############################
############################################################
print "This is JuliaGen, version 1.0 <24/07/1997>.\n";
print "Enter parameters for Julia Set computation\n";
print "------------------------------------------\n\n";
print "Enter value for `Cx'\n? ";
chomp($Cx=<STDIN>);
print "Enter value for `Cy'\n? ";
chomp($Cy=<STDIN>);
if($rectangle)
{
    print "\n";
    print "Enter \"drawing\" area:\n";
    print "-----------------------\n";
    print "Enter value for `xmin'\n? ";
    chomp($xmin=<STDIN>);
    print "Enter value for `xmax'\n? ";
    chomp($xmax=<STDIN>);
    print "Enter value for `ymin'\n? ";
    chomp($ymin=<STDIN>);
    print "Enter value for `ymax'\n? ";
    chomp($ymax=<STDIN>);
}
else
{
    $xmin=$ymin=-1.6;
    $xmax=$ymax=1.6;
}
############################################################
############### Initilizations #############################
############################################################
CASE:{
    if($size == 0) { $Xscreen = $Yscreen = 64;   last CASE;}   
    if($size == 1) { $Xscreen = $Yscreen = 128;  last CASE;}   
    if($size == 2) { $Xscreen = $Yscreen = 256;  last CASE;}   
    if($size == 3) { $Xscreen = $Yscreen = 512;  last CASE;}   
    if($size == 4) { $Xscreen = $Yscreen = 1024; last CASE;}   
    if($size == 5) { $Xscreen = $Yscreen = 2048; last CASE;}
}   
$Max       = 200;
$maxiter   = $Max;
$overflow  = 1.0e100;
$huge      = 20000.0;
#$Cx         = -0.74543; #these values for p & q create a
#$Cy         = 0.11301;  #beautiful picture.
$threshold=0.008;
$DeltaX=($xmax-$xmin)/($Xscreen-1);
$delta=$threshold*$DeltaX;
$DeltaY=($ymax-$ymin)/($Yscreen-1);
#############################################################
################## Print Headers ############################
#############################################################
&emitPSheader if ($type == 1 || $type == 2);
&emitPDFheader if($type == 3);
#############################################################
################### Start computation #######################
#############################################################
open(JL,">>$FN")||die "Something wrong has happened: I can't open file $FN\n";
$newlines=0;
$octet=0;
$hdc=1;
$HEX=0;
for($ny=0; $ny < $Yscreen; $ny++)
{
   $y0=$ymin+$ny*$DeltaY;
   for ($nx=0; $nx < $Xscreen; $nx++)
   {
      $x0=$xmin+$nx*$DeltaX;
      $distance=&JSetDist($x0,$y0);
      $sc=($nx+1) % 8; 
      if($sc == 0)
      {
	  $quotient=int($octet / 16);
          $remainder=$octet % 16;
          $hexnum= &hexdigit($quotient) . &hexdigit($remainder);
          print JL $hexnum;
          if(($hdc % 30)==0) { $newlines++; print JL "\n"; }
          $hdc++;
          $octet=0;
      } 
      if ($distance >= $delta)
      {
          $octet += 2**(7-$sc);
      }
   }
}
close JL; 
#############################################################
############### Print trailers ##############################
#############################################################
&emitPStrailer if($type == 1 || $type == 2);
&emitPDFtrailer if($type == 3);
__END__