#!/usr/bin/perl -w

use DBI;
use Getopt::Long;
use Term::ReadKey;
use Geo::Shapelib qw/:all/;
use Math::Trig;

$dbase    = "HiRISE";
$user     = "";
$pass     = "";
$host     = "pirldb.lpl.arizona.edu";
$ogrph    = 0;
$outrt    = "";
$help     = 0;

$opt = GetOptions ("dbase|d=s"         => \$dbase,
                   "host=s"            => \$host,
                   "user|u=s"          => \$user,
                   "pass|p=s"          => \$pass,
		   "graphic|g"         => \$ogrph,
		   "output|o=s"        => \$outrt,
                   "help|h"            => \$help);

if (!$opt || $help) {
 print "\n";
 print "obs_geom.pl -dbase -host -user -pass -help -output \n";
 print "\n";
 print "  dbase|d    database to connect to (default is HiRISE)\n";
 print "  host       hostname where the database is stored (default is pirldb.lpl.arizona.edu)\n";
 print "  user|u     username to connect to database (not username associated with suggestions)\n";
 print "  pass|p     password to connect to database (program will prompt you if not provided here)\n";
 print "  graphic|g  Export shapefile in 'ographic latitudes\n";
 print "  output|o   Root name for output e.g. 'xx' will produce files xx_yyyymmdd.shp etc...\n";
 print "  help|h     print this help message\n";
 print "\n";
 print "This program produces an ESRI shapefile containing HiRISE planned observations.\n";
 print "It does this by reading the HiCAT database.";
 print "Questions? Contact: Shane Byrne - shane\@lpl.arizona.edu\n";
 die("\n");
}

if ($user eq '') {
 print "MySQL username: ";
 $user = <STDIN>;
 chomp $user;
}

if ($pass eq '') {
 print "MySQL password: ";
 ReadMode('noecho');
 $pass = ReadLine(0);
 chomp $pass;
 ReadMode('normal');
}

if ($outrt eq '') {$outrt = "plan_geom"};

$dbh = DBI->connect("dbi:mysql:database=".$dbase.";host=".$host, $user, $pass);
$dbh->{RaiseError} = 1;
$dbh->{AutoCommit} = 1;

$sqlcmd1 = "SELECT
 Planned_Observations.OBSERVATION_ID
,Planned_Observations.CENTER_PLANETOCENTRIC_LATITUDE
,Planned_Observations.CENTER_LONGITUDE
,Planned_Observations.STATUS
,Planned_Observations.COMMENT
,Planned_Observations.ESTIMATED_EMISSION_ANGLE
,Planned_Observations.ESTIMATED_INCIDENCE_ANGLE
,Planned_Observations.ESTIMATED_PHASE_ANGLE
,Planned_Observations.PREDICT_TIME
,Planned_CCD_Parameters.BINNING
,Planned_CCD_Parameters.IMAGE_LINES
 FROM Planned_Observations
 LEFT JOIN 
 Observation_Geometry ON Planned_Observations.OBSERVATION_ID = Observation_Geometry.OBSERVATION_ID 
 LEFT JOIN 
 Planned_CCD_Parameters ON Planned_Observations.ID = Planned_CCD_Parameters.PLANNED_OBSERVATIONS_ID ";


$selection = " WHERE 
 Planned_Observations.STATUS != 'DEFUNCT' AND 
 Planned_Observations.TARGET_NAME = 'MARS' AND 
 Observation_Geometry.ID IS NULL AND
 Planned_CCD_Parameters.CCD_NAME = 'RED4' AND
 Planned_Observations.Observation_ID LIKE 'PSP%'";


@fldnames = qw/OBS_ID CEN_LAT CEN_LON STATUS COMMENT EMISSION INCIDENC PHASE  TIME   BINNING /;
@fldtypes = qw/String Double  Double  String String  Double   Double   Double String Double  /;

$ff   = ((3396190.0/3376200.0)*(3396190.0/3376200.0));    # Flattening to convert ocentric lats to ographic
$r2d  = 57.295779513082323;                               # Convert radians to degrees (180/pi)

($tm1, $tm2, $tm3) = (localtime)[3,4,5];
$tm2 = $tm2 + 1;
if (length($tm1) == 1) { $tm1 = "0".$tm1 };
if (length($tm2) == 1) { $tm2 = "0".$tm2 };
$nm = $outrt."_".($tm3+1900).$tm2.$tm1;

$shpfile = new Geo::Shapelib {
   Name => $nm,
   Shapetype => POLYGON,
   FieldNames => [@fldnames],
   FieldTypes => [@fldtypes]
};

print "Performing MySQL queries...\n";
@res = @{  $dbh->selectall_arrayref($sqlcmd1.$selection.";") };


print "Processing records...\n";
if ($ogrph) {print "Converting latitudes to planetographic\n"};

#############################
## Setup various parameters #
#############################

$incc = 93.0;                 # Inclination of MRO's orbit     -- CHECK THE EXACT VALUE OF THIS NUMBER --
$dpi  = 3.141592653589793116; # Double prec. PI
$d2r  = $dpi/180.0;           # Degrees to radians
$r2d  = 180.0/$dpi;           # Radians to degrees
$inc  = $incc*$d2r;           # Convert MRO inclination to radians
$ra   = 3396.1900;            # Mars equatorial radius
$rc   = 3376.2000;            # Mars polar radius

foreach $j (@res) { 
 @res2 = @{$j};

 $himg = $res2[9]*$res2[10]*0.001*(0.25 + 0.1*($res2[1]+90.0)/180.0);    # Height of a HiRISE footprint
 $wimg = 5.0 + 1.0*($res2[1]+90.0)/180.0;                                # Width  of a HiRISE footprint

 $clat = $res2[1]*$d2r; # Example image center in radians
 $clon = $res2[2]*$d2r; # Example image center in radians


$inc  = $incc*$d2r;           # Convert MRO inclination to radians
if (abs($clat) > ($dpi-$inc)) { $inc = $dpi-abs($clat) };

$rmars = ($ra*$rc) / sqrt( $rc**2*cos($clat)**2 + $ra**2*sin($clat)**2 ); # Mars radius at center latitude
$ah    = $himg/(2.0*$rmars);  # Angular half-height of HiRISE image 
$aw    = $wimg/(2.0*$rmars);  # Angular half-width  of HiRISE image 

#######################
# Do the calculations #
#######################

$x2[0]    = asin( sin($clat)/ sin($dpi-$inc) ) - $ah;
$x2[1]    = asin( sin($clat)/ sin($dpi-$inc) ) + $ah; # Angular position of image end points along MRO ground track

# lon2,lat2 are center coordinates of image ends behind and ahead of the image center respectively  
$lat2[0]  = asin(sin($dpi-$inc)*sin($x2[0])),
$lat2[1]  = asin(sin($dpi-$inc)*sin($x2[1]));
$lon2[0]  = $clon + acos( ( cos($ah) - sin($clat)*sin($lat2[0])) / (cos($clat)*cos($lat2[0])) );
$lon2[1]  = $clon - acos( ( cos($ah) - sin($clat)*sin($lat2[1])) / (cos($clat)*cos($lat2[1])) );

# lon3,lat3 are coordinates of the eastern image corners (sample zero)
$lat3[0]  = asin(  cos($aw)*sin($lat2[0]) + sin($aw)*cos($lat2[0])*sqrt(1.0-(tan($lat2[0])/tan($x2[0]))**2) );
$lat3[1]  = asin(  cos($aw)*sin($lat2[1]) + sin($aw)*cos($lat2[1])*sqrt(1.0-(tan($lat2[1])/tan($x2[1]))**2) );
$lon3[0]  = $lon2[0] + acos( (cos($aw) - sin($lat2[0])*sin($lat3[0])) / (cos($lat2[0])*cos($lat3[0])) );
$lon3[1]  = $lon2[1] + acos( (cos($aw) - sin($lat2[1])*sin($lat3[1])) / (cos($lat2[1])*cos($lat3[1])) );


# lon4,lat4 are coordinates of the western image corners (sample 20K)
$lat4[0]  = asin(  cos($aw)*sin($lat2[0]) - sin($aw)*cos($lat2[0])*sqrt(1.0-(tan($lat2[0])/tan($x2[0]))**2) );
$lat4[1]  = asin(  cos($aw)*sin($lat2[1]) - sin($aw)*cos($lat2[1])*sqrt(1.0-(tan($lat2[1])/tan($x2[1]))**2) );
$lon4[0]  = $lon2[0] - acos( (cos($aw) - sin($lat2[0])*sin($lat4[0])) / (cos($lat2[0])*cos($lat4[0])) );
$lon4[1]  = $lon2[1] - acos( (cos($aw) - sin($lat2[1])*sin($lat4[1])) / (cos($lat2[1])*cos($lat4[1])) );


# Compile these coordinates into single variables - make sure longitudes are in the 0-360 range
@cnr_lat = ($lat3[0],$lat3[1],$lat4[1],$lat4[0]);
@cnr_lon = ($lon3[0],$lon3[1],$lon4[1],$lon4[0]);

@nroi = ($cnr_lon[0]*$r2d, $cnr_lat[0]*$r2d);
for ($i = 1; $i <= 3; $i++) { push @nroi,($cnr_lon[$i]*$r2d, $cnr_lat[$i]*$r2d) };


  $ele  = ($#nroi + 1)/2.0;

  if ($ogrph) {
   for ($i = 1; $i <= $ele; $i++) { $nroi[2.0*($i-1)+1] = atan2(($ff*tan($nroi[2.0*($i-1)+1]/$r2d)),1.0)*$r2d };
  }
  
  for ($i = 1; $i <= $ele; $i++) { 
   if (($nroi[2.0*($i-1)] - 180.0) > 0) { 
    $nroi[2.0*($i-1)] = $nroi[2.0*($i-1)]-360.0;
   }
  }
  
  if ($nroi[1] - $nroi[7]) { 
   @newroi = ( $nroi[2.0*($ele-1)], $nroi[2.0*($ele-1)+1] );
   for ($i = $ele-1; $i > 0; $i--) { push @newroi,( $nroi[2.0*($i-1)], $nroi[2.0*($i-1)+1] ) };  
   @nroi = @newroi;
  }


  @vert=[( $nroi[0], $nroi[1], 0.0, 0.0)];
  for ($i = 2; $i <= $ele; $i++) { push @vert,[( $nroi[2.0*($i-1)], $nroi[2.0*($i-1)+1], 0.0, 0.0)] };
  push @vert,[( $nroi[0], $nroi[1], 0.0, 0.0)];

  foreach $i (@res2) { 
   if ((defined $i) eq "") { $i  = 0.0 };
  }  

  push @{$shpfile->{ShapeRecords}}, [@res2[0..9]];
  push @{$shpfile->{Shapes}}, {Vertices=>[@vert]};
# } # End condition of abs($clat) being less thane 180-$inc 

} # End loop over features 

$fout = open PRJ, "> $nm.prj";
print PRJ 'GEOGCS["GCS_Mars_2000",DATUM["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.8944472236118]],PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]]\n';

print "Closing files and database connections...\n";
if ($fout) {close PRJ};
$dbh->disconnect;
$shpfile->save();