bk_baseline.hh

Go to the documentation of this file.

/*
 * Copyright (c) 1999-2005 Bert Kampes
 * Copyright (c) 1999-2005 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/bk_baseline.hh,v $ *
 * $Revision: 1.3 $ *
 * $Date: 2005/04/11 13:48:45 $ *
 * $Author: kampes $ *
 *
 * definition of BASELINE class *
 * 
 * During computations these should be used with care, they
 * may be a bit approximate.  Better is to do what Doris does,
 * i.e., evaluate positions and relate that to timing.
 * But i guess it is OK.  Use this for ultra fast geocoding...
 *
 #%// Bert Kampes, 06-Apr-2005
 ****************************************************************/

#ifndef BK_BASELINE_H
#define BK_BASELINE_H
using namespace std;

#include <iostream>    // cout
#include <cstdlib>     // exit()

#include "constants.hh"   // types etc.
#include "bk_messages.hh" // info etc.
#include "slcimage.hh"    // my slc image class
#include "productinfo.hh" // my products class
#include "orbitbk.hh"     // my orbit class
#include "conversion.hh"  // rad2deg()
#include "utilities.hh"   // eye()


// BASELINE is a new type (class) that is initialized using orbits
// and then either models the baseline parameters such as Bperp
// or can give them exact.
// Usage in the programs is something like in main do:
// BASELINE baseline;
// baseline.init(orbit1,orbit2,product?master?);
// and pass baseline to subprograms, there use baseline.get_bperp(x,y) etc.
// Bert Kampes, 31-Mar-2005


// ---------------------------------------------------------------- 
// ---------------------------------------------------------------- 
// ---------------------------------------------------------------- 
class BASELINE
  {
  // ______ Private data of BASELINE class ______
  private:
    bool initialized;//          
    real8 master_wavelength;//   tmp for now used for h_amb
    real8 nearrange;// range=nearrange+drange_dp*pixel
    real8 drange_dp;// range=nearrange+drange_dp*pixel
    real8 orbit_convergence;//   tmp for now constant
    real8 orbit_heading;//       tmp for now NOT USED
    real8 H_min;//            height at which parameters are modeled
    real8 H_max;//    to model phase=f(line,pix,hei) and hei=g(line,pix,phase)
    // --- B(l,p,h) = a000 + 
    //                a100*l   + a010*p   + a001*h   +
    //                a110*l*p + a101*l*h + a011*p*h +
    //                a200*l^2 + a020*p^2 + a002*h^2
    int32 N_coeffs;//  ==10 degree of 3D-poly to model B(l,p,h)
    // --- Coefficients ---
    matrix<real8> BPERP_cf;// perpendicular baseline
    matrix<real8> BPAR_cf;// parallel baseline
    matrix<real8> THETA_cf;// viewing angle
    matrix<real8> THETA_INC_cf;// incidence angle to satellite
    //real8 avg_height_ambiguity;//center height ambiguity

    // copy constructor; copy matrices etc as public.
    BASELINE(const BASELINE& A)
      {};// prevent copy constructor usage by privatizing
    //orbit *master;
    //orbit *slave;

  /****************************************************************
    // --- B(l,p,h) = a000 + 
    //                a100*l   + a010*p   + a001*h   +
    //                a110*l*p + a101*l*h + a011*p*h +
    //                a200*l^2 + a020*p^2 + a002*h^2
   ****************************************************************/
  real8 polyval(
    const matrix<real8> &C,
    const real8 line, 
    const real8 pixel,
    const real8 height) const
    {
    #ifdef __DEBUG
    TRACE_FUNCTION("baseline::polyval (BK 05-Mar-2005)")
    #endif
    if (C.size()!=10) throw("error");
    return
      C(0,0) +
      C(1,0)*line       + C(2,0)*pixel       + C(3,0)*height +
      C(4,0)*line*pixel + C(5,0)*line*height + C(6,0)*pixel*height +
      C(7,0)*sqr(line)  + C(8,0)*sqr(pixel)  + C(9,0)*sqr(height);
    }; // END polyval


  /****************************************************************
   ****************************************************************/
  void BBparBperpTheta(real8 &B, real8 &Bpar, real8 &Bperp, real8 &theta,
             const cn Master, const cn Point, const cn Slave) const
    {
    #ifdef __DEBUG
    TRACE_FUNCTION("BBparBperp (BK 05-Mar-2005)")
    #endif
    B                  = Master.dist(Slave);// baseline. abs. value (in plane M,P,S)
    const real8 range1 = Master.dist(Point);
    const real8 range2 = Slave.dist(Point);
    Bpar               = range1-range2;                 // parallel baseline, sign ok
    const cn r1        = Master.min(Point);// points from P to M
    const cn r2        = Slave.min(Point);
    theta              = Master.angle(r1);// viewing angle
    Bperp              = sqr(B)-sqr(Bpar);
    if (Bperp < 0.0) Bperp=0.0;
    else Bperp = (theta > Master.angle(r2)) ?     // perpendicular baseline, sign ok
       sqrt(Bperp) : -sqrt(Bperp);
    }; // END BBparBperpTheta



  /****************************************************************
   * returns incidence angle in radians based on coordinate of    *
   * point P on ellips and point M in orbit                       *
   #%// Bert Kampes, 06-Apr-2005
   ****************************************************************/
  real8 IncidenceAngle(const cn Master, const cn Point) const
    {
    #ifdef __DEBUG
    TRACE_FUNCTION("IncidenceAngle (BK 05-Mar-2005)")
    #endif
    const cn r1     = Master.min(Point);// points from P to M
    const real8 inc = Point.angle(r1);// incidence angle (assume P on ellips)
    return inc;
    }; // END BBparBperpTheta


  // ______ Public function of bk_message class ______
  public:
    // ___Constructor/Destructor___
    BASELINE()
      {
      TRACE_FUNCTION("BASELINE() (BK 06-Mar-2005)");
      initialized = false;//          
      N_coeffs    = 10;
      H_min       =    0.0;//      height at which baseline is computed.
      H_max       = 5000.0;//      height at which baseline is computed.
      master_wavelength = 0.0;
      orbit_convergence = 0.0;//   tmp for now constant
      orbit_heading     = 0.0;//       tmp for now NOT USED
      }

    // ___Constructor/Destructor___
    ~BASELINE()                 //{;}// nothing to destruct
      {
      TRACE_FUNCTION("~BASELINE() (BK 06-Mar-2005)");
      ;// nothing to destruct ?
      }// dealloc


    // ___Helper functions___
    void model_parameters(
        const slcimage &master,
        const slcimage &slave,
        orbit &masterorbit,
        orbit &slaveorbit, 
        //const productinfo &ifg,  <--- includes pointer to master,slave, and they to orbits?
        const input_ell &ellips)
      {
      TRACE_FUNCTION("model_parameters (BK 06-Mar-2005)");
      if (masterorbit.is_initialized() == false)
        {
        DEBUG.print("Baseline cannot be computed, master orbit not initialized.");
        return;
        }
      if (slaveorbit.is_initialized()  == false)
        {
        DEBUG.print("Baseline cannot be computed, slave orbit not initialized.");
        return;
        }

      // --- Get on with it ---------------------------------------
      initialized = true;//          
      master_wavelength = master.wavelength;//
      // ______ Model r=nearrange+drange_dp*p, p starts at 1 ______
      nearrange  = master.pix2range(1.0);
      drange_dp  = master.pix2range(2.0)-master.pix2range(1.0);
      nearrange -= drange_dp;// (p starts at 1) 


      // ______ Loop counters ______
      register int32 cnt      = 0;// matrix index
      const int N_pointsL     = 10;// every 10km in azimuth
      const int N_pointsP     = 10;// every 10km ground range
      const int N_heights     = 4;// one more level than required for poly
      const real8 deltapixels = master.currentwindow.pixels() / N_pointsP;
      const real8 deltalines  = master.currentwindow.lines()  / N_pointsL;
      const real8 deltaheight = (H_max-H_min)                 / N_heights;
      
      // ______ Matrices for modeling Bperp (BK 21-mar-01) ______
      matrix<real8> BPERP(N_pointsL*N_pointsP*N_heights,1);// perpendicular baseline
      matrix<real8> BPAR(N_pointsL*N_pointsP*N_heights,1);// parallel baseline
      matrix<real8> THETA(N_pointsL*N_pointsP*N_heights,1);// viewing angle
      matrix<real8> THETA_INC(N_pointsL*N_pointsP*N_heights,1);// inc. angle
      matrix<real8> AMATRIX(N_pointsL*N_pointsP*N_heights,N_coeffs);// design matrix

      // ______ Loop over heights, lines, pixels to compute baseline param. ______
      for (register int32 k=0; k<N_heights; ++k) // height levels
        {
        const real8 HEIGHT = H_min + k*deltaheight;
        input_ell ELLIPS;
        ELLIPS.a   = ellips.a + HEIGHT;
        ELLIPS.b   = ellips.b + HEIGHT;
        ELLIPS.ecc1st_sqr();
        ELLIPS.ecc2nd_sqr();

        for (register int32 i=0; i<N_pointsL; ++i) // azimuthlines
          {
          const real8 line = master.currentwindow.linelo + i*deltalines;
          cn P;          // point, returned by lp2xyz
          real8 s_tazi;  // returned by xyz2t
          real8 s_trange;// returned by xyz2t
            const int32 MAXITER   = 10;
          const real8 CRITERPOS = 1e-6;
          const real8 CRITERTIM = 1e-10;
          // ______ Azimuth time for this line ______
          const real8 m_tazi = master.line2ta(line);
          // ______ xyz for master satellite from time ______
          const cn M         = masterorbit.getxyz(m_tazi);
          // ______ Loop over a pixels to compute baseline param. ______
          for (register int32 j=0; j<N_pointsP; ++j) // rangepixels
            {
            const real8 pixel = master.currentwindow.pixlo + j*deltapixels;
            // ______ Range time for this pixel ______
            //const real8 m_trange = master.pix2tr(pixel);
            lp2xyz(line,pixel,ELLIPS,master,masterorbit,P,MAXITER,CRITERPOS);
            // ______ Compute xyz for slave satellite from P ______
            xyz2t(s_tazi,s_trange,slave,slaveorbit,P,MAXITER,CRITERTIM);
            // ______ Slave position ______
            const cn S              = slaveorbit.getxyz(s_tazi);
            // ______ Compute angle between near parallel orbits ______
            const cn Mdot           = masterorbit.getxyzdot(m_tazi);
            const cn Sdot           = slaveorbit.getxyzdot(s_tazi);
            const real8 angleorbits = Mdot.angle(Sdot);
            DEBUG << "Angle between orbits master-slave (at l,p="
                 << line << "," << pixel << ") = "
                 << rad2deg(angleorbits) << " [deg]";
            DEBUG.print();
            orbit_convergence = angleorbits;// assume constant; store in member
            //const real8 heading = angle(Mdot,[1 0 0])?
            //orbit_heading = 0.0;// not yet used

            // ====== The baseline parameters, derived from the positions (x,y,z) ======
            // ______ alpha is angle counterclockwize(B, vlak met normal=rho1=rho2)
            // ______ theta is angle counterclockwize(rho1=M, r1=M-P, r2=S-P)
            real8 B,Bpar,Bperp,theta;
            BBparBperpTheta(B,Bpar,Bperp,theta,M,P,S);// return B etc.
            const real8 theta_inc = IncidenceAngle(M,P);// [rad]

            // ______ Modelling of Bperp(l,p) = a00 + a10*l + a01*p ______
            BPERP(cnt,0)     = Bperp;
            BPAR(cnt,0)      = Bpar;
            THETA(cnt,0)     = theta;
            THETA_INC(cnt,0) = theta_inc;
            // --- B(l,p,h) = a000 + 
            //                a100*l   + a010*p   + a001*h   +
            //                a110*l*p + a101*l*h + a011*p*h +
            //                a200*l^2 + a020*p^2 + a002*h^2
            AMATRIX(cnt,0)   = 1.0;
            AMATRIX(cnt,1)   = line;
            AMATRIX(cnt,2)   = pixel;
            AMATRIX(cnt,3)   = HEIGHT;
            AMATRIX(cnt,4)   = line*pixel;
            AMATRIX(cnt,5)   = line*HEIGHT;
            AMATRIX(cnt,6)   = pixel*HEIGHT;
            AMATRIX(cnt,7)   = sqr(line);
            AMATRIX(cnt,8)   = sqr(pixel);
            AMATRIX(cnt,9)   = sqr(HEIGHT);
            cnt++;
            // ______ B/alpha representation of baseline ______
            const real8 alpha  = (Bpar==0 && Bperp==0) ?
              NaN :
              theta - atan2(Bpar,Bperp);            // sign ok atan2
            // ______ hor/vert representation of baseline ______
            const real8 Bh    = B * cos(alpha);                        // sign ok
            const real8 Bv    = B * sin(alpha);                        // sign ok
            // ______ Height ambiguity: [h] = -lambda/4pi * (r1sin(theta)/Bperp) * phi==2pi ______
            const real8 hambiguity = (Bperp==0) ? 
              Inf :
              -master.wavelength*(M.min(P)).norm()*sin(theta)/(2.0*Bperp);
      
            // ______ Some extra info if in debug mode ______
            // ====== Write output to screen as DEBUG ======
            DEBUG << "The baseline parameters for (l,p,h) = " 
                  << line << ", " << pixel << ", " << HEIGHT;
            DEBUG.print();
            DEBUG << "\talpha (deg), baseline: \t" << rad2deg(alpha) << " \t" << B;
            DEBUG.print();
            DEBUG << "\tBpar, Bperp:      \t" << Bpar << " \t" << Bperp;
            DEBUG.print();
            DEBUG << "\tBh, Bv:           \t" << Bh << " \t" << Bv;
            DEBUG.print();
            DEBUG << "\tHeight ambiguity: \t" << hambiguity;
            DEBUG.print();
            DEBUG << "\ttheta (deg):      \t" << rad2deg(theta);
            DEBUG.print();
            DEBUG << "\ttheta_inc (deg):  \t" << rad2deg(theta_inc);
            DEBUG.print();
            DEBUG.precision(10);
            DEBUG.width(11);
            DEBUG.rewind();
            DEBUG << "\tM (x,y,z) = " << M.x << ", " << M.y << ", " << M.z;
            DEBUG.print();
            DEBUG << "\tS (x,y,z) = " << S.x << ", " << S.y << ", " << S.z;
            DEBUG.print();
            DEBUG << "\tP (x,y,z) = " << P.x << ", " << P.y << ", " << P.z;
            DEBUG.print();
            DEBUG.reset();// reset width/precision/pointer
            } // loop pixels
          } // loop lines
        } // loop heights

      // ====== Model the Bperp as 2d polynomial of degree 1 ======
      matrix<real8> N        = matTxmat(AMATRIX,AMATRIX);
      matrix<real8> rhsBperp = matTxmat(AMATRIX,BPERP);
      matrix<real8> rhsBpar  = matTxmat(AMATRIX,BPAR);
      matrix<real8> rhsT     = matTxmat(AMATRIX,THETA);
      matrix<real8> rhsT_INC = matTxmat(AMATRIX,THETA_INC);
      matrix<real8> Qx_hat   = N;
      choles(Qx_hat);               // Cholesky factorisation normalmatrix
      solvechol(Qx_hat,rhsBperp);   // Solution Bperp coefficients in rhsB
      solvechol(Qx_hat,rhsBpar);    // Solution Theta coefficients in rhsT
      solvechol(Qx_hat,rhsT);       // Solution Theta coefficients in rhsT
      solvechol(Qx_hat,rhsT_INC);   // Solution Theta_inc coefficients in rhsT_INC
      invertchol(Qx_hat);           // Covariance matrix of unknowns

      ;// === Copy estimated coefficients to private members ===
      BPERP_cf     = rhsBperp;//
      BPAR_cf      = rhsBpar;//
      THETA_cf     = rhsT;//
      THETA_INC_cf = rhsT_INC;//
    
      // ______Test inverse______
      for (int32 i=0; i<Qx_hat.lines(); i++)
        for (int32 j=0; j<i; j++)
          Qx_hat(j,i) = Qx_hat(i,j);// repair matrix
      const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
      DEBUG << "baseline: max(abs(N*inv(N)-I)) = " << maxdev;
      DEBUG.print();
      /*
      matrix<real8> unity = N * Qx_hat;
      real8 maxdev = abs(unity(0,0)-1);
      for (int32 i=1; i<unity.lines(); i++)
        {
        if (abs(unity(i,i)-1) > maxdev)
          maxdev = abs(unity(i,i)-1);
        for (int32 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)
        {
        WARNING << "baseline: max. deviation N*inv(N) from unity = "
           << maxdev << ". This is larger than .01: do not use this.";
        WARNING.print();// no error, no problem
        }
      else if (maxdev > .001)
        {
        WARNING << "baseline: max. deviation N*inv(N) from unity = "
           << maxdev << ". This is between 0.01 and 0.001 (do not use it)";
        WARNING.print();
        }
      // ______Some other stuff, scale is ok______
      //matrix<real8> Qy_hat = AMATRIX * (matxmatT(Qx_hat,AMATRIX));
      matrix<real8> y_hatBperp  = AMATRIX * rhsBperp;
      matrix<real8> e_hatBperp  = BPERP - y_hatBperp;
      //matrix<real8> Qe_hat  = Qy - Qy_hat;
      //matrix<real8> y_hatT  = AMATRIX * rhsT;
      //matrix<real8> e_hatT  = THETA - y_hatT;
    
    
      // ______ Output solution and give max error ______
      // --- B(l,p,h) = a000 + 
      //                a100*l   + a010*p   + a001*h   +
      //                a110*l*p + a101*l*h + a011*p*h +
      //                a200*l^2 + a020*p^2 + a002*h^2
      DEBUG.print("--------------------");
      DEBUG.print("Result of modeling: Bperp(l,p) = a000 + a100*l + a010*p + a001*h + ");
      DEBUG.print(" a110*l*p + a101*l*h + a011*p*h + a200*l^2 + a020*p^2 + a002*h^2");
      DEBUG.print("l,p,h in unnormalized, absolute, coordinates (1:N).");
      DEBUG << "Bperp_a000 = " << rhsBperp(0,0); DEBUG.print();
      DEBUG << "Bperp_a100 = " << rhsBperp(1,0); DEBUG.print();
      DEBUG << "Bperp_a010 = " << rhsBperp(2,0); DEBUG.print();
      DEBUG << "Bperp_a001 = " << rhsBperp(3,0); DEBUG.print();
      DEBUG << "Bperp_a110 = " << rhsBperp(4,0); DEBUG.print();
      DEBUG << "Bperp_a101 = " << rhsBperp(5,0); DEBUG.print();
      DEBUG << "Bperp_a011 = " << rhsBperp(6,0); DEBUG.print();
      DEBUG << "Bperp_a200 = " << rhsBperp(7,0); DEBUG.print();
      DEBUG << "Bperp_a020 = " << rhsBperp(8,0); DEBUG.print();
      DEBUG << "Bperp_a002 = " << rhsBperp(9,0); DEBUG.print();
      real8 maxerr = max(abs(e_hatBperp));
      if (maxerr > 2.00)// 
        {
        WARNING << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
        WARNING.print();
        }
      else
        {
        INFO << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
        INFO.print();
        }
      DEBUG.print();
      DEBUG.print("--------------------");
      DEBUG.print("Range: r(p) = r0 + dr*p");
      DEBUG.print("l and p in unnormalized, absolute, coordinates (1:N).");
      const real8 range1    = master.pix2range(1.0);
      const real8 range5000 = master.pix2range(5000.0);
      const real8 drange    = (range5000-range1)/5000;
      DEBUG << "range = " << range1-drange << " + " << drange << "*p";
      DEBUG.print();
      // ====== Tidy up ======
      }; // END model_parameters()

    // ___ Return RANGE to user ___
    inline real8 get_range(const real8 pixel) const
      {
      return nearrange + drange_dp*pixel;
      };// END get_bperp()

    // === modeled quantities ===
    // --- B(l,p,h) = a000 + 
    //                a100*l   + a010*p   + a001*h   +
    //                a110*l*p + a101*l*h + a011*p*h +
    //                a200*l^2 + a020*p^2 + a002*h^2
    // ___ Return BPERP to user ___
    inline real8 get_bperp(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      return polyval(BPERP_cf, line, pixel, height);
      };
    // ___ Return BPAR to user ___
    inline real8 get_bpar(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      return polyval(BPAR_cf, line, pixel, height);
      };
    // ___ Return THETA to user ___
    inline real8 get_theta(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      return polyval(THETA_cf, line, pixel, height);
      };
    // ___ Return THETA_INC to user ___
    inline real8 get_theta_inc(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      return polyval(THETA_INC_cf, line, pixel, height);
      };

    // === derived quantities ===
    // ___ Return B to user ___
    inline real8 get_b(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      return sqrt(sqr(get_bpar(line,pixel,height))+sqr(get_bperp(line,pixel,height)));
      };// END get_b()

    // ___ Return alpha baseline orientation to user ___
    inline real8 get_alpha(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      const real8 Bperp     = get_bperp(line,pixel,height);
      const real8 Bpar      = get_bpar(line,pixel,height);
      const real8 B         = get_b(line,pixel,height);
      const real8 theta     = get_theta(line,pixel,height);
      const real8 alpha     = (Bpar==0 && Bperp==0) ?
        NaN :
        theta-atan2(Bpar,Bperp);            // sign ok atan2
      return alpha;// sign ok
      };// END get_bhor()

    // ___ Return Bh to user ___
    inline real8 get_bhor(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      const real8 Bperp     = get_bperp(line,pixel,height);
      const real8 Bpar      = get_bpar(line,pixel,height);
      const real8 B         = get_b(line,pixel,height);
      const real8 theta     = get_theta(line,pixel,height);
      const real8 alpha     = get_alpha(line,pixel,height);
      return B*cos(alpha);// sign ok
      };// END get_bhor()

    // ___ Return Bv to user ___
    inline real8 get_bvert(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      const real8 Bperp     = get_bperp(line,pixel,height);
      const real8 Bpar      = get_bpar(line,pixel,height);
      const real8 B         = get_b(line,pixel,height);
      const real8 theta     = get_theta(line,pixel,height);
      const real8 alpha     = get_alpha(line,pixel,height);
      return B*sin(alpha);// sign ok
      };// END get_bvert()

    // ___ Return Height ambiguity to user ___
    inline real8 get_hamb(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      //const real8 theta     =  get_theta(line,pixel,height);
      const real8 theta_inc =  get_theta_inc(line,pixel,height);
      const real8 Bperp     =  get_bperp(line,pixel,height);
      const real8 range_MP  =  get_range(pixel);// >
      const real8 h_amb     = (Bperp==0) ?
         Inf : // inf
        -master_wavelength*range_MP*sin(theta_inc)/(2.0*Bperp);// this is wrt local
        //-master_wavelength*range_MP*sin(theta)/(2.0*Bperp);// this is wrt
      return h_amb;
      };// END get_hamb()

    // ___ Return orbit convergence to user ___
    inline real8 get_orb_conv(const real8 line, const real8 pixel, const real8 height=0.0) const
      {
      // do not use l,p..
      return orbit_convergence;
      };// END get_orb_conv()


    // --- Dump overview of all ---
    void dump(const real8 line, const real8 pixel, const real8 height=0.0)
      {
      if (initialized==false)
        {
        DEBUG.print("Exiting dumpbaseline, not initialized.");
        return;
        }
        // ______ Modeled quantities ______
        const real8 Bperp     = get_bperp(line,pixel,height);
        const real8 Bpar      = get_bpar(line,pixel,height);
        const real8 theta     = get_theta(line,pixel,height);
        const real8 theta_inc = get_theta_inc(line,pixel,height);
        // ______ Derived quantities ______
        const real8 B         = get_b(line,pixel,height);
        const real8 alpha     = get_alpha(line,pixel,height);
        const real8 Bh        = get_bhor(line,pixel,height);
        const real8 Bv        = get_bvert(line,pixel,height);
        const real8 h_amb     = get_hamb(line,pixel,height);
        // ______ Height ambiguity: [h] = -lambda/4pi * (r1sin(theta)/Bperp) * phi==2pi ______
        // ====== Write output to screen as INFO ======
        INFO << "The baseline parameters for (l,p,h) = " 
             << line << ", " << pixel << ", " << height;
        INFO.print();
        INFO << "\tBpar, Bperp:      \t" << Bpar << " \t" << Bperp;
        INFO.print();
        DEBUG << "\tB, alpha (deg):  \t" << B << " \t" << rad2deg(alpha);
        DEBUG.print();
        DEBUG << "\tBh, Bv:          \t" << Bh << " \t" << Bv;
        DEBUG.print();
        INFO << "\tHeight ambiguity: \t" << h_amb;
        INFO.print();
        INFO << "\tLook angle (deg): \t" << rad2deg(theta);
        INFO.print();
        DEBUG << "\tIncidence angle (deg): \t" << rad2deg(theta_inc);
        DEBUG.print();
      }// END dump()

  };// eof BASELINE class

#endif // BK_BASELINE_H

---