Doris: utilities.cc Source File

00001 /*
00002  * Copyright (c) 1999-2003 Bert Kampes
00003  * Copyright (c) 1999-2003 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/utilities.cc,v $  *
00032  * $Revision: 3.14 $                                            *
00033  * $Date: 2005/04/06 15:29:03 $                                 *
00034  * $Author: kampes $                                            *
00035  *                                                              *
00036  * implementation of utility routines.                          *
00037  * - math utils                                                 * 
00038  * - solving small linear system                                *
00039  * - polynomial evaluation                                      *
00040  ****************************************************************/
00041 
00042 
00043 #include "matrixbk.hh"                  // my matrix class
00044 #include "slcimage.hh"                  // my slc image class
00045 #include "constants.hh"                 // global constants
00046 #include "refsystems.hh"                // wgsa etc ellipsoid
00047 #include "utilities.hh"                 // header file
00048 #include "ioroutines.hh"                // error messages
00049 #include "conversion.hh"                // ?
00050 #include "exceptions.hh"                 // my exceptions class
00051 
00052 #include <strstream>                    // for memory stream
00053 #include <iomanip>                      // setw
00054 #include <cstring>                      // for strcat etc.
00055 #include <cstdlib>                      // system
00056 #include <cmath>                        // sqrt etc.
00057 #include <cstdio>                       // some compilers, remove function
00058 #include <ctime>                        // time functions
00059 char *strptime(const char *s, const char  *format,  struct tm *tm);        
00060 
00061 
00062 
00063 /****************************************************************
00064  *    getorb                                                    *
00065  *                                                              *
00066  * Calls getorb program to generate precise orbit file          *
00067  *  scratchorbit.                                               *
00068  * Datapoints with INTERVAL seconds in between.                 *
00069  * time1-BEFOREFIRST to timeN+BEFOREFIRST                       *
00070  * time1 in format as: 15-AUG-1996 09:49:16.913                 *
00071  *                                                              *
00072  *    Bert Kampes, 11-Dec-1998                                  *
00073  ****************************************************************/
00074 void getorb(
00075         const slcimage &image,
00076         const input_pr_orbits &inputorb,
00077         int16 ID)
00078   {
00079   TRACE_FUNCTION("getorb (BK 11-Dec-1998)")
00080   const int32 INTERVAL    = inputorb.timeinterval;
00081   const int32 BEFOREFIRST = inputorb.timebefore;
00082   char  orbdir[EIGHTY];
00083   char  dummyline[ONE27];
00084   char  startt[13];
00085   char  endt[13];
00086   char  *pstartt = startt;
00087   char  *pendt = endt;
00088   if (ID==MASTERID)
00089     strcpy(orbdir,inputorb.m_orbdir);
00090   else if (ID==SLAVEID)
00091     strcpy(orbdir,inputorb.s_orbdir);
00092   else
00093     {
00094     PRINT_ERROR("panic, impossible.")
00095     throw(unhandled_case_error);
00096     }
00097 
00098   // ______ Convert times to struct tm ______
00099   struct tm tijdstart;
00100   strptime(image.utc1,"%d-%b-%Y %T",&tijdstart);
00101   DEBUG << "value of image.utc1: " << image.utc1;
00102   DEBUG.print();
00103 
00104   // ______ Adjust interval ______
00105   int32 t1seconds = tijdstart.tm_sec + 60*tijdstart.tm_min + 3600*tijdstart.tm_hour;
00106   int32 t2seconds = t1seconds + BEFOREFIRST +
00107     int32(.5+(image.originalwindow.linehi-image.originalwindow.linelo+1)/image.prf);
00108   t1seconds -= BEFOREFIRST;                             // adjust;
00109   if (t1seconds < 0)                                    // check 0:00'00" utc
00110     {
00111     WARNING.print("getorb: couldn't adjust time lower bound, using old value.");
00112     t1seconds += BEFOREFIRST;
00113     }
00114   tijdstart.tm_hour =  t1seconds / 3600;
00115   tijdstart.tm_min  = (t1seconds-3600*tijdstart.tm_hour) / 60;
00116   tijdstart.tm_sec  = (t1seconds-3600*tijdstart.tm_hour) % 60;
00117   DEBUG << "start time: tm struct: hour=" << tijdstart.tm_hour 
00118         << "; min=" << tijdstart.tm_min << "; sec=" << tijdstart.tm_sec;
00119   DEBUG.print();
00120 
00121   if (t2seconds / 3600 > 23)                            // hours + 1 day
00122     {
00123     WARNING.print("getorb: couldn't adjust time upper bound, using old value.");
00124     t2seconds -= BEFOREFIRST;
00125     }
00126   struct tm tijdend = tijdstart;                // copy year, month, day
00127   tijdend.tm_hour =  t2seconds / 3600;
00128   tijdend.tm_min  = (t2seconds-3600*tijdend.tm_hour) / 60;
00129   tijdend.tm_sec  = (t2seconds-3600*tijdend.tm_hour) % 60;
00130   DEBUG << "end time: tm struct: hour=" << tijdend.tm_hour 
00131         << "; min=" << tijdend.tm_min << "; sec=" << tijdend.tm_sec;
00132   DEBUG.print();
00133 
00134   // ______ Check timedifference ______
00135   //if (tijdend.tm_min - tijdstart.tm_min != 0 &&
00136   //    tijdend.tm_min - tijdstart.tm_min != 1 )
00137   if (tijdend.tm_min-tijdstart.tm_min >  10 ||
00138       tijdend.tm_min-tijdstart.tm_min < -50)
00139     WARNING.print("getorb: > 10 minutes of orbit ephemerides?");
00140 
00141   // ______convert time to strings[12]  YYMMDDHHMMSS(.S)______
00142   strftime(pstartt,13,"%y%m%d%H%M%S",&tijdstart);       // including \0
00143   strftime(pendt,13,"%y%m%d%H%M%S",&tijdend);           // including \0
00144 
00145 
00146   // ====== Get precise orbits ======
00147   // ______ Call getorb (unix) ______
00148   // ______Make string starttime YYMMDDHHMMSS(.S)______
00149   char strgetorb[ONE27];                                // unix system call
00150   ostrstream omemgetorb(strgetorb,ONE27);               // combine types in string
00151   omemgetorb.seekp(0);
00152   omemgetorb << "getorb " << startt << "," << endt << "," << INTERVAL << " "
00153              << orbdir << ">scratchorbit" << ends;
00154   // BK 4 july 2000: assume ODR files are already untarred
00155   //omemgetorb << "getorb " << startt << "," << endt << "," << INTERVAL << " ";
00156   //if (removeuntarred == false)        // .ODR file is in archive directory
00157   //  omemgetorb << orbdir << ">scratchorbit" << ends;
00158   //else                                // file is untarred to current dir
00159   //  omemgetorb << ". > scratchorbit" << ends;
00160 
00161 
00162   int32 status = 3*NaN;
00163   for (register int32 i=0;i<3;i++)                      // try a few times.
00164     {
00165     DEBUG.print(strgetorb);
00166     status=(system(strgetorb));                         // run getorb
00167     DEBUG << "status of getorb: " << status;
00168     DEBUG.print();
00169     if (status == 0)
00170       {
00171       // ______Check really ok______
00172       DEBUG.print("Checking getorb result for errors.");
00173       //ifstream scratchorbit("scratchorbit",ios::in | ios::nocreate); // by getorb
00174       ifstream scratchorbit("scratchorbit",ios::in); // by getorb
00175       bk_assert(scratchorbit,"getorb: scratchorbit",__FILE__,__LINE__);
00176 
00177       int32 err=0;
00178       real8 dsec;
00179       bool getorbok=true;
00180       while(scratchorbit)
00181         {
00182         scratchorbit >> dsec >> err;
00183         DEBUG << "time: " << dsec << " err: " << err;
00184         DEBUG.print();
00185         if (err)
00186           {
00187           getorbok=false;
00188           WARNING.print("seems a problem with getorb, can you repair?");
00189           getanswer();
00190           break; // file?
00191           }
00192         scratchorbit.getline(dummyline,ONE27,'\n');
00193         }
00194       scratchorbit.close();
00195       if (getorbok) 
00196         {
00197         PROGRESS.print("getorb: program finished ok.");
00198         break;  // break try few times
00199         }
00200       }
00201     else if (status == -1 || status==-256)
00202       {
00203       WARNING.print("getorb: exit(-1): time within ODR, but outside precise part.");
00204       break;
00205       }
00206     else if (status ==  1 || status==256)
00207       {
00208       WARNING.print("getorb: exit(1): error opening ODR file.\n\tPerhaps file is tarred?");
00209       getanswer();
00210       }
00211     else if (status ==  2 || status==512)
00212       {
00213       PRINT_ERROR("code 905: getorb: exit(2): no ODR for requested epoch.")
00214       throw(some_error);
00215       }
00216     else if (status ==  3 || status==768)
00217       {
00218       PRINT_ERROR("code 905: getorb: exit(3): too many records.")
00219       throw(some_error);
00220       }
00221     else 
00222       {
00223       ERROR << "code 905: utilities.c: getorb: unknown exit code: " << status;
00224       PRINT_ERROR(ERROR.get_str())
00225       throw(some_error);
00226       }
00227     }
00228 
00229 // ______Tidy up______
00230 // BK 4 july 2000: commented to remove untar and getodr script
00231 // assume that orbit dir contains ODR files.
00232 //  if (removeuntarred)
00233 //    {
00234 //    if (remove(ODRfile))
00235 //      WARNING.print("code 101: could not remove untarred file.");
00236 //    DEBUG.print("untarred file removed.");
00237 //    }
00238   } // END getorb
00239 
00240 
00241 /****************************************************************
00242  *    convertgetorbout                                          *
00243  *                                                              *
00244  * converts output of getorb (in "scratchorbit") to format:     *
00245  *  secsofday X Y Z                                             *
00246  * input:                                                       *
00247  *  - scratchorbit                                              *
00248  * output:                                                      *
00249  *  - scratchdatapoints, this file                              *
00250  *    might already exist due to readleader                     *
00251  *  - file is removed                                           *
00252  *                                                              *
00253  *    Bert Kampes, 11-Dec-1998                                  *
00254  ****************************************************************/
00255 void convertgetorbout(
00256         int16 FILEID,
00257         const char* file)                               // default=scratchorbit
00258   {
00259   TRACE_FUNCTION("convertgetorbout (BK 11-Dec-1998)")
00260   if (FILEID != MASTERID && FILEID!=SLAVEID) 
00261     {
00262     PRINT_ERROR("fileid not ok.")
00263     throw(some_error);
00264     }
00265   char                  dummyline[ONE27];
00266 
00267   // ______ Open files ______
00268   //ifstream infile(file, ios::in | ios::nocreate);
00269   ifstream infile(file, ios::in);
00270   bk_assert(infile,file,__FILE__,__LINE__);
00271 
00272   ofstream outfile("scratchdatapoints", ios::out | ios::trunc); // do replace
00273   bk_assert(outfile,"convertgetorbout: scratchdatapoints",__FILE__,__LINE__);
00274 
00275   // ______Read line: sec85 err fie lambda h x y z______
00276   int32 err;
00277   real8 sec;
00278   real8 fie;
00279   real8 lam;
00280   real8 hei;
00281   char x[25];                                           // to be sure no round off
00282   char y[25];                                           // to be sure no round off
00283   char z[25];                                           // to be sure no round off
00284 
00285 // ______ Get number of points ______
00286   int32 numdatapoints = -1;
00287   int32 numerr        =  0;
00288   DEBUG.print("Counting number of lines in tmp file with orbit; checking err parameter");
00289   while(infile)
00290     {
00291     infile >> sec >> err >> fie >> lam >> hei >> x >> y >> z;
00292     infile.getline(dummyline,ONE27,'\n');
00293     DEBUG  << "scratchfile: sec=" << sec << " err=" << err 
00294           << " x=" << x << " y=" << y << " z=" << z;
00295     DEBUG.print();
00296     if (err!=0) numerr++;
00297     if (err==0) numdatapoints++;
00298     }
00299   infile.clear();// g++v3.2 has a problem if you don't to read it again
00300   //infile.close();// seems problem with g++v3.2 if stream is closed and opened?
00301   INFO << "Number of orbit ephemerides found: " << numdatapoints << ends;
00302   INFO.print();
00303   INFO << "Number of erroneous orbit ephemerides found: " << numerr << ends;
00304   if (numerr==0) INFO.print();
00305   if (numerr>0)  WARNING.print(INFO.get_str());
00306   INFO.reset();// rewind
00307 
00308   // ___ Write results in parameter file via 2nd tmpfile ___
00309   outfile << "\n\n*******************************************************************"
00310           << "\n*_Start_";
00311           //<< "\n*_Start_precise_datapoints"
00312   if (FILEID==MASTERID)
00313     outfile << processcontrol[pr_m_porbits];
00314   if (FILEID==SLAVEID)
00315     outfile << processcontrol[pr_s_porbits];
00316   outfile << "\n*******************************************************************"
00317           << "\n\tt(s)\tX(m)\tY(m)\tZ(m)"
00318           << "\nNUMBER_OF_DATAPOINTS: \t\t\t"
00319           <<  numdatapoints
00320           << endl;
00321 
00322 // ______Set outfile in correct digits______
00323   int32 secint;                                         // integer part
00324   real8 secfrac;                                        // fractional part
00325   outfile.setf(ios::fixed);
00326   outfile.setf(ios::showpoint);
00327   //infile.open(file);
00328   infile.seekg(0, ios::beg);// rewind file to read again
00329   DEBUG.print("Rewinding file and reading orbit ephemerides again");
00330   for (register int32 i=0;i<numdatapoints+numerr;i++)// total number of lines
00331     {
00332     infile >> sec >> err >> fie >> lam >> hei >> x >> y >> z;
00333     infile.getline(dummyline,ONE27,'\n');               // go to next line
00334     DEBUG  << "scratchfile2: sec=" << sec << " err=" << err 
00335           << " x=" << x << " y=" << y << " z=" << z;
00336     DEBUG.print();
00337 
00338     // ______Convert time sec85 to sec of day______
00339     secint  = int32(sec)%86400;                         // 60*60*24 integer part
00340     secfrac = sec - int32(sec);                         // fractional part 
00341     if (err!=0)
00342       WARNING.print("err!=0: something went wrong in getorb? (skipping point)");
00343     if (err==0)
00344       outfile << secint+secfrac << "\t" << x << "\t" << y << "\t" << z << endl;
00345     }
00346 
00347   outfile << "\n*******************************************************************"
00348           //<< "\n* End_precise_datapoints:_NORMAL"
00349           << "\n* End_";
00350   if (FILEID==MASTERID)
00351     outfile << processcontrol[pr_m_porbits];
00352   if (FILEID==SLAVEID)
00353     outfile << processcontrol[pr_s_porbits];
00354   outfile << "_NORMAL"
00355           << "\n*******************************************************************\n";
00356 
00357 // ______Tidy up______
00358   infile.close();
00359   outfile.close();
00360   if (remove(file))                                     // remove file
00361     WARNING.print("code 101: could not remove file.");
00362   } // END convertgetorbout
00363 
00364 
00365 
00366 /****************************************************************
00367  *    solve33                                                   *
00368  *                                                              *
00369  * Solves setof 3 equations by straightforward (no pivotting) LU*
00370  * y=Ax (unknown x)                                             *
00371  *                                                              *
00372  * input:                                                       *
00373  *  - matrix<real8> righthandside 3x1 (y)                       *
00374  *  - matrix<real8> partials 3x3 (A)                            *
00375  *                                                              *
00376  * output:                                                      *
00377  *  - matrix<real8> result 3x1 unknown                          *
00378  *                                                              *
00379  *    Bert Kampes, 22-Dec-1998                                  *
00380  ****************************************************************/
00381 void solve33(
00382         matrix<real8>       &RESULT,
00383         const matrix<real8> &rhs,
00384         const matrix<real8> &A)
00385   {
00386   TRACE_FUNCTION("solve33 (BK 22-Dec-1998)")
00387 #ifdef __DEBUG
00388   if (RESULT.lines() != 3 || RESULT.pixels() != 1)
00389     {
00390     PRINT_ERROR("solve33: input: size RES not 3x1.")
00391     throw(input_error);
00392     }
00393   if (A.lines() != 3 || A.pixels() != 3)
00394     {
00395     PRINT_ERROR("solve33: input: size of A not 33.")
00396     throw(input_error);
00397     }
00398   if (rhs.lines() != 3 || rhs.pixels() != 1)
00399     {
00400     PRINT_ERROR("solve33: input: size rhs not 3x1.")
00401     throw(input_error);
00402     }
00403 #endif
00404    
00405 ///  real8 L10, L20, L21;                       // used lower matrix elements
00406 ///  real8 U11, U12, U22;                       // used upper matrix elements
00407 ///  real8 b0,  b1,  b2;                        // used Ux=b
00408   const real8 L10 = A(1,0)/A(0,0);
00409   const real8 L20 = A(2,0)/A(0,0);
00410   const real8 U11 = A(1,1)-L10*A(0,1);
00411   const real8 L21 = (A(2,1)-(A(0,1)*L20))/U11;
00412   const real8 U12 = A(1,2)-L10*A(0,2);
00413   const real8 U22 = A(2,2)-L20*A(0,2)-L21*U12;
00414    
00415 // ______Solution: forward substitution______
00416   const real8 b0  = rhs(0,0);
00417   const real8 b1  = rhs(1,0)-b0*L10;
00418   const real8 b2  = rhs(2,0)-b0*L20-b1*L21;
00419    
00420 // ______Solution: backwards substitution______
00421   RESULT(2,0) = b2/U22;
00422   RESULT(1,0) = (b1-U12*RESULT(2,0))/U11;
00423   RESULT(0,0) = (b0-A(0,1)*RESULT(1,0)-A(0,2)*RESULT(2,0))/A(0,0);
00424   } // END solve33
00425 
00426 
00427 
00428 /****************************************************************
00429  *    solve22                                                   *
00430  *                                                              *
00431  * Solves setof 2 equations by straightforward substitution     *
00432  * y=Ax (unknown x)                                             *
00433  *                                                              *
00434  * input:                                                       *
00435  *  - matrix<real8> righthandside 2x1 (y)                       *
00436  *  - matrix<real8> partials 2x2 (A)                            *
00437  * output:                                                      *
00438  *  - matrix<real8> result 2x1 unknown dx,dy,dz                 *
00439  *                                                              *
00440  * Future: LaTeX software doc.                                  *
00441  * Future: this has been tested                                 *
00442  *                                                              *
00443  *    Bert Kampes, 22-Dec-1998                                  *
00444  ****************************************************************/
00445 matrix<real8> solve22(
00446         const matrix<real8> &y,
00447         const matrix<real8> &A)
00448   {
00449   TRACE_FUNCTION("solve22 (BK 22-Dec-1998)")
00450 #ifdef __DEBUG
00451   // ______ Check input ______
00452   if (A.lines() != 2 || A.pixels() != 2)
00453     {
00454     PRINT_ERROR("error...")
00455     throw(input_error);
00456     }
00457   if (y.lines() != 2 || y.pixels() != 1)
00458     {
00459     PRINT_ERROR("error...")
00460     throw(input_error);
00461     }
00462 #endif
00463 
00464   matrix<real8> Result(2,1);
00465   Result(1,0)= (y(0,0) - ((A(0,0)/A(1,0))*y(1,0))) / (A(0,1) - ((A(0,0)*A(1,1))/A(1,0)));
00466   Result(0,0)= (y(0,0) - A(0,1)*Result(1,0)) / A(0,0);
00467   return Result;
00468   } // END solve22
00469 
00470 
00471 /****************************************************************
00472  *    nextpow2                                                  *
00473  *                                                              *
00474  * Returns 32 for 31.3 etc.                                     *
00475  *                                                              *
00476  *    Bert Kampes, 03-Feb-1999                                  *
00477  * AFter bug report from Bruno who discover bug.  b-1 instead   *
00478  * of b.  Did not                                               *
00479  * have effect on Doris since this function was effectively not *
00480  * used.                                                        *
00481  #%// BK 16-Apr-2002                                            *
00482  ****************************************************************/
00483 uint nextpow2(
00484         real8 w)        
00485   {
00486   TRACE_FUNCTION("nextpow2 (BK 03-Feb-1999)")
00487   int32 b    = 0; 
00488   int32 *pnt = &b;
00489   real8 f    = frexp(w,pnt);
00490   if (f==.5) 
00491     return uint(w);
00492   // return (uint(pow(real8(2.),b-1)));
00493   return (uint(pow(real8(2.),b)));
00494   } // END nextpow2
00495 
00496 
00497 /****************************************************************
00498  *    polyval                                                   *
00499  *                                                              *
00500  * evaluates polynomial at (x,y)                                *
00501  * input:                                                       *
00502  *  - x,y: point                                                *
00503  *  - polynomial coefficients                                   *
00504  * (a00 | 01 10 | 20 11 02 | 30 21 12 03 | ...)                 *
00505  *                                                              *
00506  * output:                                                      *
00507  *  - real8 functional value at (x,y).                          *
00508  *                                                              *
00509  *    Bert Kampes, 15-Mar-1999                                  *
00510  ****************************************************************/
00511 real8 polyval(
00512         real8 x,
00513         real8 y,
00514         const matrix<real8> &coeff)
00515   {
00516   const int32 degreee = degree(coeff.size());
00517 #ifdef __DEBUG
00518   TRACE_FUNCTION("polyval (BK 15-Mar-1999)")
00519   if (coeff.size() != coeff.lines())
00520     {
00521     PRINT_ERROR("polyval: require standing data vector.")
00522     throw(input_error);
00523     }
00524   if (degreee<0 || degreee > 1000)
00525     {
00526     PRINT_ERROR("polyval: degreee < 0")
00527     throw(input_error);
00528     }
00529 #endif
00530 
00531 // ______Check default arguments______
00532 //  if (degreee==-1)                    // default applies
00533 //    const int32 degreee = degree(coeff.size());
00534 
00535 // ______Speed up______
00536 // ______Evaluate polynomial______
00537 
00538   real8 sum = coeff(0,0);
00539 
00540   if (degreee == 1)
00541     {
00542     sum += (  coeff(1,0) * x
00543             + coeff(2,0) * y );
00544     }
00545 
00546   else if (degreee == 2)
00547     {
00548     sum += (  coeff(1,0) * x
00549             + coeff(2,0) * y
00550             + coeff(3,0) * sqr(x)
00551             + coeff(4,0) * x*y
00552             + coeff(5,0) * sqr(y) );
00553     }
00554 
00555   else if (degreee == 3)
00556     {
00557     register real8 xx = sqr(x);
00558     register real8 yy = sqr(y);
00559     sum += (  coeff(1,0) * x
00560           + coeff(2,0) * y
00561           + coeff(3,0) * xx
00562           + coeff(4,0) * x*y
00563           + coeff(5,0) * yy
00564           + coeff(6,0) * xx*x
00565           + coeff(7,0) * xx*y
00566           + coeff(8,0) * x*yy
00567           + coeff(9,0) * yy*y );
00568     }
00569 
00570   else if (degreee == 4)
00571     {
00572     register real8 xx  = sqr(x);
00573     register real8 xxx = xx*x;
00574     register real8 yy  = sqr(y);
00575     register real8 yyy = yy*y;
00576     sum += (  coeff(1,0) * x
00577             + coeff(2,0) * y
00578             + coeff(3,0) * xx
00579             + coeff(4,0) * x*y
00580             + coeff(5,0) * yy
00581             + coeff(6,0) * xxx
00582             + coeff(7,0) * xx*y
00583             + coeff(8,0) * x*yy
00584             + coeff(9,0) * yyy
00585             + coeff(10,0)* xx*xx
00586             + coeff(11,0)* xxx*y
00587             + coeff(12,0)* xx*yy
00588             + coeff(13,0)* x*yyy
00589             + coeff(14,0)* yy*yy );
00590     }
00591 
00592   else if (degreee == 5)
00593     {
00594     register real8 xx   = sqr(x);
00595     register real8 xxx  = xx*x;
00596     register real8 xxxx = xxx*x;
00597     register real8 yy   = sqr(y);
00598     register real8 yyy  = yy*y;
00599     register real8 yyyy = yyy*y;
00600     sum += (  coeff(1,0) * x
00601             + coeff(2,0) * y
00602             + coeff(3,0) * xx
00603             + coeff(4,0) * x*y
00604             + coeff(5,0) * yy
00605             + coeff(6,0) * xxx
00606             + coeff(7,0) * xx*y
00607             + coeff(8,0) * x*yy
00608             + coeff(9,0) * yyy
00609             + coeff(10,0)* xxxx
00610             + coeff(11,0)* xxx*y
00611             + coeff(12,0)* xx*yy
00612             + coeff(13,0)* x*yyy
00613             + coeff(14,0)* yyyy
00614             + coeff(15,0)* xxxx*x
00615             + coeff(16,0)* xxxx*y
00616             + coeff(17,0)* xxx*yy
00617             + coeff(18,0)* xx*yyy
00618             + coeff(19,0)* x*yyyy
00619             + coeff(20,0)* yyyy*y );
00620     }
00621 
00622   else if (degreee != 0)                // degreee > 5
00623     {
00624     sum = 0.;
00625     register int32 coeffindex = 0;
00626     for (register int32 l=0; l<=degreee; l++)
00627       {
00628       for (register int32 k=0; k<=l ;k++)
00629         {
00630         sum += coeff(coeffindex,0) * pow(x,real8(l-k)) * pow(y,real8(k));
00631         coeffindex++;
00632         }
00633       }
00634     }
00635 
00636   return sum;  
00637   } // END polyval
00638 
00639 
00640 /****************************************************************
00641  *    polyval                                                   *
00642  *                                                              *
00643  * evaluates polynomial at (x,y)                                *
00644  * input:                                                       *
00645  *  - x,y: point                                                *
00646  *  - polynomial coefficients                                   *
00647  * (a00 | 01 10 | 20 11 02 | 30 21 12 03 | ...)                 *
00648  *                                                              *
00649  * output:                                                      *
00650  *  - real8 functional value at (x,y).                          *
00651  *                                                              *
00652  *    Bert Kampes, 15-Mar-1999                                  *
00653  ****************************************************************/
00654 real8 polyval(
00655         real8 x,
00656         real8 y,
00657         const matrix<real8> &coeff,
00658         int32 degreee)
00659   {
00660 #ifdef __DEBUG
00661   TRACE_FUNCTION("polyval (BK 15-Mar-1999)")
00662   if (coeff.size() != coeff.lines())
00663     {
00664     PRINT_ERROR("polyval: require standing data vector.")
00665     throw(input_error);
00666     }
00667   if (degreee<0 || degreee > 1000)
00668     {
00669     PRINT_ERROR("polyval: degreee < 0")
00670     throw(input_error);
00671     }
00672 #endif
00673 
00674 // ______Check default arguments______
00675 //  if (degreee==-1)                    // default applies
00676 //    degreee = degree(coeff.size());
00677 
00678 // ______Speed up______
00679 // ______Evaluate polynomial______
00680 
00681   real8 sum = coeff(0,0);
00682 
00683   if (degreee == 1)
00684     {
00685     sum += (  coeff(1,0) * x
00686             + coeff(2,0) * y );
00687     }
00688 
00689   else if (degreee == 2)
00690     {
00691     sum += (  coeff(1,0) * x
00692             + coeff(2,0) * y
00693             + coeff(3,0) * sqr(x)
00694             + coeff(4,0) * x*y
00695             + coeff(5,0) * sqr(y) );
00696     }
00697 
00698   else if (degreee == 3)
00699     {
00700     register real8 xx = sqr(x);
00701     register real8 yy = sqr(y);
00702     sum += (  coeff(1,0) * x
00703           + coeff(2,0) * y
00704           + coeff(3,0) * xx
00705           + coeff(4,0) * x*y
00706           + coeff(5,0) * yy
00707           + coeff(6,0) * xx*x
00708           + coeff(7,0) * xx*y
00709           + coeff(8,0) * x*yy
00710           + coeff(9,0) * yy*y );
00711     }
00712 
00713   else if (degreee == 4)
00714     {
00715     register real8 xx  = sqr(x);
00716     register real8 xxx = xx*x;
00717     register real8 yy  = sqr(y);
00718     register real8 yyy = yy*y;
00719     sum += (  coeff(1,0) * x
00720             + coeff(2,0) * y
00721             + coeff(3,0) * xx
00722             + coeff(4,0) * x*y
00723             + coeff(5,0) * yy
00724             + coeff(6,0) * xxx
00725             + coeff(7,0) * xx*y
00726             + coeff(8,0) * x*yy
00727             + coeff(9,0) * yyy
00728             + coeff(10,0)* xx*xx
00729             + coeff(11,0)* xxx*y
00730             + coeff(12,0)* xx*yy
00731             + coeff(13,0)* x*yyy
00732             + coeff(14,0)* yy*yy );
00733     }
00734 
00735   else if (degreee == 5)
00736     {
00737     register real8 xx   = sqr(x);
00738     register real8 xxx  = xx*x;
00739     register real8 xxxx = xxx*x;
00740     register real8 yy   = sqr(y);
00741     register real8 yyy  = yy*y;
00742     register real8 yyyy = yyy*y;
00743     sum += (  coeff(1,0) * x
00744             + coeff(2,0) * y
00745             + coeff(3,0) * xx
00746             + coeff(4,0) * x*y
00747             + coeff(5,0) * yy
00748             + coeff(6,0) * xxx
00749             + coeff(7,0) * xx*y
00750             + coeff(8,0) * x*yy
00751             + coeff(9,0) * yyy
00752             + coeff(10,0)* xxxx
00753             + coeff(11,0)* xxx*y
00754             + coeff(12,0)* xx*yy
00755             + coeff(13,0)* x*yyy
00756             + coeff(14,0)* yyyy
00757             + coeff(15,0)* xxxx*x
00758             + coeff(16,0)* xxxx*y
00759             + coeff(17,0)* xxx*yy
00760             + coeff(18,0)* xx*yyy
00761             + coeff(19,0)* x*yyyy
00762             + coeff(20,0)* yyyy*y );
00763     }
00764 
00765   else if (degreee != 0)                // degreee > 5
00766     {
00767     sum = 0.;
00768     register int32 coeffindex = 0;
00769     for (register int32 l=0; l<=degreee; l++)
00770       {
00771       for (register int32 k=0; k<=l ;k++)
00772         {
00773         sum += coeff(coeffindex,0) * pow(x,real8(l-k)) * pow(y,real8(k));
00774         coeffindex++;
00775         }
00776       }
00777     }
00778   return sum;  
00779   } // END polyval
00780 
00781 
00782 /****************************************************************
00783  *    polyval                                                   *
00784  *                                                              *
00785  * Evaluate 2d polynomial at regular grid.                              *
00786  * input:                                                       *
00787  *  - x,y: (standing) data axis                                 *
00788  *  - polynomial coefficients                                   *
00789  * (a00 | 01 10 | 20 11 02 | 30 21 12 03 | ...)                 *
00790  *                                                              *
00791  * output:                                                      *
00792  * ( - file ascii "TEST" with grid.)                            *
00793  *  - matrix<real4> with evaluated polynomial                   *
00794  *                                                              *
00795  * note optimixed for more y than x                             *
00796  *    Bert Kampes, 12-Mar-1999                                  *
00797  ****************************************************************/
00798 matrix<real4> polyval(
00799         const matrix<real4> &x,
00800         const matrix<real4> &y,
00801         const matrix<real8> &coeff,
00802         int32 degreee)                  // optional
00803   {
00804   TRACE_FUNCTION("polyval (BK 12-Mar-1999)")
00805 #ifdef __DEBUG
00806   if (x.size() != x.lines())
00807     {
00808     PRINT_ERROR("polyval: require standing x vector.")
00809     throw(input_error);
00810     }
00811   if (y.size() != y.lines())
00812     {
00813     PRINT_ERROR("polyval: require standing y vector.")
00814     throw(input_error);
00815     }
00816   if (coeff.size() != coeff.lines())
00817     {
00818     PRINT_ERROR("polyval: require standing coeff. vector.")
00819     throw(input_error);
00820     }
00821   if (degreee<-1)
00822     {
00823     PRINT_ERROR("polyval: degreee < -1")
00824     throw(input_error);
00825     }
00826   if (x.size()>y.size())
00827     DEBUG.print("x larger than y, while optimized for y larger x");
00828 #endif
00829 
00830 // ______Check default arguments______
00831   if (degreee==-1)                              // default applies
00832     degreee = degree(coeff.size());
00833 
00834 // ______Evaluate polynomial______
00835   matrix<real4> Result(x.size(),y.size());
00836   register int32 i,j;
00837 
00838   register real4 c00;
00839   register real4 c10;
00840   register real4 c01;
00841   register real4 c20;
00842   register real4 c11;
00843   register real4 c02;
00844   register real4 c30;
00845   register real4 c21;
00846   register real4 c12;
00847   register real4 c03;
00848   register real4 c40;
00849   register real4 c31;
00850   register real4 c22;
00851   register real4 c13;
00852   register real4 c04;
00853   register real4 c50;
00854   register real4 c41;
00855   register real4 c32;
00856   register real4 c23;
00857   register real4 c14;
00858   register real4 c05;
00859   switch (degreee)
00860     {
00861     case 0:                                     // constant
00862       Result.setdata(real4(coeff(0,0)));
00863       break;
00864 
00865     case 1:
00866       c00 = coeff(0,0);
00867       c10 = coeff(1,0);
00868       c01 = coeff(2,0);
00869       for (j=0; j<Result.pixels(); j++)
00870         {
00871         real4 c00pc01y1 = c00 + c01*y(j,0);
00872         for (i=0; i<Result.lines(); i++)
00873           {
00874           Result(i,j) =  c00pc01y1 + c10 * x(i,0);
00875           }
00876         }
00877       break;
00878 
00879     case 2:
00880         c00 = coeff(0,0);
00881         c10 = coeff(1,0);
00882         c01 = coeff(2,0);
00883         c20 = coeff(3,0);
00884         c11 = coeff(4,0);
00885         c02 = coeff(5,0);
00886         for (j=0; j<Result.pixels(); j++)
00887           {
00888           real4 y1 = y(j,0);
00889           real4 c00pc01y1 = c00 + c01 * y1;
00890           real4 c02y2 = c02 * sqr(y1);
00891           real4 c11y1 = c11 * y1;
00892           for (i=0; i<Result.lines(); i++)
00893             {
00894             real4 x1 = x(i,0);
00895             Result(i,j) =  c00pc01y1 
00896                          + c10   * x1 
00897                          + c20   * sqr(x1) 
00898                          + c11y1 * x1 
00899                          + c02y2;
00900             }
00901           }
00902         break;
00903 
00904     case 3:
00905         c00 = coeff(0,0);
00906         c10 = coeff(1,0);
00907         c01 = coeff(2,0);
00908         c20 = coeff(3,0);
00909         c11 = coeff(4,0);
00910         c02 = coeff(5,0);
00911         c30 = coeff(6,0);
00912         c21 = coeff(7,0);
00913         c12 = coeff(8,0);
00914         c03 = coeff(9,0);
00915         for (j=0; j<Result.pixels(); j++)
00916           {
00917           real4 y1 = y(j,0);
00918           real4 y2 = sqr(y1);
00919           real4 c00pc01y1 = c00 + c01 * y1;
00920           real4 c02y2 = c02 * y2;
00921           real4 c11y1 = c11 * y1;
00922           real4 c21y1 = c21 * y1;
00923           real4 c12y2 = c12 * y2;
00924           real4 c03y3 = c03 * y1 * y2;
00925           for (i=0; i<Result.lines(); i++)
00926             {
00927             real4 x1 = x(i,0);
00928             real4 x2 = sqr(x1);
00929             Result(i,j) =  c00pc01y1 
00930                          + c10   * x1 
00931                          + c20   * x2 
00932                          + c11y1 * x1 
00933                          + c02y2
00934                          + c30   * x1 * x2
00935                          + c21y1 * x2
00936                          + c12y2 * x1
00937                          + c03y3;
00938             }
00939           }
00940         break;
00941 
00942     case 4:
00943         c00 = coeff(0,0);
00944         c10 = coeff(1,0);
00945         c01 = coeff(2,0);
00946         c20 = coeff(3,0);
00947         c11 = coeff(4,0);
00948         c02 = coeff(5,0);
00949         c30 = coeff(6,0);
00950         c21 = coeff(7,0);
00951         c12 = coeff(8,0);
00952         c03 = coeff(9,0);
00953         c40 = coeff(10,0);
00954         c31 = coeff(11,0);
00955         c22 = coeff(12,0);
00956         c13 = coeff(13,0);
00957         c04 = coeff(14,0);
00958         for (j=0; j<Result.pixels(); j++)
00959           {
00960           real4 y1 = y(j,0);
00961           real4 y2 = sqr(y1);
00962           real4 c00pc01y1 = c00 + c01 * y1;
00963           real4 c02y2 = c02 * y2;
00964           real4 c11y1 = c11 * y1;
00965           real4 c21y1 = c21 * y1;
00966           real4 c12y2 = c12 * y2;
00967           real4 c03y3 = c03 * y1 * y2;
00968           real4 c31y1 = c31 * y1;
00969           real4 c22y2 = c22 * y2;
00970           real4 c13y3 = c13 * y2 * y1;
00971           real4 c04y4 = c04 * y2 * y2;
00972           for (i=0; i<Result.lines(); i++)
00973             {
00974             real4 x1 = x(i,0);
00975             real4 x2 = sqr(x1);
00976             Result(i,j) =  c00pc01y1 
00977                          + c10   * x1 
00978                          + c20   * x2 
00979                          + c11y1 * x1 
00980                          + c02y2
00981                          + c30   * x1 * x2
00982                          + c21y1 * x2
00983                          + c12y2 * x1
00984                          + c03y3
00985                          + c40   * x2 * x2
00986                          + c31y1 * x2 * x1
00987                          + c22y2 * x2
00988                          + c13y3 * x1
00989                          + c04y4;
00990             }
00991           }
00992         break;
00993 
00994     case 5:
00995         c00 = coeff(0,0);
00996         c10 = coeff(1,0);
00997         c01 = coeff(2,0);
00998         c20 = coeff(3,0);
00999         c11 = coeff(4,0);
01000         c02 = coeff(5,0);
01001         c30 = coeff(6,0);
01002         c21 = coeff(7,0);
01003         c12 = coeff(8,0);
01004         c03 = coeff(9,0);
01005         c40 = coeff(10,0);
01006         c31 = coeff(11,0);
01007         c22 = coeff(12,0);
01008         c13 = coeff(13,0);
01009         c04 = coeff(14,0);
01010         c50 = coeff(15,0);
01011         c41 = coeff(16,0);
01012         c32 = coeff(17,0);
01013         c23 = coeff(18,0);
01014         c14 = coeff(19,0);
01015         c05 = coeff(20,0);
01016         for (j=0; j<Result.pixels(); j++)
01017           {
01018           real4 y1 = y(j,0);
01019           real4 y2 = sqr(y1);
01020           real4 y3 = y2*y1;
01021           real4 c00pc01y1 = c00 + c01 * y1;
01022           real4 c02y2 = c02 * y2;
01023           real4 c11y1 = c11 * y1;
01024           real4 c21y1 = c21 * y1;
01025           real4 c12y2 = c12 * y2;
01026           real4 c03y3 = c03 * y3;
01027           real4 c31y1 = c31 * y1;
01028           real4 c22y2 = c22 * y2;
01029           real4 c13y3 = c13 * y3;
01030           real4 c04y4 = c04 * y2 * y2;
01031           real4 c41y1 = c41 * y1;
01032           real4 c32y2 = c32 * y2;
01033           real4 c23y3 = c23 * y3;
01034           real4 c14y4 = c14 * y2 * y2;
01035           real4 c05y5 = c05 * y3 * y2;
01036           for (i=0; i<Result.lines(); i++)
01037             {
01038             real4 x1 = x(i,0);
01039             real4 x2 = sqr(x1);
01040             real4 x3 = x1*x2;
01041             Result(i,j) =  c00pc01y1 
01042                          + c10   * x1 
01043                          + c20   * x2 
01044                          + c11y1 * x1 
01045                          + c02y2
01046                          + c30   * x3
01047                          + c21y1 * x2
01048                          + c12y2 * x1
01049                          + c03y3
01050                          + c40   * x2 * x2
01051                          + c31y1 * x3
01052                          + c22y2 * x2
01053                          + c13y3 * x1
01054                          + c04y4
01055                          + c50   * x3 * x2
01056                          + c41y1 * x2 * x2
01057                          + c32y2 * x3
01058                          + c23y3 * x2
01059                          + c14y4 * x1
01060                          + c05y5;
01061             }
01062           }
01063         break;
01064 
01065 // ______ solve up to 5 efficiently, do rest in loop ______
01066       default:
01067         c00 = coeff(0,0);
01068         c10 = coeff(1,0);
01069         c01 = coeff(2,0);
01070         c20 = coeff(3,0);
01071         c11 = coeff(4,0);
01072         c02 = coeff(5,0);
01073         c30 = coeff(6,0);
01074         c21 = coeff(7,0);
01075         c12 = coeff(8,0);
01076         c03 = coeff(9,0);
01077         c40 = coeff(10,0);
01078         c31 = coeff(11,0);
01079         c22 = coeff(12,0);
01080         c13 = coeff(13,0);
01081         c04 = coeff(14,0);
01082         c50 = coeff(15,0);
01083         c41 = coeff(16,0);
01084         c32 = coeff(17,0);
01085         c23 = coeff(18,0);
01086         c14 = coeff(19,0);
01087         c05 = coeff(20,0);
01088         for (j=0; j<Result.pixels(); j++)
01089           {
01090           real4 y1 = y(j,0);
01091           real4 y2 = sqr(y1);
01092           real4 y3 = y2*y1;
01093           real4 c00pc01y1 = c00 + c01 * y1;
01094           real4 c02y2 = c02 * y2;
01095           real4 c11y1 = c11 * y1;
01096           real4 c21y1 = c21 * y1;
01097           real4 c12y2 = c12 * y2;
01098           real4 c03y3 = c03 * y3;
01099           real4 c31y1 = c31 * y1;
01100           real4 c22y2 = c22 * y2;
01101           real4 c13y3 = c13 * y3;
01102           real4 c04y4 = c04 * y2 * y2;
01103           real4 c41y1 = c41 * y1;
01104           real4 c32y2 = c32 * y2;
01105           real4 c23y3 = c23 * y3;
01106           real4 c14y4 = c14 * y2 * y2;
01107           real4 c05y5 = c05 * y3 * y2;
01108           for (i=0; i<Result.lines(); i++)
01109             {
01110             real4 x1 = x(i,0);
01111             real4 x2 = sqr(x1);
01112             real4 x3 = x1*x2;
01113             Result(i,j) =  c00pc01y1 
01114                          + c10   * x1 
01115                          + c20   * x2 
01116                          + c11y1 * x1 
01117                          + c02y2
01118                          + c30   * x3
01119                          + c21y1 * x2
01120                          + c12y2 * x1
01121                          + c03y3
01122                          + c40   * x2 * x2
01123                          + c31y1 * x3
01124                          + c22y2 * x2
01125                          + c13y3 * x1
01126                          + c04y4
01127                          + c50   * x3 * x2
01128                          + c41y1 * x2 * x2
01129                          + c32y2 * x3
01130                          + c23y3 * x2
01131                          + c14y4 * x1
01132                          + c05y5;
01133             }
01134           }
01135 
01136 // ______ do rest (>5) in loop ______
01137 WARNING.print("this seems to be wrong ?? BK 9 feb.00");
01138 WARNING.print("this seems to be wrong ?? BK 9 feb.00");
01139 WARNING.print("this seems to be wrong ?? BK 9 feb.00");
01140 WARNING.print("this seems to be wrong ?? BK 9 feb.00");
01141 WARNING.print("this seems to be wrong ?? BK 9 feb.00");
01142     const int32 STARTDEGREE = 6;
01143     const int32 STARTCOEFF  = Ncoeffs(STARTDEGREE-1);   // 5-> 21 6->28 7->36 etc.
01144     for (j=0; j<Result.pixels(); j++)
01145       {
01146       real4 yy=y(j,0);
01147       for (i=0; i<Result.lines(); i++)
01148         {
01149         real4 xx=x(i,0);        // ??? this seems to be wrong (BK 9-feb-00)
01150         real4 sum = 0.;
01151         register int32 coeffindex = STARTCOEFF;
01152         for (register int32 l=STARTDEGREE; l<=degreee; l++)
01153           {
01154           for (register int32 k=0; k<=l ;k++)
01155             {
01156             sum += coeff(coeffindex,0) * pow(xx,real4(l-k)) * pow(yy,real4(k));
01157             coeffindex++;
01158             }
01159           }
01160         Result(i,j) += sum;  
01161         }
01162       }
01163     } // switch degreee
01164 
01165   return Result;
01166   } // END polyval
01167 
01168 
01169 /****************************************************************
01170  *    polyval1d                                                 *
01171  *                                                              *
01172  * evaluates polynomial at (x)                                  *
01173  * input:                                                       *
01174  *  - x: point                                                  *
01175  *  - polynomial coefficients in standing vector (a0,a1,a2,..). *
01176  *  - degree of polynomial.                                     *
01177  *                                                              *
01178  * output:                                                      *
01179  *  - real8 functional value y = f(x).                          *
01180  *                                                              *
01181  *    Bert Kampes, 03-Jun-1999                                  *
01182  ****************************************************************/
01183 real8 polyval1d(
01184         real8 x,
01185         const matrix<real8> &coeff)
01186   {
01187 #ifdef __DEBUG
01188   TRACE_FUNCTION("polyval1d (BK 03-Jun-1999)")
01189   if (coeff.size() != coeff.lines())
01190     {
01191     PRINT_ERROR("polyval: require standing data vector.")
01192     throw(input_error);
01193     }
01194 #endif
01195   real8 sum = 0.0;
01196   for (register int32 d=coeff.size()-1;d>=0;--d)
01197     {
01198     sum *= x;
01199     sum += coeff(d,0);
01200     }
01201   return sum;
01202   } // END polyval1d
01203 
01204 
01205 /****************************************************************
01206  *    normalize                                                 *
01207  *                                                              *
01208  * rescale data[min,max] to interval [-2,2]                     *
01209  *                                                              *
01210  * input:                                                       *
01211  *  - data                                                      *
01212  *  - min, max                                                  *
01213  * output:                                                      *
01214  *  - data matrix by reference                                  *
01215  *                                                              *
01216  *    Bert Kampes, 21-Jun-1999                                  *
01217  ****************************************************************/
01218 void normalize(
01219         matrix<real4> &data,
01220         real8 min,
01221         real8 max
01222         )
01223   {
01224   TRACE_FUNCTION("normalize (BK 21-Jun-1999)")
01225   // ______ zero mean, rescale ______
01226   data -= real4(.5*(min+max));                          // [.5a-.5b, -.5a+.5b]
01227   data /= real4(.25*(max-min));                         // [-2 2]
01228   } // END normalize
01229 
01230 
01231 /****************************************************************
01232  *    normalize                                                 *
01233  *                                                              *
01234  * rescale data to interval [-1,1]                              *
01235  * data -> X(data-min)/(max-min) ; data+SS                      *
01236  * X=1--1==EE-SS; SS=-1, EE=1                                   *
01237  *                                                              *
01238  * input:                                                       *
01239  *  - data                                                      *
01240  *  - min, max                                                  *
01241  * output:                                                      *
01242  *  - by reference                                              *
01243  *                                                              *
01244  *    Bert Kampes, 21-Jun-1999                                  *
01245  ****************************************************************/
01246 void normalize(
01247         matrix<real8> &data,
01248         real8 min,
01249         real8 max
01250         )
01251   {
01252   TRACE_FUNCTION("ROUTINE: normalize (BK 21-Jun-1999)")
01253   // ______ zero mean, rescale ______
01254   data -= .5*(min+max);                         // [.5a-.5b, -.5a+.5b]
01255   data /= .25*(max-min);                                // [-2 2]
01256   } // END normalize
01257 
01258 
01259 
01260 /****************************************************************
01261  *    BBparBperp                                                *
01262  * get Baseline, Bpar and Bperp based on 3 positions,           *
01263  * not efficient! See documentation for definitions             *
01264  *                                                              *
01265  #%// BK 22-Sep-2000                                            *
01266  ****************************************************************/
01267 void BBparBperp(real8 &B, real8 &Bpar, real8 &Bperp,
01268            const cn Master, const cn Point, const cn Slave)
01269   {
01270   #ifdef __DEBUG
01271   TRACE_FUNCTION("BBparBperp (BK 22-Sep-2000)")
01272   #endif
01273   B = Master.dist(Slave);               // baseline. abs. value
01274   const real8 range1 = Master.dist(Point);
01275   const real8 range2 = Slave.dist(Point);
01276   Bpar = range1-range2;                 // parallel baseline, sign ok
01277   const cn r1 = Master.min(Point);
01278   const cn r2 = Slave.min(Point);
01279   const real8 theta = Master.angle(r1);
01280   Bperp = sqr(B)-sqr(Bpar);
01281   if (Bperp < 0.0) Bperp=0.0;
01282   else Bperp = (theta > Master.angle(r2)) ?     // perpendicular baseline, sign ok
01283      sqrt(Bperp) : -sqrt(Bperp);
01284   } // END BBparBperp
01285 
01286 
01287 
01288 /****************************************************************
01289  *    BBhBv                                                     *
01290  * get Baseline, Bh and Bv based on 3 positions,                *
01291  * not efficient! See documentation for definitions             *
01292  * cannot get sign of Bh w/o Point P (?)                        *
01293  *                                                              *
01294  #%// BK 19-Oct-2000                                            *
01295  ****************************************************************/
01296 void BBhBv(
01297         real8 &B,
01298         real8 &Bh,
01299         real8 &Bv,
01300         const cn Master,
01301         const cn Slave)
01302   {
01303   TRACE_FUNCTION("BBhBv (BK 19-Oct-2000)")
01304   WARNING.print("sign Bh not computed...");
01305   B                = Master.dist(Slave);// baseline. abs. value
01306   const real8 rho1 = Master.norm();
01307   const real8 rho2 = Slave.norm();
01308   Bv               = rho2-rho1;                 // Bv, sign ok
01309   Bh               = sqr(B)-sqr(Bv);
01310   Bh               = (Bh<0.0) ? 0.0 : sqrt(Bh);
01311   } // END BBhBv
01312 
01313 
01314 
01315 /****************************************************************
01316  *    Btemp                                                     *
01317  * returns temporal baseline (days) based on 2 strings with UTC *
01318  * utc as: "30-AUG-1995 09:49:20.453"                           *
01319  *         "%d-%b-%Y %T"                                        *
01320  * (Btemp>0) if (Tslave later than Tmaster)                     *
01321  #%// BK 19-Oct-2000                                            *
01322  ****************************************************************/
01323 int32 Btemp(
01324         const char *utc_master,
01325         const char *utc_slave)
01326   {
01327   TRACE_FUNCTION("Btemp (BK 19-Oct-2000)")
01328   struct tm tm_ref;     // reference time
01329   struct tm tm_master;
01330   struct tm tm_slave;
01331   char utc_ref[ONE27] = "01-JAN-1985 00:00:01.000";
01332   strptime(utc_ref,"%d-%b-%Y %T",&tm_ref);
01333   strptime(utc_master,"%d-%b-%Y %T",&tm_master);
01334   strptime(utc_slave, "%d-%b-%Y %T",&tm_slave);
01335 
01336   time_t time_ref      = mktime(&tm_ref);
01337   time_t time_master   = mktime(&tm_master);
01338   time_t time_slave    = mktime(&tm_slave);
01339   real8 master_since85 = difftime(time_master,time_ref);        // 1>2
01340   real8 slave_since85  = difftime(time_slave,time_ref);         // 1>2
01341   int32 Btemp = int32(floor(
01342     ((slave_since85-master_since85)/(60.0*60.0*24.0))+0.5));
01343   return Btemp;
01344   } // END Btemp
01345 
01346 
01347 
01348 /****************************************************************
01349  *    BalphaBhBvBparBperpTheta                                  *
01350  * get Baseline, Bh and Bv based on 3 positions,                *
01351  * not efficient? See documentation for definitions             *
01352  * better use velocity vector, not P?                           *
01353  * input:                                                       *
01354  *  - coordinates x,y,z of Master, Slave en Point on ground     *
01355  * output:                                                      *
01356  *  - abs(B), alpha[-pi,pi], theta[approx23deg], Bh, Bv, Bpar   *
01357  *    Bperp                                                     *
01358  *                                                              *
01359  * theta  is the looking angle to point P.                      *
01360  * theta2 is the looking angle to point P for slave satellite.  *
01361  * cos theta  = (rho1^2+r1^2-p^2) / (2*rho1*r1)                 *
01362  * approx.                                                      *
01363  * cos theta2 = (rho2^2+r2^2-p^2) / (2*rho2*r2)                 *
01364  * if (theta1>theta2) then (Bperp>0)                            *
01365  * if (costheta1>costheta2) then (Bperp<0)                      *
01366  *                                                              *
01367  #%// BK 19-Oct-2000                                            *
01368  ****************************************************************/
01369 void BalphaBhBvBparBperpTheta(
01370         real8 &B,
01371         real8 &alpha,
01372         real8 &Bh,
01373         real8 &Bv,
01374         real8 &Bpar,
01375         real8 &Bperp,
01376         real8 &theta,
01377         const cn     M,
01378         const cn     P,
01379         const cn     S)
01380   {
01381   TRACE_FUNCTION("BalphaBhBvBparBperpTheta (BK 19-Oct-2000)")
01382   const cn    r1        = P.min(M);
01383   const real8 P_2       = P.norm2();
01384   const real8 M_2       = M.norm2();
01385   const real8 r1_2      = r1.norm2();
01386   const real8 costheta1 = (M_2+r1_2-P_2)/(2*sqrt(M_2*r1_2));
01387 
01388   const real8 S_2       = S.norm2();
01389   const cn    r2        = P.min(S);
01390   const real8 r2_2      = r2.norm2();
01391   const real8 costheta2 = (S_2+r2_2-P_2)/(2*sqrt(S_2*r2_2));
01392 
01393   const cn    vecB      = S.min(M);             // (vector B points to S)
01394   const real8 B_2       = vecB.norm2();
01395   Bpar  = sqrt(r1_2)-sqrt(r2_2);                // parallel baseline, sign ok
01396   Bperp = (costheta1>costheta2) ?               // perpendicular bln, sign ok
01397     -sqrt(B_2-sqr(Bpar)) :
01398      sqrt(B_2-sqr(Bpar));
01399     
01400   theta = acos(costheta1);                      // <0,pi/2>, ok
01401   alpha = theta-atan2(Bpar,Bperp);              // <-pi,pi>
01402   B     = sqrt(B_2);                            // abs. value
01403   Bv    = B * sin(alpha);                       // Bv, sign ok
01404   Bh    = B * cos(alpha);                       // Bh, sign ok
01405   } // END BalphaBhBvBparBperpTheta
01406 
01407 
01408 
01409 
01410 /****************************************************************
01411  *    shiftazispectrum                                          *
01412  * Shift spectrum of input matrix data either from fDC to zero, *
01413  * or from zero to fDC (Doppler centroid frequency).            *
01414  * Slcimage gives required polynomial fDC(column).              *
01415  * Abs(shift) gives the first range column for this polynomial. *
01416  *  First column is 1, so use this, not 0!                      *
01417  * (shift<0) indicates to shift towards zero, and               *
01418  * (shift>0) indicates to shift spectrum to fDC (back)          *
01419  * Spectrum is shifted by multiplication by a trend e^{iphase}  *
01420  * in the space domain.                                         *
01421  *                                                              *
01422  * If fDC is more or less equal for all columns, then an        *
01423  * approximation is used. If fDC is smaller then 2 percent of   *
01424  * PRF then no shifting is performed.                           *
01425  * memory intensive cause matrix is used.                       *
01426  *                                                              * 
01427  #%// BK 09-Nov-2000                                            *
01428  ****************************************************************/
01429 void shiftazispectrum(
01430         matrix<complr4> &data,  // slcdata sin space domain
01431         const slcimage  &slave, // polynomial fDC
01432         const real4      shift) // abs(shift)==firstpix in slave system)
01433                                 // if (shift>0) then shift from zero to fDC
01434   {
01435   // ______ Really debug this by defining: ______
01436   //#define REALLYDEBUG
01437   TRACE_FUNCTION("shiftazispectrum (BK 09-Nov-2000)")
01438 
01439   // ______ Evaluate fdc for all columns ______
01440   if (slave.f_DC_a0==0 && slave.f_DC_a1==0 && slave.f_DC_a2==0) // no shift at all
01441     return;
01442   
01443   // ______ Some constants ______
01444   const int32 NLINES = data.lines();            // compute e^ix
01445   const int32 NCOLS  = data.pixels();           // width
01446   const int32 SIGN   = (shift<0) ? -1:1;        // e^{SIGN*i*x}
01447   // -1
01448   const real4 PIXLO  = abs(shift)-1;            // first column for fDC polynomial
01449 
01450   // ______ Compute fDC for all columns ______
01451   // ______ Create axis to evaluate fDC polynomial ______
01452   // ______ fDC(column) = fdc_a0 + fDC_a1*(col/RSR) + fDC_a2*(col/RSR)^2 ______
01453   // ______ offset slave,master defined as: cols=colm+offsetP ______
01454   matrix<real8> FDC(1,NCOLS);
01455   FDC = slave.f_DC_a0;                                  // constant term
01456   if (slave.f_DC_a1!=0 || slave.f_DC_a2!=0)             // check to save time
01457     {
01458     matrix<real8> xaxis(1,NCOLS);                       // lying
01459     for (register int32 ii=0; ii<NCOLS; ++ii)
01460       xaxis(0,ii) = real8(ii+PIXLO);
01461     xaxis /= (slave.rsr2x/2.);
01462     FDC += (slave.f_DC_a1 * xaxis);                     // linear term
01463     FDC += (slave.f_DC_a2 * sqr(xaxis));                // cubic term
01464     #ifdef REALLYDEBUG
01465     cout << "REALLYDEBUG: matrix xaxis dumped.\n";
01466     cout << "REALLYDEBUG: SIGN   = " << SIGN << "\n";
01467     cout << "REALLYDEBUG: fda0   = " << slave.f_DC_a0 << "\n";
01468     cout << "REALLYDEBUG: PIXLO  = " << PIXLO << "\n";
01469     cout << "REALLYDEBUG: shift  = " << shift << "\n";
01470     cout << "REALLYDEBUG: fda1   = " << slave.f_DC_a1 << "\n";
01471     cout << "REALLYDEBUG: fda2   = " << slave.f_DC_a2 << "\n";
01472     cout << "REALLYDEBUG: FDC(1) = " << FDC(0,0) << "\n";
01473     dumpasc("xaxis",xaxis);
01474     cout << "REALLYDEBUG: matrix FDC with doppler centroid dumped.\n";
01475     dumpasc("FDC",FDC);
01476     #endif
01477     }
01478 
01479 
01480   // ====== Actually shift the azimuth spectrum ======
01481   matrix<real8> trend(NLINES,1);
01482   for (register int32 ii=0; ii<NLINES; ++ii)
01483     trend(ii,0) = (real8(2*ii)*PI) / slave.prf;
01484   matrix<real8> P = trend * FDC;
01485   matrix<complr4> TREND =                               // forward or backward
01486     (SIGN==-1) ? mat2cr4(cos(P),sin(-P)) :
01487                  mat2cr4(cos(P),sin(P));
01488 
01489   // free memory, delete P, by resizing...
01490   P.resize(1,1);
01491 
01492   // ______ check if this goes ok with matlab ______
01493   #ifdef REALLYDEBUG
01494   static int32 si_buffernum = 0;
01495   si_buffernum++;
01496   char OFILE[ONE27];
01497   ostrstream tmpomem(OFILE,ONE27);
01498   tmpomem.seekp(0);
01499   tmpomem << "datain." << si_buffernum << ends;
01500   cerr << "REALLYDEBUG: matrix data dumped to file " << OFILE << " (cr4).\n";
01501   cerr << "#lines, #pixs: " << data.lines() << ", " << data.pixels() << endl;
01502   ofstream of;
01503   of.open(OFILE); 
01504   of << data;
01505   of.close();
01506   tmpomem.seekp(0);
01507   tmpomem << "TREND." << si_buffernum << ends;
01508   cerr << "REALLYDEBUG: trend matrix dumped to file " << OFILE << " (cr4).\n";
01509   of.open(OFILE); 
01510   of << TREND;
01511   of.close();
01512   #endif
01513 
01514   data *= TREND;                                        // return by reference
01515 
01516   #ifdef REALLYDEBUG
01517   tmpomem.seekp(0);
01518   tmpomem << "dataout." << si_buffernum << ends;
01519   cerr << "REALLYDEBUG: matrix data dumped to file " << OFILE << " (cr4).\n";
01520   of.open(OFILE); 
01521   of << data;
01522   of.close();
01523   //cerr << "move files to other name?\n";
01524   //getanswer();        //wait
01525   #endif
01526   } // END shiftazispectrum
01527 
01528