Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

bk_baseline.hh

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 1999-2005 Bert Kampes
00003  * Copyright (c) 1999-2005 Delft University of Technology, The Netherlands
00004  *
00005  * This file is part of Doris, the Delft o-o radar interferometric software.
00006  *
00007  * Doris program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * Doris is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00020  *
00021  * Publications that contain results produced by the Doris software should
00022  * contain an acknowledgment. (For example: The interferometric processing
00023  * was performed using the freely available Doris software package developed
00024  * by the Delft Institute for Earth-Oriented Space Research (DEOS), Delft
00025  * University of Technology, or include a reference to: Bert Kampes and
00026  * Stefania Usai. \"Doris: The Delft Object-oriented Radar Interferometric
00027  * software.\" In: proceedings 2nd ITC ORS symposium, August 1999. (cdrom)).
00028  *
00029  */
00030 /****************************************************************
00031  * $Source: /users/kampes/DEVELOP/DORIS/doris/src/RCS/bk_baseline.hh,v $ *
00032  * $Revision: 1.3 $ *
00033  * $Date: 2005/04/11 13:48:45 $ *
00034  * $Author: kampes $ *
00035  *
00036  * definition of BASELINE class *
00037  * 
00038  * During computations these should be used with care, they
00039  * may be a bit approximate.  Better is to do what Doris does,
00040  * i.e., evaluate positions and relate that to timing.
00041  * But i guess it is OK.  Use this for ultra fast geocoding...
00042  *
00043  #%// Bert Kampes, 06-Apr-2005
00044  ****************************************************************/
00045 
00046 #ifndef BK_BASELINE_H
00047 #define BK_BASELINE_H
00048 using namespace std;
00049 
00050 #include <iostream>    // cout
00051 #include <cstdlib>     // exit()
00052 
00053 #include "constants.hh"   // types etc.
00054 #include "bk_messages.hh" // info etc.
00055 #include "slcimage.hh"    // my slc image class
00056 #include "productinfo.hh" // my products class
00057 #include "orbitbk.hh"     // my orbit class
00058 #include "conversion.hh"  // rad2deg()
00059 #include "utilities.hh"   // eye()
00060 
00061 
00062 // BASELINE is a new type (class) that is initialized using orbits
00063 // and then either models the baseline parameters such as Bperp
00064 // or can give them exact.
00065 // Usage in the programs is something like in main do:
00066 // BASELINE baseline;
00067 // baseline.init(orbit1,orbit2,product?master?);
00068 // and pass baseline to subprograms, there use baseline.get_bperp(x,y) etc.
00069 // Bert Kampes, 31-Mar-2005
00070 
00071 
00072 // ---------------------------------------------------------------- 
00073 // ---------------------------------------------------------------- 
00074 // ---------------------------------------------------------------- 
00075 class BASELINE
00076   {
00077   // ______ Private data of BASELINE class ______
00078   private:
00079     bool initialized;//          
00080     real8 master_wavelength;//   tmp for now used for h_amb
00081     real8 nearrange;// range=nearrange+drange_dp*pixel
00082     real8 drange_dp;// range=nearrange+drange_dp*pixel
00083     real8 orbit_convergence;//   tmp for now constant
00084     real8 orbit_heading;//       tmp for now NOT USED
00085     real8 H_min;//            height at which parameters are modeled
00086     real8 H_max;//    to model phase=f(line,pix,hei) and hei=g(line,pix,phase)
00087     // --- B(l,p,h) = a000 + 
00088     //                a100*l   + a010*p   + a001*h   +
00089     //                a110*l*p + a101*l*h + a011*p*h +
00090     //                a200*l^2 + a020*p^2 + a002*h^2
00091     int32 N_coeffs;//  ==10 degree of 3D-poly to model B(l,p,h)
00092     // --- Coefficients ---
00093     matrix<real8> BPERP_cf;// perpendicular baseline
00094     matrix<real8> BPAR_cf;// parallel baseline
00095     matrix<real8> THETA_cf;// viewing angle
00096     matrix<real8> THETA_INC_cf;// incidence angle to satellite
00097     //real8 avg_height_ambiguity;//center height ambiguity
00098 
00099     // copy constructor; copy matrices etc as public.
00100     BASELINE(const BASELINE& A)
00101       {};// prevent copy constructor usage by privatizing
00102     //orbit *master;
00103     //orbit *slave;
00104 
00105   /****************************************************************
00106     // --- B(l,p,h) = a000 + 
00107     //                a100*l   + a010*p   + a001*h   +
00108     //                a110*l*p + a101*l*h + a011*p*h +
00109     //                a200*l^2 + a020*p^2 + a002*h^2
00110    ****************************************************************/
00111   real8 polyval(
00112     const matrix<real8> &C,
00113     const real8 line, 
00114     const real8 pixel,
00115     const real8 height) const
00116     {
00117     #ifdef __DEBUG
00118     TRACE_FUNCTION("baseline::polyval (BK 05-Mar-2005)")
00119     #endif
00120     if (C.size()!=10) throw("error");
00121     return
00122       C(0,0) +
00123       C(1,0)*line       + C(2,0)*pixel       + C(3,0)*height +
00124       C(4,0)*line*pixel + C(5,0)*line*height + C(6,0)*pixel*height +
00125       C(7,0)*sqr(line)  + C(8,0)*sqr(pixel)  + C(9,0)*sqr(height);
00126     }; // END polyval
00127 
00128 
00129   /****************************************************************
00130    ****************************************************************/
00131   void BBparBperpTheta(real8 &B, real8 &Bpar, real8 &Bperp, real8 &theta,
00132              const cn Master, const cn Point, const cn Slave) const
00133     {
00134     #ifdef __DEBUG
00135     TRACE_FUNCTION("BBparBperp (BK 05-Mar-2005)")
00136     #endif
00137     B                  = Master.dist(Slave);// baseline. abs. value (in plane M,P,S)
00138     const real8 range1 = Master.dist(Point);
00139     const real8 range2 = Slave.dist(Point);
00140     Bpar               = range1-range2;                 // parallel baseline, sign ok
00141     const cn r1        = Master.min(Point);// points from P to M
00142     const cn r2        = Slave.min(Point);
00143     theta              = Master.angle(r1);// viewing angle
00144     Bperp              = sqr(B)-sqr(Bpar);
00145     if (Bperp < 0.0) Bperp=0.0;
00146     else Bperp = (theta > Master.angle(r2)) ?     // perpendicular baseline, sign ok
00147        sqrt(Bperp) : -sqrt(Bperp);
00148     }; // END BBparBperpTheta
00149 
00150 
00151 
00152   /****************************************************************
00153    * returns incidence angle in radians based on coordinate of    *
00154    * point P on ellips and point M in orbit                       *
00155    #%// Bert Kampes, 06-Apr-2005
00156    ****************************************************************/
00157   real8 IncidenceAngle(const cn Master, const cn Point) const
00158     {
00159     #ifdef __DEBUG
00160     TRACE_FUNCTION("IncidenceAngle (BK 05-Mar-2005)")
00161     #endif
00162     const cn r1     = Master.min(Point);// points from P to M
00163     const real8 inc = Point.angle(r1);// incidence angle (assume P on ellips)
00164     return inc;
00165     }; // END BBparBperpTheta
00166 
00167 
00168   // ______ Public function of bk_message class ______
00169   public:
00170     // ___Constructor/Destructor___
00171     BASELINE()
00172       {
00173       TRACE_FUNCTION("BASELINE() (BK 06-Mar-2005)");
00174       initialized = false;//          
00175       N_coeffs    = 10;
00176       H_min       =    0.0;//      height at which baseline is computed.
00177       H_max       = 5000.0;//      height at which baseline is computed.
00178       master_wavelength = 0.0;
00179       orbit_convergence = 0.0;//   tmp for now constant
00180       orbit_heading     = 0.0;//       tmp for now NOT USED
00181       }
00182 
00183     // ___Constructor/Destructor___
00184     ~BASELINE()                 //{;}// nothing to destruct
00185       {
00186       TRACE_FUNCTION("~BASELINE() (BK 06-Mar-2005)");
00187       ;// nothing to destruct ?
00188       }// dealloc
00189 
00190 
00191     // ___Helper functions___
00192     void model_parameters(
00193         const slcimage &master,
00194         const slcimage &slave,
00195         orbit &masterorbit,
00196         orbit &slaveorbit, 
00197         //const productinfo &ifg,  <--- includes pointer to master,slave, and they to orbits?
00198         const input_ell &ellips)
00199       {
00200       TRACE_FUNCTION("model_parameters (BK 06-Mar-2005)");
00201       if (masterorbit.is_initialized() == false)
00202         {
00203         DEBUG.print("Baseline cannot be computed, master orbit not initialized.");
00204         return;
00205         }
00206       if (slaveorbit.is_initialized()  == false)
00207         {
00208         DEBUG.print("Baseline cannot be computed, slave orbit not initialized.");
00209         return;
00210         }
00211 
00212       // --- Get on with it ---------------------------------------
00213       initialized = true;//          
00214       master_wavelength = master.wavelength;//
00215       // ______ Model r=nearrange+drange_dp*p, p starts at 1 ______
00216       nearrange  = master.pix2range(1.0);
00217       drange_dp  = master.pix2range(2.0)-master.pix2range(1.0);
00218       nearrange -= drange_dp;// (p starts at 1) 
00219 
00220 
00221       // ______ Loop counters ______
00222       register int32 cnt      = 0;// matrix index
00223       const int N_pointsL     = 10;// every 10km in azimuth
00224       const int N_pointsP     = 10;// every 10km ground range
00225       const int N_heights     = 4;// one more level than required for poly
00226       const real8 deltapixels = master.currentwindow.pixels() / N_pointsP;
00227       const real8 deltalines  = master.currentwindow.lines()  / N_pointsL;
00228       const real8 deltaheight = (H_max-H_min)                 / N_heights;
00229       
00230       // ______ Matrices for modeling Bperp (BK 21-mar-01) ______
00231       matrix<real8> BPERP(N_pointsL*N_pointsP*N_heights,1);// perpendicular baseline
00232       matrix<real8> BPAR(N_pointsL*N_pointsP*N_heights,1);// parallel baseline
00233       matrix<real8> THETA(N_pointsL*N_pointsP*N_heights,1);// viewing angle
00234       matrix<real8> THETA_INC(N_pointsL*N_pointsP*N_heights,1);// inc. angle
00235       matrix<real8> AMATRIX(N_pointsL*N_pointsP*N_heights,N_coeffs);// design matrix
00236 
00237       // ______ Loop over heights, lines, pixels to compute baseline param. ______
00238       for (register int32 k=0; k<N_heights; ++k) // height levels
00239         {
00240         const real8 HEIGHT = H_min + k*deltaheight;
00241         input_ell ELLIPS;
00242         ELLIPS.a   = ellips.a + HEIGHT;
00243         ELLIPS.b   = ellips.b + HEIGHT;
00244         ELLIPS.ecc1st_sqr();
00245         ELLIPS.ecc2nd_sqr();
00246 
00247         for (register int32 i=0; i<N_pointsL; ++i) // azimuthlines
00248           {
00249           const real8 line = master.currentwindow.linelo + i*deltalines;
00250           cn P;          // point, returned by lp2xyz
00251           real8 s_tazi;  // returned by xyz2t
00252           real8 s_trange;// returned by xyz2t
00253             const int32 MAXITER   = 10;
00254           const real8 CRITERPOS = 1e-6;
00255           const real8 CRITERTIM = 1e-10;
00256           // ______ Azimuth time for this line ______
00257           const real8 m_tazi = master.line2ta(line);
00258           // ______ xyz for master satellite from time ______
00259           const cn M         = masterorbit.getxyz(m_tazi);
00260           // ______ Loop over a pixels to compute baseline param. ______
00261           for (register int32 j=0; j<N_pointsP; ++j) // rangepixels
00262             {
00263             const real8 pixel = master.currentwindow.pixlo + j*deltapixels;
00264             // ______ Range time for this pixel ______
00265             //const real8 m_trange = master.pix2tr(pixel);
00266             lp2xyz(line,pixel,ELLIPS,master,masterorbit,P,MAXITER,CRITERPOS);
00267             // ______ Compute xyz for slave satellite from P ______
00268             xyz2t(s_tazi,s_trange,slave,slaveorbit,P,MAXITER,CRITERTIM);
00269             // ______ Slave position ______
00270             const cn S              = slaveorbit.getxyz(s_tazi);
00271             // ______ Compute angle between near parallel orbits ______
00272             const cn Mdot           = masterorbit.getxyzdot(m_tazi);
00273             const cn Sdot           = slaveorbit.getxyzdot(s_tazi);
00274             const real8 angleorbits = Mdot.angle(Sdot);
00275             DEBUG << "Angle between orbits master-slave (at l,p="
00276                  << line << "," << pixel << ") = "
00277                  << rad2deg(angleorbits) << " [deg]";
00278             DEBUG.print();
00279             orbit_convergence = angleorbits;// assume constant; store in member
00280             //const real8 heading = angle(Mdot,[1 0 0])?
00281             //orbit_heading = 0.0;// not yet used
00282 
00283             // ====== The baseline parameters, derived from the positions (x,y,z) ======
00284             // ______ alpha is angle counterclockwize(B, vlak met normal=rho1=rho2)
00285             // ______ theta is angle counterclockwize(rho1=M, r1=M-P, r2=S-P)
00286             real8 B,Bpar,Bperp,theta;
00287             BBparBperpTheta(B,Bpar,Bperp,theta,M,P,S);// return B etc.
00288             const real8 theta_inc = IncidenceAngle(M,P);// [rad]
00289 
00290             // ______ Modelling of Bperp(l,p) = a00 + a10*l + a01*p ______
00291             BPERP(cnt,0)     = Bperp;
00292             BPAR(cnt,0)      = Bpar;
00293             THETA(cnt,0)     = theta;
00294             THETA_INC(cnt,0) = theta_inc;
00295             // --- B(l,p,h) = a000 + 
00296             //                a100*l   + a010*p   + a001*h   +
00297             //                a110*l*p + a101*l*h + a011*p*h +
00298             //                a200*l^2 + a020*p^2 + a002*h^2
00299             AMATRIX(cnt,0)   = 1.0;
00300             AMATRIX(cnt,1)   = line;
00301             AMATRIX(cnt,2)   = pixel;
00302             AMATRIX(cnt,3)   = HEIGHT;
00303             AMATRIX(cnt,4)   = line*pixel;
00304             AMATRIX(cnt,5)   = line*HEIGHT;
00305             AMATRIX(cnt,6)   = pixel*HEIGHT;
00306             AMATRIX(cnt,7)   = sqr(line);
00307             AMATRIX(cnt,8)   = sqr(pixel);
00308             AMATRIX(cnt,9)   = sqr(HEIGHT);
00309             cnt++;
00310             // ______ B/alpha representation of baseline ______
00311             const real8 alpha  = (Bpar==0 && Bperp==0) ?
00312               NaN :
00313               theta - atan2(Bpar,Bperp);            // sign ok atan2
00314             // ______ hor/vert representation of baseline ______
00315             const real8 Bh    = B * cos(alpha);                        // sign ok
00316             const real8 Bv    = B * sin(alpha);                        // sign ok
00317             // ______ Height ambiguity: [h] = -lambda/4pi * (r1sin(theta)/Bperp) * phi==2pi ______
00318             const real8 hambiguity = (Bperp==0) ? 
00319               Inf :
00320               -master.wavelength*(M.min(P)).norm()*sin(theta)/(2.0*Bperp);
00321       
00322             // ______ Some extra info if in debug mode ______
00323             // ====== Write output to screen as DEBUG ======
00324             DEBUG << "The baseline parameters for (l,p,h) = " 
00325                   << line << ", " << pixel << ", " << HEIGHT;
00326             DEBUG.print();
00327             DEBUG << "\talpha (deg), baseline: \t" << rad2deg(alpha) << " \t" << B;
00328             DEBUG.print();
00329             DEBUG << "\tBpar, Bperp:      \t" << Bpar << " \t" << Bperp;
00330             DEBUG.print();
00331             DEBUG << "\tBh, Bv:           \t" << Bh << " \t" << Bv;
00332             DEBUG.print();
00333             DEBUG << "\tHeight ambiguity: \t" << hambiguity;
00334             DEBUG.print();
00335             DEBUG << "\ttheta (deg):      \t" << rad2deg(theta);
00336             DEBUG.print();
00337             DEBUG << "\ttheta_inc (deg):  \t" << rad2deg(theta_inc);
00338             DEBUG.print();
00339             DEBUG.precision(10);
00340             DEBUG.width(11);
00341             DEBUG.rewind();
00342             DEBUG << "\tM (x,y,z) = " << M.x << ", " << M.y << ", " << M.z;
00343             DEBUG.print();
00344             DEBUG << "\tS (x,y,z) = " << S.x << ", " << S.y << ", " << S.z;
00345             DEBUG.print();
00346             DEBUG << "\tP (x,y,z) = " << P.x << ", " << P.y << ", " << P.z;
00347             DEBUG.print();
00348             DEBUG.reset();// reset width/precision/pointer
00349             } // loop pixels
00350           } // loop lines
00351         } // loop heights
00352 
00353       // ====== Model the Bperp as 2d polynomial of degree 1 ======
00354       matrix<real8> N        = matTxmat(AMATRIX,AMATRIX);
00355       matrix<real8> rhsBperp = matTxmat(AMATRIX,BPERP);
00356       matrix<real8> rhsBpar  = matTxmat(AMATRIX,BPAR);
00357       matrix<real8> rhsT     = matTxmat(AMATRIX,THETA);
00358       matrix<real8> rhsT_INC = matTxmat(AMATRIX,THETA_INC);
00359       matrix<real8> Qx_hat   = N;
00360       choles(Qx_hat);               // Cholesky factorisation normalmatrix
00361       solvechol(Qx_hat,rhsBperp);   // Solution Bperp coefficients in rhsB
00362       solvechol(Qx_hat,rhsBpar);    // Solution Theta coefficients in rhsT
00363       solvechol(Qx_hat,rhsT);       // Solution Theta coefficients in rhsT
00364       solvechol(Qx_hat,rhsT_INC);   // Solution Theta_inc coefficients in rhsT_INC
00365       invertchol(Qx_hat);           // Covariance matrix of unknowns
00366 
00367       ;// === Copy estimated coefficients to private members ===
00368       BPERP_cf     = rhsBperp;//
00369       BPAR_cf      = rhsBpar;//
00370       THETA_cf     = rhsT;//
00371       THETA_INC_cf = rhsT_INC;//
00372     
00373       // ______Test inverse______
00374       for (int32 i=0; i<Qx_hat.lines(); i++)
00375         for (int32 j=0; j<i; j++)
00376           Qx_hat(j,i) = Qx_hat(i,j);// repair matrix
00377       const real8 maxdev = max(abs(N*Qx_hat-eye(real8(Qx_hat.lines()))));
00378       DEBUG << "baseline: max(abs(N*inv(N)-I)) = " << maxdev;
00379       DEBUG.print();
00380       /*
00381       matrix<real8> unity = N * Qx_hat;
00382       real8 maxdev = abs(unity(0,0)-1);
00383       for (int32 i=1; i<unity.lines(); i++)
00384         {
00385         if (abs(unity(i,i)-1) > maxdev)
00386           maxdev = abs(unity(i,i)-1);
00387         for (int32 j=0; j<i-1; j++)
00388           {
00389           if (abs(unity(i,j)) > maxdev)
00390             maxdev = abs(unity(i,j));
00391           if (abs(unity(j,i)) > maxdev)
00392             maxdev = abs(unity(j,i));
00393           }
00394         }
00395       */
00396       if (maxdev > .01)
00397         {
00398         WARNING << "baseline: max. deviation N*inv(N) from unity = "
00399            << maxdev << ". This is larger than .01: do not use this.";
00400         WARNING.print();// no error, no problem
00401         }
00402       else if (maxdev > .001)
00403         {
00404         WARNING << "baseline: max. deviation N*inv(N) from unity = "
00405            << maxdev << ". This is between 0.01 and 0.001 (do not use it)";
00406         WARNING.print();
00407         }
00408       // ______Some other stuff, scale is ok______
00409       //matrix<real8> Qy_hat = AMATRIX * (matxmatT(Qx_hat,AMATRIX));
00410       matrix<real8> y_hatBperp  = AMATRIX * rhsBperp;
00411       matrix<real8> e_hatBperp  = BPERP - y_hatBperp;
00412       //matrix<real8> Qe_hat  = Qy - Qy_hat;
00413       //matrix<real8> y_hatT  = AMATRIX * rhsT;
00414       //matrix<real8> e_hatT  = THETA - y_hatT;
00415     
00416     
00417       // ______ Output solution and give max error ______
00418       // --- B(l,p,h) = a000 + 
00419       //                a100*l   + a010*p   + a001*h   +
00420       //                a110*l*p + a101*l*h + a011*p*h +
00421       //                a200*l^2 + a020*p^2 + a002*h^2
00422       DEBUG.print("--------------------");
00423       DEBUG.print("Result of modeling: Bperp(l,p) = a000 + a100*l + a010*p + a001*h + ");
00424       DEBUG.print(" a110*l*p + a101*l*h + a011*p*h + a200*l^2 + a020*p^2 + a002*h^2");
00425       DEBUG.print("l,p,h in unnormalized, absolute, coordinates (1:N).");
00426       DEBUG << "Bperp_a000 = " << rhsBperp(0,0); DEBUG.print();
00427       DEBUG << "Bperp_a100 = " << rhsBperp(1,0); DEBUG.print();
00428       DEBUG << "Bperp_a010 = " << rhsBperp(2,0); DEBUG.print();
00429       DEBUG << "Bperp_a001 = " << rhsBperp(3,0); DEBUG.print();
00430       DEBUG << "Bperp_a110 = " << rhsBperp(4,0); DEBUG.print();
00431       DEBUG << "Bperp_a101 = " << rhsBperp(5,0); DEBUG.print();
00432       DEBUG << "Bperp_a011 = " << rhsBperp(6,0); DEBUG.print();
00433       DEBUG << "Bperp_a200 = " << rhsBperp(7,0); DEBUG.print();
00434       DEBUG << "Bperp_a020 = " << rhsBperp(8,0); DEBUG.print();
00435       DEBUG << "Bperp_a002 = " << rhsBperp(9,0); DEBUG.print();
00436       real8 maxerr = max(abs(e_hatBperp));
00437       if (maxerr > 2.00)// 
00438         {
00439         WARNING << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
00440         WARNING.print();
00441         }
00442       else
00443         {
00444         INFO << "Max. error bperp modeling at 3D datapoints: " << maxerr << "m";
00445         INFO.print();
00446         }
00447       DEBUG.print();
00448       DEBUG.print("--------------------");
00449       DEBUG.print("Range: r(p) = r0 + dr*p");
00450       DEBUG.print("l and p in unnormalized, absolute, coordinates (1:N).");
00451       const real8 range1    = master.pix2range(1.0);
00452       const real8 range5000 = master.pix2range(5000.0);
00453       const real8 drange    = (range5000-range1)/5000;
00454       DEBUG << "range = " << range1-drange << " + " << drange << "*p";
00455       DEBUG.print();
00456       // ====== Tidy up ======
00457       }; // END model_parameters()
00458 
00459     // ___ Return RANGE to user ___
00460     inline real8 get_range(const real8 pixel) const
00461       {
00462       return nearrange + drange_dp*pixel;
00463       };// END get_bperp()
00464 
00465     // === modeled quantities ===
00466     // --- B(l,p,h) = a000 + 
00467     //                a100*l   + a010*p   + a001*h   +
00468     //                a110*l*p + a101*l*h + a011*p*h +
00469     //                a200*l^2 + a020*p^2 + a002*h^2
00470     // ___ Return BPERP to user ___
00471     inline real8 get_bperp(const real8 line, const real8 pixel, const real8 height=0.0) const
00472       {
00473       return polyval(BPERP_cf, line, pixel, height);
00474       };
00475     // ___ Return BPAR to user ___
00476     inline real8 get_bpar(const real8 line, const real8 pixel, const real8 height=0.0) const
00477       {
00478       return polyval(BPAR_cf, line, pixel, height);
00479       };
00480     // ___ Return THETA to user ___
00481     inline real8 get_theta(const real8 line, const real8 pixel, const real8 height=0.0) const
00482       {
00483       return polyval(THETA_cf, line, pixel, height);
00484       };
00485     // ___ Return THETA_INC to user ___
00486     inline real8 get_theta_inc(const real8 line, const real8 pixel, const real8 height=0.0) const
00487       {
00488       return polyval(THETA_INC_cf, line, pixel, height);
00489       };
00490 
00491     // === derived quantities ===
00492     // ___ Return B to user ___
00493     inline real8 get_b(const real8 line, const real8 pixel, const real8 height=0.0) const
00494       {
00495       return sqrt(sqr(get_bpar(line,pixel,height))+sqr(get_bperp(line,pixel,height)));
00496       };// END get_b()
00497 
00498     // ___ Return alpha baseline orientation to user ___
00499     inline real8 get_alpha(const real8 line, const real8 pixel, const real8 height=0.0) const
00500       {
00501       const real8 Bperp     = get_bperp(line,pixel,height);
00502       const real8 Bpar      = get_bpar(line,pixel,height);
00503       const real8 B         = get_b(line,pixel,height);
00504       const real8 theta     = get_theta(line,pixel,height);
00505       const real8 alpha     = (Bpar==0 && Bperp==0) ?
00506         NaN :
00507         theta-atan2(Bpar,Bperp);            // sign ok atan2
00508       return alpha;// sign ok
00509       };// END get_bhor()
00510 
00511     // ___ Return Bh to user ___
00512     inline real8 get_bhor(const real8 line, const real8 pixel, const real8 height=0.0) const
00513       {
00514       const real8 Bperp     = get_bperp(line,pixel,height);
00515       const real8 Bpar      = get_bpar(line,pixel,height);
00516       const real8 B         = get_b(line,pixel,height);
00517       const real8 theta     = get_theta(line,pixel,height);
00518       const real8 alpha     = get_alpha(line,pixel,height);
00519       return B*cos(alpha);// sign ok
00520       };// END get_bhor()
00521 
00522     // ___ Return Bv to user ___
00523     inline real8 get_bvert(const real8 line, const real8 pixel, const real8 height=0.0) const
00524       {
00525       const real8 Bperp     = get_bperp(line,pixel,height);
00526       const real8 Bpar      = get_bpar(line,pixel,height);
00527       const real8 B         = get_b(line,pixel,height);
00528       const real8 theta     = get_theta(line,pixel,height);
00529       const real8 alpha     = get_alpha(line,pixel,height);
00530       return B*sin(alpha);// sign ok
00531       };// END get_bvert()
00532 
00533     // ___ Return Height ambiguity to user ___
00534     inline real8 get_hamb(const real8 line, const real8 pixel, const real8 height=0.0) const
00535       {
00536       //const real8 theta     =  get_theta(line,pixel,height);
00537       const real8 theta_inc =  get_theta_inc(line,pixel,height);
00538       const real8 Bperp     =  get_bperp(line,pixel,height);
00539       const real8 range_MP  =  get_range(pixel);// >
00540       const real8 h_amb     = (Bperp==0) ?
00541          Inf : // inf
00542         -master_wavelength*range_MP*sin(theta_inc)/(2.0*Bperp);// this is wrt local
00543         //-master_wavelength*range_MP*sin(theta)/(2.0*Bperp);// this is wrt
00544       return h_amb;
00545       };// END get_hamb()
00546 
00547     // ___ Return orbit convergence to user ___
00548     inline real8 get_orb_conv(const real8 line, const real8 pixel, const real8 height=0.0) const
00549       {
00550       // do not use l,p..
00551       return orbit_convergence;
00552       };// END get_orb_conv()
00553 
00554 
00555     // --- Dump overview of all ---
00556     void dump(const real8 line, const real8 pixel, const real8 height=0.0)
00557       {
00558       if (initialized==false)
00559         {
00560         DEBUG.print("Exiting dumpbaseline, not initialized.");
00561         return;
00562         }
00563         // ______ Modeled quantities ______
00564         const real8 Bperp     = get_bperp(line,pixel,height);
00565         const real8 Bpar      = get_bpar(line,pixel,height);
00566         const real8 theta     = get_theta(line,pixel,height);
00567         const real8 theta_inc = get_theta_inc(line,pixel,height);
00568         // ______ Derived quantities ______
00569         const real8 B         = get_b(line,pixel,height);
00570         const real8 alpha     = get_alpha(line,pixel,height);
00571         const real8 Bh        = get_bhor(line,pixel,height);
00572         const real8 Bv        = get_bvert(line,pixel,height);
00573         const real8 h_amb     = get_hamb(line,pixel,height);
00574         // ______ Height ambiguity: [h] = -lambda/4pi * (r1sin(theta)/Bperp) * phi==2pi ______
00575         // ====== Write output to screen as INFO ======
00576         INFO << "The baseline parameters for (l,p,h) = " 
00577              << line << ", " << pixel << ", " << height;
00578         INFO.print();
00579         INFO << "\tBpar, Bperp:      \t" << Bpar << " \t" << Bperp;
00580         INFO.print();
00581         DEBUG << "\tB, alpha (deg):  \t" << B << " \t" << rad2deg(alpha);
00582         DEBUG.print();
00583         DEBUG << "\tBh, Bv:          \t" << Bh << " \t" << Bv;
00584         DEBUG.print();
00585         INFO << "\tHeight ambiguity: \t" << h_amb;
00586         INFO.print();
00587         INFO << "\tLook angle (deg): \t" << rad2deg(theta);
00588         INFO.print();
00589         DEBUG << "\tIncidence angle (deg): \t" << rad2deg(theta_inc);
00590         DEBUG.print();
00591       }// END dump()
00592 
00593   };// eof BASELINE class
00594 
00595 #endif // BK_BASELINE_H
00596 

Generated on Fri Apr 22 15:57:48 2005 for Doris by doxygen 1.3.6