00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #include "matrixbk.hh"
00044 #include "slcimage.hh"
00045 #include "constants.hh"
00046 #include "refsystems.hh"
00047 #include "utilities.hh"
00048 #include "ioroutines.hh"
00049 #include "conversion.hh"
00050 #include "exceptions.hh"
00051
00052 #include <strstream>
00053 #include <iomanip>
00054 #include <cstring>
00055 #include <cstdlib>
00056 #include <cmath>
00057 #include <cstdio>
00058 #include <ctime>
00059 char *strptime(const char *s, const char *format, struct tm *tm);
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
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
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
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;
00109 if (t1seconds < 0)
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)
00122 {
00123 WARNING.print("getorb: couldn't adjust time upper bound, using old value.");
00124 t2seconds -= BEFOREFIRST;
00125 }
00126 struct tm tijdend = tijdstart;
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
00135
00136
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
00142 strftime(pstartt,13,"%y%m%d%H%M%S",&tijdstart);
00143 strftime(pendt,13,"%y%m%d%H%M%S",&tijdend);
00144
00145
00146
00147
00148
00149 char strgetorb[ONE27];
00150 ostrstream omemgetorb(strgetorb,ONE27);
00151 omemgetorb.seekp(0);
00152 omemgetorb << "getorb " << startt << "," << endt << "," << INTERVAL << " "
00153 << orbdir << ">scratchorbit" << ends;
00154
00155
00156
00157
00158
00159
00160
00161
00162 int32 status = 3*NaN;
00163 for (register int32 i=0;i<3;i++)
00164 {
00165 DEBUG.print(strgetorb);
00166 status=(system(strgetorb));
00167 DEBUG << "status of getorb: " << status;
00168 DEBUG.print();
00169 if (status == 0)
00170 {
00171
00172 DEBUG.print("Checking getorb result for errors.");
00173
00174 ifstream scratchorbit("scratchorbit",ios::in);
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;
00191 }
00192 scratchorbit.getline(dummyline,ONE27,'\n');
00193 }
00194 scratchorbit.close();
00195 if (getorbok)
00196 {
00197 PROGRESS.print("getorb: program finished ok.");
00198 break;
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
00230
00231
00232
00233
00234
00235
00236
00237
00238 }
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 void convertgetorbout(
00256 int16 FILEID,
00257 const char* file)
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
00268
00269 ifstream infile(file, ios::in);
00270 bk_assert(infile,file,__FILE__,__LINE__);
00271
00272 ofstream outfile("scratchdatapoints", ios::out | ios::trunc);
00273 bk_assert(outfile,"convertgetorbout: scratchdatapoints",__FILE__,__LINE__);
00274
00275
00276 int32 err;
00277 real8 sec;
00278 real8 fie;
00279 real8 lam;
00280 real8 hei;
00281 char x[25];
00282 char y[25];
00283 char z[25];
00284
00285
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();
00300
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();
00307
00308
00309 outfile << "\n\n*******************************************************************"
00310 << "\n*_Start_";
00311
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
00323 int32 secint;
00324 real8 secfrac;
00325 outfile.setf(ios::fixed);
00326 outfile.setf(ios::showpoint);
00327
00328 infile.seekg(0, ios::beg);
00329 DEBUG.print("Rewinding file and reading orbit ephemerides again");
00330 for (register int32 i=0;i<numdatapoints+numerr;i++)
00331 {
00332 infile >> sec >> err >> fie >> lam >> hei >> x >> y >> z;
00333 infile.getline(dummyline,ONE27,'\n');
00334 DEBUG << "scratchfile2: sec=" << sec << " err=" << err
00335 << " x=" << x << " y=" << y << " z=" << z;
00336 DEBUG.print();
00337
00338
00339 secint = int32(sec)%86400;
00340 secfrac = sec - int32(sec);
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
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
00358 infile.close();
00359 outfile.close();
00360 if (remove(file))
00361 WARNING.print("code 101: could not remove file.");
00362 }
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
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
00406
00407
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
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
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 }
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
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
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 }
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
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
00493 return (uint(pow(real8(2.),b)));
00494 }
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
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
00532
00533
00534
00535
00536
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)
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 }
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
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
00675
00676
00677
00678
00679
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)
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 }
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798 matrix<real4> polyval(
00799 const matrix<real4> &x,
00800 const matrix<real4> &y,
00801 const matrix<real8> &coeff,
00802 int32 degreee)
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
00831 if (degreee==-1)
00832 degreee = degree(coeff.size());
00833
00834
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:
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
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
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);
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);
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 }
01164
01165 return Result;
01166 }
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
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 }
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
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
01226 data -= real4(.5*(min+max));
01227 data /= real4(.25*(max-min));
01228 }
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
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
01254 data -= .5*(min+max);
01255 data /= .25*(max-min);
01256 }
01257
01258
01259
01260
01261
01262
01263
01264
01265
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);
01274 const real8 range1 = Master.dist(Point);
01275 const real8 range2 = Slave.dist(Point);
01276 Bpar = range1-range2;
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)) ?
01283 sqrt(Bperp) : -sqrt(Bperp);
01284 }
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
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);
01306 const real8 rho1 = Master.norm();
01307 const real8 rho2 = Slave.norm();
01308 Bv = rho2-rho1;
01309 Bh = sqr(B)-sqr(Bv);
01310 Bh = (Bh<0.0) ? 0.0 : sqrt(Bh);
01311 }
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
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;
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);
01340 real8 slave_since85 = difftime(time_slave,time_ref);
01341 int32 Btemp = int32(floor(
01342 ((slave_since85-master_since85)/(60.0*60.0*24.0))+0.5));
01343 return Btemp;
01344 }
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
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);
01394 const real8 B_2 = vecB.norm2();
01395 Bpar = sqrt(r1_2)-sqrt(r2_2);
01396 Bperp = (costheta1>costheta2) ?
01397 -sqrt(B_2-sqr(Bpar)) :
01398 sqrt(B_2-sqr(Bpar));
01399
01400 theta = acos(costheta1);
01401 alpha = theta-atan2(Bpar,Bperp);
01402 B = sqrt(B_2);
01403 Bv = B * sin(alpha);
01404 Bh = B * cos(alpha);
01405 }
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429 void shiftazispectrum(
01430 matrix<complr4> &data,
01431 const slcimage &slave,
01432 const real4 shift)
01433
01434 {
01435
01436
01437 TRACE_FUNCTION("shiftazispectrum (BK 09-Nov-2000)")
01438
01439
01440 if (slave.f_DC_a0==0 && slave.f_DC_a1==0 && slave.f_DC_a2==0)
01441 return;
01442
01443
01444 const int32 NLINES = data.lines();
01445 const int32 NCOLS = data.pixels();
01446 const int32 SIGN = (shift<0) ? -1:1;
01447
01448 const real4 PIXLO = abs(shift)-1;
01449
01450
01451
01452
01453
01454 matrix<real8> FDC(1,NCOLS);
01455 FDC = slave.f_DC_a0;
01456 if (slave.f_DC_a1!=0 || slave.f_DC_a2!=0)
01457 {
01458 matrix<real8> xaxis(1,NCOLS);
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);
01463 FDC += (slave.f_DC_a2 * sqr(xaxis));
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
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 =
01486 (SIGN==-1) ? mat2cr4(cos(P),sin(-P)) :
01487 mat2cr4(cos(P),sin(P));
01488
01489
01490 P.resize(1,1);
01491
01492
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;
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
01524
01525 #endif
01526 }
01527
01528