#! /usr/bin/perl -w
# This program displays data from the WDS.

BEGIN { push @INC, "../" }

use astroUtils;
use CGI qw(:standard); # The Perl CGI processing class.
use DBI;
use Math::Trig;
use strict;

open CT, "zTestCount" or die "Can't open zTestCount.\n";
my $ct = readline CT;
close CT;
chomp $ct;
$ct++;
open CT, ">zTestCount" or die "Can't write zTestCount.\n";
print CT "$ct\n";
close CT;

# open T1, ">TestFile" or die "Can't write TestFile.\n"; #TEST

my $pi = 3.1415926535897932384626433;
my $d2r = $pi / 180;
my $h2r = $pi / 12;
my $tpi = $pi * 2;

my $east = param('east');
my $west = param('west');
my $north = param('north');
my $south = param('south');

my $err = "";

unless($east =~ /\d/) { $err = "<H2>Eastern boundary not defined!</H2>"; }
unless($west =~ /\d/) { $err .= "<H2>Western boundary not defined!</H2>"; }
unless($north =~ /\d/) { $err .= "<H2>Northern boundary not defined!</H2>"; }
unless($south =~ /\d/) { $err .= "<H2>Southern boundary not defined!</H2>"; }

$north = ckVal($north);
$south = ckVal($south);
$east = ckVal($east);
$west = ckVal($west);

if ($err) {
  print header, start_html("Undefined boundaries");
  print $err;
  print "<P>Click your back button to continue.</P>";
  print "</BODY></HTML>\n";
  exit;
}

my $eRa = getDeg($east);
$eRa *= $h2r;
if (($eRa < 0) or ($eRa > $tpi)) {
  $err .= "<P>$east. 0 <= Eastern boundary <= 24</P>";
}

my $wRa = getDeg($west);
$wRa *= $h2r;
if (($wRa < 0) or ($wRa > $tpi)) {
  $err .= "<P>$west. 0 <= Western boundary <= 24</P>";
}

my $nDec = getDeg($north);
$nDec *= $d2r;
if (($nDec < (-$pi / 2)) or ($nDec > ($pi / 2))) {
  $err .= "<P>$north. -90 <= Northern boundary <= 90</P>";
}

my $sDec = getDeg($south);

$sDec *= $d2r;
if (($sDec < (-$pi / 2)) or ($sDec > ($pi / 2))) {
  $err .= "<P>$south. -90 <= Southern boundary <= 90</P>";
}

if (($nDec - $sDec) < 0) { 
  $err .= "<P>The southen boundary is north of the northern boundary.</P>";
}

$ct = 0;
my $maxCt = 500; # No more matches than this.
my $rCt = 0;     # Catalog list row count.

my $minMv  = param('mv');
my $deltaM = param('deltaM');
my $minSep = param('minSep');
my $maxSep = param('maxSep');

if (($minMv < 0) or ($minMv > 19)) {
  $err .= "<P>$minMv. 1 <= mv <= 19.</P>";
}

if (($deltaM < 0) or ($deltaM > 99)) {
  $err .= "<P>$deltaM. 0 <= magnitude difference <= 99.</P>";
}

if (($minSep < 0) or ($minSep > 3600)) {
  $err .= "<P>$minSep. 0 <= minimum separation <= 3600.</P>";
}

if (($maxSep < 0) or ($maxSep > 3600)) {
  $err .= "<P>$maxSep. 0 <= maximum separation <= 3600.</P>";
}

if (($maxSep - $minSep) < 0) {
  $err .= "<P>Minimum separation must be less than maximum separation.</P>";
}

if ($err) {
  print header, start_html("Out of range");
  print $err;
  print "<P>Click your back button to continue.</P>";
  print "</BODY></HTML>\n";
  exit;
}

my $list = "";
my $wdsNotes = "<BR></CENTER></P>\n"; # Text of the WDS note count.
my $wdsNtCt = 0;   # Count of the number of WDS notes in the list.

# Data from the WDS and other USNO files:
open WDS, "data4WDS_Search.txt" or die "Can't open data4WDS_Search.txt.\n";

if (($eRa - $wRa) < 0) {
  #Zero hours is in the window.
  search(0, $eRa);
  search($wRa, $tpi);
} else {
  search($wRa, $eRa);
}

$list = "<BODY BGCOLOR=navy TEXT=white LINK=cyan VLINK=#00AA00><CENTER>\n" .
"<H3>Binaries in RA: $west to $east, dec: $north to $south.</H3>\n<P><B>" .
"<H4>Min Primary mv: $minMv, Max Delta mv: $deltaM, " .
"Separation: $minSep\" - $maxSep\".</H3>\n" .  
"<A HREF=\"newExp.html\" target = \"_blank\">" .
"Explanation</A></B></P>\n<TABLE BORDER=8><TR>" .
"<TD><center>J2000 coordinates</center></TD>" .    
"<TD><center>WDS Id</center></TD>" .
"<TD><center>Dsc Id</center></TD>" .    
"<TD><center>A mv</center></TD>" .    
"<TD><center>B mv</center></TD>" .    
"<TD><center>Spectrum</center></TD>" .    
"<TD><center>first &theta;</center></TD>" .    
"<TD><center>last &theta;</center></TD>" .    
"<TD><center>first &rho;</center></TD>" .    
"<TD><center>last &rho;</center></TD>" .    
"<TD><center>First Obs</center></TD>" .    
"<TD><center>Last Obs</center></TD>" .    
"<TD><center>AKA</center></TD>" .    
"<TD><center>Comments</center></TD>" .
"<TD><center>Notes</center></TD>\n" . $list;   

print header, start_html("WDS binary stars");
print $list;

if ($rCt > $maxCt) {
  print "<P>More than $maxCt pairs found. ",
        "You might want to make your area smaller.</P>";
} else {
  print "<P>There were $rCt pairs found.</P>";
}

print "</TABLE></BODY></HTML>\n";
print $wdsNotes;

# close T1; #TEST
close WDS;

# open TST, ">test.html" or die "Can't write test.html.\n"; #Debug
# print TST "$list.\n$wdsNotes"; #Debug
# close TST; #Debug

# Search for stars within the requested parameters.
# Each data entry must have all 5 fields present, even if they are null:
# The field delimeter is a vertical bar: "|".
# 0) Right ascension (J2000)
# 1) Declination (J2000)
# 2) WDS Id
# 3) Discoverer ID 
# 4) Apparent visual magnitude of primary.
# 5) Apparent visual magnitude of secondary.
# 6) Spectrum
# 7) First observed Position angle (theta) in degrees.
# 8) Last observed Position angle (theta) in degrees.
# 9) First observed separation, in arc seconds
# 10) First observed separation, in arc seconds
# 11) First observed date
# 12) Last observed date
# 13) Other names for the star.
# 14) N = Neglected, L = Low precision coordinates only. 
#     K = K band magnitude.
# 15) WDS notes.
#
# If an orbit is found, then:
# 16) Semi Major Axis
# 17) PeriastronDate
# 18) Eccentricity
# 19) Sky angle
# 20) Ascending node angle
# 21) Period
# 22) Longitude of periastron
# 23) Quality
# 24) Theta for previous year.
# 25) Rho for previous year.
# 26) Theta for current year.
# 27) Rho for current year.
# 28) Theta for next year.
# 29) Rho for next year.
# 30) Theta for year 2.
# 31) Rho for year 2.
# 32) Orbit's quality grade.
# 33) URLs for orbit's diagrams.
#

sub search {
  my $start = $_[0]; # Westernmost RA.
  my $seek = ($start / $tpi) * 9e6;
  $seek -= $seek / 5;
  $seek = int $seek;
  my $stop = $_[1];  # Easternmost RA.

  seek(WDS, $seek, 0);
  readline WDS;

  while(<WDS>) {
    chomp $_;
    my @f = split /\|/, $_;

    if ($f[0] < $start) { next; }
    if (($f[0] > $stop) and (($f[0] - $start) < 4)) { last; }

    if (($f[1] > $sDec) and ($f[1] < $nDec)) {
      if (($f[4] > $minMv) or
          (abs($f[5] - $f[4]) > $deltaM) or
          ($f[10] > $maxSep ) or
          ($f[10] < $minSep )) {
        next;
      }

      # Write format:
      # Ra Dec wdsId discoverer mva mvb spectrum firsObs last Obs
      # 0   1    2       3       4   5     6        7       8
      # first PA, lastPa firstSep lastSep AKA Comments Notes
      #   9        10      11      12     13    14      15

      my $row = "<TR>";

      my $color = "lime";
      if (($f[12] !~ /\d/) or ($f[12] < 1970)) {
        $color = "red";
      } elsif (($f[14] =~ /F/) and ($f[14] ne "None")) {
        $color = "yellow";
      }

      $row .= "<TD><FONT COLOR=$color><center>" . r2aladinRa($f[0]) . " " .
              r2aladinDec($f[1]) . "</center></FONT></TD>\n";

      # If there is an orbit associated with this star, get the true-of-date
      # values for rho ($f[10]) and theta ($f[8]).
      my $orb = 0; # Set if the star has an orbit.
      my ($rho, $theta); # True of date values.
      if ($#f > 20) {
        ($rho, $theta) = getTrueOfDate($_);
        $orb = 1;
      }

      else { next; } #TEST This will only show data for stars with known orbits.

      for (my $i = 2; $i < 13; $i++) {
        unless (defined($f[$i])) { $f[$i] = "-"; }

        if (($orb) and (($i == 8) or ($i == 10))) {
          if ($i == 8) {
            $row .= "<TD><FONT COLOR=white><center>$f[$i]<BR>$theta" .
                    "</center></FONT></TD>\n";
          } elsif($i == 10) {
            $row .= "<TD><FONT COLOR=white><center>$f[$i]<BR>$rho" .
                    "</center></FONT></TD>\n";
          }
        } else {
          $row .= "<TD><FONT COLOR=$color><center>$f[$i]</center></FONT></TD>\n";
        }
      }

      if ($f[13] eq "None") {
        $row .= "<TD><FONT COLOR=$color><center>-</center></FONT></TD>\n";
      } else {
        $row .= "<TD><FONT COLOR=$color><center>$f[13]</center></FONT></TD>\n";
      }

      if (($f[14] eq "None") and ($orb == 0)) {

        $row .= "<TD><FONT COLOR=$color><center>-</center></FONT></TD>\n";
      } else {
        my $c = "";
        if ($f[14] =~ "B") { $c .= "Blue mags<BR>"; } 
        if ($f[14] =~ "E") { $c .= "<P>Uncertian<BR>Position</P>"; } 
        if ($f[14] =~ "F") { $c .= "Neglected<BR>"; } 
        if ($f[14] =~ "K") { $c .= "IR mags<BR>"; } 
        if ($f[14] =~ "L") { $c .= "Linear<BR>"; } 
        if ($orb) {
          $c .= "<FONT COLOR=white>Last Obs<BR>ToD ($f[23])</FONT><BR>";

          if ((defined($f[33])) and ($f[33] =~ /wds/)) {
            # There are orbital diagrams associated with this one.
            my @dgm = split /\s+/, $f[33];
            for (my $i = 0; $i <= $#dgm; $i++) {
              $c .= "<A HREF=orbitDiagrams/$dgm[$i]>" . ($i + 1) . " </A>";
            }
          }
        }
        if ($f[14] =~ "R") { $c .= "Red mags<BR>"; } 
        if ($f[14] =~ /X/) { $c .= "Dubious<BR>"; }
        if ($f[14] =~ /I/) { $c .= "Uncertian<BR>"; }
        if (($f[14] =~ /S/) or ($f[14]  =~ /U/) or ($f[14]  =~ /Y/)) {
          $c .= "Optical<BR>";
        } 
        $c =~ s/<BR>$//;

        unless ($c) {
          $row .= "<TD><FONT COLOR=$color><center>-</center></FONT></TD>\n";
        } else {
          $row .= "<TD><FONT COLOR=$color><center>$c</center></FONT></TD>\n";
        }
      }

      if ($f[15] eq "None") {
       $row .= "<TD><FONT COLOR=silver><center>-</center></FONT></TD>\n";
      } else {
       $row .= "<TD><CENTER><FONT COLOR=cyan><A HREF=\"#nt$wdsNtCt\">" .
               "Note$wdsNtCt</A></FONT></CENTER></TD>\n";

       $wdsNotes .= "<P><A ID=\"nt" . $wdsNtCt . "\"><BR>" .
                    "<FONT COLOR=cyan>Note $wdsNtCt.</FONT></P>\n" .
                    "<P>$f[15]</P>\n";
       $wdsNtCt++;
      }

      $row .= "</TR>\n\n";
      $list .= $row;

      if ($rCt++ > $maxCt) { last; }
    }
  }
}

# Return the arctangent of $_[0].
sub atan1 {
  my $val = $_[0];
  my $hpi = 3.14159265358979323846 / 2;

  while ($val > $hpi) { $val -= $hpi; }
  while ($val < -$hpi) { $val += $hpi; }

  my $r = 1;
  my $result = $val;
  my $sign = -1;
  for (my $i = 1; $i < 7; $i++) {
    my $x *= $val * $val * $sign;
    $result += $x;
    $sign *= -1;
  }
  return $result;
}

# Ensure that the boundary value can be read by getDeg.
sub ckVal {
  my $val = $_[0];

  if ($val =~ /\d+:\d+:\d+/) {
    # It's OK, do nothing.
  } elsif ($val =~ /\d+:\d+/) {
    $val = $val . ":00";
  } elsif ($val =~ /\d+/) {
    $val = $val . ":00:00";
  }

  return $val;
}

# Calculate decimal degrees from the input.
sub getDeg {

  my $raw = $_[0];
  my $deg = $raw;

  my $sign = 1;
  if ($raw =~ /\s*-/) { $sign = -1; }

  $raw =~ /(\d+):(\d+):(\d+)/;
  if (defined($1)) { 
    $deg = (($1 + ($2 / 60) + ($3 / 3600)) * $sign); 
  } else {

    $raw =~ /(\d+):(\d+)/;
    if (defined($1)) {
      $deg = (($1 + ($2 / 60)) * $sign); 
    }
  }

  return($deg);
}

# Given the orbital elements of a my, calculate the true of date position
# angle (theta) and separation (rho) of the pair.
sub getTrueOfDate {
  my ($rho, $theta); # Separation and position angle this calculates.

  my @f = split /\|/, $_[0];
  my $a = $f[16]; # Semi Major Axis.
  my $d = $f[17]; # PeriastronDate.
  my $e = $f[18]; # Eccentricity.
  my $i = $f[19]; # Sky angle.
  my $n = $f[20]; # Ascending node angle.
  my $p = $f[21]; # Period in years.
  my $w = $f[22]; # Longitude of periastron.
  my $quality = $_[23]; # Quality.
  my $tc = $f[26]; # Theta for current year.
  my $rc = $f[27]; # Rho for current year.
  my $tn = $f[28]; # Theta for next year.
  my $rn = $f[29]; # Rho for next year.
  my $pi = 3.14159265358979323846;

  # What time is it?
  my ($sec, $min, $hour, $mday, $mon, $year, $wday, $diy, $isdst) = 
  localtime(time);
  my $yrFraction = $diy / 365.25;
  my $now = $year + 1900 +  $yrFraction;

  # Interpolate values for rho and theta.
  my $iTheta = (((($tc - $tn) * $yrFraction) / 365.25) + $tc) * ($pi / 180);
  my $iRho = ((($rc - $rn) * $yrFraction) / 365.25) + $rc;

  # Use the interpolated rho and theta from the ephemerides or calculate them?
  # Interpolated values should be fine if $yrFraction is under 10 percent
  # of the period.
  if ($yrFraction > 0.5) { $yrFraction = 1 - $yrFraction; }
  if (($yrFraction / $p) > 0.1) {
    my $dt = $now - $d;
    while ($dt > $p) { $dt -= $p; }

    my $ma = ($dt / $p) * (2 * $pi);  
    my $ee = solveKepler($e, $ma);
    my $a1 = sqrt((1 + $e) / (1 - $e)) * (sin($ee / 2) /  cos($ee / 2)); 
    my $nu = 2 * atan($a1);
    my $r = $a * (1 - ($e * cos($ee)));
    my $x = cos($nu + $w);
    my $y = sin($nu + $w) * cos($i);
    $theta = resolveAtan($y, $x);

    $theta += $n;
    if ($theta > 2 * $pi) { $theta -= 2 * $pi; } 

    $rho = ($r * cos($nu + $w)) / cos($theta - $n);

    # The method above, from "Practical Astronomy", Duffett-Smith, Zwart,
    # 4th edition, is often incorrect. Compare its values with those of
    # the USNO Ephemerides table.
    if ((abs($theta - $iTheta) / (2 * $pi)) > 0.05) { $theta = $iTheta; }
  } else {
    $rho = $iRho;
    $theta = $iTheta;
  }

  my $bkupRho = $rho;
  $rho = round2($rho);
  unless ($rho =~ /\d+\.\d{2}/) {
    my $fraction = $rho - int($rho);
    if ($fraction == 0) {
      $rho = $rho . ".00";
    } elsif (length($fraction) == 1) {
      print T1 "One digit Rho: $rho. "; #TEST
      $rho .= "0";
    } elsif (length($fraction) == 3) {
      $rho .= "0";
    } elsif ($rho == 0) {
      $rho = round3($bkupRho);
    } elsif (length($fraction) > 15) {
      if ($rho =~ /\d+\.\d{1}/) {
        $rho .= "0";
      } elsif ($rho =~ /\d+\.\d{2}0/) {
        $rho =~ s/0$//;
      }
    }
  }

  $theta = round($theta * 180 / $pi);
  return ($rho, $theta);
}

# Calculate the arc tangent of x/y, and then resolve the result into the
# correct positive angle, 0 <= a <= 2*pi.  Return the angle.
sub resolveAtan {
	my $y = $_[0];
  my $x = $_[1];

  my $pi = 3.14159265358979323846;
  my $th = atan2($y, $x);

  if ($y < 0) { $th += $pi; } 
  elsif (($x < 0) and ($y > 0)) { $th += 2 * $pi; } 

  return $th;
}

# Given eccentricity e and mean anomaly angle ma, iteratively solve for ee
# in Kepler's equation.
sub solveKepler {
  my $e = $_[0];
  my $ma = $_[1];

  my $epsilon = 5e-7;  # Once $ee is less than this, consider it solved.
  my $delta = 0.7;
  my $pi = 3.14159265358979323846;
    
  # Start with an approximation.
  if ($ma < $pi) {
    if ($e < 0.1) { $delta = 0; }
    elsif ($e < 0.2) { $delta = 0.15; }
    elsif ($e < 0.4) { $delta = 0.3; }
    elsif ($e < 0.6) { $delta = 0.4; }
    elsif ($e < 0.8) { $delta = 0.6; }
  } else {
    if ($e < 0.1) { $delta = 0; }
    elsif ($e < 0.2) { $delta = -0.15; }
    elsif ($e < 0.4) { $delta = -0.3; }
    elsif ($e < 0.6) { $delta = -0.4; }
    elsif ($e < 0.8) { $delta = -0.6; }
  }

  my $dlt = 99;
  my $ee = $ma + $delta;

  while ($dlt > $epsilon) {
    $dlt = $ee - ($e * sin($ee)) - $ma;
    if ($dlt > $epsilon) {
      $delta = $dlt / (1 - ($e * cos($ee))); 
      $ee = $ee - $delta;
    }
  }

  return $ee;
}
