#!/usr/bin/perl

use Getopt::Long;
use GD::Simple;

use Math::Complex;
use Math::Round::Var;
$rnd = Math::Round::Var->new(0.00001);
$rnddisp = Math::Round::Var->new(0.001);


$PATTON = 1.86;
$RAYLEIGH = 2.00;

$Zindex = 0;
$Zindexmax = 13; #odd number (from 0..12 = 13) will have solid printed center.
$focallength = 1000; # in millimeters, so sue me.
$wavelength = 560; # nanometers, 560 is ~daylight average
$imagesize = '5000';
$scalingconstant = 10; # multiplier of answer in millimeters to a value useful for monitor viewing (arbitary)
GetOptions(
       'wave:s' => \$wavelength,
       'foc:s' => \$focallength,
       'zones:s' => \$Zindexmax,
       'px:s' => \$imagesize,
       'o:s' => \$filename,
       'scale:s' => \$scalingconstant,
       'help' => \$help
);
unless ($filename) {
	$filename = "zoneplate" . "_" . $wavelength . "_" . $focallength . "_" . $Zindexmax . ".png";
}

if ($help) {
	print "Example: ./zoneplate.pl -wave 560 -foc 200 -zones 13 -px 5000 -o zoneplate01.png\nThis would calculate the width of a zone for a zone plate for light of a wavelegnth of 560 nm"
	. ", the lens would have a focal length of 200 millimeters and 13 alternating zones (counting the middle SOLID as 1). The generated png/image named zoneplate01.png would be 5000x5000 pixels.\n"
	. "(http://superkuh.ath.cx/ , superkuh\@gmail.com)\n";
	exit;
}

# correct for arbitray units to millimeters
$wavelength *= 1e-6;  # 1 nanometer = 1.0 × 10-6 millimeters
#$wavelength *= 0.000001;  # 1 nanometer = 1.0 × 10-6 millimeters

@indices = '';
$diff = 0;
$previous = 0; # the center

print "In MILLIMETERS! Wavelength: $wavelength- Zones: $Zindexmax-\n";
print "Zone Index - Zone Width - Diff from last zone\n\n";
foreach $Zindex (1..$Zindexmax) { #radial distance out from center to block (center blocked if uneven)
	my $zonewidth = $PATTON * sqrt($focallength * $Zindex * $wavelength);
	# Zi = PATTON x squareroot(focal length X index# X wavelength)
	my $round = $rnd->round($zonewidth);
	
	#such a crap way, but it works.
	$firstzonewidth = $round if $Zindex == 1;

	if ($Zindex > 1) {
		
		if ($Zindex == 2) {
			$diff = $round - $firstzonewidth;
			$previous = $round; #set for next loop
			#print "\nCrap\n";
		} else {
			$diff = $round - $previous;
			$diff = $rnd->round($diff);
			$previous = $round; #set for next loop
		}
	}
	

	#$diff = sprintf("%.3f", $diff); # just truncation; not accurate.
	push (@indices, $round);
	print "$Zindex - $round mm - $diff mm\n";

}
# get highest value for pixel size constraints
#@sorted = (sort { $b <=> $a } @indices)[0];
#$maximagesize = pop(@sorted);
#$maximagesize *= $scalingconstant;


print "-------\nIndices: " . scalar(@indices) . "\n\n";
$guessedindexbecauseimlazy = scalar(@indices);
foreach (reverse @indices) {
		
	$guessedindexbecauseimlazy += -1;
	$radius = $_ * $scalingconstant; 
	
	$stop = 0;
	while ($stop == 0) {
		if ($imagesize < $radius) {
			print "\nDue to arbitrary scaling everything won't fit in your specified imagesize ($imagesize pixels)."
			. "\nDo you want to scale the image up to at least $radius pixels (yes)? ";
			my $answer = <STDIN>;
			if ($answer =~ /^y(es)?/i) {
				print "\n\nAnswer: $answer";
				$imagesize = $radius * 1.1;
				# create a new image
				$img = GD::Simple->new($imagesize,$imagesize); # $imagesize is pixels long
				$stop = 1;
			} else {
				print "shit\n";				
			}
		}
	}
	
	print "\nRadius:$radius\nIndex:$guessedindexbecauseimlazy";
	
	if ($guessedindexbecauseimlazy % 2 == 0) {
		print "\nEven.\n";
		drawcircle($radius,'black','white');
	} else {
		print "\nOdd.\n";
		drawcircle($radius,'black','black');
	}
}


# convert into png data
open (PLATEIMAGE,">","$filename");
print PLATEIMAGE $img->png;
close MPLATEIMAGE;
	
	
sub drawcircle {
	
	my ($zonewidthtodraw, $zonelinecolor, $zonecolor) = @_;
	
    $img->bgcolor($zonecolor);
    
    $img->fgcolor($zonelinecolor);
    
    #go to the center of the image
    $img->moveTo(($imagesize/2),($imagesize/2));
    #print a circle
	$img->ellipse($zonewidthtodraw,$zonewidthtodraw);

}

sub removezoneCantorDust {
	my $index = shift;
	# Let us start by considering the irradiance at a given point on 
	# the optical axis that is provided by an optical system with a 
	# rotationally invariant pupil function described by p(r). Within 
	# the Fresnel approximation, this magnitude is given, as a function
	# of the axial distance R from the pupil plane, as:
	# I(R) = ((2pi)/(wavelength*R))^2 | Integral[0,alpha]( p(rsub(0))exp[-i * (pi/wavelength*R) * (rsub(0))^2]*rsub(0)*drsub(0) |^2
	# s = (rsub(0)/alpha)^2 - 0.5
	# q(s) = rect(s)rect{mod[s + (p - 1)/2, p] / p}
	# visual example: take 1d line of S=3 and rotate it about the origin 2pi drawing each small angle change after each S.
	# generate S=0, draw the line where 0='white pin-hole' from (0,0)..(0,pi/2)..(0,pi)..(0,3pi/2)..(0,0)
	# generate S=1, draw the line where 0='white pin-hole' from (0,0)..(0,pi/2)..(0,pi)..(0,3pi/2)..(0,0)
	# generate S=2, draw the line where 0='white pin-hole' from (0,0)..(0,pi/2)..(0,pi)..(0,3pi/2)..(0,0)
	# generate S=3, draw the line where 0='white pin-hole' from (0,0)..(0,pi/2)..(0,pi)..(0,3pi/2)..(0,0)	
	# ...................................................... S=0
	# 000000000000000000..................000000000000000000 S=1
	# 000000......000000..................000000......000000 S=2
	# 00..00......00..00..................00..00......00..00 S=3

	# A Fractal ZP can be understood as a fresnel ZP qsub(ZP)[s, p(N,S)] with missing clear zones.
	# p(N,S) = 2 / (2N - 1)^S
	
	# check if number of zones is divisible by 3; 
	#until ($Zindexmax % 3) {
		# if not, ask user if it's okay to change. 
		# change
	#}
	
	# all zones printer for now/testing
	if (0) {
		return 1; # if this zone is ommitted
	}
	return 0; # if this zone is printed
}






sub PGPlotcircle {
		
	set_begin({
  		  file => "pgplot_test.ps",
 	});

	set_environment({
     		x_min   =>  0,
    		x_max   =>  50,
    		y_min   =>  0,
     		y_max   =>  10,
	});

	draw_circle({ x      =>              	$x,    # Required
               y      =>                	$y,    # Required
               radius =>             		$rad,  # Required
               color  =>          'Orange',
               width  =>                 10,
               style  =>  'DotDashDotDash',
               fill   =>    'CrossHatched',
 });
	
	set_end;
}



&makesvg if (0);
sub makesvg {

	 # create an SVG object
   	 my $svg= SVG->new(width=>200,height=>200);

    	 # use explicit element constructor to generate a group element
   	 my $y=$svg->group(
      	  id    => 'group_y',
      	  style => { stroke=>'red', fill=>'green' }
   	 );

    	# add a circle to the group
    	$y->circle(cx=>100, cy=>100, r=>50, id=>'circle_in_group_y');

	open(ZONEPLATE, ">zoneplate.svg") or die $!;
	# create and add a circle using the generic 'tag' method
	$z->tag('circle', cx=>50, cy=>50, r=>100, id=>'circle_in_group_z');
	print ZONEPLATE $svg->xmlify;
}

# someone elses data for checking output. mine doesn't exactly match up
# I think this might be mostly different pinhole constants and a little
# bit of rounding errors (probably on my part.
#
#ZONE	Focal Length in mm
#INDEX	25	50	75	100	125	150	175	200
#1	0.23	0.33	0.41	0.47	0.52	0.57	0.62	0.66
#2	0.33	0.47	0.57	0.66	0.74	0.81	0.88	0.94
#3	0.41	0.57	0.70	0.81	0.91	0.99	1.07	1.15
#4	0.47	0.66	0.81	0.94	1.05	1.15	1.24	1.33
#5	0.52	0.74	0.91	1.05	1.17	1.28	1.39	1.48
#6	0.57	0.81	0.99	1.15	1.28	1.41	1.52	1.62
#7	0.62	0.88	1.07	1.24	1.39	1.52	1.64	1.75
#8	0.66	0.94	1.15	1.33	1.48	1.62	1.75	1.88
#9	0.70	0.99	1.22	1.41	1.57	1.72	1.86	1.99
#10	0.74	1.05	1.28	1.48	1.66	1.82	1.96	2.10
#11	0.78	1.10	1.35	1.56	1.74	1.91	2.06	2.20