package Period; # Some::Module;  # assumes Some/Module.pm
    use strict;
    use warnings;
    BEGIN {
        use Exporter   ();
        our ($VERSION, @ISA, @EXPORT, @EXPORT_OK, %EXPORT_TAGS);
        # set the version for version checking
        $VERSION     = 1.00;
        @ISA         = qw(Exporter);
        @EXPORT      = qw(&Analysis &HighStep &MediumStep &LowStep &PlotAnalysis &Phase &FreqFormat &ShowData);
        %EXPORT_TAGS = ( );     # eg: TAG => [ qw!name1 name2! ],
        # your exported package globals go here,
        # as well as any optionally exported functions
        @EXPORT_OK   = qw(@ObsJD @ObsMag $N $GNUPlot $Montage $OUTDIR $HOME $RM $PERL $Extrema &Timespan &Brightest &Faintest);
    }
    our @EXPORT_OK;
    # exported package globals go here
    our @ObsJD;
    our @ObsMag;
    our $N;  # Number of observations

    our $GNUPlot;
    our $Montage;
    our $OUTDIR;
    our $HOME;
    our $Extrema;
    our $PERL;
    our $RM;

    # non-exported package globals go here

    our $Nb;
    our $Nc;
    our $compact;


    # initialize package globals, first exported ones

    $N = 0;

    $GNUPlot = "$ENV{PeriodGNU}pgnuplot";
    $Montage = "$ENV{PeriodIMM}montage";
    $OUTDIR = $ENV{PeriodOUT};
    $HOME = $ENV{PeriodDIR};
    $PERL = "$ENV{PeriodPERL}perl";
    $Extrema = "$PERL -s $ENV{PeriodDIR}Extrema.pl";
    $RM = $ENV{PeriodRM};

    # then the others (which are still accessible as $Some::Module::stuff)

    $Nb = 5;
    $Nc = 2;
    $compact = 1;


    # all file-scoped lexicals must be created before
    # the functions below that use them.
    # file-private lexicals go here

    # here's a file-private function as a closure,
    # callable as &$priv_func;  it cannot be prototyped.

    # make all your functions, whether exported or not;
    # remember to put something interesting in the {} stubs

    # this one isn't exported, but could be called!

    END { }       # module clean-up code here (global destructor)
    ## YOUR CODE GOES HERE


#procedure ecliptical;
#var x, y : real;
#begin
#x := cos(dec * rad) * cos(RA * 15 * rad);
#y := 0.917 * cos(dec * rad) * sin(RA * 15 * rad) + 0.398 * sin(dec * rad);
#cosbeta := sqrt(x * x + y * y);
#if y = 0 then lambda := 0 else lambda := 2 * arctan((cosbeta - x)/y)
#end;
#function helio_correction(jd : real) : real;
#(* accurate to better than 0.0001d for about a century *)
#var M : real;
#begin
#jd := jd + epoch - 2451545.0;
#M := 6.24 + 0.0172020 * jd;
#helio_correction := 0.00578 * cosbeta * (1 - 0.017 * cos(M))
#                    * cos(4.883 + 0.033 * sin(M) + 0.01720212 * jd - lambda)
#end;


# Get the average JD
sub AvgJD ()
{
	my $JD = 0;

	for (my $i=1;$i<=$N;$i++) 
	{
		$JD += $ObsJD[$i];
	}
	if ($N > 0)
	{
		return $JD/$N;
	}
	else
	{
		return 0;
	}
}


# Get the total time span of the observations
sub Timespan ()
{
	my $JD0 = $ObsJD[1];
	my $JD1 = $ObsJD[1];

	for (my $i=1;$i<=$N;$i++) 
	{
		if ($ObsJD[$i] < $JD0)
		{
			$JD0 = $ObsJD[$i];
		}
		if ($ObsJD[$i] > $JD1)
		{
			$JD1 = $ObsJD[$i];
		}
	}
	return $JD1 - $JD0;
}


# Determine the high step rate (for frequency search)
sub HighStep ()
{
	return 1/Timespan/20;
}

# Determine the medium stop rate (for frequency search)
sub MediumStep ()
{
	return 1/Timespan/15;
}

# Determine the low step rate (for frequency search)
sub LowStep ()
{
	return 1/Timespan/10;
}


# Get the brightest observation
sub Brightest ()
{
	my $mag = $ObsMag[1];

	for (my $i=1;$i<=$N;$i++) 
	{
		if ($ObsMag[$i] < $mag)
		{
			$mag = $ObsMag[$i];
		}
	}
	return $mag;
}

# Get the faintest observation
sub Faintest ()
{
	my $mag = $ObsMag[1];

	for (my $i=1;$i<=$N;$i++) 
	{
		if ($ObsMag[$i] > $mag)
		{
			$mag = $ObsMag[$i];
		}
	}
	return $mag;
}


# Show all the data
sub ShowData ($)
{
	my ($delimiter) = @_;
	if (!$delimiter)
	{
		$delimiter = " ";
	}
	for (my $i=1;$i<=$N;$i++) 
	{
		print "$ObsJD[$i]$delimiter$ObsMag[$i]\n";
	}
}


# Make a GNUPlot command file to plot the analysis
sub PlotAnalysis ($$$$)
{
	# Create PlotAnalysis command file
	my ($GifFile, $FreqFile, $analysis, $out) = @_;

#print "$out\n";
	open (GNUPLOT, $out);
	print GNUPLOT "set term gif\n";
	print GNUPLOT "set xlabel 'Frequency (c/d)'\n";
	#print GNUPLOT "set xrange [$F0:$F1]\n";
	if ($analysis eq "pdm")
	{
		print GNUPLOT "set title 'PDM Analysis'\n";
		print GNUPLOT "set ylabel 'Theta (PDM)'\n";
		print GNUPLOT "set yrange [0:1.05]\n";
		print GNUPLOT "set ytics 0.1\n";
		print GNUPLOT "set format y '%.1f'\n";
		print GNUPLOT "set mytics 2\n";
	}
	else
	{
		print GNUPLOT "set title 'Lomb-Scargle periodogram'\n";
		print GNUPLOT "set ylabel 'Amplitude'\n";
	}
	print GNUPLOT "set output '$GifFile.gif'\n";
	print GNUPLOT "plot '$FreqFile' notitle with lines\n";
	print GNUPLOT "set output\n";
	print GNUPLOT "exit\n";
	close GNUPLOT;
}


# Make a GNUPlot command file to plot the phase diagram
sub PlotPhase ($$$$)
{
	# Create PlotPhase command file
	my ($PhaseFile, $graph, $title, $out) = @_;
	my $faint = Faintest;
	my $bright = Brightest;
	my $margin = ($faint - $bright)/10;
	$faint += $margin;
	$bright -= $margin;

#print "$out\n";
	open (GNUPLOT, $out);
	print GNUPLOT "set term gif\n";
	print GNUPLOT "set nokey\n";
	print GNUPLOT "set title '$title'\n";
	print GNUPLOT "set bmargin 5\n";
	print GNUPLOT "set lmargin 10\n";
	print GNUPLOT "set rmargin 10\n";
	print GNUPLOT "set xlabel 'Phase'\n";
	print GNUPLOT "set xrange [-0.2:1.2]\n";
	print GNUPLOT "set xtics 0.2\n";
	print GNUPLOT "set mxtics 2\n";
	print GNUPLOT "set format x '%.1f'\n";
	print GNUPLOT "set ylabel 'Mag'\n";
	print GNUPLOT "set yrange [$faint:$bright]\n";
	print GNUPLOT "set output '$graph.gif'\n";
	print GNUPLOT "set multiplot\n";
	print GNUPLOT "plot '$PhaseFile' using 1:2 notitle with points pointtype 7\n";
	print GNUPLOT "set xrange [0.8:2.2]\n";
	print GNUPLOT "set noxtics\n";
	print GNUPLOT "set nomxtics\n";
	print GNUPLOT "set xlabel ''\n";
	print GNUPLOT "plot '$PhaseFile' using 1:2 notitle with points pointtype 7\n";
	print GNUPLOT "set xrange [-1.2:0.2]\n";
	print GNUPLOT "plot '$PhaseFile' using 1:2 notitle with points pointtype 7\n";
	print GNUPLOT "set nomultiplot\n";
	print GNUPLOT "set output\n";
	print GNUPLOT "exit\n";
	close GNUPLOT;
}


# Calculate the phase and output
sub Phase ($$)
{
	my ($freq, $out) = @_;
	my $phase;
	my $JDt = AvgJD;

#print "$out\n";
	open (PHASE, $out);

	# Loop through all observations
	for (my $i=1;$i<=$N;$i++) 
	{
		$phase = ($ObsJD[$i] - $JDt) * $freq;
		$phase = $phase - int($phase);
		if ($phase < 0)
		{
			$phase += 1;
		}
		printf PHASE "%.3f %s\n", $phase ,$ObsMag[$i];
	}
	close PHASE;
}


# Format the frequency (correct number of decimals depending on the step size)
sub FreqFormat ($)
{
	my ($Fstep) = @_;
	my $decimals = abs(int(log($Fstep)/log(10))) + 2;
	return sprintf "%%%d.%df", $decimals + 3,$decimals;
}


# Do the PDM analysis
sub Analysis ($$$$)
{
	my ($F0, $F1, $dF, $FreqFile) = @_;
# $F0: First frequency
# $F1: Last frequency
# $dF: Frequency step
# $FreqFile: Output destination (should contain ">" or "|")

	my $phase;
	my $bin;
	my @x;
	my @n;

	if ($N == 0)
	{
		return;
	}

	my $X = 0;
	my $X2 = 0;
	
	my $i;
	my $j;
	my $F;
	my $s2;
	my $theta = 1;
	my $rising = 0;
	my $JDt = AvgJD;

	for ($i=1;$i<=$N;$i++) 
	{
		$X += $ObsMag[$i];
		$X2 += $ObsMag[$i] * $ObsMag[$i];
	}

	my $sigma2 = ($X2 - $X * $X/$N)/($N - 1);

	my $fmt = sprintf "%s %%.3f\n", FreqFormat($dF);
#print "$N--$F0--$F1--$dF--$FreqFile--$fmt\n";

	open THETA, "$FreqFile";
	for ($F = $F0; $F <= $F1; $F += $dF)
	{
		# For each frequency
		# Initialize arrays
		for ($j = 0; $j < $Nc; $j++) 
		{
			for ($i = 0; $i < $Nb; $i++) 
			{
				$x[$j]->[$i] = 0;  # ${$x[$j]}[i] = 0;
				$n[$j]->[$i] = 0;
			}
		}

		# Loop through all observations
		for ($i=1;$i<=$N;$i++) 
		{
			for ($j = 0; $j < $Nc; $j++) 
			{
				$phase = ($ObsJD[$i] - $JDt) * $F + $j/$Nb/$Nc;
				$phase = $phase - int($phase);
				$bin = int($Nb * $phase);
				$x[$j]->[$bin] += $ObsMag[$i];
				$n[$j]->[$bin] += 1;
			}
		}

		$s2 = 0;
		for ($j = 0; $j < $Nc; $j++) 
		{
			for ($i = 0; $i < $Nb; $i++) 
			{
				if ($n[$j]->[$i] > 0) 
				{
					$s2 += $x[$j]->[$i] * $x[$j]->[$i] / $n[$j]->[$i];
				}
			}
		}
		$s2 = ($Nc * $X2 - $s2)/$Nc/($N -$Nb);

	        if ($compact)
	        {# Show only extrema of the statistic
        	  if ($s2 < $theta)
	          { # statistic descends
        	    if ($rising > 0)
  	            {# statistic was rising, so has reached minimum
        	      printf THETA $fmt, $F - $dF, $theta/$sigma2;
	            }
        	    $rising = -1;
	          }
        	  elsif ($s2 > $theta)
	          {# statistic rises
        	    if ($rising < 0)
	            {# statistic was descending, so has reached maximum
        	      printf THETA $fmt, $F - $dF, $theta/$sigma2;
	            }
        	    $rising = 1;
	          }
        	  else
		  {
	            $rising = 0;
		  }
	
        	  $theta = $s2;
	        }
	        else
		{
        		printf THETA $fmt, $F, $s2/$sigma2;
		}

# End for each frequency
	}
	close THETA;

# End Analysis
}

    1;  # don't forget to return a true value from the file