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 #include "matrixbk.hh"
00042
00043 #include "constants.hh"
00044 #include "readinput.hh"
00045 #include "orbitbk.hh"
00046 #include "slcimage.hh"
00047 #include "exceptions.hh"
00048
00049 #include "conversion.hh"
00050 #include "ioroutines.hh"
00051 #include "utilities.hh"
00052 #include "refsystems.hh"
00053
00054 #include <cmath>
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 void pol2xyz(
00069 cn &xyz,
00070 real8 phi,
00071 real8 lon,
00072 real8 height)
00073 {
00074 TRACE_FUNCTION("pol2xyz (BK 21-Dec-1998)");
00075 real8 Rph = RADIUS + height;
00076 xyz.x = Rph *cos(phi)*cos(lon);
00077 xyz.y = Rph *cos(phi)*sin(lon);
00078 xyz.z = Rph *sin(phi);
00079 }
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093 void xyz2pol(
00094 const cn &xyz,
00095 real8 &phi,
00096 real8 &lambda,
00097 real8 &height)
00098 {
00099 TRACE_FUNCTION("xyz2pol (BK 21-Dec-1998)");
00100 real8 r = sqrt(sqr(xyz.x)+sqr(xyz.y));
00101 phi = atan2(xyz.z,r);
00102 lambda = atan2(xyz.y,xyz.x);
00103 height = (r/cos(phi)) - RADIUS;
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void xyz2ell(
00123 const input_ell &ell,
00124 const cn &xyz,
00125 real8 &phi,
00126 real8 &lambda,
00127 real8 &height)
00128 {
00129 TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999)");
00130 real8 r = sqrt(sqr(xyz.x)+sqr(xyz.y));
00131 real8 nu = atan2((xyz.z*ell.a),(r*ell.b));
00132 real8 sin3 = pow(sin(nu),3);
00133 real8 cos3 = pow(cos(nu),3);
00134
00135 phi = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
00136 lambda = atan2(xyz.y,xyz.x);
00137 real8 N = ell.a / sqrt(1-ell.e2*sqr(sin(phi)));
00138 height = (r/cos(phi)) - N;
00139 }
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 void xyz2ell(
00158 const input_ell &ell,
00159 const cn &xyz,
00160 real8 &phi,
00161 real8 &lambda)
00162 {
00163 TRACE_FUNCTION("xyz2ell (BK 05-Jan-1999");
00164 real8 r = sqrt(sqr(xyz.x)+sqr(xyz.y));
00165 real8 nu = atan2((xyz.z*ell.a),(r*ell.b));
00166 real8 sin3 = pow(sin(nu),3);
00167 real8 cos3 = pow(cos(nu),3);
00168
00169 phi = atan2((xyz.z+ell.e2b*ell.b*sin3),(r-ell.e2*ell.a*cos3));
00170 lambda = atan2(xyz.y,xyz.x);
00171
00172
00173 }
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 void ell2xyz(
00190 const input_ell &ell,
00191 cn &xyz,
00192 real8 phi,
00193 real8 lambda,
00194 real8 height)
00195 {
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215 const real8 sinphi = sin(phi);
00216 const real8 N = ell.a / sqrt(1.0-ell.e2*sqr(sinphi));
00217 const real8 Nph = N + height;
00218 const real8 Nph_cosphi = Nph*cos(phi);
00219
00220 xyz.x = Nph_cosphi * cos(lambda);
00221 xyz.y = Nph_cosphi * sin(lambda);
00222 xyz.z = (Nph - ell.e2*N) * sinphi;
00223 }
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243 void tiepoint(
00244 const input_gen &generalinput,
00245 const slcimage &master,
00246 const slcimage &slave,
00247 orbit &masterorbit,
00248 orbit &slaveorbit,
00249 const input_ell &ellips)
00250 {
00251 TRACE_FUNCTION("tiepoint (Bert Kampes 14-May-2004)")
00252 if (abs(generalinput.tiepoint.x) < 0.001 &&
00253 abs(generalinput.tiepoint.y) < 0.001 &&
00254 abs(generalinput.tiepoint.z) < 0.001)
00255 {
00256 INFO.print("No tiepoint given.");
00257 return;
00258 }
00259 if (masterorbit.npoints()==0)
00260 {
00261 INFO.print("Exiting tiepoint, no orbit data master available.");
00262 return;
00263 }
00264 DEBUG.precision(30);
00265 DEBUG.width(14);
00266 INFO.precision(13);
00267 INFO.width(14);
00268 const int32 MAXITER = 10;
00269 const real8 CRITERPOS = 1e-6;
00270 const real8 CRITERTIM = 1e-10;
00271 const real8 m_minpi4cdivlam = (-4.0*PI*SOL)/master.wavelength;
00272 const real8 s_minpi4cdivlam = (-4.0*PI*SOL)/slave.wavelength;
00273 DEBUG << "Master wavelength = " << master.wavelength;
00274 DEBUG.print();
00275 DEBUG << "Slave wavelength = " << slave.wavelength;
00276 DEBUG.print();
00277
00278
00279 ;
00280 INFO.print("Computing radar coordinates and unwrapped phase of tie point");
00281 real8 tiepointphi = deg2rad(generalinput.tiepoint.x);
00282 real8 tiepointlambda = deg2rad(generalinput.tiepoint.y);
00283 real8 tiepointheight = generalinput.tiepoint.z;
00284 DEBUG << "TIEPOINT: lat/lon/hei [rad/rad/m]: "
00285 << tiepointphi << " " << tiepointlambda << " " << tiepointheight;
00286 DEBUG.print();
00287 cn tiepointpos;
00288 ell2xyz(ellips,tiepointpos,tiepointphi,tiepointlambda,tiepointheight);
00289 INFO << "TIEPOINT lat/lon/hei --> x,y,z (WGS84): "
00290 << tiepointphi << ", " << tiepointlambda << ", " << tiepointheight
00291 << " --> "
00292 << tiepointpos.x << " " << tiepointpos.y << " " << tiepointpos.z;
00293 INFO.print();
00294
00295
00296 real8 tiepointline;
00297 real8 tiepointpix;
00298 int32 n_iter = xyz2lp(tiepointline, tiepointpix,
00299 master, masterorbit, tiepointpos,
00300 MAXITER, CRITERTIM);
00301 INFO << "TIEPOINT line/pix in master: " << tiepointline << " " << tiepointpix;
00302 INFO.print();
00303
00304
00305 const real8 m_aztime = master.line2ta(tiepointline);
00306 const real8 m_trange = master.pix2tr(tiepointpix);
00307 DEBUG << "TIEPOINT master azimuth time: " << m_aztime;
00308 DEBUG.print();
00309 DEBUG << "TIEPOINT master range time: " << m_trange;
00310 DEBUG.print();
00311
00312
00313
00314 cn P_ref;
00315 lp2xyz(tiepointline,tiepointpix,
00316 ellips,master,masterorbit,P_ref,
00317 MAXITER,CRITERPOS);
00318 INFO << "Pixel of TIEPOINT on ellipsoid: x,y,z (WGS84): "
00319 << P_ref.x << " " << P_ref.y << " " << P_ref.z;
00320 INFO.print();
00321
00322
00323 if (slaveorbit.npoints()==0)
00324 {
00325 INFO.print("tiepoint: no orbit data slave available cannot compute more.");
00326 DEBUG.reset();
00327 INFO.reset();
00328 }
00329 else
00330 {
00331
00332 real8 s_aztime;
00333 real8 s_trange;
00334 xyz2t(s_aztime,s_trange,
00335 slave, slaveorbit, tiepointpos,MAXITER,CRITERTIM);
00336 DEBUG << "TIEPOINT slave azimuth time: " << s_aztime;
00337 DEBUG.print();
00338 DEBUG << "TIEPOINT slave range time: " << s_trange;
00339 DEBUG.print();
00340 INFO << "TIEPOINT line/pix in slave: "
00341 << slave.ta2line(s_aztime) << " "
00342 << slave.tr2pix (s_trange);
00343 INFO.print();
00344
00345
00346 real8 phase = m_minpi4cdivlam*m_trange - s_minpi4cdivlam*s_trange;
00347 INFO << "TIEPOINT unwrapped phase (incl. reference phase): "
00348 << phase;
00349 INFO.print();
00350
00351
00352 real8 s_tazi_ref;
00353 real8 s_trange_ref;
00354 xyz2t(s_tazi_ref,s_trange_ref,slave,
00355 slaveorbit, P_ref,
00356 MAXITER,CRITERTIM);
00357
00358
00359
00360
00361 real8 phase_ref_corrected = s_minpi4cdivlam*(s_trange_ref-s_trange);
00362 INFO << "TIEPOINT unwrapped phase (ellipsoid ref.phase corrected) [rad]: "
00363 << phase_ref_corrected;
00364 INFO.print();
00365 INFO.print("This phase can be used to correct the unwrapped interferogram");
00366 INFO.print("make sure that you compute the correct pixel position");
00367 INFO.print("the line/pix here are in original master coordinates.");
00368
00369
00370
00371
00372
00373 real8 t_az, t_rg;
00374
00375
00376 n_iter = xyz2t(t_az,t_rg,master,masterorbit,P_ref,MAXITER,CRITERTIM);
00377 const cn M = masterorbit.getxyz(t_az);
00378
00379 cn P_ref_checkM;
00380 n_iter = lp2xyz(master.ta2line(t_az), master.tr2pix(t_rg),
00381 ellips, master, masterorbit,
00382 P_ref_checkM, MAXITER, CRITERPOS);
00383
00384 n_iter = xyz2t(t_az,t_rg,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);
00385 const cn S = slaveorbit.getxyz(t_az);
00386
00387 cn P_ref_checkS;
00388 n_iter = lp2xyz(slave.ta2line(t_az), slave.tr2pix(t_rg),
00389 ellips, slave, slaveorbit,
00390 P_ref_checkS, MAXITER, CRITERPOS);
00391
00392 cn M_check;
00393 n_iter = xyz2orb(M_check,master,masterorbit,P_ref,MAXITER,CRITERTIM);
00394 cn S_check;
00395 n_iter = xyz2orb(S_check,slave, slaveorbit, P_ref,MAXITER,CRITERTIM);
00396
00397 INFO.print("Checking zero-Doppler iteration by doing forward and inverse transform:");
00398 DEBUG << "M(x,y,z) = " << M.x << "; " << M.y << "; "<< M.z;
00399 DEBUG.print();
00400 DEBUG << "M'(x,y,z)= " << M_check.x << "; " << M_check.y << "; "<< M_check.z;
00401 DEBUG.print();
00402 DEBUG << "P(x,y,z) = " << P_ref.x << "; " << P_ref.y << "; "<< P_ref.z;
00403 DEBUG.print();
00404 DEBUG << "P'(x,y,z)= " << P_ref_checkM.x << "; " << P_ref_checkM.y << "; "<< P_ref_checkM.z;
00405 DEBUG.print();
00406 DEBUG << "P''(x,y,z)=" << P_ref_checkS.x << "; " << P_ref_checkS.y << "; "<< P_ref_checkS.z;
00407 DEBUG.print();
00408 DEBUG << "S(x,y,z) = " << S.x << "; " << S.y << "; "<< S.z;
00409 DEBUG.print();
00410 DEBUG << "S'(x,y,z)= " << S_check.x << "; " << S_check.y << "; "<< S_check.z;
00411 DEBUG.print();
00412 if (abs(M.dist(M_check))>0.001)
00413 WARNING.print("check failed for M_check");
00414 if (abs(S.dist(S_check))>0.001)
00415 WARNING.print("check failed for S_check");
00416 if (abs(P_ref.dist(P_ref_checkM))>0.001)
00417 WARNING.print("check failed for P_ref_checkM");
00418 if (abs(P_ref.dist(P_ref_checkS))>0.001)
00419 WARNING.print("check failed for P_ref_checkS");
00420 INFO.print("If no warnings seen then zero-Doppler is within 1 mm accurate.");
00421 }
00422
00423
00424 PROGRESS.print("Finished tiepoint.");
00425 DEBUG.reset();
00426 INFO.reset();
00427 }
00428