#!/usr/local/bin/perl
use FileHandle;

# This program will read in the latest FUSE TLE and recreate the 
# cvz_ram_tool.html web page.

$output_filename = "/data2/violet/htdocs/support/tools/cvz_ram_tool.html";
open (OUTF, ">$output_filename") || die "Cannot open cvz_ram_tool.html file";

system("chmod ug+rw $output_filename");
system("chgrp www $output_filename");

$input_filename = "/data1/fuse/calfuse/calfiles/FUSE.TLE";
open (INP_TLE, "<$input_filename") || die "Cannot open FUSE.TLE file";

$line1 = <INP_TLE>;
$line2 = <INP_TLE>;
$line3 = <INP_TLE>;

$id2 = substr($line2, 2, 5);
$id3 = substr($line3, 2, 5);

if (($id2 != 25791) || ($id3 != 25791)) {
    print STDOUT "Error in TLE:\n";
    print STDOUT $line1;
    print STDOUT $line2;
    print STDOUT $line3;
    die "Error in TLE";
}

close (INP_TLE);

$year = substr($line2, 18, 2);
$doy = substr($line2, 20, 3);
$dayfrac = substr($line2, 23, 9);
$incl = substr($line3, 9, 8);
$raan = substr($line3, 17,8);
$eccen = substr($line3, 26, 7);
$eccen = "0." . $eccen;
$mean_mot = substr($line3, 52, 11);

if ($year > 50) {
    $year += 1900;
} else {
    $year += 2000;
}

$dayfrac = "0".$dayfrac;

&convert_doy;
&convert_dayfrac;
&calculate_a0;

$now=`date`;
$datest = substr($now,4,3)." ".substr($now,8,2).", ".substr($now,24,4);

print OUTF "<TITLE>FUSE CVZ/Ram Calculator</TITLE>\n";
print OUTF "<h1>FUSE Continuous Viewing Zone and Orbit Ram Calculator V1.3</h1>\n";
print OUTF "<h2>June 12, 1999 </h2>\n";

print OUTF "The FUSE Continuous Viewing Zone and Orbit Ram \n";
print OUTF "Calculator can determine when a given target \n";
print OUTF "direction is in the continuous viewing zone. \n";
print OUTF "It also calculates when a given target lies \n";
print OUTF "within 20 degrees of the orbital plane and \n";
print OUTF "cannot be observed due to ram avoidance\n";
print OUTF "constraints. The default orbital elements \n";
print OUTF "were last updated on $month $day, $year. \n";

print OUTF "<HR>\n";
print OUTF "<P>\n";
print OUTF "<FORM METHOD=POST ACTION=\"http://fuse.pha.jhu.edu/cgi-bin/cvz_ram_tool\">\n";

print OUTF "<INPUT TYPE=\"submit\" VALUE=\"Calculate CVZ visibilities\">\n";
print OUTF "<INPUT TYPE=\"reset\" VALUE=\"Erase the input fields\">\n";
print OUTF "<BR>\n";
  
print OUTF "<BR>\n";
print OUTF "<B>Target Parameters:</B>\n";
print OUTF "<BLOCKQUOTE>\n";
print OUTF "<BR>\n";
print OUTF "<B>Right Ascension (format HH:MM:SS.SS):</B>\n";
print OUTF "<INPUT NAME=\"right_ascension\" SIZE=15> <B>(J2000.0)</B>\n";
print OUTF "<BR>\n";
print OUTF "<B>Declination (format -DD:MM:SS.S):</B>\n";
print OUTF "<INPUT NAME=\"declination\" SIZE=15> <B>(J2000.0)</B>\n";
print OUTF "<BR>\n";
print OUTF "</BLOCKQUOTE>\n";
print OUTF "<BR>\n";
print OUTF "<B>Output table parameters:</B>\n";
print OUTF "<BLOCKQUOTE>\n";
print OUTF "<B> Minimum Earth limb angle: </B><INPUT NAME=\"limb_ang\" SIZE=4 VALUE=\"15.0\"> \n";
print OUTF "<B>(deg)</B>\n";
print OUTF "<BR>\n";
print OUTF "<B>Start Date:</B>\n";
print OUTF "<SELECT NAME=\"month\">\n";
print OUTF "<OPTION> <B>$month</B>\n";
print OUTF "<OPTION> January\n";
print OUTF "<OPTION> February\n";
print OUTF "<OPTION> March\n";
print OUTF "<OPTION> April\n";
print OUTF "<OPTION> May\n";
print OUTF "<OPTION> June\n";
print OUTF "<OPTION> July\n";
print OUTF "<OPTION> August\n";
print OUTF "<OPTION> September\n";
print OUTF "<OPTION> October\n";
print OUTF "<OPTION> November\n";
print OUTF "<OPTION> December\n";
print OUTF "</SELECT>\n";
print OUTF "<INPUT NAME=\"day\" SIZE=4 VALUE=\"1\">\n";
print OUTF "<INPUT NAME=\"year\" SIZE=8 VALUE=\"$year\">\n";
print OUTF "<BR>\n";
print OUTF "<B>Time step:<INPUT NAME=\"time_step\" SIZE=6 VALUE=\"1.0\">days</B><BR>\n";
print OUTF "<B>Number of steps:<INPUT NAME=\"number_steps\" SIZE=6 VALUE=\"10\"></B><BR>\n";
print OUTF "</BLOCKQUOTE>\n";
print OUTF "<BR>\n";
print OUTF "<B>FUSE classical Keplarian orbital elements:</B>\n";
print OUTF "<BLOCKQUOTE>\n";
print OUTF "<B>Epoch date of elements:</B>\n";
print OUTF "<SELECT NAME=\"orb_month\">\n";
print OUTF "<OPTION> $month\n";
print OUTF "<OPTION> January\n";
print OUTF "<OPTION> February\n";
print OUTF "<OPTION> March\n";
print OUTF "<OPTION> April\n";
print OUTF "<OPTION> May\n";
print OUTF "<OPTION> June\n";
print OUTF "<OPTION> July\n";
print OUTF "<OPTION> August\n";
print OUTF "<OPTION> September\n";
print OUTF "<OPTION> October\n";
print OUTF "<OPTION> November\n";
print OUTF "<OPTION> December\n";
print OUTF "</SELECT>\n";
print OUTF "<INPUT NAME=\"orb_day\" SIZE=4 VALUE=\"$day\">\n";
print OUTF "<INPUT NAME=\"orb_year\" SIZE=8 VALUE=\"$year\">\n";
print OUTF "<BR>\n";
print OUTF "<B>Epoch UT time of elements (format HH:MM:SS.SS):</B>\n";
printf OUTF "<INPUT NAME=\"orb_time\" SIZE=12 VALUE=\"%02d:%02d:%05.2f\">\n",$hour,$minute,$second;
print OUTF "<BR>\n";
print OUTF "<B>Semi-major axis (km):</B>\n";
printf OUTF "<INPUT NAME=\"orb_a\" SIZE=15 value=\"%7.3f\"> \n",$semiax;
print OUTF "<BR>\n";
print OUTF "<B>Right ascension of the ascending node (degrees):</B>\n";
print OUTF "<INPUT NAME=\"orb_raan\" SIZE=15 value=\"$raan\">\n";
print OUTF "<BR>\n";
print OUTF "<B>Orbit inclination (degrees):</B>\n";
print OUTF "<INPUT NAME=\"orb_i\" SIZE=15 value=\"$incl\">\n";
print OUTF "<BR>\n";
print OUTF "<B>Eccentricity:</B>\n";
print OUTF "<INPUT NAME=\"orb_e\" SIZE=15 value=\"$eccen\"> \n";
print OUTF "</BLOCKQUOTE>\n";
print OUTF "<BR>\n";
print OUTF "<INPUT TYPE=\"submit\" VALUE=\"Calculate CVZ visibilities\">\n";
print OUTF "<INPUT TYPE=\"reset\" VALUE=\"Erase the input fields\">\n";
print OUTF "</FORM>\n";

print OUTF "<HR>\n";
print OUTF "This page uses the <a href=\"http://star-www.rl.ac.uk/\"> STARLINK </a>\n";
print OUTF "set of astronomical subroutines.\n";
print OUTF "<BR>\n";
print OUTF "This page was automatically generated on $now";
print OUTF "<BR>\n";
print OUTF "<a href=\"mailto:emurphy\@pha.jhu.edu\"><address>emurphy\@pha.jhu.edu</address></a>\n";

close (OUTF);

### end of Perl script


sub convert_doy {
    @daytab1 = (0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
    @daytab2 = (0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31);
    @monthstr = ("", "January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December");

    if ((($year%4 == 0) && ($year%100 != 0)) || ($year%400 == 0)) {
        # $year is a leap year use daytab1
        for ($i = 1; $doy > $daytab1[$i]; $i++) {
	    $doy -= $daytab1[$i];
	}
    } else {
        # $year is not a leap year, use daytab2
        for ($i = 1; $doy > $daytab2[$i]; $i++) {
	    $doy -= $daytab2[$i];
	}
    }
    $month = $monthstr[$i];
    $day = $doy;
}

sub convert_dayfrac {
    $dayfrac *= 24.0;
    $hour = int($dayfrac);
    $minute = int(($dayfrac-$hour)*60.0);
    $second = ((($dayfrac-$hour)*60.0)-$minute)*60.0;
}

sub calculate_a0 {

$mu = 3.986005E5;
$pi = 3.14159265358979;
$radian = 0.0174532925200;

$mm=2.0*$pi*$mean_mot/(24.0*60.0);
$ke=0.74366916E-1;
$k2=5.413080E-4;
$a1=($ke/$mm)**(2.0/3.0);
$inrad = $incl * $radian;

$d1=3.0/2.0*$k2/($a1*$a1)*(3.0*cos($inrad)*cos($inrad)-1.0)/((1.0-$eccen*$eccen)**(3.0/2.0));
$a0=$a1*(1.0-$d1/3.0-$d1*$d1-134.0/81.0*$d1*$d1*$d1);
$d0=3.0/2.0*$k2/($a0*$a0)*(3.0*cos($inrad)*cos($inrad)-1.0)/((1.0-$eccen*$eccen)**(3.0/2.0));
$semiax=$a0/(1-$d0)*6378.135;
}