#!/usr/local/bin/perl -w
use Getopt::Long;
use Astro::Vex;
use Astro::Time;
use Astro::Coord;
use POSIX;
use Carp;
use JSON;
#use IO::Tee;
#
# CR: 2025-05-28 Converted old vex2atca.pl to write the bigcat json schedule
# format.

$Astro::Time::StrSep = ':';
$Astro::Time::StrZero = 2;

use strict;

sub addmode ($$$$$$);
sub modesetup ($$$\%);
sub getintent ($);
sub slewtime ($$$$);
sub convertToDate(\$\$$$);
sub checkautophase();
sub write_atca_schedule($\%\@$) ;
sub get_receiver($);
sub select_cal_source($$) ;

use constant ATCA => 'At';

#my $cycle = 5; # seconds
#my $caltime = 30; # seconds
my $useintent = 1;

#GetOptions('cycle=i'=>\$cycle, 'cal=i'=>\$caltime, 'intent!'=>\$useintent);
GetOptions('intent!'=>\$useintent);

if (@ARGV!=1) {
  die "Usage vex2atca.pl <vexfile>\n";
}

#$caltime = int (ceil $caltime/$cycle);

my $vexfile = shift;

my $vex = new Astro::Vex($vexfile);

my @intents;
if ($useintent) {
  @intents = getintent($vexfile);
}

my $globalpointing = 0;
foreach (@intents) {
  foreach (@{$_}[1..$#$_]) {
    if ($_ =~ /REFERENCE[_ ]POINTING[ _]DETERMINE/i) {
      $globalpointing = 1;
      last;
    }
  }
  last if $globalpointing;
}

#my $atca_long = str2turn('149:33:53','D');
#my $atca_lat = str2turn('-30:18:46.4','D');

my %sources;
my %modes;


my ($oldfreq1, $oldfreq2, $nfreq);

my $exper = $vex->exper->exper_name;
my $PIname = $vex->exper->pi_name;
#$PIname =~ s/\s//g;

my @scans = $vex->ant_sched(ATCA);
die "No ATCA scans in schedule\n" if (@scans==0);

my ($last_stop, $last_src);
my $calmode;
my $calstart;

sub mjd_to_utc($) {
  my ($mjd) = shift;

  my ($day, $month, $year, $ut) = mjd2cal($mjd);
  my $utstr = turn2str($ut,'H',0);

  #return sprintf('%02d/%02d/%04d %s', $day, $month, $year, $utstr);
  return sprintf('%04d-%02d-%02dT%sZ', $year, $month, $day, $utstr);
}

my $first = 1;

my (@atcascans);
foreach (@scans) {
  # extract the required scan information and store
  my $start = $_->start;
  my $scanlen = $_->stations(ATCA)->datastop;
  my $stop = $start+$scanlen->unit('day')->value;

  my $source_def = $_->source;
  if ($source_def =~ /^\d+$/) {
    die "Do not support source names with digits only (".$source_def.")";
  }
  
  my $source = $vex->source($source_def);
  if (! defined $source) {
    die "Could not find source ".$_->source;
  }
  my $src = $source->source_name();
  if (! exists($sources{$src})) {
    $sources{$src} = {};
    $sources{$src}{RA} = $source->rastr;
    $sources{$src}{RA} =~ tr/hm/:/;
    $sources{$src}{RA} =~ s/s//;
    $sources{$src}{DEC} = $source->decstr;
    $sources{$src}{DEC} =~ tr/d\'/:/;
    $sources{$src}{DEC} =~ s/\"//;
    $sources{$src}{EPOCH} = $source->ref_coord_frame;
  }

  my $modename = $_->mode;

  # Groak this mode if we have not seen it before
  if (! exists $modes{$modename}) {
    modesetup($vex, ATCA, $modename, %modes);
  }
  
  if ($first) {
    $last_stop = $start ;
    $calstart = $start ;
  } 

  my $scanlength = ($stop-$start);

  my $autophase = 0;
  my $pointing = 0;
  my $paddle = 0;
  #my $delay = undef;
  my $scanid = $_->scanid;
  my @thisintent;
  my $found = 0;
  foreach (@intents) {
    if ($scanid eq $_->[0]) {
      @thisintent = @{$_}[1..$#$_];
      $found = 1;
      last;
    }
  }
  if (!$found) {
    warn "Could not find $scanid in list of intents\n";
  } else {
    foreach (@thisintent) {
      if ($_ =~ /AUTOPHASE[_ ]DETERMINE/i) {
        $autophase = 1;
      } elsif ($_ =~ /REFERENCE[_ ]POINTING[_ ]DETERMINE/i) {
        $pointing = 1;
      } elsif ($_ =~ /ATCA:PADDLE/i) {
        $paddle = 1;
      }
    }
  }
  if ($autophase) {
    if (not checkautophase()) {
      warn "Scan $scanid too short for autophasing" ;
      $autophase = 0 ;
    }
    else {
      print "==> $scanid Phase up\n";
    }
  }

  if ($first) {
    # add a dummy scan at start, one hour early, same source and setup
    my $prestart = 1./24.;
    my (@dummyscan) = [ 
        $modename, $start-$prestart, $prestart, $src, $sources{$src}{RA},
        $sources{$src}{DEC}, 0, 0, 0, 0
      ];
    push @atcascans, @dummyscan ;
  }
  
  my @atcascan = [ 
      $modename, $start, $scanlength, $src, $sources{$src}{RA},
      $sources{$src}{DEC}, $autophase, $paddle, $pointing, $globalpointing
    ];
  push @atcascans, @atcascan ;
  $last_stop = $stop;
  $last_src = $src;
  $first = 0;
}

# write the experiment schedule
my ($schedname) = "$exper" . ".at.json";
$schedname =~ tr/A-Z/a-z/;
warn "Warning: Overwriting exisiting schedule \"$schedname\"\n" if (-e $schedname);
open(SCHED, '>', $schedname) || die "Could not open $schedname: $!\n";  
write_atca_schedule(\*SCHED, %modes, @atcascans, 0);
close(SCHED);

# write a *relative time* cal schedule for setting up the array
$schedname = "$exper" . "-cal.at.json";
$schedname =~ tr/A-Z/a-z/;
warn "Warning: Overwriting exisiting schedule \"$schedname\"\n" if (-e $schedname);
open(SCHED, '>', $schedname) || die "Could not open $schedname: $!\n";  
# get freq of one of the modes so we can select appropriate calibrator
my $freq = $modes{((keys %modes)[0])}->{FREQ1} ;
my ($cal_src, $cal_ra, $cal_dec) = select_cal_source($calstart, $freq);
my @calscan;
my @calscans;
foreach my $calmode (keys(%modes)) {
  @calscan = [
        $calmode, $calstart, 1/24., $cal_src, $cal_ra,
        $cal_dec, 0, 0, 0, 0
      ] ;
  push @calscans, @calscan ;
}
write_atca_schedule(\*SCHED, %modes, @calscans, 1);
close(SCHED);


sub write_atca_schedule ($\%\@$) {
  # convert scan info to json format for output
  my ($schedfile, $modes, $atcascans, $relative) = @_ ;

  # $relative sets schedule to be relative, otherwise utc

  # set the scan parameters for each scan
  my @output_scans ;
  foreach (@$atcascans) {
    my %scan ; 
    #for (my $i=0; $i<9; $i++) {
    #  print "$i ", $_->[$i], "\n";
    #}
    $scan{'correlator_configuration'} = $_->[0];
    if (not $relative) {
      $scan{'utc_start_date'} = mjd_to_utc($_->[1]) ;
    }
    $scan{'duration'} = ceil($_->[2] *24*3600);
    $scan{'name'} = $_->[3];
    $scan{'ra'} = $_->[4];
    $scan{'dec'} = $_->[5];
    if ($_->[6]) {
      $scan{'intents'} = ["CALIBRATE_PHASE"];
    }
    if ($_->[7]) {
      $scan{'scanType'} = "Paddle";
    } elsif ($_->[8]) {
      $scan{'scanType'} = "Point";
    } else {
      $scan{'scanType'} = "Normal";
    }
    if ($_->[9]) {
      $scan{'pointing'} = "Offset";
    } else {
      $scan{'pointing'} = "Global";
    }
    push @output_scans, \%scan
  }

  # schedule_info section
  my %info;
  $info{'observer'} = $PIname ;
  $info{'owner'} = $PIname ;
  $info{'proposal_code'} = $exper ;
  if ($relative) {
    $info{'schedule_type'} = "relative" ;
  } else {
    $info{'schedule_type'} = "utc" ;
  }

  # Set parameters for all correlator configurations in use
  my @correlator_configurations ;
  foreach my $modename (keys(%$modes)) {
    my %config;
    #print "$modename\n" ;
    #print $modes->{$modename};
    #print $modes->{$modename}->{FREQ1};
    $config{'name'} = $modename ;
    $config{'correlator_setting'} = "Continuum" ;
    $config{'frequency_configuration'}{'centreFreq1'} = 
        $modes->{$modename}->{FREQ1} ;
    $config{'frequency_configuration'}{'centreFreq2'} = 
        $modes->{$modename}->{FREQ2} ;
    if ($modes->{$modename}->{RECEIVER}) {
      $config{'frequency_configuration'}{'receiver'} = 
          $modes->{$modename}->{RECEIVER} ;
    }
    push @correlator_configurations, \%config;
  }

  # merge the components to the full schedule for converting to json
  my (%schedule) ;
  $schedule{'schedule'} = {};
  $schedule{'schedule'}{'document'} = "ATCA BIGCAT Scheduler file";
  $schedule{'schedule'}{'version'} = "1.0" ;
  @{$schedule{'schedule'}{'scans'}} = @output_scans ;
  @{$schedule{'correlator_configurations'}} = @correlator_configurations;
  $schedule{'schedule'}{'schedule_info'} = \%info ;

  #my $JSON = JSON->new ;
  my $JSON = JSON->new->utf8 ;
  # sort output alphabetically by keys
  $JSON = $JSON->canonical([1]);
  #my $json_schedule =  $JSON->encode_json(%schedule) ;
  #my $json_schedule =  encode_json \%schedule ;
  # convert to json with pretty layout
  my $json_schedule =  $JSON->pretty->encode (\%schedule) ;
  print $schedfile "$json_schedule\n" ;
}


sub modesetup ($$$\%) {
  my ($vex, $ant, $modename, $modes) = @_;

  my $mode = $vex->mode($modename)->{$ant};

  my @chans = $mode->chan_def;
  my @chan_data;
  foreach (@chans) {
    my %vals = ();
    my $sideband = $_->sideband;
    my $sign;
    if ($sideband eq 'L') {
      $sideband = 'LSB';
      $sign = -1;
    } else {
      $sideband = 'USB';
      $sign = +1;
    }
    
    $vals{SIDEBAND} = $sideband;
    $vals{SIDESIGN} = $sign;

    my $pol = $_->pol;
    if ($pol eq 'R') {
      $pol = 'RCP';
    } elsif ($pol eq 'L') {
      $pol = 'LCP';
    }
    
    $vals{POL} = $pol;
    $vals{BW} = $_->bw->value;

    my $f0  = $_->freq;
    $f0->unit('MHz');

    my $f1 = $f0;
    if ($_->sideband eq 'U') {
      $f1 += $_->bw;
    } else {
      $f0 -= $_->bw;
    }
    $vals{F0} = $f0->value;
    $vals{F1} = $f1->value;

    push @chan_data, \%vals;

  }

  @chan_data = sort {$a->{F0} <=> $b->{F0}} @chan_data;

  my $freq1;
  my $freq2;
  my $nfreq;
  my $bw;
  my $nchan = scalar(@chan_data);
  # Sort out what sort of mode it is

  my $lsb=0;
  foreach (@chan_data) {
    #warn "***  ", $_->{F0}, '--', $_->{F1}, '  ', $_->{BW}, ' ', $_->{POL}, ' ', $_->{SIDEBAND}, "\n";
    if ($_->{SIDEBAND} eq 'LSB') {
      $lsb = 1;
      last;
    }
  }

  if ($nchan==0) {
    die "No frequency channels found for mode $mode\n";
  } elsif ($nchan==1) {
    # Will use a AT16s type DAS profile

    if ($lsb) {
      warn "Does not support LSB in mode $modename\n";
      addmode($modes, $modename, 0, 0, 0, 0);
      return;
    }

    $nfreq = 1;
    $freq1 = $chan_data[0]->{F0}+$chan_data[0]->{BW};
    $freq2 = $freq1;
    $bw = $chan_data[0]->{BW};
  } elsif ($nchan==2) {

    if ($chan_data[0]->{F0}==$chan_data[1]->{F0}) {

      $nfreq = 1;
      $bw = $chan_data[0]->{BW};
      if ($bw==16) {
  if ($lsb) {
    warn "Does not support LSB in mode $modename\n";
    addmode($modes, $modename, 0, 0, 0, 0);
    return;
  }
  # VSOP mode
  $freq1 = $chan_data[0]->{F0}+$chan_data[0]->{BW};

      } elsif ($chan_data[0]->{BW}==64) {
  $freq1 = ($chan_data[0]->{F0}+$chan_data[0]->{F1})/2;
      } else {
  $freq1 = $chan_data[0]->{F0}+$chan_data[0]->{BW}/2;
      }
      $freq2 = $freq1;

    } elsif ($chan_data[0]->{BW}==64 && $chan_data[1]->{BW}==64) {
      $nfreq = 2;
      $freq1 = $chan_data[0]->{F0}+$chan_data[0]->{BW}/2;
      $freq2 = $chan_data[1]->{F0}+$chan_data[1]->{BW}/2;
    } elsif (($chan_data[0]->{F0}==$chan_data[1]->{F1} ||
        $chan_data[0]->{F1}==$chan_data[1]->{F0})
      ) {
      die "No not support single pol modes if BW != 16 MHz\n"
  if ($chan_data[0]->{BW}!=16 || $chan_data[1]->{BW}!=16);
      
      # VSOP mode
      $nfreq = 1;
      $freq1 = ($chan_data[0]->{F0}+$chan_data[1]->{F1})/2;
      $freq2 = $freq1;
    } else {
      $nfreq = 2;
      $freq1 = $chan_data[0]->{F0}+$chan_data[0]->{BW}/2;
      $freq2 = $chan_data[1]->{F0}+$chan_data[1]->{BW}/2;
    }
    $bw = $chan_data[0]->{BW};
    
  } elsif ($nchan==4) {
    # Check we understand this mode

    $bw = undef;
    for (my $i=0; $i<4; $i++) {
      if (!defined $bw) {
  $bw = $chan_data[$i]->{BW};
      } else {
  if ($bw != $chan_data[$i]->{BW}) {
    die "Does not supports mixed bandwidths\n";
  }
      }
    }
    
    # Single pol??
    my $singlepol = 1;
    my $pol = $chan_data[0]->{POL};
    foreach (@chan_data) {
      if ($_->{POL} ne $pol) {
  $singlepol = 0;
  last;
      }
    }

    if ($singlepol) {
      die "Do not support 4 channel 64 MHz (mode $modename)" if ($bw==64);

      for (my $i=0; $i<4; $i+=2) {
  if ($chan_data[$i]->{F0}+$bw!=$chan_data[$i+1]->{F0}) {
    die "Do not support this mode setup $modename";
  }
      }

      $nfreq = 2;
      $freq1 = ($chan_data[0]->{F0}+$chan_data[1]->{F1})/2;
      $freq2 = ($chan_data[2]->{F0}+$chan_data[3]->{F1})/2;

    } else {

      if ($bw!=16) {
  if ($chan_data[0]->{F0}!=$chan_data[1]->{F0} ||
      $chan_data[2]->{F0}!=$chan_data[3]->{F0}) {
    die "Do not support this mode setup $modename";
        }
        $nfreq = 2;
  $freq1 = ($chan_data[0]->{F0}+$chan_data[0]->{F1})/2;
  $freq2 = ($chan_data[2]->{F0}+$chan_data[2]->{F1})/2;
      } elsif ($chan_data[0]->{F0}==$chan_data[1]->{F0} &&
    $chan_data[2]->{F0}==$chan_data[3]->{F0} &&
    $chan_data[0]->{F0}+$bw==$chan_data[3]->{F0}) {
  # Standard dual pol
  $nfreq = 1;
  $freq1 = ($chan_data[0]->{F0}+$chan_data[2]->{F1})/2;
  $freq2 = $freq1;
      } elsif ($chan_data[0]->{F0}==$chan_data[1]->{F0} &&
         $chan_data[2]->{F0}==$chan_data[3]->{F0}) {
  # Use VSOP profile - not sure this works with LSB
  $nfreq = 2;
  $freq1 = $chan_data[0]->{F1};
  $freq2 = $chan_data[2]->{F1};
      } elsif (($chan_data[0]->{F0}+$bw==$chan_data[1]->{F0}) &&
         ($chan_data[2]->{F0}+$bw==$chan_data[3]->{F0})) {
  $nfreq = 2;
  $freq1 = $chan_data[0]->{F1};
  $freq2 = $chan_data[2]->{F1};
      } else {
  die "Do not support this mode setup $modename";
      }
    }
  } elsif ($nchan==8) {
    my $bw = undef;

    # Check not mixed bandwidth
    foreach (@chan_data) {
      if (!defined $bw) {
  $bw = $_->{BW};
      } else {
  if (($bw != $_->{BW})) {
    my $bw2 = $_->{BW};
    warn "Does not support mixed bandwith observations ($bw/$bw2)";
    addmode($modes, $modename, 0, 0, 0, 0);
    return;
  }
      }
    }

    # Check pairs of frequencies
    for (my $i=0; $i<4; $i++) {

      if ($chan_data[$i*2]->{F0}!=$chan_data[$i*2+1]->{F0} ||
    $chan_data[$i*2]->{POL} eq $chan_data[$i*2+1]->{POL}) {
        die "Do not support this mode setup $modename";
      }
    }

    # Check 16 MHz pairs
    foreach my $i (0,4) {
      if ($chan_data[$i]->{F0}+$bw!=$chan_data[$i+2]->{F0}) {
  die "Do not support this mode setup $modename";
      }
    }
    
    $nfreq = 2;
    $freq1 = ($chan_data[0]->{F0}+$chan_data[2]->{F1})/2;
    $freq2 = ($chan_data[4]->{F0}+$chan_data[6]->{F1})/2;
    $bw = '16.0';
    
  } else {
    die "Do not support $nchan channels in mode $modename\n";
  }

  addmode($modes, $modename, $freq1, $freq2, $nfreq, $bw);
}


sub addmode ($$$$$$) {
  my ($modes, $modename, $freq1, $freq2, $nfreq, $bw) = @_;

  $modes->{$modename} = {};
  $modes->{$modename}->{FREQ1} = $freq1;
  $modes->{$modename}->{FREQ2} = $freq2;
  $modes->{$modename}->{NFREQ} = $nfreq;
  $modes->{$modename}->{BANDWIDTH} = $bw;
  $modes->{$modename}->{RECEIVER} = get_receiver($freq1);
}


sub get_receiver($) {
  # convert frequency to ATCA receiver
  # Loosely based on https://www.narrabri.atnf.csiro.au/observing/users_guide/html/atug.html#table_sensitivity, but with fuzzy edges
  my $freq = shift;
  my $receiver = "None" ;
  if ($freq > 1.0e3 and $freq < 3.2e3) {
    $receiver = "16cm";
  } elsif ($freq > 3.8e3 and $freq < 12.0e3) {
    $receiver = "4cm";
  } elsif ($freq > 15.9e3 and $freq < 26.0e3) {
    $receiver = "15mm";
  } elsif ($freq > 29.0e3 and $freq < 51.0e3) {
    $receiver = "7mm";
  } elsif ($freq > 82e3 and $freq < 106.0e3) {
    $receiver = "3mm";
  } else {
    $receiver = "" ;
	  warn "Receiver for $freq MHz not known - speak to scheduler";
  }
  return $receiver;
}


sub slewtime ($$$$) { 
  my ($az1, $el1, $az2, $el2) = @_;  # In turns


  my $vslewaz = deg2turn(40.0/60.0); # Turn/sec
  my $vslewel = deg2turn(20.0/60.0);
  my $accel   = deg2turn(800.0/3600.0); # Turn/sec/sec
  my $azcrit  = $vslewaz*$vslewaz/$accel;
  my $elcrit  = $vslewel*$vslewaz/$accel;
  my $aztcrit = 2.0 * $vslewaz/$accel;
  my $eltcrit = 2.0 * $vslewel/$accel;

  my $delaz = abs($az1-$az2);
  my $delel = abs($el1-$el2);
  while ($delaz > 1) {
    $delaz -= 1;

  }

  my ($taz, $tel);
  if ($delaz <= $azcrit) {
    $taz = 2.0*sqrt($delaz/$accel);
  } else {
    $taz = $aztcrit + ($delaz-$azcrit)/$vslewaz;
  }
  if ($delel <= $elcrit) {
    $tel = 2.0*sqrt($delel/$accel);
  } else {
    $tel = $eltcrit + ($delel-$elcrit)/$vslewel;
  }

  if ($taz>$tel) {
    return $taz+3;
  } else {
    return $tel+3;
  }
}


sub convertToDate(\$\$$$) {
  my ($ra, $dec, $mjd, $epoch) = @_;

  if ($epoch eq 'J2000' || $epoch eq 'j2000') {
    ($$ra, $$dec) = J2000todate($$ra, $$dec, $mjd);
  } elsif ($epoch eq 'B1950' || $epoch eq 'b1950') {
    ($$ra, $$dec) = fk4fk5($$ra, $$dec);
    ($$ra, $$dec) = J2000todate($$ra, $$dec, $mjd);
  } else {
    carp "Cannot convert ".$epoch." to date coordinates\n";
  }
}


sub checkautophase() {
  return 1;
  # see vex2atnf.pl for how to delay start until telescope on source. Hopefully
  # won't be required with BIGCAT.
}


use constant FINDSCHED => 0;
sub getintent ($) {
  my $vexfile = shift;

  my $mode=0;
  my $scan=undef;

  open(VEX, $vexfile) || croak "Could not open $vexfile: $!\n";

  my @scans;
  while (<VEX>) {
    if ($mode==FINDSCHED) {
      if (/\$SCHED\s*\;/) {
	$mode++;
      }
    } else {
      if (/scan\s+(\S+)\s*\;/) {
	$scan=$1;
	push @scans, [$scan];
      } elsif (/\*\s*intent\s*=\s*\"([^\"]*)\"/) {
	if (!defined $scan) {
	  carp "Intent \"$1\" outside scandef = ignoring\n";
	} else {
	  push @{$scans[$#scans]}, $1;
	}
      } elsif (/endscan\s*\;/) {
	$scan=undef;
      }
    }
  }
  close(VEX);
  return @scans;
}


sub select_cal_source($$) {
  my ($calstart, $freq) = @_;
  # select a suitable ATCA calibrator for the LST
  
  my $ATCAlong = str2turn('149:34:12.00', 'D');

  my $lst = mjd2lst($calstart, $ATCAlong)*24;
  $lst -= 1;  # start roughly 1 hr before experiment start
  $lst += 24 if ($lst < 1);
  my ($src, $ra, $dec);
  if ($freq<10000) {
    if ($lst>12 || $lst<3) {
      $src="1934-638";
      $ra = "19:39:25.026";
      $dec = "-63:42:45.63";
    } else {
      $src="0823-500";
      $ra = "08:25:26.869";
      $dec = "-50:10:38.49";
    }
  } else {
    if ($lst<8) {
      $src="0420-014";
      $ra = "04:23:15.800727";
      $dec = "-01:20:33.065310";
    } elsif ($lst<15) {
      $src="1253-055";
      $ra = "12:56:11.166560";
      $dec = "-05:47:21.524580";
    } else {
      $src="1921-293";
      $ra = "19:24:51.055957";
      $dec = "-29:14:30.121150";
    }
  }

  return $src, $ra, $dec
}