diff --git a/app/otbSARAltAmbig.cxx b/app/otbSARAltAmbig.cxx index 97e0b377a6788eb50545673f902c9c0399c95448..a8b2eab147ef5e58c777ce926ea92e9a9b103080 100644 --- a/app/otbSARAltAmbig.cxx +++ b/app/otbSARAltAmbig.cxx @@ -118,6 +118,7 @@ private: bool loadOk = m_SarSensorModelAdapter->LoadState(sarKWL); SarSensorModelAdapter::Point3DType demGeoPoint; + SarSensorModelAdapter::Point3DType demGeoPointRadian; SarSensorModelAdapter::Point3DType xyzCart; SarSensorModelAdapter::Point3DType satellitePosition[2]; SarSensorModelAdapter::Point3DType satelliteVelocity[2]; @@ -127,15 +128,21 @@ private: demGeoPoint[1] = lat; demGeoPoint[2] = height; + demGeoPointRadian[0] = lon*M_PI/180.; + demGeoPointRadian[1] = lat*M_PI/180.; + demGeoPointRadian[2] = height; + SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart); otbAppLogINFO(<<"Cartesian coords : " << xyzCart[0] << " " << xyzCart[1] << " " << xyzCart[2]); - + + // Master otbAppLogINFO(<<"Inmaster Image : "); m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[0], satelliteVelocity[0]); + otbAppLogINFO(<<"Satellite Position : " << satellitePosition[0][0] << " " << satellitePosition[0][1] << " " << satellitePosition[0][2]); otbAppLogINFO(<<"Satellite Velocity : " << satelliteVelocity[0][0] << " " << satelliteVelocity[0][1] << " " << satelliteVelocity[0][2]); - searchPassageCloserToGround(satellitePosition[0], satelliteVelocity[0], incidence[0], xyzCart, demGeoPoint); + searchPassageCloserToGround(satellitePosition[0], incidence[0], xyzCart, demGeoPointRadian); otbAppLogINFO(<<"Incidence : " << incidence[0]); // Calculate lambda @@ -145,49 +152,44 @@ private: float lambda = C/radarFreq; otbAppLogINFO(<<"Lambda : " << lambda); + // Target loadOk = m_SarSensorModelAdapter->LoadState(sarKWL_target); otbAppLogINFO(<<"Target Image : "); m_SarSensorModelAdapter->WorldToSatPositionAndVelocity(demGeoPoint, satellitePosition[1], satelliteVelocity[1]); otbAppLogINFO(<<"Satellite Position : " << satellitePosition[1][0] << " " << satellitePosition[1][1] << " " << satellitePosition[1][2]); otbAppLogINFO(<<"Satellite Velocity : " << satelliteVelocity[1][0] << " " << satelliteVelocity[1][1] << " " << satelliteVelocity[1][2]); - searchPassageCloserToGround(satellitePosition[1], satelliteVelocity[1], incidence[1], xyzCart, demGeoPoint); + searchPassageCloserToGround(satellitePosition[1], incidence[1], xyzCart, demGeoPointRadian); otbAppLogINFO(<<"Incidence : " << incidence[1]); computeAltAmbig(xyzCart, satellitePosition, satelliteVelocity, lambda, incidence); } - void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType& satellitePosition, SarSensorModelAdapter::Point3DType& satelliteVelocity, double& incidence, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType demGeoPoint){ - // Incidence - SarSensorModelAdapter::Point3DType R; - R[0] = xyzCart[0] - satellitePosition[0]; - R[1] = xyzCart[1] - satellitePosition[1]; - R[2] = xyzCart[2] - satellitePosition[2]; + void searchPassageCloserToGround(SarSensorModelAdapter::Point3DType satellitePosition, double& incidence, SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType demGeoPoint){ + SarSensorModelAdapter::Point3DType rm; + rm[0] = xyzCart[0] - satellitePosition[0]; + rm[1] = xyzCart[1] - satellitePosition[1]; + rm[2] = xyzCart[2] - satellitePosition[2]; + // Incidence double normeS = sqrt(satellitePosition[0]*satellitePosition[0] + satellitePosition[1]*satellitePosition[1] + satellitePosition[2]*satellitePosition[2]); double normeCible = sqrt(xyzCart[0] * xyzCart[0] + xyzCart[1] * xyzCart[1] + xyzCart[2] * xyzCart[2]); - double normeR = sqrt(R[0] * R[0] + R[1] * R[1] + R[2] * R[2]); + double normeR = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]); incidence = acos((normeS * normeS - normeR * normeR - normeCible * normeCible) / (2 * normeCible * normeR) ) * 180. / M_PI; // Calculate the local vectors - double rm[3]; - rm[0] = xyzCart[0] - satellitePosition[0]; - rm[1] = xyzCart[1] - satellitePosition[1]; - rm[2] = xyzCart[2] - satellitePosition[2]; - + SarSensorModelAdapter::Point3DType vecti; + SarSensorModelAdapter::Point3DType vectj; + SarSensorModelAdapter::Point3DType vectk; double normeRm = sqrt( rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]); - double vecti[3]; - double vectj[3]; - double vectk[3]; double lon = demGeoPoint[0]; double lat = demGeoPoint[1]; - vecti[0] = -sin(lat)*cos(lon); vecti[1] = -sin(lat)*sin(lon); vecti[2] = cos(lat); @@ -201,26 +203,29 @@ private: vectk[2] = sin(lat); // Rm projection - otbAppLogINFO(<< -(vecti[0]*rm[0]+vecti[1]*rm[1]+vecti[2]*rm[2])/normeRm); - otbAppLogINFO(<< -(vectj[0]*rm[0]+vectj[1]*rm[1]+vectj[2]*rm[2])/normeRm); - otbAppLogINFO(<< -(vectk[0]*rm[0]+vectk[1]*rm[1]+vectk[2]*rm[2])/normeRm); + double northProjection = -(vecti[0]*rm[0]+vecti[1]*rm[1]+vecti[2]*rm[2])/normeRm; + double eastProjection = -(vectj[0]*rm[0]+vectj[1]*rm[1]+vectj[2]*rm[2])/normeRm; + double verticalProjection = -(vectk[0]*rm[0]+vectk[1]*rm[1]+vectk[2]*rm[2])/normeRm; + otbAppLogINFO(<< "North Projection : " << northProjection); + otbAppLogINFO(<< "East Projection : " << eastProjection); + otbAppLogINFO(<< "Vertical Projection : " << verticalProjection); } void computeAltAmbig(SarSensorModelAdapter::Point3DType xyzCart, SarSensorModelAdapter::Point3DType* satellitePosition, SarSensorModelAdapter::Point3DType* satelliteVelocity, double lambda, double* incidence){ otbAppLogINFO(<< "Compute Alt Ambig"); - // TODO : Get the type of phase (Bistatic by default) + // TODO : Get the type of phase (monostatic by default) double factorBperp = 1.0; double factorHamb = 2.0; - double rm[3]; + SarSensorModelAdapter::Point3DType rm; rm[0] = -(xyzCart[0] - satellitePosition[0][0]); rm[1] = -(xyzCart[1] - satellitePosition[0][1]); - rm[0] = -(xyzCart[2] - satellitePosition[0][2]); + rm[2] = -(xyzCart[2] - satellitePosition[0][2]); - double re[3]; + SarSensorModelAdapter::Point3DType re; re[0] = -(xyzCart[0] - satellitePosition[1][0]); re[1] = -(xyzCart[1] - satellitePosition[1][1]); - re[0] = -(xyzCart[2] - satellitePosition[1][2]); + re[2] = -(xyzCart[2] - satellitePosition[1][2]); // project the slave orbit slant range vector // onto the master orbit zero doppler plane @@ -230,7 +235,7 @@ private: // unitary vector with the same direction as the master velocity // this vector is orthogonal to the master zero doppler plane - double vmUnit[3]; + SarSensorModelAdapter::Point3DType vmUnit; vmUnit[0] = satelliteVelocity[0][0] / vmNorm; vmUnit[1] = satelliteVelocity[0][1] / vmNorm; vmUnit[2] = satelliteVelocity[0][2] / vmNorm; @@ -245,7 +250,7 @@ private: re[0] -= reVmunit * vmUnit[0]; re[1] -= reVmunit * vmUnit[1]; re[2] -= reVmunit * vmUnit[2]; - + double normeRm = sqrt(rm[0] * rm[0] + rm[1] * rm[1] + rm[2] * rm[2]); double normeRe = sqrt(re[0] * re[0] + re[1] * re[1] + re[2] * re[2]); @@ -265,7 +270,7 @@ private: double alt_amb = lambda * sin(inc_moy) / (factorHamb * tan(alpha)); otbAppLogINFO(<< "Incidence\tALT AMBIG\tRadial\tLateral"); - otbAppLogINFO(<< inc_moy * 180. / M_PI << "\t\t" << alt_amb << "\t\t" << xe << "\t" << ye / factorBperp); + otbAppLogINFO(<< inc_moy * 180. / M_PI << "\t" << alt_amb << "\t" << xe << "\t" << ye / factorBperp); }