#!/usr/bin/perl
use FileHandle;
use List::Util qw[min max];

$imode = 0;	#set default itab imode
$plt_flg = 0;

if (($#ARGV + 1) < 5){die <<EOS;}
*** $0
*** Copyright 2015, Gamma Remote Sensing, v3.7 18-Apr-2015 clw/uw ***
*** Generate baseline plot and output file with perpendicular baselines and delta_T values ***
*** Generate interferogram table (itab) file specifying SLCs for each interferogram  ***

usage: $0 <SLC_tab> <SLC_par> <bperp_file> <itab> <itab_type> [plt_flg] [bperp_min] [bperp_max] [delta_T_min] [delta_T_max] [delta_n_max]

    SLC_tab      (input) two column list of SLC filenames and SLC parameter filenames (including paths) (text)
                   1. SLC filename  (includes path)
                   2. SLC parameter filename (includes path)
    SLC_par      (input) reference SLC parameter filename (include path)
    bperp_file   (output) list of dates, bperp and delta_T for interferogram pairs in the itab (text)
    itab         (output) interferogram table with 4 column format:
                   1. row number in SLC_tab of the reference SLC
		   2. row number in SLC_tab of SLC-2 of the interferogram
		   3. line number in the itab
		   4. flag used to indicate if this interferogram is used in IPTA processing (0:not used  1:used)
    itab_type    itab type (enter - for default):
                   0: single reference (default)
		   1: all pairs
    pltflg       bperp plotting flag (enter - for default):
                   0: none (default)
                   1: output plot in PNG format
                   2: screen output
    bperp_min    minimum magnitude of bperp (m) (default = all, enter - for default)
    bperp_max    maximum magnitude of bperp (m) (default = all, enter - for default)
    delta_T_min  minimum number of days between passes (default = 0, enter - for default)
    delta_T_max  maximum number of days between passes
    delta_n_max  maximum scene number difference between passes

EOS

open(SLC_TAB, "<$ARGV[0]") or die "\nERROR $0: table SLC images and parameter files missing: $ARGV[0]\n";
$orb1_par = $ARGV[1];
$bpf = $ARGV[2];
$bp_plot = $bpf.".png";
$itab = $ARGV[3];

$log = "base_calc.log";
$bperp_max = -1;
$bperp_min = -1;
$dt_min = 0;
$dt_max = -1;
$n_max = -1;

-e $orb1_par or die "ERROR $0: reference rslc image parameter file does not exist: $orb1_par\n";
#-x "/usr/local/bin/convdate" or die "\nERROR: program convdate not available, install as part of the DELFT getorb package\n";

open(BP,">$bpf") or die "ERROR $0: cannot open output baseline file: $bpf\n";
open(LOG,">$log") or die "ERROR $0: cannot open log file: $log\n";
$time = localtime;
print LOG "$0  @ARGV\n";
print LOG "\n$0 processing started: $time \n";
print LOG "reference SLC parameter file: $orb1_par\n";
print LOG "bperp text file:  $bpf\n";
print LOG "output ITAB file: $itab\n\n";

print "$0 processing started: $time \n";
print "SLC_tab table SLC images and parameter files: $ARGV[0]\n";
print "reference SLC parameter file: $orb1_par\n";
print "bperp textfile:    $bpf\n";
print "output itab file:  $itab\n\n";

if ($ARGV[4] ne "-"){
  $imode = $ARGV[4];
  ($imode == 0 || $imode == 1) or die "ERROR: invalid itab_mode: $imode\n";
  if($imode == 0){
    print LOG "itab generation mode: single SLC reference scene\n";
    print "itab generation mode: single SLC reference scene\n";
  }
  if($imode == 1){
    print LOG "itab_mode: all possible pairs\n";
    print "itab_mode: all possible pairs\n";
  }
}

if(($#ARGV + 1) >= 6){
  if ($ARGV[5] ne "-"){
    $plt_flg = $ARGV[5];
  }
  print "plot mode: $plt_flg\n";
}

if(($#ARGV + 1) >= 7){
  if ($ARGV[6] ne "-"){
    $bperp_min = $ARGV[6];
    print LOG "minimum bperp magnitude (m): $bperp_min\n";
    print "minimum bperp magnitude (m): $bperp_min\n";
  }
}

if(($#ARGV + 1) >= 8){
  if ($ARGV[7] ne "-"){
    $bperp_max = $ARGV[7];
    print LOG "maximum bperp magnitude (m): $bperp_max\n";
    print "maximum bperp magnitude (m): $bperp_max\n";
  }
}
if(($#ARGV + 1) >= 9){
  if ($ARGV[8] ne "-"){
    $dt_min = $ARGV[8];
    print LOG "minimum delta_T (days): $dt_min\n";
    print "minimum delta_T (days): $dt_min\n";
  }
}
if(($#ARGV + 1) >= 10){
  if ($ARGV[9] ne "-"){
    $dt_max = $ARGV[9];
    print LOG "maximum delta_T (days): $dt_max\n";
    print "maximum delta_T (days): $dt_max\n";
  }
}
if(($#ARGV + 1) >= 11){
  if ($ARGV[10] ne "-"){
    $n_max = $ARGV[10];
    print LOG "maximum delta_n: $n_max\n";
    print "maximum delta_n: $n_max\n";
  }
}

@date = extract_param($orb1_par,"date:");	   #determine image date
@stime1 = extract_param($orb1_par, "start_time:"); #determine start time of SLC 1
$st1 = $stime1[1];
print "reference frame @date\n";
$d1 = sprintf "%4d%02d%02d",$date[1],$date[2],$date[3];
@sensor = extract_param($orb1_par,"sensor:");	#determine image format
$s1 = $sensor[1];
print "reference sensor: $s1\n\n";
$mjd1 = julday($date[1],$date[2],$date[3]) + $st1/86400.;
printf "reference SLC start time (s): %.5f\n",$st1;
printf "reference SLC MJD: %.5f\n",$mjd1;
$plt_title = "\"Radar Image Acquisitions\\nreference SLC: ".$ARGV[1]."\"";

if($plt_flg > 0){
  open (GNUPLOT, "|gnuplot");
  print GNUPLOT <<EOPLOT;
set title $plt_title font 'Verdana,16'
set xdata time
set timefmt "%Y%m%d"
set format x "%d-%b-%y"
set xtics rotate by 300
set grid
set pointsize 1.1
unset key
set size 1,0.85
set tmargin at screen 0.85
set bmargin at screen 0.15
set origin 0,0
set ylabel "Perpendicular Baseline (m)" font 'Verdana,14'
set xtics font 'Verdana,12'
set ytics font 'Verdana,12'
set mytics 5
EOPLOT
}

$i = 1;
LINE: while (<SLC_TAB>) {	#read lines of processing list file
  chomp $_;			#remove new line from record
  next LINE if /^$/; 		#skip blank lines in processing list
  next LINE if /^#/; 		#skip comments in processing list
  @fields = split;		#extract the scene identifier and other parameters on the line if present

  $orb2_par = $fields[1];
  @date = extract_param($orb2_par,"date:");	#determine image 2 date
  @stime2 = extract_param($orb2_par, "start_time:"); #determine start time of image 2
  $st2 = $stime2[1];
  @sensor = extract_param($orb2_par,"sensor:");	#determine image format
  $d2[$i-1] = sprintf "%4d%02d%02d",$date[1],$date[2],$date[3];

  $mjd2 = julday($date[1],$date[2],$date[3]) + $st2/86400.;
  $delta[$i-1] = $mjd2 - $mjd1;
#  printf "SLC2 start time: %12.4f  MJD: %.5f\n", $st2, $mjd2;
  if($delta[$i-1] == 0){
    $rn = $i;
  }
  $borb = `base_orbit $orb1_par $orb2_par base.out`;
  print LOG "base_orbit $orb1_par $orb2_par base.out\n";
  print LOG $borb, "\n";
  @lines = split /^/m, $borb;	#split base_orbit output into lines

  foreach $l1 (@lines) {	#examine each line
    chomp $l1;
    @tok = split /:/, $l1;	#break line on semicolon and check if "base perpendicular" is in th line
    if ($tok[0] =~ m/baseline perpendicular/) {
      $bperp[$i-1] = $tok[1];	#save bperp values
      printf "%3d  ref.: %s  %s  Bperp:%10.4f  delta_T:%12.5f\n",$i, $d1, $d2[$i-1], $bperp[$i-1], $delta[$i-1];
    }
  }
  $i++;
}

$nim = $i-1;
print "\nnumber of SLC images: $nim\n";
$ymin = min(@bperp);
$ymax = max(@bperp);
$xmin = min(@d2);
$xmax = max(@d2);
print "starting date: $xmin  ending date: $xmax\n";
$jd_min = julday(substr($xmin,0,4),  substr($xmin,4,2), substr($xmin,6,2)) - 100;
$jd_max = julday(substr($xmax,0,4),  substr($xmax,4,2), substr($xmax,6,2)) + 100;
#print "starting julian day: $jd_min\n";
#print "ending julian day:   $jd_max\n";

($yr0, $m0, $day0) = caldat($jd_min);	#determine x axis plot range
($yr1, $m1, $day1) = caldat($jd_max);
$xmin = sprintf "%4d%02d%02d",$yr0,$m0,$day0;
$xmax = sprintf "%4d%02d%02d",$yr1,$m1,$day1;
print "starting plot date: $xmin  ending plot date: $xmax\n";

$dy = $ymax - $ymin;	#set y axis plot range
$ymin = $ymin - 0.05*$dy;
$ymax = $ymax + 0.05*$dy;

if($plt_flg == 1){
  printf GNUPLOT "set terminal png size 800,600 font 'Verdana,14'\n";
  printf GNUPLOT "set output \"$bp_plot\"\n";
  print LOG "bperp plot file: $bp_plot \n";
  print "bperp plot file:  $bp_plot\n";
}

if ($plt_flg == 2){	#plot to screen
  my $os = $^O;
  print "computert OS: $os\n";
  my $gterm = $ENV{'GNUTERM'};
  if (undefined $gterm) {
    $gterm = 'wxt';
  }
  print "Gnuplot terminal: $gterm\n";
  printf GNUPLOT "set terminal $gterm size 800 600 font 'Verdana,12'\n";
}

if($plt_flg >= 1){
  printf GNUPLOT "set multiplot\n";
  printf GNUPLOT "set xrange [\"$xmin\":\"$xmax\"]\n";
  printf GNUPLOT "set yrange [$ymin:$ymax]\n";
  printf GNUPLOT "plot '-' using 1:2:3 with labels right offset 0,.2 point tc lt 3 font 'Verdana,10'\n";
  for ($j=0; $j < $nim; $j++) {
    printf GNUPLOT "%s  %10.4f  %3d\n", $d2[$j], $bperp[$j], $j+1;
  }
  printf GNUPLOT "e\n";
  printf GNUPLOT "set xrange [GPVAL_X_MIN:GPVAL_X_MAX]\n";
  printf GNUPLOT "set yrange [GPVAL_Y_MIN:GPVAL_Y_MAX]\n";
  printf GNUPLOT "set label \"Bperp    min: %5.1f  max: %6.1f\" at graph .55,.10 tc lt 1 font 'Verdana,12'\n",$bperp_min, $bperp_max;
  printf GNUPLOT "set label \"Delta_T min: %5.1f  max: %6.1f\" at graph .55,.05 tc lt 1 font 'Verdana,12'\n",$dt_min, $dt_max;
}

$nc = 0;
if($imode == 0){	#single reference
  $bpmag_ave = 0;
  $bp_ave = 0;
  $nbp = 0;
  print "single reference mode\n";

  open(ITAB,">$itab") or die "ERROR $0: cannot output itab file: $itab\n";
  print "reference orbit number: $rn\n";
  if($plt_flg >= 1){
    printf GNUPLOT "plot '-' using 1:2:3 with labels right offset 0,.2 point lt 1 tc lt 1 font 'Verdana,10'\n";
  }
  for($j=0; $j < $nim; $j++){
    $jj = $j+1;
    if((($bperp_min == -1) || ($bperp_min <= abs($bperp[$j]))) &&
       (($bperp_max == -1) || ($bperp_max >= abs($bperp[$j]))) &&
       (($dt_max == -1) || ($dt_max >= abs($delta[$j]))) &&
       (($n_max == -1) || ($n_max >= abs($jj-$rn))) &&
       ($dt_min <= abs($delta[$j]))) {
      $nc++;
      printf ITAB "%4d %4d %4d  1\n",$rn,$jj,$nc;
      printf BP "%4d  $d2[$rn-1] $d2[$j]   %10.4f %10.5f %10.4f %10.4f\n",$nc, $bperp[$j], $delta[$j], $bperp[$j], $bperp[$k];
      printf "%4d  ref: $d2[$rn-1] $d2[$j]  Bperp:%10.4f  delta_T:%10.5f  MJD1:%11.4f  MJD2:%11.4f  Bperp1:%10.4f  Bperp2:%10.4f\n",$nc,$bperp[$j],$delta[$j],0.0,$delta[$j], $bperp[$k], $bperp[$j];
      if($plt_flg >= 1){
        printf GNUPLOT "%s  %10.4f  %3d\n", $d2[$j], $bperp[$j], $j+1;
      }
      if(abs($bperp[$j]) > 0.0){
        $bpmag_ave += abs($bperp[$j]);
        $bp_ave += $bperp[$j];
        $nbp++;
      }
    }
  }

  if($plt_flg >= 1){	#plot lines from reference to each image that meets the criterea
    printf GNUPLOT "e\n";
    printf GNUPLOT "plot '-' using 1:2 with lines lw 0.5 lt 1\n";

    for($j=0; $j < $nim; $j++){
      $jj = $j+1;
      if((($bperp_min == -1) || ($bperp_min <= abs($bperp[$j]))) &&
	 (($bperp_max == -1) || ($bperp_max >= abs($bperp[$j]))) &&
	 (($dt_max == -1) || ($dt_max >= abs($delta[$j]))) &&
	 (($n_max == -1) || ($n_max >= abs($jj - $rn))) &&
	 ($dt_min <= abs($delta[$j]))) {
	printf GNUPLOT "%s  %10.4f\n", $d1, 0.0;
        printf GNUPLOT "%s  %10.4f\n\n", $d2[$j], $bperp[$j];
      }
    }
    printf GNUPLOT "e\n";
    printf GNUPLOT "unset multiplot\n";
  }

  printf "\nnumber of baselines in the ITAB with length > 0: $nbp\n";
  if($nbp > 0){
    $bpmag_ave /= $nbp;
    $bp_ave /= $nbp;
    printf "average bperp magnitude (m): %8.3f \n",$bpmag_ave;
    printf "average bperp (m): %8.3f \n",$bp_ave;
  }
}

if($imode == 1){			#all unique pairs
  print "*** checking all SLC pairs ***\n";
  $bp_ave = 0.0;
  $nbp = 0;
  open(ITAB,">$itab") or die "ERROR $0: cannot output ITAB file: $itab\n";
  if($plt_flg >= 1){
    printf GNUPLOT "plot '-' using 1:2:3 with labels right offset 0,.2 point lt 1 tc lt 1 font 'Verdana,10'\n";
  }
  for($k = 0; $k < $nim-1; $k++){
    $kk = $k + 1;

    for($j = $k + 1; $j < $nim; $j++){		#don't check with self
      $jj = $j + 1;
      $bp = $bperp[$j] - $bperp[$k];
      $dt = $delta[$j] - $delta[$k];

      if((($bperp_min == -1) || ($bperp_min <= abs($bp))) &&
         (($bperp_max == -1) || ($bperp_max >= abs($bp))) &&
         (($n_max == -1) || ($n_max >= abs($jj-$kk))) &&
         (($dt_max == -1) || ($dt_max >= abs($dt))) &&
	  ($dt_min <= abs($dt))){

        $nc++;
        printf ITAB "%4d %4d %4d  1\n",$kk,$jj,$nc;
        printf BP "%4d  $d2[$k] $d2[$j]  %11.5f %11.4f %11.5f %11.5f %10.4f %10.4f\n",$nc,$bperp[$j] - $bperp[$k], $delta[$j] - $delta[$k], $delta[$k], $delta[$j], $bperp[$k], $bperp[$j];
        printf "%4d  ref.: $d2[$k] $d2[$j] %3d  %3d Bperp:%10.4f  delta_T:%10.5f  MJD1:%11.4f  MJD2:%11.4f  Bperp1:%10.4f  Bperp2:%10.4f\n",$nc, $kk, $jj, $bperp[$j]-$bperp[$k], $delta[$j]-$delta[$k], $delta[$k], $delta[$j], $bperp[$k], $bperp[$j];
        if($plt_flg >= 1){
	  printf GNUPLOT "%s  %10.4f  %3d\n", $d2[$k], $bperp[$k], $kk;
	  printf GNUPLOT "%s  %10.4f  %3d\n", $d2[$j], $bperp[$j], $jj;
        }

        if(abs($bperp[$j] - $bperp[$k]) > 0.0){
	  $bp_ave += abs($bperp[$j] - $bperp[$k]);
	  $nbp++;
	}
      }
    }
  }

  if($plt_flg >= 1){	#plot lines from reference to each image that meets the criterea
    printf GNUPLOT "e\n";
    printf GNUPLOT "plot '-' using 1:2 with lines lw 0.5 lt 1\n";

    for($k = 0; $k < $nim-1; $k++){
      $kk = $k + 1;

      for($j = $k + 1; $j < $nim; $j++){		#don't check with self
	$jj = $j + 1;
	$bp = $bperp[$j] - $bperp[$k];
	$dt = $delta[$j] - $delta[$k];

	if((($bperp_min == -1) || ($bperp_min <= abs($bp))) &&
	   (($bperp_max == -1) || ($bperp_max >= abs($bp))) &&
	   (($n_max == -1) || ($n_max >= abs($jj-$kk))) &&
	   (($dt_max == -1) || ($dt_max >= abs($dt))) &&
	    ($dt_min <= abs($dt))){
	  printf GNUPLOT "%s  %10.4f\n", $d2[$k], $bperp[$k];
	  printf GNUPLOT "%s  %10.4f\n\n", $d2[$j], $bperp[$j];
	}
      }
    }
    printf GNUPLOT "e\n";
    printf GNUPLOT "unset multiplot\n";
  }
 
  printf "\nnumber of baselines in the ITAB with length > 0: $nbp\n";
  if($nbp > 0){
    $bp_ave /= $nbp;
    printf "average bperp magnitude (m): %8.3f \n",$bp_ave;
  }
}

print"\nnumber of SLC pairs that meet the bperp_max and delta_T_max criteria: $nc\n";
if ($plt_flg == 1){
  print "bperp plot file:  $bp_plot\n";
}

printf "output bperp_file: $bpf\n";

$time = localtime;
print LOG "\n$0 processing completed: $time \n";
print "\n$0 processing completed: $time \n";
exit 0;

sub execute{
  my ($command, $log) = @_;
  if (-e $log){open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";}
  else {open(LOG,">$log") or die "ERROR $0: cannot open log file: $log\n";}
  LOG->autoflush;
  print "$command\n";
  print LOG ("\n${command}\n");
  close LOG;
  $exit = system("$command 1>> $log");
  $exit == 0 or die "ERROR $0: non-zero exit status: $command\n"
}

sub extract_param{
  my ($infile,$keyword) = @_;
  open(PAR_IN,$infile) || die "ERROR $0: cannot open parameter file: $infile\n";

  while(<PAR_IN>){
    chomp;
    @tokens = split;
    if($tokens[0] eq $keyword){close PAR_IN; return @tokens;}
  }
  close PAR_IN;
  die "ERROR $0: keyword $keyword not found in file: $infile\n\n";
}

sub julday {
##################################
# Julian Day from Gregorian Date #
##################################

use integer;
my ( $y, $m, $d) = @_;
#print "y,m,d: $y $m $d\n";
my $jd = (1461 * ($y + 4800 +
($m - 14) / 12)) / 4 + (367 *
($m - 2 - 12 * (($m - 14) / 12)))
/ 12 - (3 * (($y + 4900 + ($m - 14)
/ 12) / 100)) / 4 + $d - 32075;
#print "jd: $jd\n";
return $jd;
};

sub caldat {
##################################
# Gregorian Date from Julian Day #
##################################

use integer;
my ($d,$i,$j,$jd,$l,$m,$n,$y);
$jd = shift;
$l = $jd + 68569;
$n = ( 4 * $l ) / 146097;
$l = $l - ( 146097 * $n + 3 ) / 4;
$i = ( 4000 * ( $l + 1 ) ) / 1461001;
$l = $l - ( 1461 * $i ) / 4 + 31;
$j = ( 80 * $l ) / 2447;
$d = $l - ( 2447 * $j ) / 80;
$l = $j / 11;
$m = $j + 2 - ( 12 * $l );
$y = 100 * ( $n - 49 ) + $i + $l;
return ($y,$m,$d);
};

=pod

Algorithms taken from Peter Meyers web site
hermetic.magnet.ch/cal_stud/cal_art.htm

Test numbers:
Given $y = 1; $m = 1; $d = 1;   $jd == 1721426;
Given $y = -4713; $m = 11; $d = 24;  $jd == 0;
Given $y = 1858; $m = 11; $d = 16;   $jd == 2400000;

$jd = julday(1858,11,16);
printf "julian day 16-Nov-1858: %d\n", $jd;
printf "date of julian day %d: %d %d %d\n\n",$jd, caldat($jd);

=cut