astro

#!/usr/bin/perl -w
# Calculate Sun and Moon position, calculate EME conditions.
# License: GPLv3

# A portion is derived from the WSJT/WSJT-X project.
#
# Copyright (C) 2001 - 2017 by Joe Taylor, K1JT.
# Copyright (c) 2018, 2019, 2020, 2022 molo1134@github

use strict;
use utf8;
use feature 'unicode_strings';
use Astro::Coord::ECI;
use Astro::Coord::ECI::Moon;
use Astro::Coord::ECI::Sun;
#use Astro::Coord::ECI::Utils qw{deg2rad rad2deg};
use URI::Escape;
use Math::Round;
use Math::Trig;
use POSIX qw(floor strftime);
binmode(STDOUT, ":utf8");

use File::Basename;
use Cwd 'realpath';
use lib dirname(realpath(__FILE__));
use Colors;
use Location;
use Util;

my $username = $ENV{'USER'} || $ENV{'USERNAME'} || getpwuid($<);

# change this to 0 to always exit with success
our $exitnonzeroonerror = 1;
$exitnonzeroonerror = 0 if $username eq getEggdropUID();

if ($#ARGV < 0 || length($ARGV[0]) == 0) {
  print "usage: $0 [--moon] [--eme] [--sun] [--graves] \\\n\t<grid>|<qth> [<elev>ft|<elev>m]\n";
  exit $exitnonzeroonerror;
}

my $doEme = 0;
my $doMoon = 0;
my $doSun = 0;
my $doQuarter = 0;
my $doGraves = 0;

my $now = time();
#$now = 1487653200; ## high degredation
#$now = 1488053200; ## new moon
#$now = 1515572994; ## GRAVES illuminated

my ($lat, $lon) = (undef, undef);

# eggdrop doesn't split args on spaces
@ARGV=split(' ', join(' ',@ARGV));
my $query = "";
my $elev_m = "0";
my $i = 0;
while ($i <= $#ARGV) {
  if ( $ARGV[$i] =~ /^([0-9]+)(ft|m)/ ) {
    if ($2 eq "m") {
      $elev_m = $1;
    } elsif ($2 eq "ft") {
      $elev_m = $1 * 0.3048;  # feet to meters
    }
  } elsif ($ARGV[$i] =~ /^--moon$/) {
    $doMoon = 1;
  } elsif ($ARGV[$i] =~ /^--sun$/) {
    $doSun = 1;
  } elsif ($ARGV[$i] =~ /^--(quarter|solstice|equinox)$/) {
    $doQuarter = 1;
  } elsif ($ARGV[$i] =~ /^--graves$/) {
    $doGraves = 1;
  } elsif ($ARGV[$i] =~ /^--eme$/) {
    $doMoon = 1;
    $doEme = 1;
  } elsif ($ARGV[$i] =~ /^--date$/) {
    # undocumented -- unix timestamp
    $now = $ARGV[++$i];
  } elsif ($ARGV[$i] =~ /^--geo$/) {
    ($lat, $lon) = split(/,/, $ARGV[++$i]);
  } else {
    $query .= " " unless $query eq "";
    $query .= $ARGV[$i];
  }
  $i++;
}
#print "query: $query\n";
#print "elev_m: $elev_m\n";

if ($doMoon == 0 && $doSun == 0 && $doGraves == 0 && $doQuarter == 0) {
  # no data was selected, so do all.
  $doMoon = 1;
  $doEme = 1;
  $doSun = 1;
  $doGraves = 1;
  $doQuarter = 1;
}

if ($query ne "") {
  my $result;
  $result = argToCoords($query) unless $query eq "" and $doGraves == 1;

  if (!defined $result and $query ne "") {
    print "error: \"$query\" not found\n";
    exit $exitnonzeroonerror;
  }
  ($lat, $lon) = split(',', $result) if defined $result;
}

#print "$lat, $lon\n";

if ($username eq getEggdropUID() and not defined $lat and not defined $lon and ($doMoon == 1 || $doSun == 1)) {
  print "usage: !eme <qth>|<grid>\n" if $doEme == 1;
  print "usage: !sun <qth>|<grid>\n" if $doSun == 1;
  print "usage: !moon <qth>|<grid>\n" if $doMoon == 1 && $doEme == 0;
  exit $exitnonzeroonerror;
}

if ((not defined $lat or not defined $lon) and ($doEme == 1 or $doSun == 1 or $doMoon == 1)) {
  print "error: no position specified\n";
  exit $exitnonzeroonerror;
}

# 1600 Pennsylvania Ave, Washington DC USA
# latitude 38.899 N, longitude 77.038 W,
# altitude 16.68 meters above sea level
#my $lat = deg2rad (38.899);    # Radians
#my $long = deg2rad (-77.038);  # Radians
#my $alt = 16.68 / 1000;        # Kilometers

my $moon = Astro::Coord::ECI::Moon->new();
my $sun = Astro::Coord::ECI::Sun->new();
my $sta = Astro::Coord::ECI->universal($now)->
  geodetic(deg2rad($lat), deg2rad($lon), $elev_m / 1000) if defined $lat and defined $lon;  #kilometers
$sta->set(refraction => 0) if defined $sta;
$moon->set(refraction => 0);
$sun->set(refraction => 0);
$sta->universal($now) if defined $sta;
$moon->universal($now);
$sun->universal($now);

if ($doMoon) {

  my ($moon_az, $moon_el, $moon_range) = getAzElRng($sta, $moon);

  my $earthradius = 6371;
  my ($minrange, $maxrange) = (356375 - $earthradius, 406720 - $earthradius); # perigee, apogee
  my $loss2m = emePathLoss($moon_range, 144.1);
  my $minloss2m = emePathLoss($minrange, 144.1); # loss at minimum perigee
  #my $maxloss2m = emePathLoss($maxrange - $earthradius, 144.1); # loss at maximum apogee
  my $percentRange = round(($moon_range - $minrange) / ($maxrange - $minrange) * 100);
  my $lossdiff = nearest(0.1, $loss2m - $minloss2m);

  my ($phase, $illum) = $moon->phase($now);
  $phase = rad2deg($phase);
  my ($phasetext, $phaseunicode) = phaseText($phase);
  $phasetext = red($phasetext)
      if ($doEme == 1 && $phasetext eq "new moon");

  $illum = round($illum * 100);

  if ($doMoon == 1) {
    if ($moon_el < 0) {
      print bold("Moon"), " is ", red("set"), ";";
    } else {
      print bold("Moon"), " az/el ", bold("$moon_az°/$moon_el°"), "; dist $moon_range km;";
    }
    print " phase: $phasetext $phaseunicode, $illum% illum.;";
  }

  # rise/set
  my ($eventtime, $rise) = $sta->next_elevation ($moon);
  my ($h, $m) = getTimeRemaining($now, $eventtime);
  print " Moon will @{[$rise ? 'rise' : 'set']} at ",
	strftime("%Y-%m-%d %H:%Mz", gmtime($eventtime)), " (in ${h}h${m}m);"
	if $doMoon == 1;

  my ($moon_ra, $moon_dec, $moon_ra_range) = $moon->equatorial();
  $moon_ra = rad2deg($moon_ra);
  my $moon_h = ($moon_ra / 360) * 24;
  my $moon_m = ($moon_h - floor($moon_h)) * 60;
  my $moon_s = ($moon_m - floor($moon_m)) * 60;
  $moon_h = floor($moon_h);
  $moon_m = floor($moon_m);
  $moon_s = round($moon_s);
  $moon_dec = rad2deg($moon_dec);
  my $moon_dd = $moon_dec;
  my $moon_dm = ($moon_dd - floor($moon_dd)) * 60;
  my $moon_ds = ($moon_dm - floor($moon_dm)) * 60;
  $moon_dd = floor($moon_dd);
  $moon_dm = floor($moon_dm);
  $moon_ds = round($moon_ds);
  $moon_ra_range = round($moon_ra_range);

  #print "\nMoon ra, dec, range: $moon_ra, $moon_dec, $moon_ra_range\n";
  #print "\nMoon ra/dec: ${moon_h}h${moon_m}m${moon_s}s  $moon_dd°$moon_dm'$moon_ds\"  range $moon_ra_range km\n";

  my $degd = degredation($moon_ra, $moon_dec, 144.100);
  print " EME degrd(2m): $degd dB (" . degredationQuality($degd) . ");"
    if $doEme == 1;
    print " EME path loss(2m): $loss2m dB (", lossDiffColor($lossdiff), ", $percentRange%);" if $doEme == 1 and $moon_el >= 0;
  print "\n" if $doMoon == 1;
}

if ($doSun == 1) {
  my ($sun_az, $sun_el, $sun_range) = getAzElRng($sta, $sun);
  #my ($az, $el, $range) = $sta->azel($sun);
  if ($sun_el < 0) {
    print bold("Sun"), " is set;";
  } else {
    print bold("Sun"), " az/el ", bold("$sun_az°/$sun_el°"), ";";
  }
  my ($time, $rise) = $sta->next_elevation ($sun);
  my $h = floor(($time - $now) / (3600));
  my $m = floor(((($time - $now) / 3600) - $h) * 60);
  print " Sun will @{[$rise ? 'rise' : 'set']} at ",
	strftime("%Y-%m-%d %H:%Mz", gmtime($time)), " (in ${h}h${m}m)\n";
}

if ($doQuarter == 1) {
  my %qtrDesc = (
    0 => "March equinox",
    1 => "June solstice",
    2 => "September equinox",
    3 => "December solstice" );

  my ($time, $qtr, $desc) = $sun->next_quarter();
  $desc = $qtrDesc{$qtr};

  my $dur = int($time) - $now;
  my $m = floor($dur / (30 * 24 * 60 * 60));
  my $d = floor($dur % (30 * 24 * 60 * 60) / (24 * 60 * 60));
  my $h = floor($dur % (24 * 60 * 60) / (60 * 60));
  my $min = ($dur / 60) % 60;
  my $sec = $dur % 60;

  print "Time until $desc: ";
  print "$m month, " if $m == 1;
  print "$m months, " if $m > 1;
  print "$d day, " if $d == 1;
  print "$d days, " if $d > 1;
  print "$h hours, $min min, $sec sec\n";
}

# GRAVES radar -- see https://en.wikipedia.org/wiki/Graves_(system)
# see http://www.itr-datanet.com/~pe1itr/pdf/The%20143.050MHz%20Graves%20Radar%20a%20VHF%20Beacon.pdf
# illuminates southern sky in France
if ($doGraves == 1) {
  $lat = 47.3480;
  $lon = 4.5151;
  $elev_m = 204;
  my $graves = Astro::Coord::ECI->
    universal ($now)->
    geodetic (deg2rad($lat), deg2rad($lon), $elev_m / 1000);  #kilometers
  $graves->set(refraction => 0);
  $graves->universal($now);
  $moon->universal($now);
  my ($graves_az, $graves_el, $graves_range) = getAzElRng($graves, $moon);
  my $loss = emePathLoss($graves_range, 143.05);
  print bold("GRAVES"), " Radar: Moon ";
  if ($graves_az >= 90 && $graves_az <= 270 && $graves_el >= 15 && $graves_el <= 40) {
    print green("Illuminated"), ", 143.050 MHz: az/el $graves_az°/$graves_el°; path loss $loss dB\n";
  } else {
    if ($graves_el < 0) {
      print "is ", red("set"), "\n";
    } else {
      print red("outside of beam"), "; az/el $graves_az°/$graves_el°\n";
    }
  }
}

exit 0;


sub log10 {
  my $n = shift;
  return log($n)/log(10);
}

# emePathLoss( dist_km, freq_MHz )
# from wikipedia:
#   https://en.wikipedia.org/w/index.php?title=Earth%E2%80%93Moon%E2%80%93Earth_communication&oldid=759977305
sub emePathLoss {
  my $d = shift;
  my $f = shift;
  my $loss = ((20 * log10($f)) + (40 * log10($d)) - 14.49);
  return nearest(0.1, $loss);
}

# From wsjt astro.f90 -- translated from Fortran
sub degredation {
  my $RAMoon = shift;
  my $DecMoon = shift;
  my $freq = shift;

  # data nt144/                                             &
  #        234, 246, 257, 267, 275, 280, 283, 286, 291, 298,  &
  #        305, 313, 322, 331, 341, 351, 361, 369, 376, 381,  &
  #        383, 382, 379, 374, 370, 366, 363, 361, 363, 368,  &
  #        376, 388, 401, 415, 428, 440, 453, 467, 487, 512,  &
  #        544, 579, 607, 618, 609, 588, 563, 539, 512, 482,  &
  #        450, 422, 398, 379, 363, 349, 334, 319, 302, 282,  &
  #        262, 242, 226, 213, 205, 200, 198, 197, 196, 197,  &
  #        200, 202, 204, 205, 204, 203, 202, 201, 203, 206,  &
  #        212, 218, 223, 227, 231, 236, 240, 243, 247, 257,  &
  #        276, 301, 324, 339, 346, 344, 339, 331, 323, 316,  &
  #        312, 310, 312, 317, 327, 341, 358, 375, 392, 407,  &
  #        422, 437, 451, 466, 480, 494, 511, 530, 552, 579,  &
  #        612, 653, 702, 768, 863,1008,1232,1557,1966,2385,  &
  #       2719,2924,3018,3038,2986,2836,2570,2213,1823,1461,  &
  #       1163, 939, 783, 677, 602, 543, 494, 452, 419, 392,  &
  #        373, 360, 353, 350, 350, 350, 350, 350, 350, 348,  &
  #        344, 337, 329, 319, 307, 295, 284, 276, 272, 272,  &
  #        273, 274, 274, 271, 266, 260, 252, 245, 238, 231/

  my @nt144 = ( 234, 246, 257, 267, 275, 280, 283, 286, 291, 298,
	        305, 313, 322, 331, 341, 351, 361, 369, 376, 381,
	        383, 382, 379, 374, 370, 366, 363, 361, 363, 368,
	        376, 388, 401, 415, 428, 440, 453, 467, 487, 512,
	        544, 579, 607, 618, 609, 588, 563, 539, 512, 482,
	        450, 422, 398, 379, 363, 349, 334, 319, 302, 282,
	        262, 242, 226, 213, 205, 200, 198, 197, 196, 197,
	        200, 202, 204, 205, 204, 203, 202, 201, 203, 206,
	        212, 218, 223, 227, 231, 236, 240, 243, 247, 257,
	        276, 301, 324, 339, 346, 344, 339, 331, 323, 316,
	        312, 310, 312, 317, 327, 341, 358, 375, 392, 407,
	        422, 437, 451, 466, 480, 494, 511, 530, 552, 579,
	        612, 653, 702, 768, 863,1008,1232,1557,1966,2385,
	       2719,2924,3018,3038,2986,2836,2570,2213,1823,1461,
	       1163, 939, 783, 677, 602, 543, 494, 452, 419, 392,
	        373, 360, 353, 350, 350, 350, 350, 350, 350, 348,
	        344, 337, 329, 319, 307, 295, 284, 276, 272, 272,
	        273, 274, 274, 271, 266, 260, 252, 245, 238, 231);

  # call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb)
  my ($el, $eb) = coord(0.0, 0.0, -1.570796, 1.161639,
			deg2rad($RAMoon), deg2rad($DecMoon));

  # longecl_half=nint(rad*el/2.0)
  my $longecl_half = round(rad2deg($el)/2.0);

  # if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180
  if ($longecl_half < 1.0 || $longecl_half > 180.0) {
    $longecl_half = 180;
  }

  # t144=nt144(longecl_half)
  my $t144 = $nt144[$longecl_half - 1];

  # tsky=(t144-2.7)*(144.0d6/freq8)**2.6 + 2.7      !Tsky for obs freq
  my $tsky = ($t144 - 2.7) * (144.0 / $freq) ** 2.6 + 2.7;

  # tr=80.0                              !Good preamp
  my $tr = 80.0;

  # tskymin=13.0*(408.0d6/freq8)**2.6      !Cold sky temperature
  my $tskymin = 13.0 * (408.0 / $freq) ** 2.6;

  # tsysmin=tskymin+tr
  my $tsysmin = $tskymin + $tr;

  # tsys=tsky+tr
  my $tsys = $tsky + $tr;

  # dgrd=-10.0*log10(tsys/tsysmin) + dbMoon
  my $dgrd = 10.0 * log10($tsys/$tsysmin); # XXX dbMoon?

  return nearest(0.1, $dgrd);
}

sub degredationQuality {
  my $degd = shift;
  if ($degd <= 1.5) {
    return green("Excellent");
  }
  if ($degd > 1.5 && $degd <= 2.5) {
    return green("Good");
  }
  if ($degd > 2.5 && $degd <= 4.0) {
    return yellow("Fair");
  }
  if ($degd > 4.0 && $degd <= 5.5) {
    return red("Poor");
  }
  if ($degd > 5.5) {
    return darkRed("Very Poor");
  }
}

# from wsjt coord.f90 -- translated from Fortran
sub coord {
  my $a0 = shift;
  my $b0 = shift;
  my $ap = shift;
  my $bp = shift;
  my $a1 = shift;
  my $b1 = shift;
  my ($a2, $b2);

  # SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2)
  #
  # ! Examples:
  # ! 1. From ha,dec to az,el:
  # !      call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
  # ! 2. From az,el to ha,dec:
  # !      call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
  # ! 3. From ra,dec to l,b
  # !      call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
  # !        ra,dec,l,b)
  # ! 4. From l,b to ra,dec
  # !      call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
  # !        0.478220215d0,l,b,ra,dec)
  # ! 5. From ra,dec to ecliptic latitude (eb) and longitude (el):
  # !      call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
  # ! 6. From ecliptic latitude (eb) and longitude (el) to ra,dec:
  # !      call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec)
  #
  #   SB0=sin(B0)
  my $sb0 = sin($b0);

  #   CB0=cos(B0)
  my $cb0 = cos($b0);

  #   SBP=sin(BP)
  my $sbp = sin($bp);

  #   CBP=cos(BP)
  my $cbp = cos($bp);

  #   SB1=sin(B1)
  my $sb1 = sin($b1);

  #   CB1=cos(B1)
  my $cb1 = cos($b1);

  #   SB2=SBP*SB1 + CBP*CB1*cos(AP-A1)
  my $sb2 = $sbp*$sb1 + $cbp*$cb1*cos($ap-$a1);

  #   CB2=SQRT(1.e0-SB2**2)
  my $cb2 = sqrt(1.0 - $sb2**2);

  #   B2=atan(SB2/CB2)
  $b2 = atan($sb2/$cb2);

  #   SAA=sin(AP-A1)*CB1/CB2
  my $saa = sin($ap-$a1)*$cb1/$cb2;

  #   CAA=(SB1-SB2*SBP)/(CB2*CBP)
  my $caa = ($sb1 - $sb2 * $sbp)/($cb2*$cbp);

  #   CBB=SB0/CBP
  my $cbb = $sb0/$cbp;

  #   SBB=sin(AP-A0)*CB0
  my $sbb = sin($ap-$a0)*$cb0;

  #   SA2=SAA*CBB-CAA*SBB
  my $sa2 = $saa*$cbb-$caa*$sbb;

  #   CA2=CAA*CBB+SAA*SBB
  my $ca2 = $caa * $cbb + $saa * $sbb;

  #   TA2O2=0.0 !Shut up compiler warnings. -db
  my $ta2o2 = 0.0;

  #   IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2
  if ($ca2 <= 0.0) {
    $ta2o2=(1.0-$ca2)/$sa2;
  }

  #   IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2)
  if ($ca2 > 0.0) {
    $ta2o2 = $sa2 / (1.0 + $ca2);
  }

  #   A2=2.e0*atan(TA2O2)
  $a2 = 2.0 * atan($ta2o2);
  #   IF(A2.LT.0.e0) A2=A2+6.2831853
  if ($a2 < 0.0) {
    $a2 = $a2 + 6.2831853
  }

  #   RETURN
  # END SUBROUTINE COORD

  return ($a2, $b2);
}

sub getAzElRng {
  my $sta = shift;
  my $obj = shift;
  my ($az, $el, $rng) = $sta->azel($obj);
  $az = nearest(0.1, rad2deg($az));
  $el = nearest(0.1, rad2deg($el));
  $rng = round($rng);
  return ($az, $el, $rng);
}

sub getTimeRemaining {
  my $now = shift;
  my $event = shift;
  my $h = floor(($event - $now) / (3600));
  my $m = floor(((($event - $now) / 3600) - $h) * 60);
  return ($h, $m);
}

sub phaseText {
  my $phase = shift;  # phase in degrees; 0 = new moon; 180 = full moon
  my $text = undef;
  my $unicode = undef;

  if ($phase < 22.5 || $phase > 337.5) {
    $text = "new moon";
    $unicode = "🌑";
  } elsif (22.5 <= $phase && $phase < 67.5) {
    $text = "wax. cres.";
    $unicode = "🌒";
  } elsif (67.5 <= $phase && $phase < 112.5) {
    $text = "1st qrtr.";
    $unicode = "🌓";
  } elsif (112.5 <= $phase && $phase < 157.5) {
    $text = "wax. gibb.";
    $unicode = "🌔";
    #$unicode = "\xF0\x9F\x8C\x94"; # 🌔  -- U+1F314 is outside of tcl range
  } elsif (157.5 <= $phase && $phase < 202.5) {
    $text = "full moon";
    $unicode = "🌕";
  } elsif (202.5 <= $phase && $phase < 247.5) {
    $text = "wan. gibb.";
    $unicode = "🌖";
  } elsif (247.5 <= $phase && $phase < 292.5) {
    $text = "3rd qrtr.";
    $unicode = "🌗";
  } elsif (292.5 <= $phase && $phase < 337.5) {
    $text = "wan. cres.";
    $unicode = "🌘";
  }
  return ($text, $unicode);
}

sub lossDiffColor {
  my $lossdiff = shift;
  my $text = "+$lossdiff dB";
  if ($lossdiff <= 1) {
    return green($text);
  } elsif ($lossdiff <= 2) {
    return yellow($text);
  } else {
    return red($text);
  }
}