#!/usr/bin/perl -w

use strict;

use IO::Select;
use IO::Socket;
use Getopt::Long;
use POSIX qw/floor fmod/;
use File::Basename;
use CGI qw/:standard *table *Tr/;
use Cwd qw(cwd);

use Astro::Time;
use Astro::Coord;
use Astro::Vex;
use DIFX::Input;
use RtFC;

sub Initialize ();
sub Shutdown ();
sub child_reaper;
sub process_obs_command ($$$$);
sub process_obsdata_command ($$);
sub process_control_command ($);
sub setup_correlation ($);
sub correlate ($$$$$$$@);
sub launch_corr ($$$$);
sub plotfringe($$$$$$$$\%$@);
sub plotit ($$$$$$$$$$$$$);
sub findvex ($;$);
sub getEOP ($);
sub arraygrep($@);
sub doPing (@);

#my $centreZoom = 1;

$Astro::Time::StrZero = 2;

use constant DEBUG => 0;
use constant WWWROOT => '/var/www/RtFC';
use constant DEFAULT_EXECUTABLE => "/home/vlbi/difx/bin/mpifxcorr";

my $AIPSPROC = 0;

my $PGDEV = 'png';

my $control_host = defined($ENV{RTFC_CONTROL}) ?
                                         $ENV{RTFC_CONTROL}: 'localhost';

my $corrhost = defined($ENV{RTFC_CORRHOST}) ? $ENV{RTFC_CORRHOST}: 'localhost';

my $nthread = defined($ENV{RTFC_NTHREAD}) ? $ENV{RTFC_NTHREAD}: 1;

my $difx_executable = defined($ENV{RTFC_DIFX}) ? $ENV{RTFC_DIFX}: DEFAULT_EXECUTABLE;

my $wwwroot = defined($ENV{RTFC_WWWROOT}) ? $ENV{RTFC_WWWROOT}: WWWROOT;

my $nocorr = 0;
my $fringesearch = 0;

GetOptions('control=s'=>\$control_host, 'nocorr'=>\$nocorr,
	   'search'=>\$fringesearch, 
	   'nthread=i'=>$nthread, 'corr=s'=>\$corrhost,
           'difxpath=s'=>\$difx_executable);

if (!@ARGV && ! $nocorr) {
  die "Usage: corr.pl <vexfile> [<vexfile> ...]\n";
}


my @corrhosts = split /:/, $corrhost;

# Open the vexfile. This really should be a changeable value from the UI.

my @vex;
foreach (@ARGV) {
  push @vex, new Astro::Vex($_);

  if ((!defined $vex[$#vex])) {
    die "Error reading $_\n";
  }
}

$SIG{INT} = \&Shutdown;
$SIG{CHLD} = \&child_reaper;

my ($sock_obs, $sock_control) = Initialize();

# Loops for ever, waiting for data on the various sockets. We might
# also want to run periodic commands

my $sel = new IO::Select($sock_obs, $sock_control);

my @ready;

# Lists of observatory OBJECTS which are currently connected
my @Obs = ();
my @ObsData = ();

# Hash of outstanding (ie still downloading) correlation objects
my %Corrs = ();

# Get the timeranges for all the vex files
my @vextimes = ();

foreach (@vex) {
  my @scans = $_->sched();

  my $start = $scans[0]->start;
  my $finish = $start+$scans[0]->maxscan/(60*60*24);

  foreach (@scans) {
    if ($_->start < $start) {
      $start = $_->start;
    }
    if ($_->start + $_->maxscan/(60*60*24) > $finish) {
      $finish = $_->start + $_->maxscan/(60*60*24);
    }
  }

  push @vextimes, [$start, $finish, $_->file];

  my ($day, $month, $year, $ut) = mjd2cal($start);
  my $start_str = sprintf("%02d/%02d/%04d %s", $day, $month, $year,
			  turn2str($ut, 'H', 0));
  ($day, $month, $year, $ut) = mjd2cal($finish);
  my $stop_ut = turn2str($ut, 'H', 0);

  printf("%s: $start_str  --  $stop_ut\n", $_->exper->exper_name);

}

# Set the xterm titlebar
print "\033]0;RtFC corr.pl\007\n";
$| = 1;  # Flush stdout

my $slept = 0;
my $count = 0;
while (1) {
  @ready = $sel->can_read(60);
  if (!@ready) {
    $count++;
    if ($count>=10) {
      print '.';
      $slept = 1;
    }
    doPing(@Obs);
    next;
  }
  if ($slept) {
    print "\n";
    $slept = 0;
  }

  foreach my $fh (@ready) {
    my $handled = 0;
    if (defined $sock_obs && $fh == $sock_obs) {
      # Observatory connection request
      $handled = 1;
      
      print "Got a new observatory socket connection...\n";
      # Create a new socket
      my $newobs = $sock_obs->accept();
      push @Obs, new RtFC::ObsLink($newobs);
      $sel->add($newobs);
    } elsif ($fh == $sock_control) {
      #print "PARENT: Got message from control socket\n";
      # Control process message
      $handled = 1;
      
      my $command = recv_command($sock_control,$sel,'control');

      if (ref $command) { # Must be a command
	#print "PARENT: Got message from control socket\n";
	process_control_command($command);
      } else {
	next if (! defined $command); # There was a read error
	if ($command==-1) { # Socket was closed remotely
	  warn "Control socket gone. Shutting down...\n";
	  Shutdown();
	}
      }
    } else {
      #print "PARENT: Got message from obs\n";
      # Must have data from one of the client obs sockets...
      
      my @killedsock = ();
      my $i=0;
      foreach my $obs (@Obs) {
	my $sock = $obs->socket;
	if ($fh==$sock) { # Yep, this is a obs socket with some data
	  $handled = 1;
	  
	  my $command = recv_command($sock,$sel,'observatory',$i,
				       \@killedsock);
	  
	  if (ref $command) { # Must be a command
	    process_obs_command($obs, $command, $i, \@killedsock);
	  } else {
	    #warn "PARENT: Received message which was not a command!\n";
	    # We can ignore these errors at the moment
	  }
	}
	$i++;
      }
      # Remove any obs sockets which were closed
      foreach (@killedsock) {
	splice @Obs, $_, 1; # Remove the i'th element
      }
      
      @killedsock = ();
      $i=0;
      foreach my $obs (@ObsData) { # Loop through the child connections
	my $sock = $obs->socket;
	if ($fh==$sock) { # Yep, this is a obs socket with some data
	  $handled = 1;
	  
	  #print "PARENT: Got message from obs data socket\n";
	  my $command = recv_command($sock,$sel,'observatory data',$i,
				     \@killedsock);
	  
	  if (ref $command) { # Must be a command
	    process_obsdata_command($obs, $command);
	  } else {
	    # We can ignore these errors at the moment
	  }
	}
	$i++;
      }
      # Remove any parent sockets which were closed
      foreach (@killedsock) {
	splice @ObsData, $_, 1; # Remove the i'th element
      }
    }
    if (!$handled) {
      warn "Did not handle message from $fh\n";
    }
  }
}

sub process_obs_command ($$$$) {
  my ($obs, $command, $i, $killedsock) = @_;

  if ($command->type==ANTENNA_NAME_COMMAND) {
    printf "Connection is from %s\n", $command->antenna_name;
    $obs->name($command->antenna_name);

  } elsif ($command->type==ANTENNA_ID_COMMAND) {
    $obs->code($command->antenna_id);

  } elsif ($command->type==REDEFINE_COMMAND) {
    if ($command->redef==OBSDATA) {
      # This is the secondary pipe for actually dumping the data
      # First, remove from @Obs list
      unshift @{$killedsock}, $i; # Store largest index first
      my $obsdata = create_data_process($obs);
      if (defined $obsdata) {
	push @ObsData, $obsdata;
	$sel->add($obsdata->socket);
      }
    } else {
      die "Cannot redefine connection to $command->redef at ";
    }
  } elsif ($command->type==PING) {
    #print "Got PING\n";
  } else {
    printf STDERR "Unknown command (%d) received\n", $command->type;
  }
}

sub process_obsdata_command ($$) {
  my ($obs, $command) = @_;

  if ($command->type==DATA_RECEIVED_COMMAND) {
    #print "PARENT: Got data from Child\n";
    # Got some data. See what other data has arrived so we can correlate it

    my $code = $command->code;
    my $mjd = $command->mjd; # Should check wrt $corr->mjd
    my $dir = $command->dir;
    my $file = $command->file;
    my $channels = $command->channels;
    my $id = $obs->code;

    printf "Try and correlate data from %s (%s)\n", $id, "$dir/$file";
    if (exists $Corrs{$code}) {
      # Check this data is expected and not already received
      my $corr = $Corrs{$code};
      my @expected_obs = @{$corr->observatories};
      my $obsfiles = $corr->obsfiles;

      # First check we expect data from this telescope
      my $found = 0;
      foreach (@expected_obs) {
	if ($_ eq $id) {
	  $found=1;
	  last;
	}
      }
      if ($found) {
	$obsfiles->{$id} = "$dir/$file"; # Save this for later

	correlate($mjd, '', $corr, $obsfiles, $dir, $id, $channels, @expected_obs);

      } else {
	warn "Not expecting data from $id!! (@expected_obs)\n";
      }

    } else {
      warn "Received data from unknown request!!\n";
    }
  } elsif ($command->type==PING) {
    #print "Got PING\n";
  } else {
    printf STDERR "Unknown command (%d) received\n", $command->type;
  }
}

sub process_control_command ($) {
  my ($command) = @_;

  if ($command->type==ANTENNA_NAME_COMMAND) {
    printf "Connection is from %s\n", $command->antenna_name;
  } elsif ($command->type==ERROR_COMMAND) {
    if ($command->error==DUPLICATE_CONNECTION) {
      # Cannot connect as there is another connection already
      warn "Correlator connection already present. Aborting...\n";
      Shutdown();
    } else {
      printf STDERR "Unknown error code: %d\n", $command->error;
    }

  } elsif ($command->type==CORRGRAB_COMMAND) {
    setup_correlation($command);

  } elsif ($command->type==PING) {
    #print "Got PING\n";
  } else {
    printf STDERR "Unknown command (%d) received\n", $command->type;
  }
}

sub create_data_process ($) {
  my $obs = shift;

  # We want to fork and create a local loop polling on the data pipe

  my $child = new IO::Socket;
  my $parent = new IO::Socket;

  # Create a socket connection between the parent and child
  socketpair($child, $parent, AF_UNIX, SOCK_STREAM, PF_UNSPEC)
    or  die "Could not create socketpair: $! at ";

  #print "Send PING message\n";
  my $ping = new RtFC::Command::Ping;
  my $status = $ping->send($parent, 'Ping test');
  if (!$status) {
    print "Error sending PING command";
  }

  my $pid;

  if ($pid = fork) { # Parent fork
    $sel->remove($obs->socket);
    $obs->socket($child);
    $obs->pid($pid);
    return $obs;

  } else {
    # Child fork
    die "cannot fork: $! at " unless defined $pid;
    #$child->close;
    my $childmsg = sprintf 'Child(%s):', $obs->code;

    my ($buf);

    # Set up own own private loop accepting data from the observatory
    my @ready;
    my $data_pending = 0;
    my $last_msg = undef;
    my $data_code = undef;
    my $datadir = undef;
    my $datafile = undef;
    my $mjd = undef;
    my $channels = undef;

    my $sel = new IO::Select($parent, $obs->socket);

    while (1) {
      while(@ready = $sel->can_read) {
	foreach my $fh (@ready) {
	  if ($fh == $obs->socket) {
	    # Message from observatory
	    if ($data_pending) { # A previous message said we would now be
	                         # getting data
	      my $status = $fh->recv($buf, $data_pending);
	      $last_msg = 1e999; # Very large

	      if (! defined $status) {
		warn "$childmsg Error reading data from observatory socket\n";
		$data_pending=0; # What else to do...
		$data_code=undef;
		$datadir = undef;
		$datafile = undef;
		$mjd = undef;
		$channels = undef;
	      } elsif (length $buf == 0) {
		# Socket must have closed
		print "Remote obsdata socket closed\n";
		exit;
	      } else { # Got data OK do somethign with it

		if ($data_pending<length($buf)) {
		  warn "$childmsg Something went really wrong with that data ".
		    "transfer\n";
		  warn "$childmsg More data that expected was received\n";
		  $data_pending=0;
		  $data_code = undef;
		  $datadir = undef;
		  $datafile = undef;
		  $mjd = undef;
		  $channels = undef;
		} else {

		  # Temporary - development
		  if ($last_msg-$data_pending>100000) {
		    #printf "$childmsg Still got %d kbytes to go\n", 
		    #  $data_pending/1024;
		    $last_msg = $data_pending;
		  }

		  $status = syswrite(OBSDATA, $buf);
		  if (! $status) {
		    warn "Error writing data file: $!\n";
		    Shutdown();
		  } elsif ($status != length($buf)) {
		    warn "Error writing data file\n";
		    Shutdown();
		  }

		  $data_pending -= length($buf);


		  if ($data_pending==0) { # Now got all data

		    printf("$childmsg Got all data\n");
		    close(OBSDATA);

		    my $message = new RtFC::Command::Data_Received($data_code,
								   $mjd,
								   $datadir,
								   $datafile,
								   $channels);
		    #print "Sent data received back to parent\n";

		    $status = $message->send($parent, 'Parent Process');
		    #print "done\n";
		    $data_code = undef;
		    $datadir = undef;
		    $datafile = undef;
		    $mjd = undef;
		    $channels = undef;
		  }
		}
	      }
	    } else {
	      my $command = recv_command($fh, undef,'observatory2');
	      if (ref $command) { # Must be a command
		if ($command->type==TRANSFER_COMMAND) {
		  if ($command->transfer_type==DATA_TRANSFER) {
		    # This is what we are after
		    $data_code = $command->code;
		    $mjd = $command->mjd;
		    $channels = $command->channels;
		    $data_pending = $command->datasize;

		    my ($day, $month, $year, $ut) = mjd2cal($mjd);
		    my ($sign, $hour, $minute, $second) =
		                             time2hms($ut, 'H', 0);

		    $datafile = sprintf("%s.data", $obs->code);
		    $datadir = sprintf("%04d%02d%02d-%02d%02d%02d",
				       $year, $month, $day, $hour,
				       $minute, $second);

		    printf "Got %d bytes to read\n", $data_pending;

		    #my @chans;
		    #if ($channels==0) {
		    #  @chans = ('*');
		    #} else {
		    #  @chans = mask2chan($channels);
		    #}

		    #my $chanfile = "$datadir/$datafile";
		    #$chanfile =~ s/\.data$/\.chan/;
		    #open(OBSCHAN, '>', "$chanfile")
		    #  || die "Could not open $$chanfile: $!\n";

		    #print OBSCHAN "@chans\n";
		    #close(OBSCHAN);

		    open(OBSDATA, '>', "$datadir/$datafile")
		              || die "Could not open $datadir/$datafile: $!\n";
		  } else {
		    warn "Received unexpected data transfer request - ".
		      "aborting\n";
		    exit;
		  }
		} else {
		  printf STDERR "$childmsg Unknown command (%d) received\n",
		    $command->type;
		}
	      } elsif (! defined $command) { # There was a read error
		warn "$childmsg Observatory data connection lost\n";
		exit;
	      } else {
		warn "$childmsg Unexpected mesage received!!! ";
	      }
	    }
	  } elsif ($fh == $parent) {
	    # Message from parent
	    my $command = recv_command($fh, undef, 'parent');

	    if (ref $command) { # Must be a command
	      print "$childmsg: Got message from parent\n";

	      printf " of type %s\n", $command->name;

	    } else {
	      next if (! defined $command); # There was a read error
	      if ($command==-1) { # Socket was closed remotely
		warn "$childmsg Parent data connection lost\n";
		exit;
	      } else {
		warn "$childmsg Unknown error ($command) on Parent ".
		  "connection\n";
	      }
	    }

	  } else {
	    warn "$childmsg Got data from unknown socket!!\n"
	  }
	}
      }
    }
  }
}

sub setup_correlation ($) {
  my ($command) = @_;

  my ($day, $month, $year, $ut) = mjd2cal($command->mjd);
  my ($sign, $hour, $minute, $second) = time2hms($ut, 'H', 0);

  my $newcorr = $command->newcorr;
  
  my $msg;
  if ($newcorr) {
    $msg = 'Grab';
  } else {
    $msg = 'Re-correlation';
  }

  printf("$msg requested for %02d/%02d/%04d %s %.1f\n", $day, $month, $year,
	 turn2str($ut, 'H', 0), $command->integration);

  my $dir = sprintf("%04d%02d%02d-%02d%02d%02d", $year, $month, $day,
		    $hour, $minute, $second);
  if ($newcorr) {
   
    # Fill in obsfile field
    my %obsfile = ();
    foreach ($command->observatories) {
      $obsfile{$_} = undef;
    }
    $command->obsfiles(\%obsfile);

    $Corrs{$command->code} = $command;

    # Return to control process channels needed

    my $vex = findvex($command->mjd, $command->vex);
    if (!defined $vex) {
      my $epoch = mjd2timestr($command->mjd);
      if ((defined $command->vex) && $command->vex) {
	warn "Could not find named schedule (".$command->vex.") corresponding to $epoch\n";
      } else {
	warn "Could not find schedule corresponding to $epoch\n";
      }
      return;
    }

    # Check there is a valid scan associated with this

    my $scan = findscan($vex, $command->mjd);
    if (!$scan) {
      my $epoch = mjd2timestr($command->mjd);
      warn "Could not find a scan corresponding to $epoch in ".$vex->file."\n";
      return;
    }
  
    my $modename = $scan->mode;
    my $mode = $vex->mode($modename);
    #print "DEBUG: $modename has stations ".join(' ',keys(%$mode))."\n";
    
    # Create a directory to store the data
    if (-d $dir) {
      #warn "Requested to create existing directory $dir\n";
    } else {
      mkdir $dir || die "Could not create $dir: $!\n";
    }

    # Create Obsgrab commands and send back to control process to send to observatory
    
    my $grabcommand = new RtFC::Command::ObsGrab($command->code, 
						 $command->mjd,
						 $command->integration,
						 0, 0, '', 'LBA');
    my @vexstations = $vex->stationlist;

    foreach my $obs (@{$command->observatories}) {

      next if (! arraygrep($obs, @vexstations));
      if (exists $command->secoffsets->{lc($obs)}) {
	$grabcommand->secoffset($command->secoffsets->{lc($obs)});
      } else {
	$grabcommand->secoffset(0);
      }
      if (exists $command->chans->{lc($obs)}) {
	$grabcommand->channels($command->chans->{lc($obs)});
      } else {
	$grabcommand->channels(0);
      }
      $grabcommand->antid($obs);

      ## Check if mark5 format
      my $das = $vex->das($obs);

      my $dataformat;
      if ($das->electronics_rack_type eq 'LBA') {
	  #print "DEBUG: $obs is LBADAS\n";
	$dataformat = 'LBA';
      } else {
	my $stationmode = $mode->{$obs};
	my $nbits = 2;  # Assume 2 for now
	$nbits = 1 if ($obs eq 'Wz');
	my $nchan = $stationmode->chan_def;
	my $samplerate = $stationmode->sample_rate->unit('Ms/sec')->value;
	$nbits = 8 if ($samplerate==1024);
	my $rate = $nbits*$nchan * $samplerate;
	#my $fanout = 4;  # Need to set
	my $fanout = 2;  # Need to set

	my $framesize;
	my $tformat = uc $stationmode->tracks->track_frame_format;

	if ($tformat eq 'MARK4') {
	  $dataformat = "MKIV1_${fanout}-${rate}-${nchan}-$nbits";
	  print "dataformat($obs) = $dataformat\n";
	} elsif ($tformat eq 'MARK5B') {
	  $dataformat = "Mark5B-${rate}-${nchan}-$nbits";
	} elsif ($tformat eq 'VLBA') {
	  my $fanout = 4;  # Need to set
	  $dataformat = "MKIV1_${fanout}-${rate}-${nchan}-$nbits";
	} elsif ($tformat =~ /^(VDIF[DCL]?)\/(\d+)\/(\d+)/ ||
		 $tformat =~ /^(VDIF[DCL]?)\/(\d+)/ ||	
		 $tformat =~ /^(VDIF[DCL]?)(\d+)/ ||
		 $tformat =~ /^(VDIF[DCL]?)/) {
	  if (defined $2) {
	    $framesize=$2-32;
          } else {
	    $framesize = '8192';
	  }
	  my $format = $1;
	  $framesize += 16 if ($format eq 'VDIFL');
	  $format = 'VDIFC' if ($format eq 'VDIFD');
	  $nbits = $3 if (defined $3);
	  $rate = $nbits*$nchan * $samplerate;
	  $dataformat = "${format}_${framesize}-${rate}-${nchan}-$nbits";
	}
      }
      $grabcommand->mode($dataformat);

      my $status = $grabcommand->send($sock_control, $obs);

      if (!defined $status) { # Ignore other error states for the moment
	warn "Failed to send Observatory Grab command to control for $obs\n";
      }
    }


  } else {
    # Does the directory exist
    if (! -d $dir) {
      warn "Directory $dir does not exists. Cannot continue re-correlation\n";
      #delete $command->value->[0]; ##Not sure what this was meant to do...
      return;
    }

    do_recorr($command, $dir);
  }
}

sub do_recorr($$) {
  my ($corr, $dir) = @_;

  my @obs = @{$corr->observatories};
  
  # Fill in obsfile field
  my %obsfile = ();
  foreach (@obs) {
    my $file = "$dir/${_}.data";
    if (! -f $file) {
      warn "$file does not exist. Skipping $_\n";
      next;
    }
    $obsfile{$_} = $file;
  }
  $corr->obsfiles(\%obsfile);

  correlate($corr->mjd, $corr->vex, $corr, \%obsfile, $dir, undef, undef, @obs);
}

sub correlate($$$$$$$@) {
  my ($mjd, $vexname, $corr, $obsfile, $dir, $newant, $channels, @expected_obs) = @_;
  my $pwd = `pwd`;
  chomp $pwd;
  print "Grabbed pwd to be $pwd\n";

  # Determine the values to write to the config file

  my ($day, $month, $year, $ut) = mjd2cal($mjd);
  my ($sign, $hour, $minute, $second) = time2hms($ut, 'H', 0);
  my $epoch = sprintf('%04d%02d%02d:%02d%02d%02d', $year, $month,
		      $day, $hour, $minute, $second);
  my $time = sprintf('%02d:%02d:%02d', $hour, $minute, $second);
  my $session = sprintf("%s%d", month2str($month), $year);
  my $integration = $corr->integration;

  # Find which vexfile this mjd corresponds to

  my $vex = findvex($mjd, $vexname);
  if (!defined $vex) {
    if ((defined $vexname) && $vexname) {
      warn "Could not find named schedule ($vexname) corresponding to $epoch\n";
    } else {
      warn "Could not find schedule corresponding to $epoch\n";
    }
    return;
  }

  # Read the vexfile for experiment id, source name and position
  # and observing frequency
  my $exper_name = $vex->exper->exper_name;

  printf("Correlating data for %s\n", $exper_name);

  my $found = 0;
  my ($sourcename, $modename, $rate, @scanstations);

  my $scan = findscan($vex, $mjd);

  if ($scan) {
    printf("Found: %s start=%s source=%-8s\n", $scan->scanid,  $scan->startepoch, 
	   $scan->source);
    $sourcename = $scan->source;
    $modename = $scan->mode;
    my @stations = $scan->stations;
    foreach (@stations) {
      push @scanstations, $vex->station($_->station)->site_id;
    }
    @scanstations = sort @scanstations;
  } else {
    print "No scan found\n";
    return;
  }

  my $source = $vex->source($sourcename);
  $sourcename = $source->source_name;

  # Get a list of all antenna we can correlate (ie the valid ones)
  
  my @obsgot;
  foreach my $o (@expected_obs) {
    if (defined $obsfile->{$o}) {
      push @obsgot, $o;
    }
  }
  return if (scalar(@obsgot)<2); # Only one antenna present
  @obsgot = sort @obsgot;

  my $npoint = $corr->npoint;
  my $pack = $corr->pack;

  my $offsets = $corr->offsets;
  my $secoffsets = $corr->secoffsets;
  my $clockrates = $corr->clockrates;

  my @offsets;
  # Fill in default values
  foreach (@obsgot) {
    print "Doing defaults for $_\n";
    if (!exists $offsets->{lc($_)}) {
      $offsets->{lc($_)} = 0.0;
      print "Setting $_ clock delay to 0.0\n";
    } else {
      print "$_ delay ", $offsets->{lc($_)}, "\n";
    }
    
    push @offsets, $offsets->{lc($_)};

    if (!exists $clockrates->{lc($_)}) {
      $clockrates->{lc($_)} = 0.0;
    }
  }

  my $ntel = scalar(@obsgot);

  # Need to work out how many channels in data file. First count the number of x's
  #my @packmode = reverse(split //, $packmode);

  my $vexfile = $vex->file;
  
  my $newvex = fileparse($vexfile);
  print "Copying $vexfile to $newvex\n";
  system("cp -f $vexfile $dir/$newvex");

  #my $mjdStart = sprintf('%.9f', $mjd);
  my $mjdStart = mjd2vextime($mjd);
  my $mjdStop;
  if ($integration<1) {
    $mjdStop = mjd2vextime($mjd+1.0/60.0/60.0/24.0);
    #$mjdStop = sprintf('%.12f', $mjd+$integration/60.0/60.0/24.0);
  } else {
    $mjdStop = mjd2vextime($mjd+$integration/60/60/24);
    #$mjdStop = sprintf('%.9f', $mjd+$integration/60/60/24);
  }

  # Determine if zoombands are needed
  my $mode = $vex->mode($modename);

  my $zband = undef;
  my %bandwidths = ();
  foreach my $ant (@obsgot) {
    my $stationmode = $mode->{$ant};
    if (! defined $stationmode) {
      warn "WARNING: $ant not defined in $modename\n";
      next;
    }
    my @chans = $stationmode->chan_def;
    foreach (@chans) {
      my $bw = $_->bandwidth->unit('MHz')->value;
      if (!exists($bandwidths{$bw})) {
	$bandwidths{$bw} = 1;
      }
    }
  }
  my @bandwidths = keys(%bandwidths);
  if (@bandwidths>1) {
    @bandwidths = sort {$a<=>$b} @bandwidths;
    print "Got multiple bandwidths => @bandwidths\n";

    # This is fairly crude but should work for all our needs
    $zband = $bandwidths[0];
    print "Will try zoom band of $zband\n";
  } 

  # Create the config file
  my $fconfig = "$dir/$exper_name.v2d";
  my $finput = "$dir/$exper_name.input";
  print "Writing $fconfig\n";
  my $roundInt = sprintf("%.3f", $integration);
  $roundInt = '1.0' if ($AIPSPROC);
  open(V2D, '>', $fconfig) || die "Could not open $fconfig: $!\n";

  my $outputFormat = 'ASCII';
  $outputFormat = 'SWIN' if ($AIPSPROC);
  print V2D <<EOF;

vex=$newvex

maxLength = 86400
maxGap = 86400
maxSize = 100000
mjdStart=$mjdStart
mjdStop=$mjdStop
minLength=0.0001
startSeries = 0
threadsFile = threads
outputFormat = $outputFormat
tweakIntTime = True
#maxReadSize = 400000000

SETUP default
{
  tInt = $roundInt
  guardNS = 2000
EOF

  my $corrchan = $npoint;
  if (defined $zband) {
    # Determine spectral resolution. Assume # points was for narrowest band
    my $specres = $bandwidths[0]/$npoint;
    printf V2D ("  specRes = %.16f\n", $specres);
    $corrchan /= 2;
  } else {
    print V2D "  nChan=$npoint\n";
  }
  if ($corrchan > 128) {
    # Choose an XMAC length
    my $xmac = 512;
    while ($xmac > 0) {
      if (($corrchan % $xmac)==0) {
	last;
      }
      $xmac--;
    }
  }
  print V2D "}\n";

  # Correlate all antennas present so far
    
  print V2D "antennas = ", join(', ', @obsgot), "\n";

  my $nfreq = undef;
  my $first = 1;
  my %plotvals = ();

  my %obschan = ();
  foreach my $ant (@obsgot) {

    my $station = $vex->station($ant);
    my $das = $vex->das($ant);
    my $site_id = uc($station->site_id);
    my $thisoffset = $offsets->{lc($ant)};
    if (exists $secoffsets->{lc($ant)} && $das->electronics_rack_type ne 'LBA') {
	$thisoffset += $secoffsets->{lc($ant)}*1e6;
    }
    my $clockdelay = $thisoffset;
    my $clockrate = $clockrates->{lc($ant)};
    my $datafile = $obsfile->{$ant};
    my $mode = $vex->mode($modename)->{$ant};
    my @chans = $mode->chan_def;

    print V2D <<EOF;

ANTENNA $site_id {
  phaseCalInt = 0
  clockOffset = $clockdelay
#CLOCK RATE (us/s):  $clockrate
  file=$pwd/$datafile
EOF
    if ($ant eq 'Ka') {
      print V2D "  format = VDIF544\n";
    } elsif ($ant eq 'Xx') {
      print V2D<<EOF;
#  format = VDIF8224
  sampling = COMPLEX_DSB
EOF
    } elsif ($ant eq 'Md') {
      print "*******SWAPPING MEDUSA POLS*******\n";
      print  V2D "  polSwap = True\n";
    } elsif ($ant eq 'Ak') {
      #print "*******SWAPPING ASKAP POLS*******\n";
      #print  V2D "  polSwap = True\n";
    } elsif ($ant eq 'offWa') {
      #print "*******SWAPPING Warkworth 30m POLS*******\n";
      #print  V2D "  polSwap = True\n";
    }

    if ($das->electronics_rack_type eq 'LBA') {
      # Assume all channels per ant are the same (pretty safe)
      my $bw = $chans[0]->bw;
      if ($bw==64) {
	print V2D "  format = LBASTD\n";
      } else {
	print V2D "  format = LBAVSOP\n";
      }
    }

    if (defined $zband) {
      my @customZoom = ();
      if (-e 'customZoom') {
	open(ZOOM, '<', 'customZoom') || die "Could not open 'customZoom': $!\n";
	chomp (@customZoom = <ZOOM>);
	close(ZOOM);
	warn("Using custom Zoom frequencies\n");
      }

      
      my %uniqfreq = ();
      foreach (@chans) {
	my $f0 = $_->freq;
	if ($_->sideband ne 'U') {
	  $f0 -= $_->bandwidth;
	}
	$f0 = $f0->unit('MHz')->value;

	my $bw = $_->bandwidth->unit('MHz')->value;
	next if ($bw == $zband);

	my $fstart = $f0;
	for (my $f=$fstart; $f+$zband<=$f0+$bw; $f += $zband) {
	  $uniqfreq{$f} = 1 if (!exists $uniqfreq{$f});
	}
      }

      my @uniqfreq = sort {$a <=> $b} keys(%uniqfreq);
      @uniqfreq = @customZoom if (@uniqfreq>0 and @customZoom>0);
      foreach (@uniqfreq) {
	printf V2D "  addZoomFreq = freq@%.3f/bw@%.2f\n", $_, $zband;
      }
    }
    print V2D "}\n";
  }
  
  print "DEBUG: getting EOPs\n";
  my $status = getEOP($mjd);
  close(V2D) || die "Error closing $fconfig\n";
  return if ($status);

  # Run vex2difx
  printf("Running vex2difx\n");
  system("cd $dir && vex2difx $exper_name.v2d");
  if ($?) {
    if ($?==-1) {
      warn "Failed to run vex2difx: $!\n";
    } elsif ($? & 127) {
      my $sig = $? & 127;
      warn "Error: vex2difx died with signal $sig\n";
    } else {
      my $ret =$? >> 8;
      warn "vex2difx returned $ret. Terminating correlation\n";
    }
    return;
  }
  #system("cd $dir && emacs $exper_name.input");

  # Run calcdifx
  printf("Running Calc\n");
  system("cd $dir && difxcalc -f ${exper_name}.calc"); # TODO Check error status
  #system("cd $dir && calcif2 ${exper_name}"); # TODO Check error status
  if ($?) {
    if ($?==-1) {
      warn "Failed to run difxcalc: $!\n";
    } elsif ($? & 127) {
      my $sig = $? & 127;
      warn "Error: difxcalc died with signal $sig\n";
    } else {
      my $ret =$? >> 8;
      warn "difxcalc returned $ret. Terminating correlation\n";
    }
    return;
  }

  # Create machine file
  if (! open(MACHINES, '>', "$dir/machines")) {
    warn "Could not open $dir/machines: $!\n";
    return;
  }
  print MACHINES "localhost\n"x(scalar(@obsgot)+1);

  foreach (@corrhosts) {
    print MACHINES "$_\n";
  }
  close(MACHINES);

  # Create threads file
  if (! -f "$dir/threads") {
    if (! open(THREADS, '>', "$dir/threads")) {
      warn "Could not open $dir/threads: $!\n";
      return;
    }
    print THREADS <<EOF;
NUMBER OF CORES:    1
$nthread
EOF
    
    close(THREADS);
  }

  if (!$nocorr) {
    
    launch_corr("${exper_name}.input", $pwd, $dir,
		scalar(@obsgot)+scalar(@corrhosts)+1);

    if (!$AIPSPROC) {
      my $exportpath = plotfringe($vex, $nfreq, $dir, $sourcename, $time,
				  $session, $corr->postfix, $finput,
				%obschan, $offsets, @obsgot);
      
      if (defined $exportpath) {
	makehtml($exportpath, $exper_name, @scanstations);
      }
    }
    print "Done!\n";
  }
}

sub Initialize () {
  # Create a listening socket for the observatory connection
  # to the corr process
  # Connect to the central "server"

  my $sock_obs = IO::Socket::INET->new(LocalPort => CORRPORT,
				       Listen    => 10,
				       ReuseAddr => 1)
    || die "Couldn't create server socket for Observatories: $@\n";

  my $sock_control = IO::Socket::INET->new(PeerAddr  => $control_host,
					   PeerPort  => OBSPORT,
					   ReuseAddr => 1
					  )
    || die "Couldn't create socket connection to control process: $@\n";

  # Register this as a correlator connection
  my $command = new RtFC::Command::Antenna_Name('CORRELATOR');
  my $status = $command->send($sock_control, 'control process');
  # Ignore return value for now

  return ($sock_obs, $sock_control);

}

sub Shutdown () {
  warn "Calling shutdown routine\n";

  foreach my $obs (@Obs) {
    $obs->socket->shutdown(2);
    warn sprintf("Closed obs %s\n", $obs->name);
  }

  foreach my $obs (@ObsData) {
    $obs->socket->shutdown(2);
    warn sprintf("Closed obsdata %s\n", $obs->name);
  }

  $sock_obs->shutdown(2) if (defined $sock_obs);
  $sock_control->shutdown(2) if (defined $sock_control);

  exit(0);
}

sub child_reaper {
  #print "Child reaper was called @_\n";
}

sub launch_corr ($$$$) {
  my ($input, $pwd, $dir, $np) = @_;

  # Remove old .output files, they cause trouble
  my @outputfiles = glob "$dir/*.output";
  unlink @outputfiles if (@outputfiles);

#  my $exec = "cd $pwd/$dir \&\& mpirun -machinefile machines -np $np --map-by node $difx_executable $input";
  my $exec = "cd $pwd/$dir \&\& mpirun -machinefile machines -np $np $difx_executable $input";

  print "Launching:   $exec\n";
  system $exec;

  printf "mpifxcorr exited with value %d\n", $? >> 8;
}

sub plotfringe ($$$$$$$$\%$@) {
  my ($vex, $nfreq, $dir, $source, $time, $session, $postfix, $finput, $obschan, $offsets, @obs) = @_;

  my $exper = $vex->exper->exper_name;

  my @antname;
  my @antid;
  foreach (@obs) {
    my $station = $vex->station($_);
    push @antname, $station->site_name;
    push @antid, $station->site_id;
  }

  #  umask 0000;
  umask 0002;

  my $exportpath = $wwwroot;

  if (! -d $exportpath) {
    warn "$exportpath does not exists!. Saving locally\n";
    $exportpath = ".";
  }

  $exportpath .= "/$session";
  if (! -d $exportpath) {
   my $status = mkdir $exportpath;
   if (! $status) {
     warn "Could not create directory $exportpath: $!\n";
     return undef;
   }
 }
  
  $exportpath .= "/$exper";
  if (! -d $exportpath) {
   umask 0002;
   my $status = mkdir $exportpath;
   if (! $status) {
     warn "Could not create directory $exportpath: $!\n";
     return undef;
   }
 }

  $exportpath .= "/$time";
  if (defined $postfix && $postfix) {
    $exportpath .= "-$postfix";
  }
  if (! -d $exportpath) {
    umask 0002;
    my $status = mkdir $exportpath;
    if (! $status) {
      warn "Could not create directory $exportpath: $!\n";
      return undef;
    }
  }

  my $input = new DIFX::Input($finput); # TODO CHECK RETURN

  my @telescopes = $input->telescope;
  my @datastream = $input->datastream->datastreams;
  my @freqs = $input->freq;

  my $startmjd = $input->common->startmjd;
  my $startseconds = $input->common->startseconds;
  $startseconds += floor($input->common->executetime/2);
  if ($startseconds>24*60*60) {
    $startseconds -= 24*60*60;
    $startmjd += 1;
  }

  my $timestr = sprintf("%d/%02d:%02d:%02d",
			$input->common->startmjd,
			int $startseconds/3600,
			int $startseconds/60 % 60,
			int $startseconds %60);
  # Do the cross correlations
  my $b=0;
  foreach ($input->baseline)  {  # Loop over baselines
    my $dsA = $_->dstreamA;
    my $dsB = $_->dstreamB;
    my $ant1 = $telescopes[$dsA]->name;
    my $ant2 = $telescopes[$dsB]->name;
    my $outfile = "$ant1-$ant2.$PGDEV";

    #print "DEBUG: Sorting out $ant1-$ant2\n";

    my $freqpol;
    my @sideband;
    my @bandwidth; # Assume same for all
    foreach ($_->freqs) {        # Loop over frequencies
      my $baselineindexA = $_->[0]->bandA;
      my $nrecbandA = $datastream[$dsA]->nrecband;
      my $nrecbandB = $datastream[$dsB]->nrecband;
      my $freqindex;
      if ($baselineindexA<$nrecbandA) {
	my $localfreqindex = $datastream[$dsA]->recbands($baselineindexA)->index;
	#print "DEBUG:  localfreqindex = $localfreqindex\n";
	$freqindex = $datastream[$dsA]->freqs($localfreqindex)->index;
      } else {
	my $zoomfreqindex = $baselineindexA-$nrecbandA;
	my $localzoomfreqindex = $datastream[$dsA]->zoomband($zoomfreqindex)->index;
	#print "DEBUG:  localzoomfreqindex = $localzoomfreqindex\n";
	$freqindex = $datastream[$dsA]->zoomfreq($localzoomfreqindex)->index;
      }
      #print "DEBUG:  freqindex = $freqindex\n";
      my $freq = $freqs[$freqindex]->freq;
      my $sideband = $freqs[$freqindex]->sideband;
      print "DEBUG: *sideband*", $freqs[$freqindex]->sideband, "\n";
      #print "DEBUG: That means freq=$freq\n";

      push @bandwidth, $freqs[$freqindex]->bw;
      push @sideband, $sideband;
      $freq =~ s/0+$//;
      $freq =~ s/\.$//;
      my @pols;
      my ($polA, $polB);
      foreach (@$_) {
	if ($_->bandA<$nrecbandA) {
	  $polA = $datastream[$dsA]->recbands($_->bandA)->pol;
	} else {
	  $polA = $datastream[$dsA]->zoomband($_->bandA-$nrecbandA)->pol;
	}
	if ($_->bandB<$nrecbandB) {
	  $polB = $datastream[$dsB]->recbands($_->bandB)->pol;
	} else {
	  $polB = $datastream[$dsB]->zoomband($_->bandB-$nrecbandB)->pol;
	}

	push @pols, sprintf "%s%s", $polA, $polB;
      }
      my $thisfreq = "$freq-".join(':',@pols);
      if (defined $freqpol) {
	$freqpol .= "-$thisfreq";
      } else {
	$freqpol = $thisfreq;
      }
    }
    my $sideband = join('-', @sideband);
    print "DEBUG: Using sideband $sideband\n";

    my $cmd = "rtfc_plot.pl -ant1 $ant1 -ant2 $ant2 -baseline $b "
      ."-time $timestr -freqpol $freqpol -sideband $sideband "
	." -source $source -path $dir -outfile $outfile"
	  ." -outpath $exportpath -dev /$PGDEV -bandwidth $bandwidth[0]";

    if ($fringesearch) {
      my $delay1 = $offsets->{lc($ant1)};
      my $delay2 = $offsets->{lc($ant2)};
      $cmd .= " -search -delay1 $delay1 -delay2 $delay2"
    }
    print "$cmd\n";
    system $cmd;
  
    $b++;
  }

  # Do the autocorrelations
  my $d=0;
  foreach my $ds (@datastream) {
    my $freqpol;
    my @sideband;
    my $zoom = '';
    if ($ds->nzoomband>0) {
      $zoom = sprintf("-zoom %d", $ds->nrecband);
      foreach ($ds->zoomband) {
	my $freqindex = $ds->zoomfreq($_->index)->index;
	my $freq = $freqs[$freqindex]->freq;
	push @sideband, $freqs[$freqindex]->sideband;
	$freq =~ s/0+$//;
	$freq =~ s/\.$//;
	my $thisfreq = "$freq-".$_->pol;
	if (defined $freqpol) {
	  $freqpol .= "-$thisfreq";
	} else {
	  $freqpol = $thisfreq;
	}
      }

    } else {
      foreach ($ds->recbands) {
	my $freqindex = $ds->freqs($_->index)->index;
	my $freq = $freqs[$freqindex]->freq;
	push @sideband, $freqs[$freqindex]->sideband;
	$freq =~ s/0+$//;
	$freq =~ s/\.$//;
	my $thisfreq = "$freq-".$_->pol;
	if (defined $freqpol) {
	  $freqpol .= "-$thisfreq";
	} else {
	  $freqpol = $thisfreq;
	}
      }
    }
    my $sideband = join('-', @sideband);
    my $ant1 = $telescopes[$d]->name;
    my $outfile = "$ant1.$PGDEV";

    my $cmd = "rtfc_plot.pl -ant1 $ant1 -datastream $d "
      ."-time $timestr -freqpol $freqpol -sideband $sideband $zoom"
	." -source $source -path $dir -outfile $outfile"
	  ." -outpath $exportpath -dev /$PGDEV";
    print "$cmd\n";
    system $cmd;

    $d++;
  }

  return($exportpath);
}

sub findvex ($;$) {
  my $mjd = shift;
  my $vexname = shift;
  # Find which vexfile this mjd corresponds to

  print "Looking for $mjd\n";

  if (defined $vexname and $vexname =~ /\S/) {
    print "Choosing name search for $vexname\n";
    for (my $i=0; $i<@vextimes; $i++) {
      if ($vextimes[$i]->[2] eq $vexname) {
	return $vex[$i];
      }
    }
  } else {
    for (my $i=0; $i<@vextimes; $i++) {
      if ($vextimes[$i]->[0] < $mjd && $vextimes[$i]->[1] > $mjd) {
	print "Choosing ",$vextimes[$i]->[2],"\n";
	return $vex[$i];
      }
    }
  }
  return undef;
}

sub makehtml ($$@) {
  my ($exportpath, $exper, @stations) = @_;

  my $filename = "$exportpath/index.html";
  my $status = open(HTML, '>', $filename);
  if (!$status) {
    warn "Could not open $filename: $!\n";
    return;
  }

  print HTML start_html(-title=>"RtFC $exper",
			-script=>{-language=>'JavaScript',
				  -code=>'function change(imgurl) {
                                           if (document.images) {
                                              document["fringe"].src=imgurl;
                                            }
                                          }'
				 },
			-style => {-code => join('',<DATA>)},
		       );

  print HTML start_table({-border=>1}), start_Tr, td();
  foreach (@stations) {
    print HTML th("&nbsp;&nbsp;$_&nbsp;&nbsp;");
  }
  print HTML end_Tr;

  for (my $i=0; $i<@stations; $i++) {
    my $ant1 = $stations[$i];
    print HTML start_Tr;
    print HTML th("&nbsp;&nbsp;$ant1&nbsp;&nbsp;");

    for (my $j=0; $j<$i; $j++) { print HTML td()}
    for (my $j=$i; $j<@stations; $j++) {
      my $ant2 = $stations[$j];
      my $baseline;
      if (open(BASELINE, "$exportpath/".uc("$ant1-$ant2").".html")) {
	my @baseline = <BASELINE>;
	chomp @baseline;
	$baseline = join(' ', @baseline);
	close(BASELINE);
      } else {
	$baseline = '-';
      }
      print HTML td($baseline);
    }
    print HTML end_Tr;
  }

  print HTML end_table;

  print HTML img {src=>'none.png', align=>'center', name=>"fringe", height=>600};

  print HTML end_html;

  close(HTML);

  return;
}

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

  my $home = $ENV{HOME};
  my $eops = defined($ENV{RTFC_EOP}) ? $ENV{RTFC_EOP}: "$home/.eops";
  my $dutc = defined($ENV{RTFC_DUTC}) ? $ENV{RTFC_DUTC}: "$home/.ut1ls";

  if (! -f $eops) {
    warn "$eops does not exist. Need EOPS to correlate\n";
    return(1);
  }

  if (! -f $dutc) {
    warn "$dutc does not exist. Need DUTC1 values to correlate\n";
    return(1);
  }

  my $targetJD = mjd2jd($mjd);


  if (! open(DUTC, $dutc)) {
    warn "Could Not open $dutc: $!\n";
    return(1)
  };

  my $tai_utc = undef;
  while (<DUTC>) {
    if (/^\s*(\d+)\s+(\S+)\s+(\d+)\s*=\s*JD\s+(\S+)\s+TAI-UTC\s*=\s*(\S+)/) {
      my $thisjd = $4;
      my $leapsec = $5;
      last if ($thisjd > $targetJD);
      $tai_utc = $leapsec;
    } else {
      warn "Skipping $_\n";
    }
  }
  close(DUTC);
  if (!defined $tai_utc) {
    warn "Could not find TAI-UTC value\n";
    return(1);
  }
  print "TAI-UTC = $tai_utc\n";

  if (! open(EOPS, $eops)) {
    warn "Could Not open $eops: $!\n";
    return(1)
  };

  $targetJD = mjd2jd($mjd);

  my $first=1;
  my @fields;
  while (<EOPS>) {
    # Skip the first line
    if ($first) {
      $first = 0;
      next;
    }
    chomp;
    s/#.*$//; # Remove comment character
    next if (length()==0);
    
    @fields = split;
    # print an EOP line if we're within 2.5 days of the target day
    if (abs($fields[0] - $targetJD) < 2.5) {
      printf V2D "EOP %.0f { xPole=%f yPole=%f tai_utc=$tai_utc ut1_utc=%f }\n",
                jd2mjd($fields[0]), $fields[1]/10, $fields[2]/10, 
		  $tai_utc+$fields[3]/1000000;
      printf "Added EOP %s\n", $fields[0];
    }
  }

  close(EOPS);
  return(0);

}

sub findscan ($$) {
  my ($vex, $mjd) = @_;

  my @scans = $vex->sched;
  foreach (@scans) {
    return undef if ($_->start > $mjd);
    next if ($_->stop < $mjd);
    return $_;
  }
}

sub mjd2timestr ($) {
  my $mjd = shift;
  my ($day, $month, $year, $ut) = mjd2cal($mjd);
  my ($sign, $hour, $minute, $second) = time2hms($ut, 'H', 0);
  return sprintf('%04d%02d%02d:%02d%02d%02d', $year, $month,
		 $day, $hour, $minute, $second);

}

sub arraygrep($@) {
  my $search = shift;
  foreach (@_) {
    return 1 if ($search eq $_);
  }
  return 0;
}

sub doPing (@) {
  my (@Obs) = @_;
  my $status;
  my $ping = new RtFC::Command::Ping;

  foreach my $obs (@Obs) {
    $status = $ping->send($obs->socket, 'Ping');
    if (!$status) {
      print "Error sending Ping to Obs process\n";
    }
  }
}


__DATA__

th {
  background-color: #0066ff;
  color: #ffffff;
  text-align: center;
  border-style: solid;
  border-color: #0066ff;
  padding-left: 0.0em
}
td {
  text-align: center;
  padding-left: 0.5em;
  padding-right: 0.5em
}