#!/usr/bin/perl
use FileHandle;
use File::Basename;
use Getopt::Long;

$mflag = 0;
$rsf = 1;
$azsf = 0;
$ccw = 3;
$ccw_min = 3;	#min and max window sizes for cc_ad()
$ccw_max = 9;
$r_samp0 = -1;
$az_line0 = -1;
$scale = 1.0;	#display scale factor
$exp = 0.35;	#display exponent
$ccflg = 0;	#use cc_wave() to estimate correlation
$reforb = "-";	#use phase_sim_orb to estimate correlation, provide geometric reference scene
$tdxflg = 0;	#Tandem-X single-pass mode
$orbflg = 0;	#use phase_sim_orb() to simulate the topographic phase
$aflg = 0;	#add phase from residual baseline simulation

#command line options parsing, the rest of the arguments are passed to @ARGV
GetOptions ('e=f' => \$exp, 's=f' => \$scale, 'c' => \$ccflg, 'o' => \$orbflg, 'r=s' => \$reforb, 't' =>\$tdxflg, '-a' => \$aflg);

if (($#ARGV + 1) < 10){die <<EOS ;}
*** $0
*** Copyright 2015, Gamma Remote Sensing, v4.4 11-Apr-2015 clw ***
*** Calculate 2D diff. interferograms using RSLC_tab, itab, DEM and deformation rate in Range-Doppler Coordinates (RDC) ***

usage: $0 <RSLC_tab> <itab> <bflag> <DEM_rdc> <def> <mli> <mli_dir> <diff_dir> <rlks> <azlks> <cc_win> [rsflag] [azflag] [mflag] [r_samp] [az_line] 

    RSLC_tab  (input) two column list of resampled SLC filenames and SLC parameter filenames (including paths) (text)
                1. SLC filename  (includes path)
                2. SLC parameter filename (includes path)
    itab      (input) table associating interferogram stack records with pairs of SLC stack records (text)
                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 for IPTA processing (0:not used  1:used)
    bflag     baseline flag:
                0: use initial baseline in the baseline file
                1: use precision baseline in the baseline file
    DEM_rdc   (input) terrain height in radar range-Doppler coordinates (meters, float, enter - for none)
    def       (input) deformation rate (m/year) (enter - for none)
    mli       (input) reference MLI image with same rlks and azlks as the interferogram used for background
    mli_dir   directory containing MLI images of the coregistered SLCs                      
    diff_dir  differential interferograms after subtraction of simulated unwrapped phase   
    rlks      range looks for interferogram generation
    azlks     azimuth looks for interferogram generation
    cc_win    correlation estimation window size in pixels with linear weighting (default: $ccw)
    rsflag    range spectral shift filtering flag:
                0: off
                1: on (default)
    azflag    azimuth common-band filtering flag:
                0: off (default)
                1: on 
    mflag     initial baseline estimation and refinement flag:
                0: orbit state vector data (default)
                1: orbit state vector data + baseline refinement using 2-D FFT
                2: use existing baseline file
    r_samp    range pixel offset to center of FFT window for baseline refinement using 2-D FFT (default: image center)
    az_line   azimuth line offset to center of FFT window for baseline refinement using 2-D FFT (default: image center)
    -s scale  (option) set image display scale factor (default: $scale)   
    -e exp    (option) set image display exponent (default: $exp)
    -c        (option) use of cc_ad rather than cc_wave to estimate correlation
    -o        (option) simulate interferogram phase using orbit state vectors with phase_sim_orb
    -a        (option) when using orbit state vectors for simulation of the phase, add phase calculated from residual baseline obtaine by base_ls
    -r        (option) SLC parameter file from the scene used for coregistration, required by phase_sim_orb
    -t        (option) Tandem-X single-pass interferometry mode
EOS

print "*** Copyright 2015 Gamma Remote Sensing, v4.4 11-Apr-2015 clw ***\n";
print "*** Calculate 2D diff. interferograms using RSLC_tab, itab, DEM and deformation rate in Range-Doppler Coordinates (RDC) ***\n";

$itab     = $ARGV[1];
$bflag    = $ARGV[2];
($bflag >= 0 && $bflag <= 1) or die "\nERROR $0: invalid value for baseline flag: $bflag\n";
$dem      = $ARGV[3];
$def      = $ARGV[4];
$rmli     = $ARGV[5];
$rmli_dir = $ARGV[6];
$diff_dir = $ARGV[7];
$rlks     = $ARGV[8];
$azlks    = $ARGV[9];

#flip the sign of $tdxflg 
if ($tdxflg == 1){
 $tdxflg = 0;
}
else {
 $tdxflg = 1;
}

($rlks >= 1 && $rlks <= 32) or die "\nERROR $0: number of range looks is outside of valid range : $rlks\n";
($azlks >= 1 && $azlks <= 128) or die "\nERROR $0: number of azimuth looks is outside of valid range : $azlks\n";

if($#ARGV >= 10){$ccw = $ARGV[10];}
($ccw >= 1 && $ccw <= 15) or  die "\nERROR $0: invalid value for corrrelation window parameter: $ccw\n"; 

if($#ARGV >= 11){$rsf = $ARGV[11];}
($rsf <= 1 && $rsf >=0) or  die "\nERROR $0: invalid value for range spectral shift flag: $rsf\n"; 

if($#ARGV >= 12){$azf = $ARGV[12];}
($azf <= 1 && $azsf >=0) or  die "\nERROR $0: invalid value for azimuth spectral shift flag: $azsf\n"; 

if($#ARGV >= 13){
  $mflag = $ARGV[13];
  ($mflag == 0 || $mflag == 1 || $mflag== 2) or die "\nERROR $0: invalid value for mflag baseline estimation command line parameter: $mflag\n"; 
  if($#ARGV >= 14){
    $r_samp = $ARGV[14];
    $r_samp =~ /\d+$/ && $r_samp > 0 or die "ERROR $0: invalid range sample offset r_samp: $rlks\n";
  }
  if($#ARGV >= 15){
    $az_line = $ARGV[15];
    $az_line =~ /\d+$/ && $az_line > 0 or die "ERROR $0: invalid azimuth line offset az-line: $az_line\n";
  } 
}

open(SLC_TAB, "<$ARGV[0]") or die "\nERROR $0: table of resampled SLC images and parameter files does not exist: $ARGV[0]\n\n";
print "table of resampled SLC images and parameter files: $ARGV[0]\n"; 

open(ITAB, "<$itab") or die "ERROR $0: cannot open itab file: $itab\n";
-e $dem || $dem =~ /-/ or die "ERROR $0: DEM data in radar range-Doppler coordinates does not exist: $dem\n";
-e $def || $def =~ /-/ or die "ERROR $0: deformation rate data in radar range-Doppler coordinates does not exist: $def\n";
-e $rmli or die "ERROR $0: rmli reference image file does not exist: $rmli\n";
-d $rmli_dir or die "ERROR $0: directory containing rmli images does not exist: $rmli_dir\n";

#check if cc_ad available
if($ccflg == 0){	#try running cc_ad to see if it is available
  if(length(`cc_ad`) == 0){
     print "WARNING: program cc_ad adaptive correlation estimation not available\n";
     $ccflg = 1;
  }
}
if($ccflg == 1){
  print "NOTE: using cc_ad to estimate correlation\n";
}

unless (-d $diff_dir) {
  print "creating diff directory: $diff_dir\n";
  (mkdir $diff_dir) or die "ERROR $0: non-zero exit status, could not create differential interferogram directory $diff_dir\n";
}

$diff_par_in = "$diff_dir/diff_par.in";
$off_par_in = "$diff_dir/off_par.in";
$log = "$diff_dir/mk_diff_2d.log"; 
 
print "$0 log file: $log\n";
open(LOG,">$log") or die "ERROR $0: cannot open log file: $log\n";
LOG->autoflush; 
print LOG "$0 @ARGV\n\n";	#print command line that was used
$rlks =~ /\d+$/ && $rlks > 0 or die "ERROR $0: invalid number of range looks: $rlks\n";
$azlks =~ /\d+$/ && $azlks > 0 or die "ERROR $0: invalid number of azimuth looks: $azlks\n";

$time = localtime;
print "start: $time  log file: $log\n\n";
print "itab:                     $itab\n";
print "DEM data:                 $dem\n";
print "deformation data:         $def\n";
print "baseline flag:            $bflag\n";
print "baseline estimation mode: $mflag\n";
print "diff. directory:          $diff_dir\n";
print "interf. looks range: $rlks  azimuth: $azlks\n";
print "correlation window size: $ccw\n";
print "display scale factor: $scale  exponent: $exp\n";
 
if ($orbflg == 1){
  print "simulating phase from orbits with phase_sim_orb\n\n";
  print "geometric reference orbit: $reforb\n";
}
else {
  print("simulating phase using baseline model with phase_sim\n\n");
}

if ($rsf == 1){print "NOTE: applying range spectral shift filtering\n";}
else {print "NOTE: not applying range spectral shift filtering\n";}

if ($azf == 1){print "NOTE: applying azimuth common band filtering\n";}
else {print "NOTE: not applying azimuth common band filtering\n";}
print "\n";

print LOG "$0  @ARGV\n";
print LOG "start: $time  log file: $log\n";
print LOG "itab:                   $itab\n";
print LOG "DEM data:               $dem\n";
print LOG "deformation data:       $def\n";
print LOG "baseline flag:          $bflag\n";
print LOG "diff. directory:        $diff_dir\n";
print LOG "interf. looks range: $rlks  azimuth: $azlks\n";
print LOG "correlation window size: $ccw\n";
print LOG "display scale factor: $scale  exponent: $exp\n\n";

if($rsf == 1){print LOG "NOTE: applying range spectral shift filtering\n";}
else {print LOG "NOTE: not applying range spectral shift filtering\n";}

if ($azf == 1){print LOG  "NOTE: applying azimuth common band filtering\n";}
else {print LOG "NOTE: not applying azimuth common band filtering\n";}
print LOG "\n";

$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
  $slc[$i] = $fields[0];
  $slc_par[$i] = $fields[1];
  -e $slc[$i] or die "ERROR $0: SLC image file does not exist: $slc[$i]\n";     
  -e $slc_par[$i] or die "ERROR $0: ISP image parameter file does not exist: $slc_par[$i]\n";
  printf "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  printf LOG "%2d  $slc[$i]  $slc_par[$i]\n",$i;
  $i++;
} 
close LOG;

$n = 1;
LINE: while (<ITAB>) {		#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
  print "ITAB line: @fields\n";
  next LINE if ($fields[3] == 0);  #skip if the itab column flag = 0
  $slc1 = $slc[$fields[0]];
  $slc1_par = $slc_par[$fields[0]];
  print "SLC1 filename: $slc1\n";
  my ($s1, $dir1, $ext1) = fileparse($slc1, qr/\.[^.]*/); #extension is the last bit after the last period   

  $slc2 = $slc[$fields[1]];
  $slc2_par = $slc_par[$fields[1]];
  print "SLC2 filename: $slc2\n";
  my ($s2, $dir2, $ext2) = fileparse($slc2, qr/\.[^.]*/); #extension is the last bit after the last period   
 
  $fn = $diff_dir."/".$s1."_".$s2;
  $int = $fn.".int";
  $off = $fn.".off";		
  printf "\n*** %d  SLC1: $slc1  SLC2: $slc2  OFF: $off ***\n",$n;
  
  open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";
  printf LOG "\n*** %d  SLC1: $slc1  SLC2: $slc2  OFF: $off ***\n",$n;
  close LOG;

  unless(-e $off){
    printf "NOTE: $off from SLC resampling does not exist, creating new ISP offset parameter file\n";
    $off = $fn.".off";
    open(MPAR, ">$off_par_in") or die "ERROR $0: cannot create file with inputs to create_diff_par: $off_par_in\n\n";
    print MPAR "$s1 $s2\n\n\n\n\n\n\n";
    close MPAR;  
    execute("create_offset $slc1_par $slc2_par $off 1 $rlks $azlks < $off_par_in", $log);
    print "\n";
  }  
  
  $diff = $fn.".diff";
  $diff_par = $fn.".diff_par";
  $sim = $fn.".sim";
  $sim_orb =  $fn.".sim_orb";
  $sim_orb_tmp =  $fn.".sim_orb_tmp";
  $int = $fn.".int";
  $base = $fn.".base";
  $base2 = $fn.".base2";
  $mli1 = $rmli_dir."/".$s1.".rmli";
  $mli2 = $rmli_dir."/".$s2.".rmli";
  unless (-e $mli1){
    $mli1 = $rmli_dir."/".$s1.".mli";
    $mli2 = $rmli_dir."/".$s2.".mli";
    -e $mli1 or die "MLI-1 does not exist: $mli1\n";
    -e $mli2 or die "MLI-2 does not exist: $mli2\n";
  }
  $cc = $fn.".cc";
   
  @slc1_date = extract_param($slc1_par,"date:");
  @slc2_date = extract_param($slc2_par,"date:");
  $jd1 = julday($slc1_date[1], $slc1_date[2], $slc1_date[3]);
  $jd2 = julday($slc2_date[1], $slc2_date[2], $slc2_date[3]);
  print "SLC-1 date: $slc1_date   julian day: $jd1\n";
  print "SLC-2 date: $slc2_date   julian day: $jd2\n";
  $delta_t = $jd2 - $jd1;
  print "interferogram delta_t (days): $delta_t\n";
  @slc1_freq = extract_param($slc1_par,"radar_frequency:");
  @slc2_freq = extract_param($slc2_par,"radar_frequency:");
  print "SLC-1 radar frequency: @slc1_freq\n";
  print "SLC-2 radar frequency: @slc2_freq\n\n";

  $f2flg = 0;
  if($slc1_freq[1] != $slc2_freq[1]){
    printf "SLC scenes have different radar carrier frequencies!\n"; 
    $f2flg = 1
  }
 
  if ($bflag == 0 && $mflag != 2) {
    execute("base_orbit $slc1_par $slc2_par $base",$log);
  }

  open(MPAR, "> $diff_par_in") or die "ERROR $0: cannot create file with inputs to create_diff_par: $diff_par_in\n\n";
  print MPAR "$s1 $s2\n\n\n\n\n";
  close MPAR;
  execute("create_diff_par $off - $diff_par 0 < $diff_par_in", $log);
  print "\n";

  if ($orbflg == 1) {	#simulate phase using orbit state vectors with possible refinement
    execute("phase_sim_orb $slc1_par $slc2_par $off $dem $sim_orb $reforb $def $delta_t $tdxflg", $log);
    @b1 = extract_param($base, "precision_baseline(TCN):");	#baseline file must exist created when running mk_diff_2d

    if ($b1[2] != 0.0 or  $b1[3] != 0.0) {	#check if baseline is non-zero
      if($f2flg == 1){
        execute("phase_sim $slc1_par $off $base $dem $sim 0 1 $def $delta_t $slc2_par", $log);
      }
      else {
        execute("phase_sim $slc1_par $off $base $dem $sim 0 1 $def $delta_t $tdxflg", $log);  
      }
      rename $sim_orb, $sim_orb_tmp;
      execute("sub_phase $sim_orb_tmp $sim $diff_par $sim_orb 0 1",$log); #add phases together
      unlink $sim_orb_tmp;
    }
  }
  else {
    if($f2flg == 1){
      execute("phase_sim $slc1_par $off $base $dem $sim 0 $bflag $def $delta_t $slc2_par", $log);
    }
    else {
      execute("phase_sim $slc1_par $off $base $dem $sim 0 $bflag $def $delta_t $tdxflg", $log);  
    }
  }

  @b1 = extract_param($base,"initial_baseline(TCN):");
  if ($b1[2] != 0.0 or  $b1[3] != 0.0){		#check if non-zero baseline
    if ($orbflg == 0) {	#phase simulated from baseline model in $base and DEM
      execute("SLC_diff_intf $slc1 $slc2 $slc1_par $slc2_par $off $sim $diff $rlks $azlks $rsf $azf", $log);
    }
    else {	#phase simulated using orbit state vectors and DEM
      execute("SLC_diff_intf $slc1 $slc2 $slc1_par $slc2_par $off $sim_orb $diff $rlks $azlks $rsf $azf", $log);
    }
  }
  else {	#zero baseline
    execute("SLC_intf $slc1 $slc2 $slc1_par $slc2_par $off $diff $rlks $azlks 0 - $rsf $azf",$log); 
  }

# calculate interferogram using SLC_intf
# execute("SLC_intf $slc1 $slc2 $slc1_par $slc2_par $off $int $rlks $azlks 0 - $rsf $azf",$log); 
# execute("sub_phase $int $sim $diff_par $diff 1",$log);	#subtract simulated phase
  @width = extract_param($off, "interferogram_width:");
  @lines = extract_param($off, "interferogram_azimuth_lines:");

  if($mflag > 0 && $mflag < 2){				
    $nrfft = 1024;
    if ($width[1] < 1024){
      $nrfft = 256*(int $width[1]/256);
    }
    $nazfft = 1024;
    if ($lines[1] < 1024){
      $nazfft = 256*(int $lines[1]/256);
    }
    $r_samp = int ($width[1]/2);
    $az_line = int ($lines[1]/2);
    if($r_samp0 >= 0){$r_samp = $r_samp0;}	#check for command line input values of r_samp and az_line 
    if($az_line0 >= 0){$az_line = $az_line0};

    execute("base_init $slc1_par $slc2_par $off $diff $base2 4 $nrfft $nazfft $r_samp $az_line",$log); #2D FFT 

    @b1 = extract_param($base,"initial_baseline(TCN):");
    @b1r = extract_param($base,"initial_baseline_rate:");

    @b2 = extract_param($base2,"initial_baseline(TCN):");
    @b2r = extract_param($base2,"initial_baseline_rate:");

    for($j = 1; $j <= 3; $j += 1){			#add baseline correction
      $b3[$j] = $b1[$j] + $b2[$j];
      $b3r[$j] = $b1r[$j] + $b2r[$j];
    }

    $sv1 = sprintf "%14.7f %14.7f %14.7f",$b3[1], $b3[2], $b3[3]; 
    $sv2 = sprintf "%14.7f %14.7f %14.7f",$b3r[1], $b3r[2], $b3r[3]; 
    execute("set_value $base $base \"initial_baseline(TCN)\"  \"$sv1\"",$log); #write the updated values into the baseline file
    execute("set_value $base $base \"initial_baseline_rate\" \"$sv2\"",$log);
    
    if($f2flg == 1){	#resimulate phse with updated baseline!
      execute("phase_sim $slc1_par $off $base $dem $sim 0 $bflag $def $delta_t $slc2_par", $log);
    }
    else {
      execute("phase_sim $slc1_par $off $base $dem $sim 0 $bflag $def $delta_t", $log);  
    }    
    execute("SLC_diff_intf $slc1 $slc2 $slc1_par $slc2_par $off $sim $diff $rlks $azlks $rsf $azf", $log);
  } 

# generate raster file of the differential interferogram
  execute("rasmph_pwr $diff $rmli $width[1] 1 1 0 1 1 $scale $exp",$log);
# estimate correlation coefficients
  if ($ccflg == 0){  
    execute("cc_wave $diff $mli1 $mli2 $cc $width[1] $ccw $ccw 1", $log);
  } 
  else {
    eval{
      execute("cc_ad $diff $mli1 $mli2 - - $cc $width[1] $ccw_min $ccw_max 0", $log);
    } or do {
      $ccflg = 0;
      print "WARNING: program cc_ad not available, using cc_wave to estimate correlation. cc_ad is in the LAT package of the Gamma software\n";
      execute("cc_wave $diff $mli1 $mli2 $cc $width[1] $ccw $ccw 1", $log);
    };
  }
  execute("rascc $cc $rmli $width[1] 1 1 0  1 1 .1 .9 $scale $exp", $log);
  $n++;
}

$time = localtime;
print "\n $0 processing end: $time\n";
open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log\n";
print LOG "\n $0 processing end: $time\n";
exit 0;

sub execute{
  my ($command, $log) = @_;  
  if (-e $log){open(LOG,">>$log") or die "ERROR $0: cannot open log file: $log  command: $command\n";}
  else {open(LOG,">$log") or die "ERROR $0 : cannot open log file: $log  command: $command\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";
}

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