#!/usr/bin/perl -w

# Generates a rmsun11 script which can in theory remove solar interference from
# the shortest 11 spacings, assuming that remove_confsrc is in your path.
# You may want to edit and chmod +x rmsun11 before running it.
#
# Also uses <a href="../misc/Astangles.pm">Astangles.pm</a>.

use Astangles qw(dec2hms);
use DRAORed qw(survey);

####################### "main()" #######################

my $srvy = survey();

open(RMSUNS, "> rmsun11") or die "Could not open rmsun11 for writing.\n";
print RMSUNS "#!/bin/tcsh\n\n";
print RMSUNS "foreach p ( q u )\n";
print RMSUNS "  cd \$p\n";
print RMSUNS "  foreach b ( a b c d )\n";
print RMSUNS "    cd \$b\n";

open(RDTBL, "${srvy}.TBL_DAT");
while (<RDTBL>){
  if(/^0*([0-9]+)==============================/){
    $ph2dn = $1;
  }
  elsif(/^JULDATE : ([0-9.]+)/){
    $jdate = $1 + 0.25;
  }
  elsif(/^SPACINGS: *([0-9]+)/){
    $spacing = $1
  }
  elsif(/^SPCTR /){   # This is just a convenient marker for the end of a day.
    if($spacing < 12){
      my ($sunra, $sundec) = solar_position($jdate);

      # It may work a bit better with the Sun's image reflected by the equator.
      # Waiting for testing.
      # if($sundec < 0.0){
      #   $sundec = -$sundec;
      # }
      print RMSUNS "    remove_confsrc -l \"$sunra, $sundec\" -n $ph2dn ";
      print RMSUNS "-i $spacing -a $spacing --xy 1.5 SUN${spacing}_D$ph2dn\n"
    }
  }
}
close(RDTBL) or die "Error reading ${srvy}.TBL_DAT\n";

print RMSUNS "    cd ..\n";
print RMSUNS "  end\n";
print RMSUNS "  cd ..\n";
print RMSUNS "end\n";
close(RMSUNS) or die "Error writing rmsun11\n";

#############################################################

sub solar_position{
  my $jdate = $_[0];

  my $ddate = int $jdate;  # int truncates, and jdate > 0, so ddate is the
                           # start of the day that includes jdate.
  my $udate = $ddate + 1;

  # Look up the Sun's position for both $ddate and $udate.
  open(SJDD, "/home/rreid/gpsr/polscripts/sunjddir.dat") or
    die "Could not open ~/gpsr/polscripts/sunjddir.dat\n";
  while(<SJDD>){
    if(/^(${ddate}\.[0-9]+) +([^ ]+) \+?(.*)/){
      $sddate = $1;
      ($dradd, $ddecdd) = dec2hms($2, $3);
    }
    elsif(/^(${udate}\.[0-9]+) +([^ ]+) \+?(.*)/){
      $sudate = $1;
      ($uradd, $udecdd) = dec2hms($2, $3);
    }
  }
  close(SJDD) or die "Error closing ~/gpsr/polscripts/sunjddir.dat\n";

  # Interpolate to $jdate.
  my ($radd, $decdd) = ($dradd, $ddecdd);
  my $howfar = ($jdate - $sddate) / ($sudate - $sddate);

  $radd  += $howfar * ($uradd - $dradd);
  $decdd += $howfar * ($udecdd - $ddecdd);

  return dec2hms($radd, $decdd);
}

# Copyleft (C) 2004 Rob Reid