#! /bin/env perl

# RELEASE 1.5

##########################################################################
# Calibration constants.  
##########################################################################

#
# This is the point in the image that is zenith on the sky
#

my $xcen = 298.521;
my $ycen = 250.5;

##########################################################################
# Other constants
##########################################################################

my $PI = 4.0*atan2(1.0,1.0);

##########################################################################
# Extract command line arguments.
##########################################################################

if ($#ARGV + 1 != 4) {
	print "Must specify: input_file output_file alt and az\n";
	die;
}

my $input_file = $ARGV[0];
my $output_file = $ARGV[1];
my $alt = $ARGV[2];
my $az = $ARGV[3];

##########################################################################
# Begin computation
##########################################################################

#
# We use zenith distance.
#

my $zd = 90.0 - $alt;

#
# 10.175 deg offset in image orientation
#

my $theta = $az + 12.936;

if ($theta > 360.0) {
	$theta -= 360.0;
}

if ($theta < 0.0) {
	$theta += 360.0;
}

$theta *= $PI / 180.0;

#
# this is the fit to get $r in image from zenith angle. 
#

$r = 3.203*$zd + 0.002456*$zd*$zd - 0.000121*$zd*$zd*$zd;

#
# convert back to X and Y
#

$x = $xcen - $r*sin($theta);
$y = $r*cos($theta) + $ycen;

$y = 480.0 - $y;

##########################################################################
# Computation finished.
##########################################################################

my $x1 = $x - 10;
my $x2 = $x + 10;
my $y1 = $y - 10;
my $y2 = $y + 10;

system("gm convert -depth 8 -gamma 1.5"
	. " -fill orange1"
    . " -draw \'line $x1,$y $x2,$y\'"
    . " -draw \'line $x,$y1 $x,$y2\'"
    . " $input_file $output_file");