#!/usr/local/bin/perl -w
use Getopt::Long;
use Astro::Vex;
use Astro::Time;
use Astro::Coord;
use POSIX;
use Carp;
#use IO::Tee;

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

use strict;

sub addsource ($$$$$$$$$$$$$$);
sub addmode ($$$$$$);
sub modesetup ($$$\%);
sub printSchedStart ($$$$);
sub makeCalSched ($$\@\%$);
sub getintent ($);
sub slewtime ($$$$);
sub convertToDate(\$\$$$);

use constant ATCA => 'At';

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

my $refocus = 1;  # Telescopes will need refocusing between scans - ensure first scan is refocused

GetOptions('cycle=i'=>\$cycle, 'cal=i'=>\$caltime,
	   '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 || $_ =~ /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 $lastut = 100;
my @calmodes = ();
my $calstart;

sub mjdstr($) {
  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);
}


my $MAX_SCANS = 699 ;
my $OVERLAP_SCANS = 200 ;
my $schedfilenum = 0;
my $scans_end = 0; 
my $scans_start = 0; 

while ($scans_end == 0 || @scans > $scans_start+$MAX_SCANS) {
  my $first = 1;

  # limit of 699 scans per file for ATCA scheduler. So break up the scan list into groups <= 699 but with some overlap so schedule changes
  # can happen at a flexible time.
  ++$schedfilenum;
  my $schedname ;
  if (@scans <= $MAX_SCANS) {
    $schedname = "$exper.sch"; 
  }
  else {
    $schedname = "$exper" . $schedfilenum . ".sch";
  }
  $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";  

  $scans_start = $scans_end - $OVERLAP_SCANS if $scans_end > 0 ;
  $scans_end = $scans_start + $MAX_SCANS;
  $scans_end = $#scans-1 if $#scans < $scans_end ;
  #print "scans_start=$scans_start, scans_end=$scans_end \n";
  my @scans_this_file = @scans[$scans_start .. $scans_end+1] ;
  #print "nscans=", scalar(@scans), " scans_start=", $scans_start, " scans_end=", $scans_end, "\n";

  foreach (@scans_this_file) {
    # Run through the N scans (<= $MAX_SCANS) that we will write to this file
    # Round times up to nearest cycle;
    
    ##my $start = ceil($_->start*60*60*24/$cycle-0.001)/(60*60*24/$cycle);
    my $start = $_->start;
    my $scanlen = $_->stations(ATCA)->datastop;
    my $stop = $start+$scanlen->unit('day')->value;
    ##$stop = ceil($stop*60*60*24/$cycle-0.001)/(60*60*24/$cycle);
  
    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);
  
      push @calmodes, $modename;
  
    }
    
    my $freq1 = $modes{$modename}->{FREQ1};
    my $freq2 = $modes{$modename}->{FREQ2};
    my $bw = $modes{$modename}->{BANDWIDTH};
    my $scanlength = '';
  
    print SCHED "\$SCAN*V5\n";
  
    if ($first) {
      $oldfreq1 = $freq1;
      $oldfreq2 = $freq2;
      $nfreq = $modes{$modename}->{NFREQ};
  
      my $config = 'UNKNOWN';
      my $band1 = '64.0';
      my $band2 = '64.0';
  
      $calstart = $start;
      $last_stop = $start;
  
      my ($day, $month, $year, $ut) = mjd2cal($start);
      #$ut = turn2str($ut,'H',0);
      my $date = sprintf("%02d/%02d/%04d", $day, $month, $year);
  
      $freq1 = 'SETME' if ($freq1 == 0);
      $freq2 = 'SETME' if ($freq2 == 0);
  
      printSchedStart(\*SCHED, $exper, $PIname, 'UTC');
  
      # Start 1hr early
      addsource($src, $sources{$src}{RA}, $sources{$src}{DEC}, $sources{$src}{EPOCH},
  	      $start-1/24.0, 1/24.0, $freq1, $freq2, undef, $caltime, 1, 0, 0, 0);
      print SCHED "\$SCAN*V5\n";
  
    } else { # Not first
  
    #   if ($freq1==$oldfreq1) {
    #     $freq1 = undef;
    #   } else {
    #     $oldfreq1 = $freq1;
    #   }
    #   if (defined $freq2 && defined $oldfreq2 && $freq2==$oldfreq2) {
    #     $freq2 = undef;
    #   } else {
    #     $oldfreq2 = $freq2;
    #   }
    }
  
    $scanlength = ($stop-$last_stop);
  
    my $autophase = 0;
    my $pointing = 0;
    my $paddle = 0;
    my $delay = undef;
    if ($useintent) {
      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 DETERMINE/i) {
  	  $autophase = 1;
  	} elsif ($_ =~ /REFERENCE_POINTING_DETERMINE/i || $_ =~ /REFERENCE POINTING DETERMINE/i) {
  	  $pointing = 1;
  	} elsif ($_ =~ /ATCA:PADDLE/i) {
  	  $paddle = 1;
  	}
        }
      }
      if ($autophase) {
        # What if first source
        print "==> $scanid Phase up\n";
  
        my $slew;
        if (!$first) {
  	# Az/El When we start the slew
  	my $ra = str2turn($sources{$last_src}{RA},'H');
  	my $dec = str2turn($sources{$last_src}{DEC},'D');
  	convertToDate($ra, $dec, $last_stop, $sources{$last_src}{EPOCH});
  	my $lst = mjd2lst($last_stop, $atca_long);
  	my $ha = $lst - $ra;
  	my ($az1, $el1) = eqazel($ha, $dec, $atca_lat);
        
  	# Caculate Az/El of new source at time slew starts
  	$ra = str2turn($sources{$src}{RA},'H');
  	$dec = str2turn($sources{$src}{DEC},'D');
  	convertToDate($ra, $dec, $last_stop, $sources{$src}{EPOCH});
  	$ha = $lst - $ra;
  	my ($az2, $el2) = eqazel($ha, $dec, $atca_lat);
  	
  	$slew = slewtime($az1, $el1, $az2, $el2);
  	
  	# update time at end of slew and recalculate slew time
  	my $mjd2 = $last_stop + $slew/(60*60*24);
  	$lst = mjd2lst($mjd2, $atca_long);
  	$ha = $lst - $ra;
  	($az2, $el2) = eqazel($ha, $dec, $atca_lat);
  
  	$slew = slewtime($az1, $el1, $az2, $el2);
        } else {
  	$slew=0;
        }
        my $schedslew = 0;
        if ($_->stations(ATCA)->datastart>0) {
  	$schedslew = ($start-$last_stop)*60*60*24 + $_->stations(ATCA)->datastart;
        }
        $slew = $schedslew->unit('sec')->value if ($schedslew>$slew);
  	
        $delay = ceil($slew/$cycle - 0.01)+$caltime+3; # Allow 3 extra cycles and allow for the 3 to actually phase up on
        if ($delay*$cycle>$scanlength*60*60*24) {
  	my $d = int $delay*$cycle;
  	my $s = int $scanlength*60*60*24;
  	warn "Not enough time for autophase on scan $scanid ($d/$s)\n";
  	$delay = undef;
        }
      }
    }
  
    addsource($src, $sources{$src}{RA}, $sources{$src}{DEC},
  	    $sources{$src}{EPOCH}, $last_stop, $scanlength,
  	    $freq1, $freq2, $delay, $caltime, 0, $pointing, $paddle, $globalpointing);
  
    $first = 0;
    $last_stop = $stop;
    $last_src = $src;
  }
  
  close(SCHED);
}

makeCalSched($exper, $PIname, @calmodes, %modes, $calstart);

sub addsource ($$$$$$$$$$$$$$) {
  my ($source, $ra, $dec, $epoch, $start, $length, $freq1, $freq2, 
      $phasedelay, $caltime, $first, $pointing, $paddle, $globalpointing) = @_;
  my $zoom;


  my $command = "";
  if ($refocus) { # Note assmption this is only change to $command till line for pointing just below
    $command = "focus default;";
    $refocus = 0;
  }

  if ($pointing) {
    print SCHED "Pointing=Update\n";
    if ($freq1 > 80000 || $freq2 > 80000) {
      $freq1 = 17000;
      $freq2 = 17000;
      $command = "focus default;" if ($command eq "");  # This assumes nothing added to $command before this line
      $refocus = 1;
    }
  } elsif ($globalpointing) {
    print SCHED "Pointing=Offset\n";
  } else {
    print SCHED "Pointing=Global\n";
  }

  if ($pointing) {
    print SCHED "ScanType=Point\n";
  } elsif ($paddle) {
    print SCHED "ScanType=Paddle\n";
  } else {
    print SCHED "ScanType=Normal\n";
  }

  if (defined $freq1 && $freq1 == 0) {
    $freq1 = 'SETME';
  }
  if (defined $freq2 && $freq2 == 0) {
    $freq2 = 'SETME';
  }

  if (defined $freq1) {
    if ($freq1>1300&&$freq1<=1728) {
      $freq1 += 512;
      $zoom=49;
    } elsif ($freq1>3927&&$freq1<4928) {
      $freq1 += 512;
      $zoom=17;
    } elsif ($freq1>25200&&$freq1<25712) {
      $freq1 -= 512;
      $zoom=17;
    } else {
      $zoom=33;
    }
    print SCHED "Freq-1=$freq1\n";
    print SCHED "Zoom1-1=$zoom\n";
  }
  if (defined $freq2) {
    if ($freq2>1300&&$freq2<1721) {
      $freq2 += 512;
      $zoom=49;
    } elsif ($freq2>3927&&$freq2<4928) {
      $freq2 += 512;
      $zoom=17;
    } elsif ($freq2>25200&&$freq2<25712) {
      $freq2 -= 512;
      $zoom=17;
    } else {
      $zoom=33;
    }
    print SCHED "Freq-2=$freq2\n";
    print SCHED "Zoom1-2=$zoom\n";
  }

  $length = turn2str($length,'H',0);

  print SCHED<<EOF;
Source=$source
RA=$ra
Dec=$dec
Epoch=$epoch
ScanLength=$length
EOF

  if (defined $start) {
    my ($day, $month, $year, $ut) = mjd2cal($start);
    my $utstr = turn2str($ut,'H',0);
    print SCHED "Time=$utstr\n";
    printf SCHED "Date=%02d/%02d/%04d\n", $day, $month, $year if ($ut<$lastut);
    $lastut = $ut;
  }

  if ($first) {
    $command .= "corr nncal $caltime; set rise 60;";  # Code above will add 'focus default' to first scan
  }
  if (defined $phasedelay) {
    $command .= "wait $phasedelay; corr pcal;";
  }
  print SCHED "Cmd=$command\n";

  print SCHED "\$SCANEND\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;
}

sub printSchedStart ($$$$) {
  my ($fh, $exper, $PIname, $timecode) = @_;

  print $fh <<EOF;
Project=$exper
Observer=$PIname
TimeCode=$timecode
CalCode=
Valid=true
Averaging=1
Environment=0
PointingOffset1=0.00
PointingOffset2=0.00
ChnBW-1=64
ChnBW-2=64
EOF

}

sub makeCalSched ($$\@\%$) {

  my ($exper, $PIname, $calmodes, $modes, $calstart) = @_;

  my $ATCAlong = str2turn('149:34:12.00', 'D');


  my $lst = mjd2lst($calstart, $ATCAlong)*24;
  $lst -= 1;  # e start roughly 1 hr before experiment start
  $lst += 24 if ($lst < 1);
  my $freq = $modes->{$calmodes->[0]}->{FREQ1};
  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";
    }
  }



  my $schedname = "$exper-cal.sch";
  $schedname =~ tr/A-Z/a-z/;
  warn "Warning: Overwriting exisiting schedule \"$schedname\"\n"
    if (-e $schedname);

  open(SCHED, '>', $schedname) || die "Could not open cal $schedname: $!\n";

  # print start of scan block maybe

  my $first = 1;
  foreach (@$calmodes) {
    print SCHED "\$SCAN*V5\n";
    if ($first) {
      printSchedStart(\*SCHED, $exper, $PIname, 'LST');
    }

    my $freq1 = $modes->{$_}->{FREQ1};
    my $freq2 = $modes->{$_}->{FREQ2};

    addsource($src, $ra, $dec, 'J2000', undef, 1/24, $freq1, $freq2, undef, $caltime, 0, 0, 0, 0);
    $first=0;
  }

}

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 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";
  }
}