referencephase.cc

Go to the documentation of this file.

/*
 * Copyright (c) 1999-2003 Bert Kampes
 * Copyright (c) 1999-2003 Delft University of Technology, The Netherlands
 *
 * This file is part of Doris, the Delft o-o radar interferometric software.
 *
 * Doris program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * Doris is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 * Publications that contain results produced by the Doris software should
 * contain an acknowledgment. (For example: The interferometric processing
 * was performed using the freely available Doris software package developed
 * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
 * University of Technology, or include a reference to: Bert Kampes and
 * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
 * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
 *
 */
/****************************************************************
 * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/referencephase.cc,v $ *
 * $Revision: 3.18 $                                            *
 * $Date: 2005/04/16 20:47:55 $                                 *
 * $Author: kampes $                                            *
 *                                                              *
 * -computation flat earth correction.                          *
 * -computation radarcoding dem + interpolation to (1,1) grid   *
 ****************************************************************/

#include "matrixbk.hh"
#include "orbitbk.hh"
#include "slcimage.hh"                  // my slc image class
#include "productinfo.hh"               // my 'products' class
#include "constants.hh"
#include "referencephase.hh"            // proto types
#include "ioroutines.hh"                // error etc.
#include "utilities.hh"                 // isodd; ones(), etc.
#include "coregistration.hh"            // distribute points
#include "conversion.hh"                // lp2xyz etc.
#include "refsystems.hh"                // for dsk WGS84
#include "exceptions.hh"                 // my exceptions class

#include <iomanip>                      // setw only for test..
#include <cstdlib>                      // system
#include <algorithm>                    // max
#include <netinet/in.h>                 // ntohl byteorder x86-HP unix



/****************************************************************
 *    flatearth                                                 *
 *                                                              *
 * Compute polynomial model for 'flat earth' correction.        *
 *  fie(l,p) = sumj=0:d sumk=0:d Ajk l^j p^k (NOT bert 8sept99) *
 * precise orbits are used to compute delta range for Npoints   *
 * after which the polynomial model is fitted (LS).             *
 *                                                              *
 * input:                                                       *
 *  - inputoptions                                              *
 *  - info structs                                              *
 *  - platform data points                                      *
 * output:                                                      *
 *  - void (result to file "scratchresflat")                    *
 *  - coefficients normalized wrt. original window of master    *
 *                                                              *
 *    Bert Kampes, 09-Mar-1999                                  *
 *    Bert Kampes, 26-Oct-1999 normalization of coeff.,         *
 *    dump to logfile: var(unknowns) == diag(inv(AtA))          *
 ****************************************************************/
void flatearth(
        const input_comprefpha &comprefphainput,
        const input_ell        &ellips,
        const slcimage         &master,
        const slcimage         &slave,
        const productinfo      &interferogram,
        orbit                  &masterorbit,
        orbit                  &slaveorbit)
  {
  TRACE_FUNCTION("flatearth (BK 26-Oct-1999)")
  const int32 MAXITER   = 10;
  const real8 CRITERPOS = 1e-6;
  const real8 CRITERTIM = 1e-10;
  char                  dummyline[ONE27];

  INFO << "FLATEARTH: MAXITER: "   << MAXITER   << "; "
                  << "CRITERPOS: " << CRITERPOS << " m; "
                  << "CRITERTIM: " << CRITERTIM << " s";
  INFO.print();

// ______ Normalization factors for polynomial ______
  const real8 minL = master.originalwindow.linelo;
  const real8 maxL = master.originalwindow.linehi;
  const real8 minP = master.originalwindow.pixlo;
  const real8 maxP = master.originalwindow.pixhi;
  INFO << "flatearth: polynomial normalized by factors: "
       << minL << " " << maxL << " " << minP << " " << maxP
       << " to [-2,2]";
  INFO.print();

// ______Handling of input______
  const real8 m_minpi4cdivlam = (-4.0*PI*SOL)/master.wavelength;
  const real8 s_minpi4cdivlam = (-4.0*PI*SOL)/slave.wavelength;
  DEBUG << "master wavelength = " << master.wavelength;
  DEBUG.print();
  DEBUG << "slave  wavelength = " << slave.wavelength;
  DEBUG.print();
  const int32 DEGREE    = comprefphainput.degree;
  const int32 Nunk      = Ncoeffs(DEGREE);              // Number of unknowns
  bool pointsrandom = true;
  if (specified(comprefphainput.ifpositions))
    pointsrandom = false;                       // only use those points



// ______ Distribute points wel distributed over win ______
// ______ or read from ascii file ______
// ______(i,0): line, (i,1): pixel, (i,2) flagfromdisk______
  matrix<uint> Position;
  const uint  Npoints = comprefphainput.Npoints;
  register int32 i,j,k,index;

  if (pointsrandom)                             // no filename specified
    {
    Position = distributepoints(Npoints,interferogram.win);
    }
  else // read from file
    {
    Position.resize(Npoints,3);
    //ifstream ifpos(comprefphainput.ifpositions, ios::in);
    ifstream ifpos;
    openfstream(ifpos,comprefphainput.ifpositions);
    bk_assert(ifpos,comprefphainput.ifpositions,__FILE__,__LINE__);
    uint ll,pp;
    for (i=0; i<Npoints; ++i)
      {
      ifpos >> ll >> pp;
      Position(i,0) = uint(ll);
      Position(i,1) = uint(pp);
      Position(i,2) = uint(1);                // flag from file
      ifpos.getline(dummyline,ONE27,'\n');       // goto next line.
      }
    ifpos.close();

// ______ Check last point ivm. EOL after last position in file ______
    if (Position(Npoints-1,0) == Position(Npoints-2,0) &&
        Position(Npoints-1,1) == Position(Npoints-2,1))
      {
      Position(Npoints-1,0) = uint(.5*(minL + maxL) + 27);      // random
      Position(Npoints-1,1) = uint(.5*(minP + maxP) + 37);      // random
      WARNING << "refpha: there should be no EOL after last point in file: "
           << comprefphainput.ifpositions;
      WARNING.print();
      }

    // ______ Check if points are in overlap ______
    // ______ no check for uniqueness of points ______
    }

  matrix<real8>         y(Npoints,1);                   // observation
  matrix<real8>         A(Npoints,Nunk);                // designmatrix

// ______Check redundancy______
  if (Npoints < Nunk)
    {
    PRINT_ERROR("flatearth: Number of points is smaller than parameters solved for.")
    throw(input_error);
    }



// ======Compute delta r for all points======
  for (i=0; i<Npoints; ++i)
    {
    const real8 line  = Position(i,0);
    const real8 pixel = Position(i,1);
   
// ______ Compute azimuth/range time of this pixel______
    //const real8 m_trange = pix2tr(pixel,master.t_range1,master.rsr2x);
    const real8 m_trange = master.pix2tr(pixel);
   
// ______ Compute xyz of this point P from position in image ______
    cn P;                                       // point, returned by lp2xyz
    lp2xyz(line,pixel,ellips,master,masterorbit,
           P,MAXITER,CRITERPOS);
 
// ______ Compute xyz for slave satelite from P ______
    real8 s_tazi;                               // returned
    real8 s_trange;                             // returned
    xyz2t(s_tazi,s_trange,slave,
          slaveorbit,
          P,MAXITER,CRITERTIM);
   

// ______Compute delta range ~= phase______
// ______ real8 dr = dist(m_possat,pospoint) - dist(s_possat,pospoint);
// ______ real8 phase = -pi4*(dr/LAMBDA);
// ______  dr    == M-S         want if no flatearth M-S - flatearth = M-S-(M-S)=0
// ______  phase == -4pi*dr/lambda == 4pi*(S-M)/lambda
// BK: 24-9: actually defined as: phi = +pi4/lambda * (r1-r2) ???
    // real8 phase = pi4*((dist(s_possat,pospoint)-dist(m_possat,pospoint))/LAMBDA);
    //y(i,0) = pi4*((dist(s_possat,pospoint)-dist(m_possat,pospoint))/LAMBDA);
    //y(i,0) = pi4divlam*(s_possat.dist(pospoint)-m_possat.dist(pospoint));
    //y(i,0) = minpi4cdivlam * (m_trange - s_trange);
    y(i,0) = m_minpi4cdivlam*m_trange - s_minpi4cdivlam*s_trange;
    DEBUG << "l="   << line     << " p="  << pixel 
          << " t1=" << m_trange << " t2=" << s_trange
          << " fe=" << y(i,0) << " [rad]";
    DEBUG.print();


// ______Set up system of equations______
// ______Order unknowns: A00 A10 A01 A20 A11 A02 A30 A21 A12 A03 for degree=3______
// ______  normalize data [-2,2] ______
    //real8 posP                = (Position(i,1) - minP) / (maxP - minP);
    real8 posL          = normalize(line,minL,maxL);
    real8 posP          = normalize(pixel,minP,maxP);
 
    index = 0;
    for (j=0; j<=DEGREE; j++)
      {
      for (k=0; k<=j; k++)
        {
        A(i,index) = pow(posL,real8(j-k)) * pow(posP,real8(k));
        index++;
        }
      }
    }


// ======Compute polynomial for these phases (LS)======
// ______Compute Normalmatrix, rghthandside______
  matrix<real8> N   = matTxmat(A,A);
  matrix<real8> rhs = matTxmat(A,y);


// ______Compute solution______
  matrix<real8> Qx_hat = N;
  choles(Qx_hat);               // Cholesky factorisation normalmatrix
  solvechol(Qx_hat,rhs);        // Estimate of unknowns in rhs
  invertchol(Qx_hat);           // Covariance matrix


// ______Test inverse______
  for (i=0; i<Qx_hat.lines(); i++)
    for (j=0; j<i; j++)
      Qx_hat(j,i) = Qx_hat(i,j);// repair Qx_hat
  const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
  INFO << "flatearth: max(abs(N*inv(N)-I)) = " << maxdev;
  INFO.print();
  /*
  matrix<real8> unity = N * Qx_hat;
  real8 maxdev = abs(unity(0,0)-1);
  for (i=1; i<unity.lines(); i++)
    {
    if (abs(unity(i,i)-1) > maxdev)
      maxdev = abs(unity(i,i)-1);
    for (j=0; j<i-1; j++)
      {
      if (abs(unity(i,j)) > maxdev)
        maxdev = abs(unity(i,j));
      if (abs(unity(j,i)) > maxdev)
        maxdev = abs(unity(j,i));
      }
    }
*/
  if (maxdev > .01) 
    {
    ERROR << "Deviation too large. Decrease degree or number of points?";
    PRINT_ERROR(ERROR.get_str())
    throw(some_error);
    }
  else if (maxdev > .001) 
    {
    WARNING << "Deviation quite large. Decrease degree or number of points?";
    WARNING.print();
    }
  else
    {
    INFO.print("Deviation is OK.");
    }

// ______Some other stuff, scale is ok______
//  matrix<real8> Qy_hat        = A * (matxmatT(Qx_hat,A));
  matrix<real8> y_hat   = A * rhs;
  matrix<real8> e_hat   = y - y_hat;


// ______Overall model test (variance factor)______
// ... ?





// ______ Wrap offset ______ 
// BK 30/9/99 do not do this, later absolute ref. phase is used.
// in s2h rodriguez for example.
// it does not change anything for compinterfero etc.
//  rhs(0,0) = remainder(rhs(0,0),2*PI);



// ______Write results to file______
  ofstream scratchlogfile("scratchlogflat", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"flatearth: scratchlogflat",__FILE__,__LINE__);

  scratchlogfile << "\n\n*******************************************************************"
                 //<< "\n*_Start_flatearth"
                 << "\n*_Start_" << processcontrol[pr_i_comprefpha]
                 << "\nDegree_flat:\t" << DEGREE
                 << "\nEstimated coefficients:\n"
                 << "\nx_hat \tstd:\n";
  for (i=0; i<Nunk; i++)
    scratchlogfile
                 << setiosflags(ios::fixed)
                 << setiosflags(ios::showpoint)
                 << setiosflags(ios::right)
                 << setw(8) << setprecision(4) 
                 <<  rhs(i,0) << " \t" << sqrt(Qx_hat(i,i)) << endl;

  scratchlogfile << "\nCovariance matrix estimated parameters:"
                 << "\n---------------------------------------\n";
  for (i=0; i<Nunk; i++)
    {
    for (j=0; j<Nunk; j++)
      {
      scratchlogfile 
                 << setiosflags(ios::fixed)
                 << setiosflags(ios::showpoint)
                 << setiosflags(ios::right)
                 << setw(8) << setprecision(4) 
                 << Qx_hat(i,j) << " ";
      }
    scratchlogfile << endl;
    }

  scratchlogfile << "\nMaximum deviation N*inv(N):"
                 << setiosflags(ios::scientific)
                 << maxdev
                 << "\nSome more info for each observation:"
                 << "\nline \tpixel \tobs \t\tobs_hat \t\t err_hat\n";
  for (i=0; i<Npoints; i++)
    scratchlogfile <<  Position(i,0) << "\t" <<  Position(i,1) << "\t"
                   << y(i,0) << "\t" << y_hat(i,0) << "\t"
                   << setiosflags(ios::fixed)
                   << setiosflags(ios::showpoint)
                   << setiosflags(ios::right)
                   << setw(8) << setprecision(4) 
                   << e_hat(i,0) << "\n";
  scratchlogfile << "\nMaximum absolute error: \t\t\t"
                 <<  max(abs(e_hat))
                 << "\n*******************************************************************\n";
  scratchlogfile.close();

  ofstream scratchresfile("scratchresflat", ios::out | ios::trunc);
  bk_assert(scratchresfile,"flatearth: scratchresflat",__FILE__,__LINE__);

  scratchresfile.setf(ios::scientific, ios::floatfield);
  scratchresfile.setf(ios::right, ios::adjustfield);
  scratchresfile.precision(8);
  scratchresfile.width(18);

  scratchresfile << "\n\n*******************************************************************"
                 //<< "\n*_Start_flat_earth"
                 << "\n*_Start_" << processcontrol[pr_i_comprefpha]
                 << "\n*******************************************************************"
                 << "\nDegree_flat:\t" << DEGREE
                 << "\nEstimated_coefficients_flatearth:\n";
  int32 coeffL = 0;
  int32 coeffP = 0;
  for (i=0; i<Nunk; i++)
    {
    if (rhs(i,0) < 0.)
      scratchresfile <<         rhs(i,0);
    else
      scratchresfile << " " <<  rhs(i,0);

// ______ Add coefficient number behind value ______
    scratchresfile << " \t" <<  coeffL << " " << coeffP << "\n";
    coeffL--;
    coeffP++;
    if (coeffL == -1)
      {
      coeffL = coeffP;
      coeffP = 0;
      }
    }

  scratchresfile << "*******************************************************************"
                 << "\n* End_" << processcontrol[pr_i_comprefpha] << "_NORMAL"
                 << "\n*******************************************************************\n";


// ====== Compute coordinates of corners of interferogram here ======
// ______ (though better place this somewhere else ....
  real8 phi,lambda,height;                              // rad, returned
  lp2ell(interferogram.win.linelo,interferogram.win.pixlo,
         ellips,master,masterorbit,
         phi,lambda,height,MAXITER,CRITERPOS);
  INFO << "Coordinates of corner interferogram: "
       << interferogram.win.linelo << ", " << interferogram.win.pixlo
       << " = " << rad2deg(phi) << ", " << rad2deg(lambda);
  INFO.print();
  lp2ell(interferogram.win.linehi,interferogram.win.pixlo,
         ellips,master,masterorbit,
         phi,lambda,height,MAXITER,CRITERPOS);
  INFO << "Coordinates of corner interferogram: "
       << interferogram.win.linehi << ", " << interferogram.win.pixlo
       << " = " << rad2deg(phi) << ", " << rad2deg(lambda);
  INFO.print();
  lp2ell(interferogram.win.linelo,interferogram.win.pixhi,
         ellips,master,masterorbit,
         phi,lambda,height,MAXITER,CRITERPOS);
  INFO << "Coordinates of corner interferogram: "
       << interferogram.win.linelo << ", " << interferogram.win.pixhi
       << " = " << rad2deg(phi) << ", " << rad2deg(lambda);
  INFO.print();
  lp2ell(interferogram.win.linehi,interferogram.win.pixhi,
         ellips,master,masterorbit,
         phi,lambda,height,MAXITER,CRITERPOS);
  INFO << "Coordinates of corner interferogram: "
       << interferogram.win.linehi << ", " << interferogram.win.pixhi
       << " = " << rad2deg(phi) << ", " << rad2deg(lambda);
  INFO.print();

// ______Tidy up______
  scratchresfile.close();
  } // END flatearth



/****************************************************************
 *    radarcodedem                                              *
 *                                                              *
 * Compute reference phase based on DEM (gtopo30)               *
 * DEM on equiangular grid (lat/lon) assumed                    *
 * DEM seems stored from North to South                         *
 *                                                              *
 * Implementation:                                              *
 * 0) find out dphi,dlambda for multilooked interferogram       *
 * 1) get dem on equidistant grid into matrix                   *
 * 2) interpolate this DEM dense enough                         *
 * 3) convert interpolated DEM to l,p,refpha                    *
 * 3) for line=first:mlL:last                                   *
 *        pixel=first:mlP:last  (interferogram window)          *
 *      nearest neighbour to l,p from dem                       *
 *                                                              *
 * input:                                                       *
 * output:                                                      *
 *                                                              *
 *    Bert Kampes, 18-Jan-2000                                  *
 *                                                              *
 * small changes wrt formats, more verbose, error checking.     *
 #%// BK 13-Apr-2002                                            *
 * phi is latitude, lambda is longitude                         *
 #%// BK 19-Apr-2002                                            *
 ****************************************************************/
void radarcodedem(
        const input_gen        &generalinput,
        const input_ell        &ellips,
        const input_comprefdem &refdeminput,
        const slcimage         &master,
        const slcimage         &slave,
        const productinfo      &interferogram,
        orbit                  &masterorbit,
        orbit                  &slaveorbit)
  {
  TRACE_FUNCTION("radarcodedem (BK 18-Jan-2000)")
  const int32 MAXITER   = 10;
  const real8 CRITERPOS = 1e-6;
  const real8 CRITERTIM = 1e-10;

  // ====== Handle input ======
  // ______ See gtopo30.HDR file, assume equiangular grid ______
  //  // ______ Equi angular grid with spacing delatlat/lon ______
  //  // ______ Assume increasing lon per row ______
  //    NROWS         6000
  //    NCOLS         4800
  //    NBANDS        1
  //    NBITS         16
  //    NODATA        -9999
  //    ULXMAP        -19.99583333333333
  //    ULYMAP        89.99583333333334
  //    XDIM          0.00833333333333
  //    YDIM          0.00833333333333
  //    format=...
  //const real8 lat0file    = deg2rad(89.99583333333334);// first pix on disk w02090
  //const real8 DEMdeltalat = deg2rad(0.00833333333);   // in radians

  char demfilename[ONE27];
  char refdemofilename[ONE27];

  strcpy (demfilename,refdeminput.firefdem);
  strcpy (refdemofilename,refdeminput.forefdem);
  const real8 lat0file    = refdeminput.demlatleftupper;// first pix on disk w02090
  const real8 lon0file    = refdeminput.demlonleftupper;// first pix on disk
  const real8 DEMdeltalat = refdeminput.demdeltalat;    // in radians
  const real8 DEMdeltalon = refdeminput.demdeltalon;    // in radians
  const int32 numberoflonpixels = refdeminput.demcols;  // NCOLS on file
  const int32 numberoflatpixels = refdeminput.demrows;  // NROWS on file
  const real8 NODATA      =  refdeminput.demnodata;     // (BK 4 may 2001)
  bool onlyrefphasetopo   = !refdeminput.includerefpha; // true: phase DEM w.r.t. ellipsoid
  const real8 EXTRADENSE  =  refdeminput.extradense;    // account for overlay etc.
  const bool outputdem    =  specified(refdeminput.fodem);// if spec. then output
  const bool outputdemi   =  specified(refdeminput.fodemi);// if spec. then output
  // Bert Kampes, 16-Apr-2005: addition output for geocoding using
  // an external DEM; you can dump the height in meters in radarcoordinates
  // then use geocode to obtain lat/lon matrices.
  // extra output for Zbigniew Perski (04/2005)
  // if specified, then output (always do phase too)
  const bool outputrefdemhei   =  specified(refdeminput.forefdemhei);
  const bool do_nearest   = (refdeminput.method==crd_nearest) ? true : false;

  // ______ Some other variables ______
  const int32 mlL = interferogram.multilookL;   // multilookfactor interferogram
  const int32 mlP = interferogram.multilookP;   // multilookfactor interferogram
  const real8 m_min4picdivlam = (-4.0*PI*SOL)/master.wavelength;
  const real8 s_min4picdivlam = (-4.0*PI*SOL)/slave.wavelength;
  DEBUG << "master wavelength = " << master.wavelength;
  DEBUG.print();
  DEBUG << "slave  wavelength = " << slave.wavelength;
  DEBUG.print();

  const real8 latNfile = lat0file-DEMdeltalat*(numberoflatpixels-1); // upper=max. lat value
  const real8 lonNfile = lon0file+DEMdeltalon*(numberoflonpixels-1); // left=min. lon value
  real8 phi, lambda, height;                            // latitude longitude height

  // ______ Extra info ______
  INFO << "DEM input: w/e/s/n:          \t"
       << rad2deg(lon0file) << "/" << rad2deg(lonNfile) << "/"
       << rad2deg(latNfile) << "/" << rad2deg(lat0file);
  INFO.print();


  // ______ Open output file ______
  ofstream refdemofile;
  openfstream(refdemofile,refdemofilename,generalinput.overwrit);
  bk_assert(refdemofile,refdemofilename,__FILE__,__LINE__);
  // if request for height in radar coordinates l,p
  ofstream refdemheiofile;// for Zbigniew Perski
  if (outputrefdemhei==true)
    {
    openfstream(refdemheiofile,refdeminput.forefdemhei,generalinput.overwrit);
    bk_assert(refdemheiofile,refdeminput.forefdemhei,__FILE__,__LINE__);
    }

  // ______ Get corners of interferogram (approx) to select DEM ______
  // ______ in radians (if height were zero)______
  lp2ell(interferogram.win.linelo, interferogram.win.pixlo,
         ellips, master, masterorbit,
         phi, lambda, height);          // returned
  real8 phil0p0    = phi;
  real8 lambdal0p0 = lambda;

  lp2ell(interferogram.win.linehi, interferogram.win.pixlo,
         ellips, master, masterorbit,
         phi, lambda, height);          // returned
  real8 philNp0    = phi;
  real8 lambdalNp0 = lambda;

  lp2ell(interferogram.win.linehi, interferogram.win.pixhi,
         ellips, master, masterorbit,
         phi, lambda, height);          // returned
  real8 philNpN    = phi;
  real8 lambdalNpN = lambda;

  lp2ell(interferogram.win.linelo, interferogram.win.pixhi,
         ellips, master, masterorbit,
         phi, lambda, height);          // returned
  real8 phil0pN    = phi;
  real8 lambdal0pN = lambda;

  // find out spacing of interferogram, use last pixel for densest gridding.
  lp2ell(interferogram.win.linelo+mlL, interferogram.win.pixhi-mlP,
         ellips, master, masterorbit,
         phi, lambda, height);          // returned
  const real8 dphi_interferogram    = abs(phi-phil0pN);
  const real8 dlambda_interferogram = abs(lambda-lambdal0pN);

  // ______ Select DEM values based on rectangle outside l,p border ______
  real8 phimin    = min(min(min(phil0p0,philNp0),philNpN),phil0pN);
  real8 phimax    = max(max(max(phil0p0,philNp0),philNpN),phil0pN);
  real8 lambdamin = min(min(min(lambdal0p0,lambdalNp0),lambdalNpN),lambdal0pN);
  real8 lambdamax = max(max(max(lambdal0p0,lambdalNp0),lambdalNpN),lambdal0pN);

  // ______ a little bit extra at edges to be sure ______ 
  const real8 extrafarrange = (1.5*DEMdeltalon + deg2rad(1.0/25.0));
  const real8 extralatitude = (1.5*DEMdeltalon + deg2rad(1.0/25.0));
  phimin    -= extralatitude;
  phimax    += extralatitude;
  lambdamax += extrafarrange;
  lambdamin -= extrafarrange;

  // ______ Extra info ______
  INFO << "DEM input required: w/e/s/n: \t"
       << rad2deg(lambdamin) << "/" << rad2deg(lambdamax) << "/"
       << rad2deg(phimin)    << "/" << rad2deg(phimax);
  INFO.print();
  INFO << "For window (l0,lN,p0,pN):    \t"
       << interferogram.win.linelo << " "
       << interferogram.win.linehi << " "
       << interferogram.win.pixlo << " "
       << interferogram.win.pixhi;
  INFO.print();


  // ______ Check corners of DEM ______
  // check if DEM is appropriate for interferogram
  // DEM should at least partially cover IFG
  // note: phi is [90:-90]
  if (phimax <= latNfile)// DEM is more north than IFG
    {
    ERROR << "IFG outside DEM: most South latitude: " << rad2deg(latNfile)
         << " [deg]; IFG requires: " << rad2deg(phimax) 
         << " [deg]";
    PRINT_ERROR(ERROR.get_str())
    throw(some_error);
    }
  // DEM is more south than IFG
  if (phimin >= lat0file)// largest latitude at first line of file
    {
    ERROR << "IFG outside DEM: most North latitude: " << rad2deg(lat0file)
         << " [deg]; IFG requires: " << rad2deg(phimax)
         << " [deg]";
    PRINT_ERROR(ERROR.get_str())
    throw(some_error);
    }
  if (lambdamax <= lon0file)
    {
    ERROR << "IFG outside DEM: most West longitude: " << rad2deg(lon0file)
         << " [deg]; IFG window requires: " << rad2deg(lambdamax)
         << " [deg]";
    PRINT_ERROR(ERROR.get_str())
    throw(some_error);
    }
  if (lambdamin >= lonNfile)
    {
    ERROR << "IFG outside DEM: most East longitude: " << rad2deg(lonNfile)
         << " [deg]; IFG window requires: " << rad2deg(lambdamin)
         << " [deg]";
    PRINT_ERROR(ERROR.get_str())
    throw(some_error);
    }


  // ______ Line/pixel of first point in original master coordinates ______
  // ______ maybe this should be changed to be x+(ml/2) ?? but depends on
  // ______ definition of range_to_first_bin is to center or start..
  // Bert Kampes, 08-Apr-2005: chose center by adding ml/2
  const int32 numlinesml    = interferogram.win.lines()  / mlL;
  const int32 numpixelsml   = interferogram.win.pixels() / mlP;
  const real8 veryfirstline = real8(interferogram.win.linelo) + real8(mlL)/2.0;
                                //(real8(mlL)-1.0) / 2.0;
  const real8 firstpixel    = real8(interferogram.win.pixlo)  + real8(mlP)/2.0;
                                //(real8(mlP)-1.0) / 2.0;
  const real8 lastpixel     = firstpixel + real8((numpixelsml-1)*mlP);

  // ______ Compute oversamplefactor for interpolation of DEM ______
  const int32 DEMphioversamplefactor =                  // e.g. 4.2 -> 5
    int32(ceil(EXTRADENSE * DEMdeltalat / dphi_interferogram));// >=1
  const int32 DEMlambdaoversamplefactor =
    int32(ceil(EXTRADENSE * DEMdeltalon / dlambda_interferogram));// >=1
  if (DEMphioversamplefactor > 20)
    {
    WARNING << "Resolution IFG latitude much higher than DEM: "
         << rad2deg(dphi_interferogram) << " [deg] <--> "
         << rad2deg(DEMdeltalat) << " [deg].";
    WARNING.print();
    }
  if (DEMlambdaoversamplefactor > 20)
    {
    WARNING << "Resolution IFG longitude much higher than DEM: "
         << rad2deg(dlambda_interferogram) << " [deg] <--> "
         << rad2deg(DEMdeltalon) << " [deg].";
    WARNING.print();
    }

  // ______ Extra info ______
  INFO << "Approximate longitude resolution IFG: "
       << rad2deg(dlambda_interferogram) << " [deg].";
  INFO.print();
  INFO << "Approximate latitude resolution IFG:  "
       << rad2deg(dphi_interferogram) << " [deg].";
  INFO.print();
  INFO << "DEM oversampled with factors (lat/lon): "
       << DEMphioversamplefactor << " " << DEMlambdaoversamplefactor;
  INFO.print();

  // ______ Added some variables for verboser output BK 13-Apr-2002 ______
  int32 numNOneighbor   = 0;// number of no neigbor points. in buffer
  int32 totalNOneighbor = 0;// number of no neigbor points.
  int32 numvalid        = 0;// number of good values, not NODATA in buffer
  int32 numNODATA       = 0;// number of NODATA values in buffer
  real8 meancroppedDEM  = 0.0;// to detect byte order problems, formats
  real8 min_input_dem   =  100000.0;// stats
  real8 max_input_dem   = -100000.0;// stats


  // ______ Buffer variables ______
  // memory usage is approximately 4arrays*8Byte*ofsize(DEMinterp)
  // and size of DEMinterp = sizeDEM*oversamplingfactors
  // approximate this a bit by saying DEMinterp
  // the input DEM is cropped a little bit larger in far range,
  // account for this here. by multiplying with a factor.
  const real8 BUFFERMEMSIZE    = generalinput.memory;// Bytes
  const real8 numelements_line_DEMI = 
    ((lambdamax-lambdamin)/DEMdeltalat) * DEMlambdaoversamplefactor;
  // ______ 4 arrays required of this size, of 8B doubles ______
  const real8 numlines_possible_DEMi =
    BUFFERMEMSIZE / (4*8*numelements_line_DEMI);
  // ______ bufferlines are of output IFG ______
  int32 bufferlines = 
    int32(ceil(numlines_possible_DEMi/DEMphioversamplefactor));
  if (bufferlines>numlinesml) bufferlines=numlinesml;
  int32 numfullbuffers = numlinesml / bufferlines;
  int32 restlines      = numlinesml % bufferlines;
  // ______ Avoid a very small last buffer (max. 20% of normal size) ______
  if (real8(restlines)/real8(bufferlines) < 0.20)
    if (restlines!=0)
      bufferlines = int32((numlinesml/numfullbuffers)); 
  // ______ And compute the buffersizes finally ______
  numfullbuffers = numlinesml / bufferlines;
  restlines      = numlinesml % bufferlines;
  const int32 EXTRABUFFER = (restlines == 0) ? 0 : 1;   // smaller last buffer
 
  // ______ Extra info ______
  INFO << "Computatation in: " << numfullbuffers << " buffers of " << bufferlines
       << " lines and " << EXTRABUFFER << " extra buffer of "
       << restlines << " lines.";
  INFO.print();
  if (bufferlines < 50) 
    WARNING.print("Small buffersize: increase MEMORY or decrease EXTRADENSE?");

  // ______ The big loop per buffer ______
  register int32 i,j;// DEM index grid counter
  for (register int32 buffer=0; buffer<numfullbuffers+EXTRABUFFER; ++buffer)
    {
    const real8 firstline = veryfirstline + real8(buffer*bufferlines*mlL);
    const int32 blines = (buffer == numfullbuffers) ? restlines : bufferlines;
    matrix<real4> OBUFFER(blines,numpixelsml);
    if (buffer == numfullbuffers) // smaller last one, EXTRABUFFER==1
      OBUFFER.resize(restlines,numpixelsml);
    const real8 lastline  = firstline + real8((OBUFFER.lines()-1)*mlL);

    // ______ Extra info ______
    PROGRESS << "Buffer# [l0:lN, p0:pN]: " << buffer+1 << " ["
         << firstline  << ": " << lastline  << ", "
         << firstpixel << ": " << lastpixel << "]";
    PROGRESS.print();

    // ______ find out corners of this buffer in phi/lambda ______
    lp2ell(firstline, firstpixel,
           ellips, master, masterorbit,
           phi, lambda, height);                // returned
    phil0p0    = phi;
    lambdal0p0 = lambda;
    lp2ell(lastline, firstpixel,
           ellips, master, masterorbit,
           phi, lambda, height);                // returned
    philNp0    = phi;
    lambdalNp0 = lambda;
    lp2ell(lastline,lastpixel,
           ellips, master, masterorbit,
           phi, lambda, height);                // returned
    philNpN    = phi;
    lambdalNpN = lambda;
    lp2ell(firstline,lastpixel,
           ellips, master, masterorbit,
           phi, lambda, height);                // returned
    phil0pN    = phi;
    lambdal0pN = lambda;

    // ______ Corners for DEM 2b read from disk ______
    // BK 14-Feb-00
    phimin    = min(min(min(phil0p0,philNp0),philNpN),phil0pN);
    phimax    = max(max(max(phil0p0,philNp0),philNpN),phil0pN);
    lambdamin = min(min(min(lambdal0p0,lambdalNp0),lambdalNpN),lambdal0pN);
    lambdamax = max(max(max(lambdal0p0,lambdalNp0),lambdalNpN),lambdal0pN);

    // ______ Read extra large part from DEM because due to layover etc. ______
    // BK 20-Apr-2002
    //phimin    -= deg2rad(1/200.); // 1 degree is 100 km (equator)
    // longitude far range can be problem if height, say 2 km extra.)
    // because of skewed geometry of IFG over DEM North/South.
    // 1 degree is 100 km, 2 km extra should do it, if mountain is < 2 km?.
    //phimin    -= 1.5*DEMdeltalat;
    //const real8 extrafarrange = (1.5*DEMdeltalon + deg2rad(1.0/100.0));
    //const real8 extralatitude = (1.5*DEMdeltalon + deg2rad(1.0/100.0));
    phimin    -= extralatitude;
    phimax    += extralatitude;
    lambdamax += extrafarrange;
    lambdamin -= extrafarrange;
    INFO << "Using larger DEM input for far range: "
         << rad2deg(extrafarrange) << " [deg]";
    INFO.print();


    // ====== Read DEM from file ======
    // ______ Get index of first line/pixel to load from file ______
    // ______ Index boundary: [0:numberofx-1] ______
    // ______ upperleft == max_lat,min_lon ______
    int32 indexphi0DEM = int32(floor((lat0file-phimax)/DEMdeltalat));
    if (indexphi0DEM < 0) 
      {
      indexphi0DEM=0;           // default start at first
      WARNING.print("DEM does not cover entire interferogram.");
      WARNING.print("input DEM should be extended to the North.");
      }
    int32 indexphiNDEM = int32(ceil((lat0file-phimin)/DEMdeltalat));
    if (indexphiNDEM > numberoflatpixels-1)
      {
      indexphiNDEM=numberoflatpixels-1;
      WARNING.print("DEM does not cover entire interferogram.");
      WARNING.print("input DEM should be extended to the South.");
      }
    int32 indexlambda0DEM = int32(floor((lambdamin-lon0file)/DEMdeltalon));
    if (indexlambda0DEM < 0)
      {
      indexlambda0DEM=0;                // default start at first
      WARNING.print("DEM does not cover entire interferogram.");
      WARNING.print("input DEM should be extended to the West.");
      }
    int32 indexlambdaNDEM = int32(ceil((lambdamax-lon0file)/DEMdeltalon));
    if (indexlambdaNDEM > numberoflonpixels-1)
      {
      indexlambdaNDEM=numberoflonpixels-1;
      WARNING.print("DEM does not cover entire interferogram.");
      WARNING.print("input DEM should be extended to the East.");
      }

    // ______ lat/lon for first pixel in matrix read from file ______
    // ______ upper is max. latitude, left is min. longitude ______
    const real8 upperleftphi    = lat0file-indexphi0DEM*DEMdeltalat;
    const real8 upperleftlambda = lon0file+indexlambda0DEM*DEMdeltalon;

    window zerooffset  (0,0,0,0);
    window winfromfile (indexphi0DEM,indexphiNDEM,
                        indexlambda0DEM,indexlambdaNDEM);

    // BK: COMMENTED OUT 4may2001, see below.
    //      // gtopo30 (lat,lon) ... format signed short
    //      // leftupper is first
    //      // ?? assume lat 0:90, what if 0:-90 ?
    //      // check how to with other formats?
    //      matrix<int16> DEM(indexphiNDEM-indexphi0DEM+1,indexlambdaNDEM-indexlambda0DEM+1);
    //      readfile(DEM,demfilename,numberoflatpixels,winfromfile,zerooffset);
    //      // phi (latitude) is stored rowwise, starting at leftupperphi
    //      // lambda (longitude) are the columns, start=leftupperlambda.

    // ______ Read in grdfile of DEM in matrix R4 (raw data, no header) _______
    // ______ added formats (BK 4-May-2001) ______
    PROGRESS << "Reading crop of DEM for buffer: " << buffer+1;
    PROGRESS.print();
    DEBUG.print("Reading input DEM into real4 matrix (buffer).");
    matrix<real4> DEM(indexphiNDEM-indexphi0DEM+1,indexlambdaNDEM-indexlambda0DEM+1);
    switch (refdeminput.iformatflag)
      {
      // ______ Read as short BE, then convert to host order ______
      case FORMATI2_BIGENDIAN:
        {
        matrix<int16> DEMi2(indexphiNDEM-indexphi0DEM+1,indexlambdaNDEM-indexlambda0DEM+1);
        readfile(DEMi2,demfilename,numberoflatpixels,winfromfile,zerooffset);
        for (int32 iii=0; iii<DEM.lines(); ++iii)
          for (int32 jjj=0; jjj<DEM.pixels(); ++jjj)
            DEM(iii,jjj) = real4(ntohs(DEMi2(iii,jjj)));// cast to real4
        DEMi2.resize(1,1);// dealloc...
        INFO.print("Read crop of input DEM: format: SHORT SIGNED INT.");
        break;
        }

      case FORMATI2:
        {
        matrix<int16> DEMi2(indexphiNDEM-indexphi0DEM+1,indexlambdaNDEM-indexlambda0DEM+1);
        readfile(DEMi2,demfilename,numberoflatpixels,winfromfile,zerooffset);
        for (int32 iii=0; iii<DEM.lines(); ++iii)
          for (int32 jjj=0; jjj<DEM.pixels(); ++jjj)
            DEM(iii,jjj) = DEMi2(iii,jjj);// cast to real4
        DEMi2.resize(1,1);// dealloc...
        INFO.print("Read crop of input DEM: format: SHORT SIGNED INT.");
        break;
        }

      case FORMATR4:
        readfile(DEM,demfilename,numberoflatpixels,winfromfile,zerooffset);
        INFO.print("Read crop of input DEM: format: REAL4.");
        break;
      case FORMATR8:
        {
        matrix<real8> DEMr8(indexphiNDEM-indexphi0DEM+1,indexlambdaNDEM-indexlambda0DEM+1);
        readfile(DEMr8,demfilename,numberoflatpixels,winfromfile,zerooffset);
        for (int32 iii=0; iii<DEM.lines(); ++iii)
          for (int32 jjj=0; jjj<DEM.pixels(); ++jjj)
            DEM(iii,jjj) = DEMr8(iii,jjj);// cast to real4
        DEMr8.resize(1,1);// dealloc...
        INFO.print("Read crop of input DEM: format: REAL8.");
        break;
        }
      default:
        PRINT_ERROR("totally impossible, checked input.")
        throw(unhandled_case_error);
      }



    // ----- Loop over DEM for stats ------------------------
    real8 min_dem_buffer =  100000.0;
    real8 max_dem_buffer = -100000.0;
    for (i=0; i<DEM.lines(); ++i)
      {
      // ----- Loop over oversampled matrix in x ------
      for (j=0; j<DEM.pixels(); ++j)
        {
        if(DEM(i,j)!=NODATA)
          {
          numvalid++;
          meancroppedDEM += DEM(i,j);// divide by numvalid later
          if (DEM(i,j)<min_dem_buffer) min_dem_buffer=DEM(i,j);//buffer
          if (DEM(i,j)>max_dem_buffer) max_dem_buffer=DEM(i,j);// stats
          }
        else
          {
          numNODATA++;
          }
        }//loop dem for stats
      }//loop dem for stats
    min_input_dem = min(min_input_dem,min_dem_buffer);//global stats
    max_input_dem = max(max_input_dem,max_dem_buffer);//global stats


    // ______ Interpolate DEM with dense enough grid spacing ______
    // ______ Cubic conv. first point upto (not including) last one ______
    PROGRESS << "Oversampling DEM by bilinear splines for buffer: "
         << buffer+1;
    PROGRESS.print();
    matrix<real4> DEMinterpolated(
      (DEM.lines()-1)*DEMphioversamplefactor,
      (DEM.pixels()-1)*DEMlambdaoversamplefactor);

    const real8 deltalatDEMinterpolated = 
      DEMdeltalat/real8(DEMphioversamplefactor);
    const real8 deltalonDEMinterpolated = 
      DEMdeltalon/real8(DEMlambdaoversamplefactor);
    // ------ Check if they are equal, due to small extradense, -----
    // ------ high input resolution, or large multilooking ------
    if (DEMphioversamplefactor==1 && DEMlambdaoversamplefactor==1)
      {
      INFO.print("Interpolated DEM same as input DEM, making a copy");
      INFO.print("-> extradense small, e.g., method trilinear");
      INFO.print("-> or maybe input dem has already high resolution");
      INFO.print("-> or maybe large multilooking factors.");
      DEMinterpolated = DEM;// could be done with swap?
      }
    // ______ Do bilinear interpolation ______
    else
      {
      // ______ Do bilinear interpolation ______
      // f(row,col) = a00 + a01*c + a10*r +a11*r*c 
      // P1(0,0); P2(0,1); P3(1,0); P4(1,1); P[0:1,0:1];
      //
      //         cols->  
      //      r P1    P2
      //      o    .P
      //      w P3    P4
      //      s
      // ----- Loop over oversampled matrix in y ------
      for (i=0; i<DEM.lines()-1; ++i)   // bi-c: use i,i+1 -> 4 points
          {
        if ((i%100)==0)
          {
          PROGRESS << "Interpolation of DEM: line: " << i << " (" 
               << floor(.5+(100.*real8(i)/real8(DEM.lines())))
               << "%)";
          PROGRESS.print();
          }
        // ----- Loop over oversampled matrix in x ------
        for (j=0; j<DEM.pixels()-1; ++j)
          {
          // ______ 4 points required for bilinear ______
          // ______ other interpolators may need more points ______
          // ______ We assume the user can externally oversample ______
          // ______ the input DEM to high enough resolution  ______
          const real4 P1 = DEM(i,j);
          const real4 P2 = DEM(i,j+1);
          const real4 P3 = DEM(i+1,j);
          const real4 P4 = DEM(i+1,j+1);
          // ______ check if NODATA then set h=zero -> flatearth correction ______
          if (P1==NODATA || P2==NODATA || P3==NODATA || P4==NODATA)
            {
            // for (k=i*factor; k<(i+1)*factor; ++k)
            //   for (l=j*factor; l<(j-1)factor; ++l)
            //     DEMinterpolated(k,l)=0.;
            // DEMinterpolated is already 0 initialized.
            //break;
            continue;   // with next j
            }
          // ______ Else bilinear interpolation to finer grid ______
          // ______ ignore last line/row for speed. ______
          const real4 a00 = P1;
          const real4 a01 = P2-P1;
          const real4 a10 = P3-P1;
          const real4 a11 = P4-P3-P2+P1;
          const int indexrowinterpolated = i*DEMphioversamplefactor;
          const int indexcolinterpolated = j*DEMlambdaoversamplefactor;
          register real4 row = 0.0;
          for (int32 k=indexrowinterpolated;
                     k<indexrowinterpolated+DEMphioversamplefactor;
                     ++k)
            {
            register real4 col = 0.0;
            for (int32 l=indexcolinterpolated;
                       l<indexcolinterpolated+DEMlambdaoversamplefactor;
                       ++l)
              {
              DEMinterpolated(k,l) = a00 + a01*col + a10*row + a11*col*row;
              col += 1.0/DEMlambdaoversamplefactor;
              }
            row += 1.0/DEMphioversamplefactor;
            }
          } // loop DEM pixels
        } // loop DEM lines
      }// End check interpolation required


    // ====== Output debugging matrices if requested ======
    if (outputdem)
      {
      // create a new file per buffer for cropped DEM that is used.
      ofstream demofile;                                // input DEM
      char            newfilename[ONE27];
      ostrstream      qomem(newfilename,ONE27);
      qomem << refdeminput.fodem << "_buffer_" << buffer << ends;
      char filename[ONE27];
      strcpy(filename,newfilename);
      openfstream(demofile,filename,generalinput.overwrit);
      bk_assert(demofile,filename,__FILE__,__LINE__);
      // ___ Give info to debug user ___
      DEBUG << "DEM output to file: " << filename
           << "; format: real4;  size (l,p): "
           << DEM.lines() << ", " << DEM.pixels();
      DEBUG.print();
      demofile << DEM;
      demofile.close();
      }
    if (outputdemi)
      {
      ofstream demiofile;                               // interpolated input DEM
      char            newfilename[ONE27];
      ostrstream      qomem(newfilename,ONE27);
      qomem << refdeminput.fodemi << "_buffer_" << buffer << ends;
      char filename[ONE27];
      strcpy(filename,newfilename);
      openfstream(demiofile,filename,generalinput.overwrit);
      bk_assert(demiofile,filename,__FILE__,__LINE__);
      // ___ Give info to debug user ___
      DEBUG << "Interpolated DEM file: " << filename
           << "; format: real4; size (l,p): "
           << DEMinterpolated.lines() << " " << DEMinterpolated.pixels();
      DEBUG.print();
      demiofile << DEMinterpolated;
      demiofile.close();
      }
    // --- remove DEM, from now use DEMinterpolated ---
    DEBUG.print("resizing DEM to 1,1 to reduce memory requirements");
    DEM.resize(1,1);// "dealloc"

    // ====== DEMinterpolated contains now values from leftupper with ======
    // ______ spacing (deltalatDEMinterpolated,deltalonDEMinterpolated) ______
    // ______ Transform DEM to l,p,refphase; store in LPPHASE ______
    PROGRESS.print("Converting oversampled DEM to radar system for this buffer.");
    const int32 numpointsDEMinterpolated = DEMinterpolated.size();
    const int32 numpointsBUFFER = OBUFFER.size();
    // ______ Extra info ______
    INFO << "Number of points in BUFFER (IFG):    "
         << numpointsBUFFER;
    INFO.print();
    INFO << "Number of points in oversampled DEM: "
         << numpointsDEMinterpolated;
    INFO.print();
    if (numpointsBUFFER > numpointsDEMinterpolated)
      if (do_nearest==true)
        WARNING.print("N points IFG OBUFFER > N_points interpolated DEM (increase EXTRADENSE?)");
      else
        INFO.print("N points IFG OBUFFER > N_points interp. dem (is ok)");

    matrix<real8> LPPHASE(numpointsDEMinterpolated,3);// real8 is important
    if (outputrefdemhei==true)
      {
      DEBUG.print("Adding column for height in radar coordinates.");
      LPPHASE.resize(numpointsDEMinterpolated,4);// extra column for height
      }

    // --- Loop bilinearly interpolated DEM ---
    cn P;
    int32 LPPHASEindex=0;
    phi = upperleftphi;
    for (i=0; i<DEMinterpolated.lines(); ++i)
      {
      if ((i%100)==0)
        {
        // ______ Extra info ______
        PROGRESS << "Radarcoding interpolated DEM line: " << i << " ("
             << floor(.5+(100.*real8(i)/real8(DEMinterpolated.lines())))
             << "%)";
        PROGRESS.print();
        }

      lambda = upperleftlambda;
      for (j=0; j<DEMinterpolated.pixels(); ++j)
        {
        height = DEMinterpolated(i,j);
        ell2xyz(ellips,P,phi,lambda,height);            // returns P(x,y,z)
        real8 t_range_master;                           // range (fast) time
        real8 t_azi_master;                             // azimuth (slow) time
        xyz2t(t_azi_master,t_range_master,              // t returned for (l,p)
              master,masterorbit,
              P,MAXITER,CRITERTIM);
        real8 t_azi_dummy;
        real8 t_range_slave;
        xyz2t(t_azi_dummy,t_range_slave,                // t returned
              slave,slaveorbit,
              P,MAXITER,CRITERTIM);

        // ______ (line,pixel) in master system ______
        LPPHASE(LPPHASEindex,0) = master.ta2line(t_azi_master);
        LPPHASE(LPPHASEindex,1) = master.tr2pix(t_range_master);
        // only store t_range_slave here, compute ref_phase later!
        LPPHASE(LPPHASEindex,2) = t_range_slave;
        if (outputrefdemhei==true)
          LPPHASE(LPPHASEindex,3) = height;// for Zbigniew Perski

        // ______ update longitude of next pixel ______
        lambda += deltalonDEMinterpolated;
        ++LPPHASEindex;
        } // loop DEMinterpolated pixels

      // ______ update latitude of next line ______
      phi -= deltalatDEMinterpolated;           // upper left is max. value
      } // loop DEMinterpolated lines


    // ====== Nearest neighbor to get reference phase to OBUFFER ======
    PROGRESS << "Start nearest neighbor for buffer: " << buffer+1;
    PROGRESS.print();
    // tested PNT2 as a real8* matrix, no gain in speed... BK 15-feb-00
    // matrix<real8*> PNT2LPPHI(numlinesOBUFFER,numpixelsOBUFFER);// pointers
    // Store squared distances in DISTANCES, OBUFFER could be used if cleaned.
    const int32 numlinesOBUFFER  = OBUFFER.lines();
    const int32 numpixelsOBUFFER = OBUFFER.pixels();
    matrix<int32> PNT2LPPHI(numlinesOBUFFER,numpixelsOBUFFER);// 'pointers' ...
    matrix<real4> DISTANCES(numlinesOBUFFER,numpixelsOBUFFER);  // not really required..
    const real4 DEFAULTDIST      = 99999.;// must be a positive large number
    DISTANCES.setdata(DEFAULTDIST);// store here the pixel distance for nearest neighbor.
    for (i=0; i<numpointsDEMinterpolated; ++i)
      {
      if ((i%100000)==0)
        {
        // ______ Extra info ______
        PROGRESS << "Nearest neighbor point number: " << i << " ("
             << floor(.5+(100.*real8(i)/real8(numpointsDEMinterpolated)))
             << "%)";
        PROGRESS.print();
        }
      // ______ Index in output buffer ______
      // ______ firstline/pixel is center of multilooked pixel in original radar coord. ___
      // ______ distance of radarcoded point to the center is important ______
      // --- compute distance with center location of multilooked pixel
      // --- firstline is center; want correct bin in multilooked output buffer
      //real8(mlL)/2.0 was added to get center of firstline, so subtract it again
      const int32 Y_in_obuffer = 
        int32((LPPHASE(i,0)-(firstline-real8(mlL)/2.0))/real8(mlL));// bin
      const int32 X_in_obuffer = 
        int32((LPPHASE(i,1)-(firstpixel-real8(mlP)/2.0))/real8(mlP));// bin
      //
      //const real8 indexobufferLr8 = (real8(LPPHASE(i,0))-firstline) /real8(mlL);
      //const real8 indexobufferPr8 = (real8(LPPHASE(i,1))-firstpixel)/real8(mlP);
      //const int32 indexobufferL   =  int32(indexobufferLr8+0.5); // round
      //const int32 indexobufferP   =  int32(indexobufferPr8+0.5); // round
      // ______ Check if index is within boundaries ______
      // ______ It is possible that the index is negative since ______
      // ______ the whole DEM is radarcoded, also outside the IFG ______
      if (Y_in_obuffer>=0 && Y_in_obuffer<numlinesOBUFFER &&
          X_in_obuffer>=0 && X_in_obuffer<numpixelsOBUFFER  )
        {
        //const real4 squareddistance = sqr(indexobufferLr8-indexobufferL) +
        //                            sqr(indexobufferPr8-indexobufferP);
        //const real4 dist = sqrt(sqr(indexobufferLr8-real8(indexobufferL)) +
        //                      sqr(indexobufferPr8-real8(indexobufferP)));
        const real8 line_desired  = firstline  + real8(Y_in_obuffer*mlL);// 
        const real8 pixel_desired = firstpixel + real8(X_in_obuffer*mlP);
        const real4 dist = sqrt(sqr(real8(LPPHASE(i,0))-line_desired) +
                                sqr(real8(LPPHASE(i,1))-pixel_desired));
        if (dist < DISTANCES(Y_in_obuffer,X_in_obuffer)) 
          {
          DISTANCES(Y_in_obuffer,X_in_obuffer) = dist;
          PNT2LPPHI(Y_in_obuffer,X_in_obuffer) = i;// 'pointer' to range time slave
          }
        }
      } // all points in LPPHASE

 
// ______ Check if all points are assigned a value (?) ______
// ______ And compute reference phase ______
    PROGRESS << "Start reference phase computation for buffer: "
         << buffer+1;
    PROGRESS.print();
    // ======================================================================
    // ====== OLD METHOD; USES NEAREST NEIGHBOR DIRECTLY ====================
    // Bert Kampes, 08-Apr-2005
    // ======================================================================
    matrix<real4> OBUFFERHEI;// buffer for Zbigniew Perski
    if (outputrefdemhei==true) OBUFFERHEI.resize(1,numpixelsml);
    if (do_nearest==true)
      {
      WARNING.print("Using old method; please consider using TRI_LINEAR.");
      for (i=0; i<numlinesOBUFFER; ++i)
        {
        if ((i%25)==0)
          {
          // ______ Extra info ______
          PROGRESS << "Computing reference phase for nearest neighbor line: " << i 
               << " (" << floor(.5+(100.*real8(i)/real8(numlinesOBUFFER))) 
               << "%)";
          PROGRESS.print();
          }
        // --- Loop in range direction ---
        for (j=0; j<numpixelsOBUFFER; ++j)
          {
          if (DISTANCES(i,j) == DEFAULTDIST)
            {
            numNOneighbor++;
            WARNING << "No nearest neighbor for pixel at: "
                 << buffer*bufferlines+i << ", " << j
                 << "(i.e., " << firstline + real8(i*mlL) 
                 << ", " << firstpixel + real8(j*mlP) << ").";
            WARNING.print();
            continue; // with next j, leave 0 in OBUFFER.
            }
          // ______ Compute ref phase, nearest neighbor is known ______
          const int32 lpindex       = PNT2LPPHI(i,j);// index ('pointer') to LPP matrix
          const real8 line          = LPPHASE(lpindex,0);
          const real8 pixel         = LPPHASE(lpindex,1);
          const real8 t_range_slave = LPPHASE(lpindex,2);
          // ______ Compute master time for this pixel ______
          // ______ If appropriate, ref.phase = phi_topo ______
          real8 ref_phase;
          if (onlyrefphasetopo)                 // do not include flat earth phase
            {
            lp2xyz(line,pixel,
                   ellips,                      // h==0
                   master, masterorbit,
                   P,MAXITER,CRITERPOS);                        // P returned
            real8 t_range_flatearth;
            real8 t_azi_dummy;
            xyz2t(t_azi_dummy,t_range_flatearth,
                  slave,slaveorbit,
                  P,MAXITER,CRITERTIM);                 // P on h=0
            ref_phase = m_min4picdivlam*t_range_flatearth-
                        s_min4picdivlam*t_range_slave;
            }
          else // include flatearth, ref.pha = phi_topo+phi_flatearth
            {
            ref_phase = 
              m_min4picdivlam*master.pix2tr(pixel)-
              s_min4picdivlam*t_range_slave;
            }
          OBUFFER(i,j) = real4(ref_phase);

          // --- Dump height in radarcoordinates directly if requested
          if (outputrefdemhei==true)// For Zbigniew Perski
            {
            OBUFFERHEI(0,j) = real4(LPPHASE(lpindex,3));// output line buffer
            }
          } // j pixel loop OBUFFER
        // For Zbigniew Perski
        if (outputrefdemhei==true)// For Zbigniew Perski
          {
          DEBUG.print("Writing buffer line for Zbigniew Perski");
          refdemheiofile << OBUFFERHEI;
          }
        } // i line loop OBUFFER
      }
    // ======================================================================
    // ====== NEW METHOD, USES TRI-LINEAR INTERPOLATION OF NEAREST POINTS ===
    // Bert Kampes, 08-Apr-2005
    // ====== Now find 3 closest points A,B,C and do trigrid lin. interpolation
    // ====== point of interpolation is "+", local coordinate system has origin there.
    // ====== Then, cn A.z = f(x,y) = a00+a10*A.x+a01*A.y;
    // ======       cn B.z = f(x,y) = a00+a10*B.x+a01*B.y;
    // ======       cn C.z = f(x,y) = a00+a10*C.x+a01*C.y;
    // ====== and solve for a00, which is the Z value at origin, i.e., it!
    //
    //     A
    //       +   C
    //  B
    //
    // ====== To find A,B,C close to +, I first do neareast neighbor
    // ====== and then a local search from the output grid in each direction.
    // ====== Note that all interpolations are linear, i.e., the input
    // ====== DEM should contain all peaks and should be dense.
    // ====== also note that this throws out some of the points that
    // ====== were previously oversampled.
    // ====== In future, using this framework with DD,
    // ------- it is very simple to create a better interpolator...
    // Bert Kampes, 07-Apr-2005
    // ======================================================================
    else
      {
      DEBUG.print("Using new comprefdem interpolator method.");
      // --- Set up matrices that specify order for search of closeby points ---
      const int search_radius = 20;// in x,y smaller is faster, but risky?
      const real4 radius      = real4(search_radius);// for Type templates
      const real4 NX          = 2*search_radius+1;// in x,y
      const real4 NY          = NX;// square matrix; used for Type ones()
      const real4 one         = real4(1);// used for Type ones()
      // --- Compare meshgrid in matlab:// define matrix from -10:10
      const matrix<real4> XX  = ones(NX,one)*(indgen(NX)-real4(radius));// [-10:10]
      const matrix<real4> YY  = matTxmat((indgen(NY)-radius),ones(one,NY));
      // 1 pixel increase in line direction is what in multilooked?
      // XX*ml_X.^2 *ground_resolution??  or use Euclidian distance?
      DEBUG.print("Assuming that multilooked interferogram has more or less square pixels");
      DEBUG.print("i.e., going in X or Y with 1 pixel is same distance (?)");
      const matrix<real4> DD  = dotmult(XX,XX)+dotmult(YY,YY);// sqr.distances
      // --- now sort the distances and go through them until border of input matrix;
      matrix<real4> SORTED_DD(DD.size(),3);// distance,Y,X
      int32 cnt=0;
      for (int32 y=0; y<XX.lines();++y)
        {
        for (int32 x=0; x<XX.pixels();++x)
          {
          SORTED_DD(cnt,0) = DD(y,x);// distance
          SORTED_DD(cnt,1) = YY(y,x);// delta Y to go to
          SORTED_DD(cnt,2) = XX(y,x);// delta X to go to
          ++cnt;
          }
        }
      // --- qsort the distances in ascending order ---
      mysort2(SORTED_DD);// sort ascending on distance to center
      // --- Loop over output grid in azimuth direction --------------
      real8 height_A=0.0, height_B, height_C;// For Zbigniew Perski
      for (i=0; i<numlinesOBUFFER; ++i)
        {
        // non-multilooked center (close to LPP)
        const real8 line_desired = firstline + real8(i*mlL);
        if ((i%100)==0)
          {
          // ______ Extra info ______
          PROGRESS << "Computing reference phase trilinear for line: " << i 
               << " (" << floor(.5+(100.*real8(i)/real8(numlinesOBUFFER))) 
               << "%)";
          PROGRESS.print();
          }
        // --- Loop over output grid in range direction ---------------
        for (j=0; j<numpixelsOBUFFER; ++j)
          {
          // non-multilooked center (close to LPP)
          const real8 pixel_desired = firstpixel + real8(j*mlP);
          // --- Find 3 closest points that have a nearest neighbor assigned ---
          // --- Could easily find 10 if you'd like ---
          int32 n_found = 0;// number of points found closeby
          cn A,B,C;// {x,y,phase} triplet of nearest three points
          //cerr << endl;
          for (cnt=0; cnt<XX.size(); ++cnt)
            {
            //const real4 this_D = SORTED_DD(cnt,0);// distance to next closest point in matrix
            const int32 this_Y = i + int32(SORTED_DD(cnt,1));// next closest point
            const int32 this_X = j + int32(SORTED_DD(cnt,2));// next closest point
            if (this_Y>=0 && this_Y<numlinesOBUFFER &&
                this_X>=0 && this_X<numpixelsOBUFFER)
              {
              // --- Check if this is a point that has a nearest neighbor: ---
              if (DISTANCES(this_Y,this_X) != DEFAULTDIST)
                {
                // ______ Compute ref phase for this point ______
                const int32 lpindex       = PNT2LPPHI(this_Y,this_X);// index ('pointer') to LPP matrix
                const real8 line          = LPPHASE(lpindex,0);
                const real8 pixel         = LPPHASE(lpindex,1);
                const real8 t_range_slave = LPPHASE(lpindex,2);
                // ______ Compute master time for this pixel ______
                // ______ If appropriate, ref.phase = phi_topo ______
                real8 ref_phase;
                if (onlyrefphasetopo)// do not include flat earth phase
                  {
                  lp2xyz(line,pixel,
                         ellips,                        // h==0
                         master, masterorbit,
                         P,MAXITER,CRITERPOS);                  // P returned
                  real8 t_range_flatearth;
                  real8 t_azi_dummy;
                  xyz2t(t_azi_dummy,t_range_flatearth,
                        slave,slaveorbit,
                        P,MAXITER,CRITERTIM);                   // P on h=0
                  ref_phase = m_min4picdivlam*t_range_flatearth-
                              s_min4picdivlam*t_range_slave;
                  }
                else // include flatearth, ref.pha = phi_topo+phi_flatearth
                  {
                  ref_phase = 
                    m_min4picdivlam*master.pix2tr(pixel)-
                    s_min4picdivlam*t_range_slave;
                  }
                switch (n_found)
                  {
                  case 0:
                    A.x = pixel - pixel_desired;// center desired position 
                    A.y = line  - line_desired;// ordinate
                    A.z = ref_phase;// functional value
                    if (outputrefdemhei==true)// For Zbigniew Perski
                      height_A = LPPHASE(lpindex,3);
                    break;
                  case 1:
                    B.x = pixel - pixel_desired;// center desired position 
                    B.y = line  - line_desired;// ordinate
                    B.z = ref_phase;// functional value
                    if (outputrefdemhei==true)// For Zbigniew Perski
                      height_B = LPPHASE(lpindex,3);
                    break;
                  case 2:
                    C.x = pixel - pixel_desired;// center desired position 
                    C.y = line  - line_desired;// ordinate
                    C.z = ref_phase;// functional value
                    if (outputrefdemhei==true)// For Zbigniew Perski
                      height_C = LPPHASE(lpindex,3);
                    break;
                  default:
                    ;// impossible
                  }// switch
                n_found++;
                }// point found
              }//index within matrix
            if (n_found==3) break;
            }
          // --- Interpolate phase to pixel position ---
          // +++ here add tri-interpolation of closest three points
          // +++ Bert Kampes, 07-Apr-2005
          real8 interpolated_refphase = 0.0;
          if (n_found==3)
            {
            interpolated_refphase = interp_trilinear(A,B,C);// optimized inline
            }
          else if (n_found>=1)
            {
            interpolated_refphase = A.z;// first point found is closest
            }
          // === Store in output buffer ===
          OBUFFER(i,j) = real4(interpolated_refphase);
          // --- Dump height in radarcoordinates directly if requested
          // Bert Kampes, 16-Apr-2005
          if (outputrefdemhei==true)// For Zbigniew Perski
            {
            real4 closest_refphase = A.z;// tmp copy
            real8 interpolated_refheight;
            A.z = height_A;// height of first point
            B.z = height_B;// height of second point
            C.z = height_C;// height of third point
            if (n_found==3)
              {
              // could be optimized to compute a00 for hei/refpha at same time.
              interpolated_refheight = interp_trilinear(A,B,C);
              // --- do not trust result, e.g. points on line  ---
              if (abs(interpolated_refheight-A.z)>500.0 || //[m]
                  interpolated_refheight > 10000.0 ||      //[m]
                  interpolated_refheight < -1000.0)//[m]
                {
                DEBUG.print("Following interpolated value is suspect (using A.z):");
                A.showdata();
                B.showdata();
                C.showdata();
                DEBUG << "interpolated P.z = " << interpolated_refheight;
                DEBUG.print();
                interpolated_refheight = A.z;//
                OBUFFER(i,j) = closest_refphase;
                }
              }
            else
              {
              interpolated_refheight = A.z;// use nearest point
              }
            // --- write the height to file (hopefully real4?) ---
            OBUFFERHEI(0,j) = real4(interpolated_refheight);//
            }
          } // j pixel loop OBUFFER
        // For Zbigniew Perski
        if (outputrefdemhei==true)// For Zbigniew Perski
          {
          DEBUG.print("Writing buffer line for Zbigniew Perski");
          refdemheiofile << OBUFFERHEI;
          }
        } // i line loop OBUFFER
      }// method selector nn or interp


    // ___ Give statistics for this buffer ___ 
    if (numvalid<50) WARNING.print("Less than 50 valid pixels in input DEM.");
    meancroppedDEM    /= numvalid;
    real8 numpoints    = real8(numvalid+numNODATA);// in this buffer DEM
    real8 numpointsOUT = real8(OBUFFER.lines()*OBUFFER.pixels());
    INFO << "Number of points with NODATA value: " << numNODATA 
         << " (" << 100.0*numNODATA/numpoints << "%)";
    INFO.print();
    if (numNODATA/numpoints>0.5)
      WARNING.print(INFO.get_str());
    INFO << "Number of points w/o nearest neighbor: " << numNOneighbor
         << " (" << 100.0*numNOneighbor/numpointsOUT << "%)";
    INFO.print();
    if (numNOneighbor/numpointsOUT>0.5)
      WARNING.print(INFO.get_str());
    INFO << "Min. value of cropped input DEM: " << min_dem_buffer; 
    INFO.print();
    INFO << "Max. value of cropped input DEM: " << max_dem_buffer; 
    INFO.print();
    INFO << "Mean value of cropped input DEM: " << meancroppedDEM 
         << " (over " << numvalid << " valid pixels.)";
    INFO.print();
    if (meancroppedDEM>5000)
      WARNING.print("I suspect byte order problem, mean of DEM>5000");
    if (meancroppedDEM<-100)
      WARNING.print("I suspect byte order problem, mean of DEM<-100");
    if (abs(meancroppedDEM)<1)
      {
      WARNING.print("I suspect input format error, abs(mean of DEM)<1");
      WARNING.print("Or are you subtracting ellipsoid reference phase like this?");
      }
    if (numvalid==0)
      WARNING.print("NO valid pixels found in crop of input DEM.");
    // ______ Nearest neighbor distances stats ______
    int32 num_stat  = 0;
    real8 stat_dist = 0;
    for (i=0; i<numlinesOBUFFER ;++i)
      {
      for (j=0; j<numpixelsOBUFFER ;++j)
        {
        if (DISTANCES(i,j) != DEFAULTDIST)
          num_stat++;
          stat_dist += DISTANCES(i,j);// squared distance was stored.
        }
      }
    INFO << "Mean distance for " << num_stat << " nearest neighbors: "
         << stat_dist/num_stat << " [pixels]";
    INFO.print();

    // ______ Reset counters for next buffer ______
    totalNOneighbor += numNOneighbor;
    numvalid         = 0;
    numNODATA        = 0;
    numNOneighbor    = 0;
    meancroppedDEM   = 0.;

    // ====== Write OBUFFER to ofile ======
    PROGRESS << "Writing reference phase for buffer: " << buffer+1;
    PROGRESS.print();
    refdemofile << OBUFFER;
    } // buffer loop
  refdemofile.close();  // OBUFFER
  if (outputrefdemhei==true)// For Zbigniew Perski
    {
    INFO.print("Closing file for additional output");
    refdemheiofile.close();// OBUFFERHEI
    }



// ====== Write output information ======
  char croppeddem[ONE27];
  char croppeddemi[ONE27];
  strcpy(croppeddem, "NO output requested");
  strcpy(croppeddemi,"NO output requested");
  if (outputdem)  strcpy(croppeddem,refdeminput.fodem);
  if (outputdemi) strcpy(croppeddemi,refdeminput.fodemi);
  INFO << "Min. value of input DEM covering interferogram: " << min_input_dem; 
  INFO.print();
  INFO << "Max. value of input DEM covering interferogram: " << max_input_dem; 
  INFO.print();

  ofstream scratchlogfile("scratchlogcomprefdem", ios::out | ios::trunc);
  bk_assert(scratchlogfile,"comprefdem: scratchlogcomprefdem",__FILE__,__LINE__);
  scratchlogfile
    << "\n*_Start_" << processcontrol[pr_i_comprefdem]
    << "\n*******************************************************************"
    << "\n1) DEM source file:                   \t" <<  demfilename
    << "\nFormat:                               \t";
    switch (refdeminput.iformatflag)
      {
      case FORMATI2:
        {
        scratchlogfile << "SHORT SIGNED INTEGER (HOST ENDIANNESS)";
        break;
        }
      case FORMATI2_BIGENDIAN:
        {
        scratchlogfile << "SHORT SIGNED INTEGER, BIG ENDIAN";
        break;
        }
      case FORMATR4:
        {
        scratchlogfile << "REAL4 SIGNED FLOAT";
        break;
        }
      case FORMATR8:
        {
        scratchlogfile << "REAL8 SIGNED DOUBLE";
        break;
        }
      default:
        {
        scratchlogfile << "UNKNOWN? IMPOSSIBLE...";
        break;
        }
      }
  scratchlogfile
    << "\nByte order:                           \t" <<  "check yourself..."
    << "\nNumber of lines:                      \t" <<  numberoflatpixels
    << "\nNumber of pixels:                     \t" <<  numberoflonpixels
    << "\nResolution latitude:                  \t" 
    << rad2deg(DEMdeltalat) << " [deg]"
    << "\nResolution longitude:                 \t" 
    << rad2deg(DEMdeltalon) << " [deg]"
    << "\nMost West point in input DEM:         \t" << rad2deg(lon0file)
    << "\nMost East point in input DEM:         \t" << rad2deg(lonNfile)
    << "\nMost South point in input DEM:        \t" << rad2deg(latNfile)
    << "\nMost North point in input DEM:        \t" << rad2deg(lat0file)
    << "\nMin. value of input DEM covering interferogram: " << min_input_dem
    << "\nMax. value of input DEM covering interferogram: " << max_input_dem
    << "\n2) Output file cropped DEM:           \t" <<  croppeddem
    << "\nFormat:                               \t" << "REAL4"
    << "\nByte order:                           \t" << "(same as host)"
//    << "\nMean value:                           \t" <<  meancroppedDEM
    << "\n3) Output file interpolated crop DEM: \t" <<  croppeddemi
    << "\nFormat:                               \t" << "REAL4"
    << "\nByte order:                           \t" << "(same as host)"
    << "\nOversampling factor latitude:         \t" << DEMphioversamplefactor
    << "\nOversampling factor longitude:        \t" << DEMlambdaoversamplefactor
    << "\nApproximate resolution latitude:      \t" 
    << rad2deg(DEMdeltalat/DEMphioversamplefactor) << " [deg]"
    << "\nApproximate resolution longitude:     \t" 
    << rad2deg(DEMdeltalon/DEMlambdaoversamplefactor) << " [deg]"
    << "\n4) Output file synthetic phase:       \t" <<  refdemofilename
    << "\nFormat:                               \t" << "REAL4"
    << "\nByte order:                           \t" << "(same as host)"
    << "\nApproximate resolution latitude:      \t" 
    << rad2deg(dphi_interferogram) << " [deg]"
    << "\nApproximate resolution longitude:     \t" 
    << rad2deg(dlambda_interferogram) << " [deg]"
    << "\nNumber of lines (multilooked):        \t" <<  numlinesml
    << "\nNumber of pixels (multilooked):       \t" <<  numpixelsml
    << "\nNumber of points w/o nearest neighbor:\t" << totalNOneighbor
    << " (" << 100*totalNOneighbor/(numlinesml*numpixelsml) << "%)"

// this is not correct, only stats per buffer...
//    << "\n\n----- Other STATS -----"
//    << "\nTotal points in cropped DEM:          \t" << numpoints
//    << "\nNumber of valid points in DEM:        \t" << numvalid
//    << " (" << 100*numvalid/numpoints << "%)"
//    << "\nNumber of NODATA points in DEM:       \t" << numNODATA
//    << " (" << 100*numNODATA/numpoints << "%)"
//    << "\nMean height in meters at valid points:\t" << meancroppedDEM
    << "\n*******************************************************************\n\n";
  scratchlogfile.close();


  ofstream scratchresfile("scratchrescomprefdem", ios::out | ios::trunc);
  bk_assert(scratchresfile,"comprefdem: scratchrescomprefdem",__FILE__,__LINE__);
  scratchresfile
    << "\n\n*******************************************************************"
    << "\n*_Start_" << processcontrol[pr_i_comprefdem]
    << "\n*******************************************************************";
  if (do_nearest==true) scratchresfile
    << "\nMethod:                               \tNEAREST";
  else scratchresfile
    << "\nMethod:                               \tTRILINEAR";
  if (onlyrefphasetopo==true) scratchresfile
    << "\nInclude 'flat earth':                 \tNo";
  else scratchresfile
    << "\nInclude 'flat earth':                 \tYes";
  scratchresfile
    << "\nDEM source file:                      \t" <<  demfilename
    << "\nMin. of input DEM:                    \t" << min_input_dem
    << "\nMax. of input DEM:                    \t" << max_input_dem
    << "\nData_output_file:                     \t" <<  refdemofilename
    << "\nData_output_format:                   \t" << "real4"
    << "\nFirst_line (w.r.t. original_master):  \t"
    <<  interferogram.win.linelo
    << "\nLast_line (w.r.t. original_master):   \t"
    <<  interferogram.win.linehi
    << "\nFirst_pixel (w.r.t. original_master): \t"
    <<  interferogram.win.pixlo
    << "\nLast_pixel (w.r.t. original_master):  \t"
    <<  interferogram.win.pixhi
    << "\nMultilookfactor_azimuth_direction:    \t" <<  mlL
    << "\nMultilookfactor_range_direction:      \t" <<  mlP
    << "\nNumber of lines (multilooked):        \t" <<  numlinesml
    << "\nNumber of pixels (multilooked):       \t" <<  numpixelsml
    << "\n*******************************************************************"
    << "\n* End_" << processcontrol[pr_i_comprefdem] << "_NORMAL"
    << "\n*******************************************************************\n";
  scratchresfile.close();

  } // END radarcodedem

---