#!/usr/bin/perl
# 
# Terrain Geocoding Script 21-Mar-2006 clw
# v1.9 add ability to geocode to equiangular projection
# v2.0 added search radius for interpolation 
# v2.3 added parameters for offset estimation
# v2.4 modified geocode_back to use spline_log interpolation for MLI images
# v2.5 patch size for initial offset is now adjusted for small images
# v2.6 display of DEM in SAR coordinates now generated using rashgt
# v2.9 corrected error in calculation of $rad_max
# v3.0 20091217 now closing LOG file before execution of Gamma programs
# v3.1 20110223 added SCH projection
# v3.2 20110402 added scale factor and exponent for output geocoded image raster file generation
# v3.3 20110428 changed display scale factor and exponent to be command line options with the -s and -e keys
# v3.4 20120328 added capability to support different raster file types sun raster and BPM
# v3.7 20121030 added output for u,v, inc, and psi angles and ls_map, pixel area
# v3.8 20150409 remove dependence on Linux/Unix

use FileHandle;
use Getopt::Long;

$gc_map = "gc_map";	#gc_map program
$gc_map_grd = "gc_map_grd"; #gc_map_grd program
$ls_mode = 2;		#output lookup table values in regions of layover, shadow, or DEM gaps:
			#0:set to (0.,0.)  1:linear interpolation across these regions (default)  2:actual value  3:nn-thinned
$r_ovr = 16.0;		#range oversampling factor for nn-thinned layover/shadow mode
$frame = 8;		#number of DEM pixels to add around area covered by SAR image

#refinement and geocoding parameters
$interp_mode = 2; 	#geocode DEM interpolation mode: (0: 1/dist, 1: nearest neighbor, 2: SQR(1/dist), 3: const, 4: Gaussian)
$gc_n_ovr = 2;		#geocode interpolation oversampling factor
$rad_max = "-"; 	#search radius for geocode program
$rlks = 1;		#default number of range looks for initial offset estimate 
$azlks = 1;		#default number of azimuth looks for initial offset estimate
$thres = 10.0;		#threshold for accepting an initial offset estimate
$offr = 0;		#initial range offset for offset estimation
$offaz = 0;		#initial azimuth offset for offset estimation
$rpatch = 256;		#range patch size for lookup table refinement
$azpatch = 256;		#azimuth patch size for lookup table refinement
$nrp = 20;		#number of range patches
$nazp = 20;		#number of azimuth patches
$href = 100.0;		#reference height for ellipsoid geocoding
$rinit = 0;		#initial range offset for lookup table refinement
$azinit = 0;		#initial range offset for lookup table refinement
$rpos = "-";		#default azimuth and range positions for determining initial offset in lookup table refinement
$azpos = "-";
$pflg = 1;		#copy initial offsets to polynomials in call to init_offsetm
$psize = 1024;		#default patch size for initial offset
$psize_min = 64;	#minimum patch size for initial offset estimate
$hscale = 100.0;	#height scale for DEM display is 100 m/cycle
$ext = "ras";		#default image extension Sun Raster

$scd = 0.6;		#image display intensity scale factor
$expd = 0.4;		#image display exponent
#command line options parsing, the rest of the arguments are passed to @ARGV
GetOptions ('e=f' => \$expd, 's=f' => \$scd, 'b' => \$ext_flg);
if ($ext_flg == 1){
  $ext = "bmp";
}

if (($#ARGV + 1) < 10){die <<EOS ;}
*** $0
*** Copyright 2015, Gamma Remote Sensing, v3.8 9-Apr-2015 clw ***
*** Terrain geocoding of SAR images with lookup table refinement and resample DEM to SAR Range-Doppler Coordinates (RDC) ***

usage: $0 <MLI> <MLI_par> <DEM> <DEM_par> <DEM_seg> <DEM_seg_par> <GEO_dir> <scene_id> <post> <mode> [ls_mode] [r_ovr] [n_ovr] [rad_max] [rlks] [azlks] [thres] [rpos] [azpos] [roff] [azoff] [r_patch] [az_patch] [nr] [naz]

    MLI          (input) MLI SAR image (including path)
    MLI_par      (input) ISP image parameter file of the MLI image (including path)
    DEM          (input) DEM in desired output projection (including path)
    DEM_par      (input) DEM parameter file (including path)
    DEM_seg      (output) DEM segment for output image products (including path)
    DEM_seg_par  (input/output) DEM parameter file for output image products (including path), regenerated each time 
    GEO_dir      directory for output images, lookup tables and DEM products
    scene_id     scene name to identify output files
    post         output image sample spacing in meters or degrees for geographic (EQA) projection
    mode         mk_geo processing mode:    
                   0: generate initial lookup table, simulated SAR image, and DEM segment parameters
                   1: generate initial lookup table and simulated SAR image using existing DEM segment to determine image bounds
                   2: measure initial offset between simulated SAR image and actual SAR image
                   3: perform refinement of lookup table by offset measurement with respect to the simulated SAR image
                   4: update lookup table and produce terrain geocoded SAR image and DEM in SAR range-Doppler coordinates (RDC) 
                   5: ellipsoid geocoding of the SAR image without a DEM: generate new DEM segment parameters          
                   6: ellipsoid geocoding of the the SAR image without a DEM: use existing DEM segment parameters          
    ls_mode      algorithm selection in gc_map for regions of layover, shadow, or DEM gaps:
                   0: set to (0.,0.)
                   1: linear interpolation across these regions 
                   2: use actual value (default)
                   3: nearest neighbor thinned (nn-thinned)
    r_ovr        range over-sampling parameter for gc_map parameter ls_mode = 3 (nn-thinned) (default: $r_ovr)
    n_ovr        interpolation oversampling factor in geocode (default = 2.0)
    rad_max      maximum interpolation search radius in geocode (enter - for default: 4*n_ovr)
    
    NOTE: n_ovr and rad_max are parameters used by the program geocode for transformation of the simulated image 
          and DEM into SAR geometry. The parameters rlks, azlks, thres, rpos, azpos, roff, azoff are used  
	  for estimation of the initial offset of the SAR image with respected to the simulated SAR image. 
		 		 
    rlks         number of range looks for the initial offset estimate (default: $rlks)
    azlks        number of azimuth looks for the initial offset estimate (default: $azlks)
    thres        SNR threshold for offset estimates (default: $thres)
    rpos         range position for initial offset (enter - for default)
    azpos        azimuth position for initial offset (enter - for default)
    roff         initial range offset estimate (enter - for current value in DIFF_par file)
    azoff        initial azimuth offset estimate (enter - for current value in DIFF_par file)
    r_patch      range patch size for offset estimation (default: $rpatch samples)
    az_patch     azimuth patch size for offset estimation (default: $azpatch lines)
    nr           number of range patches for offset estimation (default: $nrp)
    naz          number of azimuth patches for offset estimation (default: $nazp)
    -s scale     (option) set image display scale factor (default: $scd)   
    -e exp       (option) set image display exponent (default: $expd)
    -b           (option) set raster image type to BMP (default: $ext)
EOS

$mli =      $ARGV[0];
$mli_par =  $ARGV[1];
$dem =      $ARGV[2];
$dem_par =  $ARGV[3];
$dem_seg =  $ARGV[4];
$dem_seg_par = $ARGV[5];
$geo_dir =  $ARGV[6];
$scene_id = $ARGV[7];
$post =     $ARGV[8];
$mode =     $ARGV[9];

($mode >= 0 && $mode <= 5) or die "\nERROR $0: invalid procesing mode for $0: $mode\n\n";
if($#ARGV >= 10){$ls_mode = $ARGV[10];}
if($#ARGV >= 11){$r_ovr =  $ARGV[11];}
if($#ARGV >= 12){$gc_n_ovr =  $ARGV[12];}
if($#ARGV >= 13){$rad_max =$ARGV[13];}
if($#ARGV >= 14){$rlks =   $ARGV[14];}
if($#ARGV >= 15){$azlks =  $ARGV[15];}
if($#ARGV >= 16){$thres =  $ARGV[16];}
if($#ARGV >= 17){$rpos =   $ARGV[17];}
if($#ARGV >= 18){$azpos =  $ARGV[18];}
if($#ARGV >= 19){$rinit =  $ARGV[19];}
if($#ARGV >= 20){$azinit = $ARGV[20];}
if($#ARGV >= 21){$rpatch = $ARGV[21];}
if($#ARGV >= 22){$azpatch = $ARGV[22];}
if($#ARGV >= 23){$nrp  = $ARGV[23];}
if($#ARGV >= 24){$nazp = $ARGV[24];}

if($rad_max eq "-"){
  $rad_max = 4 * $gc_n_ovr;
}

-e $mli or die "\nERROR $0: MLI image does not exist: $mli\n";
-e $mli_par or die "\nERROR $0: MLI image parameter file does not exist: $mli_par\n";

#if($mode < 4){
#  -e $dem or die "\nERROR $0: input DEM does not exist: $dem\n";
#}
-e $dem_par or die "\nERROR $0: input DEM parameter file does not exist: $dem_par\n";

$mli_ras = "$mli.$ext";

unless (-d $geo_dir){
  print "creating directory for output geocoded image products: $geo_dir\n";
  (mkdir $geo_dir) or die "ERROR $0: non-zero exit status, could not create geocoded data directory $geo_dir\n";
}

$log = "$geo_dir/mk_geo"."_$mode".".log"; 
open(LOG,">$log") or die "ERROR $0: cannot open log file: $log\n";
LOG->autoflush;  
    
$time = localtime;
print LOG "$0 @ARGV\n\n";
print LOG "*** $0 processing start: $time  mode: $mode ***\n";
print "*** $0 processing start: $time  mode: $mode ***\n";
print "log file: $log\n";
print "raster image extension: $ext\n";
@mli_width = extract_param($mli_par,"range_samples:");
@mli_lines = extract_param($mli_par,"azimuth_lines:"); 
@mli_rps = extract_param($mli_par,"range_pixel_spacing:"); 
@mli_azps = extract_param($mli_par,"azimuth_pixel_spacing:"); 

printf "\nMLI image: $mli\n";
print  "MLI raster image: $mli_ras\n";
printf "MLI parameter file: $mli_par\n";
printf "MLI width: %d   lines: %d\n",$mli_width[1],$mli_lines[1];
printf "MLI pixel spacing (m) range: $mli_rps[1]  azimuth: $mli_azps[1]\n";
print "display scale factor: $scd  exponent: $expd\n\n";
printf LOG "MLI image: $mli\n";
printf LOG "MLI raster image: $mli_ras\n";
printf LOG "MLI parameter file: $mli_par\n";
printf LOG "MLI width: %d   lines: %d\n",$mli_width[1],$mli_lines[1];
printf LOG "MLI pixel spacing (m) range: $mli_rps[1]  azimuth: $mli_azps[1]\n";
printf LOG "display scale factor: $scd  exponent: $expd\n\n";

unless ( -e $mli_ras) {
  close LOG;
  execute("raspwr $mli $mli_width[1] 1 0 1 1 $scd $expd 1 $mli_ras", $log);
  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
}

@dem_width = extract_param($dem_par,"width:");
@dem_lines = extract_param($dem_par,"nlines:"); 
@dem_proj = extract_param($dem_par,"DEM_projection:"); 

printf "input DEM: $dem\n";
printf "DEM parameter file: $dem_par\n";printf "DEM projection: $dem_proj[1]\n";
printf "DEM width: %d   lines: %d\n",$dem_width[1], $dem_lines[1];
printf "DEM projection: %s\n",$dem_proj[1];
printf LOG "\ninput DEM: $dem\n";
printf LOG "DEM parameter file: $dem_par\n";
printf LOG "DEM width: %d   lines: %d\n",$dem_width[1],$dem_lines[1];
 
if ($dem_proj[1] eq "EQA"){
  @dem_post_lat = extract_param($dem_par,"post_lat:"); 
  @dem_post_lon = extract_param($dem_par,"post_lon:");

  $sc_north = abs($dem_post_lat[1]/$post);
  $sc_east = abs($dem_post_lon[1]/$post);
    
  printf "DEM posting Northing (deg): %e\n",$dem_post_lat[1];
  printf "DEM posting Easting (deg):  %e\n",$dem_post_lon[1];
  printf "output posting (deg): %e\n",$post;
  printf LOG "DEM posting Northing (deg): %e\n",$dem_post_lat[1];
  printf LOG "DEM posting Easting (deg):  %e\n",$dem_post_lon[1];
  printf LOG "output posting (deg): %e\n",$post;
} 
else { 
  if($dem_proj[1] eq "SCH"){
    @dem_post_s = extract_param($dem_par,"SCH_post_s:"); 
    @dem_post_c = extract_param($dem_par,"SCH_post_c:");

    printf "DEM posting S (m) %f\n",$dem_post_s[1];
    printf "DEM posting C (m):  %f\n",$dem_post_c[1];
    printf "output posting (m): %f\n",$post;
    printf LOG "DEM posting S (m): %f\n",$dem_post_s[1];
    printf LOG "DEM posting C (m):  %f\n",$dem_post_c[1];
    printf LOG "output posting (m): %f\n",$post;
    $sc_north = abs($dem_post_s[1]/$post);
    $sc_east = abs($dem_post_c[1]/$post);
  }
  else {
    @dem_post_north = extract_param($dem_par,"post_north:"); 
    @dem_post_east = extract_param($dem_par,"post_east:");

    printf "DEM posting Northing (m) %f\n",$dem_post_north[1];
    printf "DEM posting Easting (m):  %f\n",$dem_post_east[1];
    printf "output posting (m): %f\n",$post;
    printf LOG "DEM posting Northing (m): %f\n",$dem_post_north[1];
    printf LOG "DEM posting Easting (m):  %f\n",$dem_post_east[1];
    printf LOG "output posting (m): %f\n",$post;
    $sc_north = abs($dem_post_north[1]/$post);
    $sc_east = abs($dem_post_east[1]/$post);
  }
}

$gc0 = "$geo_dir/$scene_id"."_0.map_to_rdc";
$gc1 = "$geo_dir/$scene_id"."_1.map_to_rdc";
$sim0 = "$geo_dir/$scene_id"."_0.sim";
$sim0_ras = "$sim0.$ext";
$sim_rdc0 = "$geo_dir/$scene_id"."_0.sim_rdc";
$sim_rdc1 = "$geo_dir/$scene_id"."_1.sim_rdc";
$sim_rdc0_ras = "$sim_rdc0.$ext";
$sim_rdc1_ras = "$sim_rdc1.$ext";
$diff_par = "$geo_dir/$scene_id".".diff_par";
$offs = "$geo_dir/$scene_id".".offs";
$snr = "$geo_dir/$scene_id".".snr";
$coffs = "$geo_dir/$scene_id".".coffs";
$dem_seg_rdc = "$geo_dir/$scene_id"."_dem.rdc";
$dem_seg_rdc_ras = "$geo_dir/$scene_id"."_dem.rdc.$ext";
$mli_map = "$geo_dir/$scene_id"."_map.mli";
$mli_map_ras = "$geo_dir/$scene_id"."_map.mli.$ext";
$mli_gec = "$geo_dir/$scene_id"."_gec.mli";
$mli_gec_ras = "$geo_dir/$scene_id"."_gec.mli.$ext";
$diff_par_in = "$diff_par".".in";

$ls_map = "$geo_dir/$scene_id"."_0.ls_map";
$inc_map = "$geo_dir/$scene_id"."_0.inc_map";
$pix_map = "$geo_dir/$scene_id"."_0.pix_map";
$u_map = "$geo_dir/$scene_id"."_0.u_map";
$v_map = "$geo_dir/$scene_id"."_0.v_map";
$psi_map = "$geo_dir/$scene_id"."_0.psi_map";

printf "\noutput DEM segment: $dem_seg\n";
printf "DEM segment parameter file: $dem_seg_par\n";
printf "DEM interpolation factor northing: $sc_north  easting: $sc_east\n";
printf "output terrain geocoded product directory: $geo_dir\n";

printf LOG "\noutput DEM segment: $dem_seg\n";
printf LOG "DEM segment parameter file: $dem_seg_par\n";
printf LOG "DEM interpolation factor northing: $sc_north  easting: $sc_east\n";
printf LOG "DEM interpolation factor northing: $sc_north  easting: $sc_east\n";
printf LOG "output terrain geocoded product directory: $geo_dir\n";
close LOG;

if ($mode == 0 || $mode == 1){
  if($mode == 0){
    printf "\n$0 mode $mode: generate initial lookup table, simulated SAR image, and DEM segment parameters\n";
    printf "delete existing DEM parameter file\n";
    execute("rm -f $dem_seg $dem_seg_par", $log);
  }  
  if($mode == 1){
    printf "\n$0 mode $mode: generate initial lookup table and simulated SAR image using existing DEM segment\n";
    -e $dem_seg_par or die "ERROR: DEM segment parameter file does not exist $dem_seg_par\n";
    printf "use existing DEM segment parameter file: $dem_seg_par\n";
  }
  
  ($ls_mode == 0 || $ls_mode == 1 || $ls_mode == 2 || $ls_mode == 3) or die "\nERROR $0: invalid ls_mode for gc_map: $ls_mode\n\n";
  if($ls_mode == 0){printf "Layover, shadow, or DEM gaps gc_map layover-shadow mode: 0 (set to 0)\n\n";} 
  if($ls_mode == 1){printf "Layover, shadow, or DEM gaps gc_map layover-shadow mode: 1 (linear interpolation)\n\n";} 
  if($ls_mode == 2){printf "Layover, shadow, or DEM gaps gc_map layover-shadow mode: 2 (actual value)\n\n";} 
  if($ls_mode == 3){
    printf "Layover, shadow, or DEM gaps region gc_map mode: 3 (nn-thinned)\n";
    printf "range oversampling factor for nn_thinned arrays: $r_ovr\n\n";
  }
  @image_geometry = extract_param($mli_par,"image_geometry:");

  if ($image_geometry[1] eq "GROUND_RANGE"){
    execute("$gc_map_grd $mli_par $dem_par $dem $dem_seg_par $dem_seg $gc0 $sc_north $sc_east $sim0 $u_map $v_map $inc_map $psi_map $pix_map $ls_map $frame $ls_mode $r_ovr", $log);
  }
  else {
    execute("$gc_map $mli_par - $dem_par $dem $dem_seg_par $dem_seg $gc0 $sc_north $sc_east $sim0 $u_map $v_map $inc_map $psi_map $pix_map $ls_map $frame $ls_mode $r_ovr", $log);
  }

#get DEM segment width and nlines
  @dem_seg_width = extract_param($dem_seg_par,"width:");
  @dem_seg_lines = extract_param($dem_seg_par,"nlines:"); 

  execute("geocode $gc0 $sim0 $dem_seg_width[1] $sim_rdc0 $mli_width[1] $mli_lines[1] $interp_mode 0 1 1 $gc_n_ovr $rad_max", $log); 
  execute("raspwr $sim_rdc0 $mli_width[1] 1 0 1 1 1. .35 1 $sim_rdc0_ras", $log);
  execute("dis2ras $sim_rdc0_ras $mli_ras&", $log);

  print "\nsimulated SAR image in map coordinates: $sim0\n";
  print "simulated SAR image in range-Doppler coordinates (RDC): $sim_rdc0\n";
  print "raster image of simulated SAR image in RDC: $sim_rdc0_ras\n";
  print "DEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n\n";
  $time = localtime;
  print "\n*** $0 mode: $mode completed  $time ***\n";

  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
  print LOG "\nsimulated SAR image in map coordinates: $sim0\n";
  print LOG "simulated SAR image in range-Doppler coordinates: $sim_rdc0\n";
  print LOG "raster image of simulated SAR image in RDC: $sim_rdc0_ras\n";
  print LOG "DEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print LOG "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n\n";
  print LOG "\n*** $0 mode: $mode completed  $time ***\n";
  close LOG;
}

if ($mode == 2){
  open(DPAR, "> $diff_par_in") or die "ERROR $0: cannot create file with inputs to create_diff_par: $diff_par_in\n\n";
  print DPAR "$scene_id\n0 0\n$nrp $nazp\n$rpatch $azpatch\n$thres\n";
  close DPAR;
  print "creating DIFF/GEO parameter file: $diff_par\n";
  if(-e $diff_par){
    unlink $diff_par;
  }
  execute("create_diff_par $mli_par - $diff_par 1 < $diff_par_in", $log);
  -e $sim_rdc0 or die "\nERROR $0: simulated SAR image does not exist: $sim_rdc0\n";
    
  if($rpos eq "-"){
    $rpos = $mli_width[1]/2;
  }
  if($azpos eq "-"){
    $azpos = $mli_lines[1]/2;
  }
  print "\n\ninitial offset patch center range: $rpos  azimuth: $azpos\n"; 
  if($rinit ne "-"){
    execute("set_value $diff_par $diff_par initial_range_offset $rinit",$log); 
  }
  if($azinit ne "-"){
    execute("set_value $diff_par $diff_par initial_azimuth_offset $azinit",$log); 
  }
  print "initial offset range (pixels): $rinit  azimuth: $azinit\n";
  print "number of looks for initial offset estimate range: $rlks  azimuth: $azlks\n\n";  

  while (($mli_lines[1] < ($azpos + $psize/2)) || ($azpos - $psize/2 < 0)){
    $psize = $psize/2;
  }
 
  #determine range patch size constraint
  while (($mli_width[1] < ($rpos + $psize/2)) || ($rpos - $psize/2 < 0)){
    $psize = $psize/2;
  } 
  print "patch size for determining the initial offset (range, azimuth): $psize, $psize\n";

  ($psize >= $psize_min) or die "\nERROR: insufficient number of lines in the MLI image: $mli_lines[1]\n";
  @offr = extract_param($diff_par,"initial_range_offset:");
  @offaz = extract_param($diff_par,"initial_azimuth_offset:");
  print ("initial range offset: $offr[1]   initial azimuth offset: $offaz[1]\n");

  execute("init_offsetm $sim_rdc0 $mli $diff_par $rlks $azlks $rpos $azpos $offr[1] $offaz[1] $thres $psize $pflg", $log); 
  @polyr = extract_param($diff_par,"range_offset_polynomial:");
  @polyaz = extract_param($diff_par,"azimuth_offset_polynomial:");

  print "\ninitial range offset (samples): $polyr[1]   azimuth offset (lines): $polyaz[1]\n";
  print "setting range and azimuth offset polynomial values to initial offsets in file: $diff_par\n";
  print "\n*** $0 mode: $mode completed  $time ***\n";
  
  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
  print LOG "creating DIFF/GEO parameter file: $diff_par\n";
  print LOG "\ninitial offset patch center range: $rpos  azimuth: $azpos\n";
  print LOG "initial offset range (pixels): $rinit  azimuth: $azinit\n";
  print LOG "number of looks for initial offset estimate range: $rlks  azimuth: $azlks\n\n";
  print LOG "\ninitial range offset (samples): $polyr[1]   azimuth offset (lines): $polyaz[1]\n";
  print LOG "setting range and azimuth offset polynomial values to initial offsets in file: $diff_par\n";
  print LOG "\n*** $0 mode: $mode completed  $time ***\n";
  close LOG;
}

if ($mode == 3){
#read initial range and azimuth offsets set in mode 2
  @offr = extract_param($diff_par,"initial_range_offset:");
  @offaz = extract_param($diff_par,"initial_azimuth_offset:");

  unless (-e $diff_par) {			#check if DIFF_par exists, if not create it 
    open(DPAR, "> $diff_par_in") or die "ERROR $0: cannot create file with inputs to create_diff_par: $diff_par_in\n\n";
    print DPAR "$scene_id\n $offr[1] $offaz[1]\n$nrp $nazp\n$rpatch $azpatch\n$thres\n";
    close DPAR;
    execute("create_diff_par $mli_par - $diff_par 1 < $diff_par_in", $log);
  }
  -e $sim_rdc0 or die "\nERROR $0: simulated SAR image does not exist: $sim_rdc0\n";
  -e $diff_par or die "\nERROR $0: DIFF/GEO parameter file does not exist, run mode 2 to create $diff_par\n";

  @nptr = extract_param($diff_par,"offset_estimation_range_samples:");
  @nptaz = extract_param($diff_par,"offset_estimation_azimuth_samples:");
  @pwidth = extract_param($diff_par,"offset_estimation_patch_width:");
  @phgt = extract_param($diff_par,"offset_estimation_patch_height:");

  print "initial offset (pixels) range: $offr[1]  azimuth: $offaz[1]\n";
  printf "number of range patches:   %3d   patch width:  %3d\n", $nptr[1], $pwidth[1];
  printf "number of azimuth patches: %3d   patch height: %3d\n\n", $nptaz[1], $phgt[1];

  execute("offset_pwrm $sim_rdc0 $mli $diff_par $offs $snr", $log);
  execute("offset_fitm $offs $snr $diff_par $coffs - $thres 3",$log);
  print "DIFF/GEO offset parameter file: $diff_par\n";
  execute("cat $diff_par",$log);
  print "\n*** $0 mode: $mode completed  $time ***\n";
  
  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
  printf LOG "number of range patches:   %3d   patch width:  %3d\n", $nptr[1], $pwidth[1];
  printf LOG "number of azimuth patches: %3d   patch height: %3d\n\n", $nptaz[1], $phgt[1];
  print LOG "DIFF/GEO offset parameter file: $diff_par\n";
  print LOG "\n*** $0 mode: $mode completed  $time ***\n";
  close LOG;
}

if ($mode == 4){
  -e $gc0 or die "\nERROR: initial lookup table does not exist: $gc0\n";
  printf "lookup table (MAP --> SAR RDC): $gc0\n\n";

  @dem_seg_width = extract_param($dem_seg_par,"width:");
  @dem_seg_lines = extract_param($dem_seg_par,"nlines:"); 

#refinement of the gecoding lookup table
  execute("gc_map_fine $gc0 $dem_seg_width[1] $diff_par $gc1 1",$log);
  execute("geocode $gc1 $sim0 $dem_seg_width[1] $sim_rdc1 $mli_width[1] $mli_lines[1] $interp_mode 0 1 1 $gc_n_ovr $rad_max", $log); 
  execute("raspwr $sim_rdc1 $mli_width[1] 1 0 1 1 0.8 0.4 1 $sim_rdc1_ras", $log);

#transform DEM into SAR RDC
  execute("geocode $gc1 $dem_seg $dem_seg_width[1] $dem_seg_rdc $mli_width[1] $mli_lines[1] $interp_mode 0 1 1 $gc_n_ovr $rad_max", $log); 
  execute("rashgt $dem_seg_rdc $sim_rdc1 $mli_width[1] 1 1 0 1 1 $hscale 1. .35 1 $dem_seg_rdc_ras", $log);

#terrain geocoded intensity image 
  execute("geocode_back $mli $mli_width[1] $gc1 $mli_map $dem_seg_width[1] $dem_seg_lines[1] 2 0",$log);
  execute("raspwr $mli_map $dem_seg_width[1] 1 0 1 1 $scd $expd 1 $mli_map_ras", $log);

  print "\n";
  execute("disras_dem_par $mli_map_ras $dem_seg_par&", $log);
  sleep 1;
  print "\n";
  execute("dis2ras $sim_rdc1_ras $mli_ras&", $log);		#confirm refinement of lookup table
  sleep 1;
  print "\n";
  execute("disras $dem_seg_rdc_ras&", $log); 
  sleep 1;
  print "\nrefined lookup table (MAP --> SAR RDC): $gc1\n";
  print "DEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n";
  print "terrain geocoded SAR image: $mli_map\n";
  print "terrain geocoded SAR image raster file: $mli_map_ras\n";
  print "terrain geocoded SAR image DEM parameter file: $dem_seg_par\n";
  print "DEM in SAR RDC coordinates: $dem_seg_rdc\n";
  print "DEM in SAR RDC coordinates raster image: $dem_seg_rdc_ras\n";   
  print "\n*** $0 mode: $mode completed  $time ***\n";
    
  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
  print LOG "refined lookup table (MAP --> SAR RDC): $gc1\n";
  print LOG "DEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print LOG "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n";
  print LOG "\nterrain geocoded SAR image: $mli_map\n";
  print LOG "terrain geocoded SAR image raster file: $mli_map_ras\n";
  print LOG "terrain geocoded SAR image DEM parameter file: $dem_seg_par\n";
  print LOG "DEM in SAR RDC coordinates: $dem_seg_rdc\n";
  print LOG "DEM in SAR RDC coordinates raster image: $dem_seg_rdc_ras\n";
  print LOG "\n*** $0 mode: $mode completed  $time ***\n";
  close LOG;
}

if ($mode == 5 || $mode == 6){
  printf "\n$0 mode $mode: ellipsoid geocoding of a SAR image\n";
  if ($mode == 5){
    execute("rm -f $dem_seg_par ", $log);
  }
  execute("gec_map $mli_par - $dem_par $href $dem_seg_par $gc0 $sc_north $sc_east", $log);

#get DEM segment windth and nlines
  @dem_seg_width = extract_param($dem_seg_par,"width:");
  @dem_seg_lines = extract_param($dem_seg_par,"nlines:"); 

  print "\nDEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n\n";

  execute("geocode_back $mli $mli_width[1] $gc0 $mli_gec $dem_seg_width[1] $dem_seg_lines[1] 2 0",$log);
  execute("raspwr $mli_gec $dem_seg_width[1] 1 0 1 1 $scd $expd 1 $mli_gec_ras", $log);
  execute("disras_dem_par $mli_gec_ras $dem_seg_par&",$log);

  print "ellipsoid geocoded SAR image: $mli_gec\n";
  print "ellipsoid geocoded SAR image raster file: $mli_gec_ras\n";
  print "ellipsoid geocoded SAR image DEM parameter file: $dem_seg_par\n";
  $time = localtime;
  print "\n*** $0 mode: $mode completed  $time ***\n";
    
  open(LOG,">>$log") or die "\nERROR $0: cannot open log file: $log\n";
  print LOG "\nDEM segment parameter file for geocoded SAR image: $dem_seg_par\n";
  print LOG "DEM segment width: $dem_seg_width[1]   lines: $dem_seg_lines[1]\n\n";  print LOG "\ngellipsoid geocoded SAR image: $mli_gec\n";
  print LOG "ellipsoid geocoded SAR image raster file: $mli_gec_ras\n";
  print LOG "ellipsoid geocoded SAR image DEM parameter file: $dem_seg_par\n"; 
  print LOG "\n*** $0 mode: $mode completed  $time ***\n";
  close LOG; 
}
exit 0;

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