#!/usr/bin/perl -w

use strict;

use Getopt::Long;
use CGI qw/:standard *table *Tr *td/;
use PGPLOT;
use Astro::Time;

sub scale ($$$);
sub subplot ($$);
sub matchchan($$$$$);
sub makeplot($$$$$$$$$$*$);
sub doplot ($@);

#my $FRINGE_CONVERT = '/home/phi196/difx-trunk/bin/fringe_find';
my $FRINGE_CONVERT =  defined($ENV{DIFXROOT}) ?
                                         "$ENV{DIFXROOT}/bin/fringe_find": 'fringe_find';

use constant CHAN     => 0;
use constant AMP      => 1;
use constant PHASE    => 2;
use constant DELAY    => 3;
use constant FRINGE   => 4;
use constant ANT1     => 5;
use constant ANT2     => 6;
use constant POL      => 7;
use constant FREQ     => 8;
use constant SIDEBAND  => 9;
use constant MAXDELAY => 10;
use constant CROSSPOL => 11;

my $dev = '/xs';
my $ant1 = undef;
my $ant2 = undef;
my $baseline = undef;
my $path = undef;
my $outfile = undef;
my $outpath = undef;
my $time = undef;
my $source = undef;
my $chans1 = undef;
my $chans2 = undef;
my $freqpol = undef;
my $sideband = undef;
my $delay1 = undef;
my $delay2 = undef;
my $search = 0;
my $bandwidth = 16;
my $zoom = 0;

use constant DEBUG => 0;

GetOptions('device=s'=>\$dev, 'ant1=s'=>\$ant1, 'ant2=s'=>\$ant2,
	   'baseline=i'=>\$baseline, 'path=s'=>\$path, 'outpath=s'=>\$outpath,
	   'outfile=s'=>\$outfile, 'time=s'=>\$time, 'source=s'=>\$source,
	   'freqpol=s'=>\$freqpol, 'delay1=f'=>\$delay1, 'delay2=f'=>\$delay2,
	   'search'=>\$search, 'datastream=i'=>\$baseline, 'zoom=i'=>\$zoom,
	   'bandwidth=f'=>\$bandwidth, 'sideband=s'=>\$sideband);

sub usage {
  my $msg = shift;

  warn "\n Error: $msg\n" if (defined $msg);
  print STDERR<<EOF;

    Usage: rtfc_plot.pl -ant1 <ant1> -ant2 <ant2> -freqpol <freq-pol>
                        -baseline <baseline> -time <time>

EOF
  exit(1);
}

usage("Missing ant1") if (!defined $ant1);
usage("Missing freqpol") if (!defined $freqpol);
usage("Missing baseline") if (!defined $baseline);
usage("Missing time") if (!defined $time);

my @freq = ();
my @pol = ();
my $fp = $freqpol;
while ($fp =~ /^([^-]+)-([^-]+)-?/) {
  push @freq, $1;
  push @pol, [split ':', $2];
  $fp = $';
}
die "Error in freqpol specification $freqpol\n" if (length($fp)>0);

my @sideband = ();
if ($sideband) {
  @sideband = split /-/,  $sideband;
}

my $auto = 0;
if (!defined $ant2 || $ant1 eq $ant2) {
  exit if $search; # No point doing autos
  # Auto correlation
  $auto=1;
  $ant2 = $ant1;
} else {

}


my ($mjd, $dayno);
if ($time =~/^(\d+)\/(\d\d:\d\d:\d\d)/) {
  $mjd = $1;
  $time = $2;
  my ($year, $ut);
  ($dayno, $year, $ut) = mjd2dayno($mjd);
} else {
  die "Invalid time $time\n";
}

my @plots = ();

if ($search) {
  open(DELAYLOG, '>>', "$ant1-$ant2.delay")
    || die "Could not open delay log: $!\n";
  if (defined $delay1 && defined $delay2) {
    print DELAYLOG "$delay1 $delay2 ";
  }
}

for (my $i=0; $i<@freq; $i++) {
  for (my $j=0; $j<@{$pol[$i]}; $j++) {

    makeplot($ant1, $ant2, $baseline, $i+$zoom, $j, $freq[$i], $pol[$i][$j],
	     $sideband[$i], $time, $search, *DELAYLOG, 0);

    if ($auto) {
      makeplot($ant1, $ant2, $baseline, $i+$zoom, $j, $freq[$i], $pol[$i][$j],
	       $sideband[$i], $time, $search, *DELAYLOG, 1);
    }
  }
}

my $nplots = scalar(@plots);
if ($nplots==0) {
  warn "No products found!!\n";
  exit(1);
}

sub sortplots {
  return -1 if ($a->[FREQ]<$b->[FREQ]);
  return  1 if ($a->[FREQ]>$b->[FREQ]);

  return -1 if ($a->[CROSSPOL]<$b->[CROSSPOL]);
  return  1 if ($a->[CROSSPOL]>$b->[CROSSPOL]);

  return $a->[POL] cmp $b->[POL];
}

@plots = sort sortplots @plots;


my @plotsdone = ();
my $p1 = 0;
my $iplot = 0;
while ($p1<$nplots) {
  my $p2 = $p1+3;  # Plotting p1:p2
  $p2 = $nplots-1 if ($p2>=$nplots);

  warn "*** Plotting $p1-$p2\n" if DEBUG;
  $iplot++ if ($nplots>4);
  doplot($iplot, @plots[$p1..$p2]) if (!$search);
  $p1 = $p2+1;
}

if ($search) {
  print DELAYLOG "\n";
  close(DELAYLOG);
} else {
  my $filename = "${ant1}-${ant2}.html";
  $filename = "$outpath/$filename" if (defined $outpath);

  my $status = open(HTML, '>', $filename);
  if (!$status) {
    die "Could not open $filename: $!\n";
    return;
  }

  my $first = 1;
  foreach (@plotsdone) {
    foreach my $freq (@{$_->[1]}) {
      if ($first) {
	$first = 0;
      } else {
	print HTML ", ";
      }
      print HTML a({-href=>$_->[0], onMouseover=>"change('$_->[0]')"},$freq);
    }
  }
  close HTML;
}

sub doplot ($@) {
  my $p1 = 0;
  my $npage = shift;
  my $nplots = scalar(@_);

  my ($nx, $ny);
  if ($nplots==1) {
    $nx=1;
    $ny=1;
  } elsif ($nplots==2) {
    $nx=1;
    $ny=2;
  } elsif ($nplots<=4) {
    $nx=2;
    $ny=2;
  } else {
    die "Cannot plot $nplots products!!\n";
  }

  my $thisdev = $dev;
  my $thisfile = $outfile;
  if (defined $thisfile) {
    if ($thisfile =~ /(.*)\.(.*)/) { # Add check for $1 not defined
      $thisfile = "$1-$npage.$2";
    } else {
      $thisfile .= "-IF$npage";
    }
  }

  if (defined $thisfile) {
    $thisdev = "$thisfile$thisdev";
    $thisdev = "$outpath/$thisdev" if (defined $outpath);
  }

  if ($auto) {
    pgbeg(0, $thisdev, $nx, $ny) || die "Error opening PGPLOT device ($thisdev)\n";
  } else {
    pgbeg(0, $thisdev, $nx, $ny*3) || die "Error opening PGPLOT device ($thisdev)\n";
  }

  my ($maxchan, $minchan, $maxamp, $minamp, $maxphase, $minphase,
      $mindelay, $maxdelay, $minfringe, $maxfringe);

  pgscr(0,1,1,1);
  pgscr(1,0,0,0);

  my %freqused = ();
  # Get max/min
  my $first=1;
  foreach (@_) {
    if ($first) {
      $first=0;

      $maxchan = $_->[0][0];
      $minchan = $maxchan;

      $maxamp = $_->[1][2];
      $minamp = $maxamp;


      $maxphase = $_->[2][2];
      $minphase = $maxphase;

      $maxdelay = $_->[3][0];
      $mindelay = $maxdelay;

      $maxfringe = $_->[4][0];
      $minfringe = $maxfringe;

    }

    foreach (@{$_->[0]}) {
      if ($_>$maxchan) {
	$maxchan = $_;
      } elsif ($_<$minchan) {
	$minchan = $_;
      }
    }
    my @a = @{$_->[1]};
    foreach (@a[1..$#a-1]) {
      if ($_>$maxamp) {
	$maxamp = $_;
      } elsif ($_<$minamp) {
	$minamp = $_;
      }
    }

    foreach (@{$_->[2]}) {
      if ($_>$maxphase) {
	$maxphase = $_;
      } elsif ($_<$minphase) {
	$minphase = $_;
      }
    }
    foreach (@{$_->[3]}) {
      if ($_>$maxdelay) {
	$maxdelay = $_;
      } elsif ($_<$mindelay) {
	$mindelay = $_;
      }
    }
    foreach (@{$_->[4]}) {
      if ($_>$maxfringe) {
	$maxfringe = $_;
      } elsif ($_<$minfringe) {
	$minfringe = $_;
      }
    }
  }

  scale($minamp, $maxamp, 0.1);
  scale($minphase, $maxphase, 0.05);
  scale($minfringe, $maxfringe, 0.1) if (!$auto);

  my ($x, $y);

  $y = 1;
  $x = 1;
  for (my $i=0; $i<@_; $i++) {
    my $data = $_[$i];

    if ($auto) {
      pgsch(1.2);
      pgpanl($x, $y);
      pgsci(1);
      pgvstd();
      #$maxamp = 5 if ($maxamp>2);
      $minamp = 0;
      pgwindow($minchan, $maxchan, $minamp, $maxamp);
      if ($data->[CROSSPOL]) {
	pgbox('BCNST', 0.0, 0, 'BNST', 0.0, 0);
      } else {
	pgbox('BCNST', 0.0, 0, 'BCN', 0.0, 0);
      }
      pgsch(1.2);
      pglab('Channel', 'Amplitude', '');
      pgsci(2);
      pgline(scalar(@{$data->[CHAN]}), $data->[CHAN], $data->[AMP]);

      pgsci(4);
      pgsch(1.5);
      pgmtxt('T', -1.5,0.97,1,$data->[POL]);

      if (exists $data->[FREQ]) {
	pgmtxt('T', -2.6,0.97,1,$data->[FREQ]." MHz");
	$freqused{$data->[FREQ]} = 1;
      }

      if (exists $data->[SIDEBAND]) {
	pgmtxt('T', -3.8,0.97,1,$data->[SIDEBAND]);
      }

      if ($y == 1) {
	my $ant = $data->[ANT1];
	if (defined $source) {
	  $ant .= " $source";
	}
	pgmtxt('T', -1.5,0.02,0,$ant);

	if (defined $time) {
	  pgmtxt('T', -1.5,0.90,1,"$dayno/$time");
	}
      }

      if ($data->[CROSSPOL]) {
	# Plot phase if crosspol
	pgsci(1);
	#pgwindow($minchan, $maxchan, $minphase, $maxphase);
	pgwindow($minchan, $maxchan, -200, 200);
	pgbox('', 0.0, 0, 'MCT', 0.0, 0);
	pgsci(3);
	pgmtxt('R', 2.5, 0.5, 0.5, 'Phase');
	pgsch(0.8);
	pgpt(scalar(@{$data->[CHAN]}), $data->[CHAN], $data->[PHASE], 17);
      }

      $x++;
      if ($x>$nx) {
	$x = 1;
	$y++;
      }
    } else { # Crosscorr

      pgsch(1);
      pgpanl($x, $y);
	pgsci(1);
      pgvstd();
      pgwindow($mindelay, $maxdelay, $minfringe, $maxfringe);
      pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0);
      pgsch(2);
      pglab('Delay', 'Amplitude', '');
      pgsci(2);
      pgline(scalar(@{$data->[DELAY]}), $data->[DELAY], $data->[FRINGE]);

      pgsci(4);
      pgsch(4.5);
      pgmtxt('T', -1.5,0.97,1,$data->[POL]);
      if ($x==1 && $y==1) {
	#pgsch(4);
	my $baseline = $data->[ANT1].'-'.$data->[ANT2];
	if (defined $source) {
	  $baseline .= " $source";
	}
	pgmtxt('T', -1.5,0.02,0,$baseline);

       	if (defined $time) {
	  #gsch(5);
	  pgmtxt('T', -1.5,0.90,1,"$dayno/$time");
	}

	if (defined $delay1 && defined $delay2) {
	  my $offset = sprintf("Offset %.2f %.2f \\gms", $delay1, $delay2);
	  pgmtxt('B', -0.7,0.03,0,$offset);
	}

      }
      if (exists $data->[FREQ]) {
	pgmtxt('T', -2.6,0.97,1,$data->[FREQ]." MHz");
	$freqused{$data->[FREQ]} = 1;
      }
      
      if ($data->[POL] eq 'RR'  || $data->[POL] eq 'LL') {
	pgmtxt('B', -0.7,0.97,1,sprintf('Delay = %.2f',$data->[MAXDELAY]));
      }
      
      $y++;
      pgsch(1);
      pgpanl($x, $y);
      pgsci(1);
      pgvstd();
      pgwindow($minchan-1, $maxchan+1, $minamp, $maxamp);
      pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0);
      pgsch(2);
      pglab('Channel', 'Amplitude', '');
      pgsci(2);
      pgline(scalar(@{$data->[CHAN]}), $data->[CHAN], $data->[AMP]);
      
      $y++;
      pgsch(1);
      pgpanl($x, $y);
      pgsci(1);
      pgvstd();
      pgwindow($minchan, $maxchan, $minphase, $maxphase);
      pgbox('BCNST', 0.0, 0, 'BCNST', 0.0, 0);
      pglab('Channel', 'Phase', '');
      pgsci(2);
      pgpt(scalar(@{$data->[CHAN]}), $data->[CHAN], $data->[PHASE], 17);

      $x++;
      if ($x>$nx) {
	$x = 1;
	$y++;
      } else {
	$y -= 2;
      }
    }
  }
  pgend();

  if (defined $outfile) {
    my @freqused = keys %freqused;
    if (@freqused>0) {
      foreach (@freqused) {
	s/(\.\d*?)0+$/$1/;   # Remove trailing 0's after decimal
	s/\.$//;             # Remvove trailing decimal
      }
      @freqused = sort @freqused;
      push @plotsdone, [$thisfile, [@freqused]];
    }
  }
}

sub scale ($$$) {
  my ($min, $max, $scale) = @_;
  
  my $delta = ($max-$min)*$scale;

  $_[0] -= $delta;
  $_[1] += $delta;
}
  
sub matchchan($$$$$) {
  my ($pol1, $pol2, $freq1, $freq2, $crosspol) = @_;

  #warn "Comparing $pol1:$freq1  to $pol2:$freq2  ($crosspol)\n";

  return 0 if ($freq1 ne $freq2);

  if ($pol1 eq $pol2) {
    return (! $crosspol);
  } else {
    return $crosspol;
  }
}

sub makeplot($$$$$$$$$$*$) {
  my ($ant1, $ant2, $baseline, $ifreq, $product, $freq, $pol, $sideband, $time, $search,
      $dlog, $crosspol) = @_;

  warn "*** Plotting $ant1-$ant2 $baseline $ifreq $product $pol  $freq MHz $sideband\n" if DEBUG;
  my $auto = 0;
  $auto =1 if ($ant1 eq $ant2);

  $ant1 = uc($ant1);
  $ant2 = uc($ant2);

  my $hms = $time;
  $hms =~ s/://g;

  my $fname;
  if ($auto) {
    $fname = "datastream_${baseline}_crosspolar_${crosspol}_product_${ifreq}_${mjd}_${hms}_*_bin_0.output";
  } else {
    $fname = "baseline_${baseline}_freq_${ifreq}_product_${product}_${mjd}_${hms}_*_source_0_bin_0.output";
  }
  if (defined $path) {
    $fname = "$path/$fname";
  }

#  printf "DEBUG: FNAME=$fname\n";
  my @fglob = glob $fname; # Expand wildcard, as we don't know the fractional second

  if (!@fglob) {
    #warn "Could not expand $fname\n";
    return;
  } else {
      #printf "DEBUG:     @fglob\n";
    $fname = $fglob[0];
  }

  if (! -f $fname) {
    warn "Cannot find $fname\n";
    return;
  }

  open(CORR, $fname) || die "Could not open $fname: $!\n";

  my @chan;
  my @amp;
  my @phase;
  my @delay;
  my @fringe;
  my $maxdelay;
  my $maxfringe;

  while (<CORR>) {
    if (/^(\S+)\s+(\S+)(\s+(\S+))?\s*/) {
      push @chan, $1;
      push @amp, $2;
      push @phase, $4/3.14*180 if defined $4;
    }
  }
  close(CORR);

  if (!$auto) {
    # Convert to lag spectrum
    my $nchan = scalar(@chan);

    my $cmd = "$FRINGE_CONVERT $fname $nchan";
    if (-e $FRINGE_CONVERT) {
	print "DEBUG: $cmd\n";
      my $fringefind = `$cmd`;
      print $fringefind;
    } else {
      warn "Cannot run: $cmd\n";
    }

    open(LAGSPACE, 'lagspace.out') || die "Could not open lagspace.out: $!\n";

    while (<LAGSPACE>) {
      if (/^(\S+)\s+(\S+)\s*/) {
	push @delay, $1;
	push @fringe, $2;
      }
    }

    close(LAGSPACE) || die "Error reading lagspace.out\n";

    $maxfringe = abs($fringe[0]);
    $maxdelay = $delay[0];
    for (my $i=0; $i<@fringe; $i++) {
      if (abs($fringe[$i])>$maxfringe) {
	$maxfringe = abs($fringe[$i]);
	$maxdelay = $delay[$i];
      }
    }

    $maxdelay *= 1/(2*$bandwidth*1e6)*1e6;
    
    if ($search) {
      printf $dlog "  %.2f  %.2f", abs($maxfringe), $maxdelay;
    }
  }
  

  if ($auto) {
    if ($crosspol) {
      if ($pol eq 'R') {
	$pol .= 'L';
      } elsif ($pol eq 'L') {
	$pol .= 'R';
      } else {
	$pol .= '?';
      }
    } else {
      $pol .= $pol;
    }
  } else {
    $crosspol = 1 if ($pol ne 'RR' && $pol ne 'LL');
  }

  push @plots, [[@chan], [@amp], [@phase], [@delay], [@fringe], $ant1, $ant2,
		$pol, $freq, $sideband, $maxdelay, $crosspol];

}