math - online - coordenadas polares y rectangulares
Cómo convertir las coordenadas esféricas de velocidad en cartesianas (2)
Tengo un vector de velocidad en altitud, longitud y altitud. Me gustaría convertirlo en coordenadas cartesianas, vx, vy, vz. El formato es del estándar WGS84.
aquí está la fórmula
//------------------------------------------------------------------------------
template <class T>
TVectorXYZ<T> WGS84::ToCartesian(T latitude, T longitude, T elevation)
//------------------------------------------------------------------------------
{
double sinlat, coslat;
double sinlon, coslon;
sincos_degree(latitude, sinlat, coslat);
sincos_degree(longitude, sinlon, coslon);
const double v = a / sqrt(1 - WGS84::ee * sinlat*sinlat);
TVectorXYZ<T> coord
(
static_cast<T>((v + elevation) * coslat * sinlon),
static_cast<T>(((1 - WGS84::ee) * v + elevation) * sinlat),
static_cast<T>((v + elevation) * coslat * coslon)
);
return coord;
}
Tome la fórmula que usa para convertir posiciones de coordenadas geográficas a coordenadas cartesianas. Ese es un vector p (λ, φ, h) ∈ ℝ³, es decir, convierte la latitud, la longitud y la altitud en un vector de tres elementos de coordenadas x, y, z. Ahora calcule las derivadas parciales de esta fórmula con respecto a los tres parámetros. Obtendrás tres vectores, que deberían ser ortogonales entre sí. La derivada con respecto a la longitud λ debe apuntar localmente hacia el este, la con respecto a la latitud φ apuntando al norte, la derivada con respecto a la altitud h hacia arriba. Multiplica estos vectores con las velocidades que tienes para obtener un vector de velocidad cartesiana.
Observe cómo las unidades coinciden: la posición está en metros, las dos primeras derivadas son metros por grado y la velocidad sería grados por segundo. O algo más, quizás millas y radianes.
Todo esto es bastante fácil para la esfera. Para el elipsoide WGS84, la fórmula de posición está un poco más involucrada, y esa complejidad se realizará a través del cálculo.
Aceptar en función de su pregunta anterior y el flujo de comentarios largos permite suponer que su entrada es:
lon [rad], lat [rad], alt [m] // WGS84 position
vlon [m/s], vlat [m/s], alt [m/s] // speed in WGS84 lon,lat,alt directions but in [m/s]
Y quiere salida:
x,y,z // Cartesian position [m/s]
vx,vy,vz // Cartesian velocity [m/s]
Y tenga una transformación válida en coordenadas cartesianas para las posiciones a su disposición, esto es mío:
void WGS84toXYZ(double &x,double &y,double &z,double lon,double lat,double alt) // [rad,rad,m] -> [m,m,m]
{
const double _earth_a=6378137.00000; // [m] WGS84 equator radius
const double _earth_b=6356752.31414; // [m] WGS84 epolar radius
const double _earth_e=8.1819190842622e-2; // WGS84 eccentricity
const double _aa=_earth_a*_earth_a;
const double _ee=_earth_e*_earth_e;
double a,b,x,y,z,h,l,c,s;
a=lon;
b=lat;
h=alt;
c=cos(b);
s=sin(b);
// WGS84 from eccentricity
l=_earth_a/sqrt(1.0-(_ee*s*s));
x=(l+h)*c*cos(a);
y=(l+h)*c*sin(a);
z=(((1.0-_ee)*l)+h)*s;
}
Y rutina para normalizar el vector al tamaño de la unidad:
void normalize(double &x,double &y,double &z)
{
double l=sqrt(x*x+y*y+z*z);
if (l>1e-6) l=1.0/l;
x*=l; y*=l; z*=l;
}
Sí, puedes intentar derivar la fórmula que sugiere @MvG, pero de tus errores de novato dudo mucho que conduzca a un resultado exitoso. En cambio, puedes hacer esto:
obtenga vectores de dirección
lon,lat,altpara su posición(x,y,z)eso es fácil, simplemente use un pequeño incremento de paso en la posición WGS84 , conviértalo en resta cartesiana y normalícelos en vectores unitarios. Llamemos a estos vectores de base de dirección
U,V,Wdouble Ux,Uy,Uz; // [m] double Vx,Vy,Vz; // [m] double Wx,Wy,Wz; // [m] double da=1.567e-7; // [rad] angular step ~ 1.0 m in lon direction double dl=1.0; // [m] altitide step 1.0 m WGS84toXYZ( x, y, z,lon ,lat,alt ); // actual position WGS84toXYZ(Ux,Uy,Uz,lon+da,lat,alt ); // lon direction Nort WGS84toXYZ(Vx,Vy,Vz,lon,lat+da,alt ); // lat direction East WGS84toXYZ(Wx,Wy,Wz,lon,lat ,alt+dl); // alt direction High/Up Ux-=x; Uy-=y; Uz-=z; Vx-=x; Vy-=y; Vz-=z; Wx-=x; Wy-=y; Wz-=z; normalize(Ux,Uy,Uz); normalize(Vx,Vy,Vz); normalize(Wx,Wy,Wz);convertir velocidad desde
lon,lat,altavx,vy,vzvx = vlon*Ux + vlat*Vx + valt*Wx; vy = vlon*Uy + vlat*Vy + valt*Wy; vz = vlon*Uz + vlat*Vz + valt*Wz;
Espero que sea lo suficientemente claro. Como de costumbre, tenga cuidado con las unidades deg/rad y m/ft/km porque las unidades importan mucho.
Los vectores base Btw U,V,W forman el marco de referencia NEH y al mismo tiempo son las derivadas de dirección que MvG menciona.
[Editar1] conversiones más precisas
//---------------------------------------------------------------------------
//--- WGS84 transformations ver: 1.00 ---------------------------------------
//---------------------------------------------------------------------------
#ifndef _WGS84_h
#define _WGS84_h
//---------------------------------------------------------------------------
// http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion
//---------------------------------------------------------------------------
// WGS84(a,b,h) = (long,lat,alt) [rad,rad,m]
// XYZ(x,y,z) [m]
//---------------------------------------------------------------------------
const double _earth_a=6378137.00000; // [m] WGS84 equator radius
const double _earth_b=6356752.31414; // [m] WGS84 epolar radius
const double _earth_e=8.1819190842622e-2; // WGS84 eccentricity
//const double _earth_e=sqrt(1.0-((_earth_b/_earth_a)*(_earth_b/_earth_a)));
const double _earth_ee=_earth_e*_earth_e;
//---------------------------------------------------------------------------
const double kmh=1.0/3.6; // [km/h] -> [m/s]
//---------------------------------------------------------------------------
void XYZtoWGS84 (double *abh ,double *xyz ); // [m,m,m] -> [rad,rad,m]
void WGS84toXYZ (double *xyz ,double *abh ); // [rad,rad,m] -> [m,m,m]
void WGS84toXYZ_posvel(double *xyzpos,double *xyzvel,double *abhpos,double *abhvel); // [rad,rad,m],[m/s,m/s,m/s] -> [m,m,m],[m/s,m/s,m/s]
void WGS84toNEH (reper &neh ,double *abh ); // [rad,rad,m] -> NEH [m]
void WGS84_m2rad (double &da,double &db,double *abh); // [rad,rad,m] -> [rad],[rad] representing 1m angle step
void XYZ_interpolate (double *pt,double *p0,double *p1,double t); // [m,m,m] pt = p0 + (p1-p0)*t in ellipsoid space t = <0,1>
//---------------------------------------------------------------------------
void XYZtoWGS84(double *abh,double *xyz)
{
int i;
double a,b,h,l,n,db,s;
a=atanxy(xyz[0],xyz[1]);
l=sqrt((xyz[0]*xyz[0])+(xyz[1]*xyz[1]));
// estimate lat
b=atanxy((1.0-_earth_ee)*l,xyz[2]);
// iterate to improve accuracy
for (i=0;i<100;i++)
{
s=sin(b); db=b;
n=divide(_earth_a,sqrt(1.0-(_earth_ee*s*s)));
h=divide(l,cos(b))-n;
b=atanxy((1.0-divide(_earth_ee*n,n+h))*l,xyz[2]);
db=fabs(db-b);
if (db<1e-12) break;
}
if (b>0.5*pi) b-=pi2;
abh[0]=a;
abh[1]=b;
abh[2]=h;
}
//---------------------------------------------------------------------------
void WGS84toXYZ(double *xyz,double *abh)
{
double a,b,h,l,c,s;
a=abh[0];
b=abh[1];
h=abh[2];
c=cos(b);
s=sin(b);
// WGS84 from eccentricity
l=_earth_a/sqrt(1.0-(_earth_ee*s*s));
xyz[0]=(l+h)*c*cos(a);
xyz[1]=(l+h)*c*sin(a);
xyz[2]=(((1.0-_earth_ee)*l)+h)*s;
}
//---------------------------------------------------------------------------
void WGS84toNEH(reper &neh,double *abh)
{
double N[3],E[3],H[3]; // [m]
double p[3],xyzpos[3];
const double da=1.567e-7; // [rad] angular step ~ 1.0 m in lon direction
const double dl=1.0; // [m] altitide step 1.0 m
vector_copy(p,abh);
// actual position
WGS84toXYZ(xyzpos,abh);
// NEH
p[0]+=da; WGS84toXYZ(N,p); p[0]-=da;
p[1]+=da; WGS84toXYZ(E,p); p[1]-=da;
p[2]+=dl; WGS84toXYZ(H,p); p[2]-=dl;
vector_sub(N,N,xyzpos);
vector_sub(E,E,xyzpos);
vector_sub(H,H,xyzpos);
vector_one(N,N);
vector_one(E,E);
vector_one(H,H);
neh._rep=1;
neh._inv=0;
// axis X
neh.rep[ 0]=N[0];
neh.rep[ 1]=N[1];
neh.rep[ 2]=N[2];
// axis Y
neh.rep[ 4]=E[0];
neh.rep[ 5]=E[1];
neh.rep[ 6]=E[2];
// axis Z
neh.rep[ 8]=H[0];
neh.rep[ 9]=H[1];
neh.rep[10]=H[2];
// gpos
neh.rep[12]=xyzpos[0];
neh.rep[13]=xyzpos[1];
neh.rep[14]=xyzpos[2];
neh.orto(1);
}
//---------------------------------------------------------------------------
void WGS84toXYZ_posvel(double *xyzpos,double *xyzvel,double *abhpos,double *abhvel)
{
reper neh;
WGS84toNEH(neh,abhpos);
neh.gpos_get(xyzpos);
neh.l2g_dir(xyzvel,abhvel);
}
//---------------------------------------------------------------------------
void WGS84_m2rad(double &da,double &db,double *abh)
{
// WGS84 from eccentricity
double p[3],rr;
WGS84toXYZ(p,abh);
rr=(p[0]*p[0])+(p[1]*p[1]);
da=divide(1.0,sqrt(rr));
rr+=p[2]*p[2];
db=divide(1.0,sqrt(rr));
}
//---------------------------------------------------------------------------
void XYZ_interpolate(double *pt,double *p0,double *p1,double t)
{
const double mz=_earth_a/_earth_b;
const double _mz=_earth_b/_earth_a;
double p[3],r,r0,r1;
// compute spherical radiuses of input points
r0=sqrt((p0[0]*p0[0])+(p0[1]*p0[1])+(p0[2]*p0[2]*mz*mz));
r1=sqrt((p1[0]*p1[0])+(p1[1]*p1[1])+(p1[2]*p1[2]*mz*mz));
// linear interpolation
r = r0 +(r1 -r0 )*t;
p[0]= p0[0]+(p1[0]-p0[0])*t;
p[1]= p0[1]+(p1[1]-p0[1])*t;
p[2]=(p0[2]+(p1[2]-p0[2])*t)*mz;
// correct radius and rescale back
r/=sqrt((p[0]*p[0])+(p[1]*p[1])+(p[2]*p[2]));
pt[0]=p[0]*r;
pt[1]=p[1]*r;
pt[2]=p[2]*r*_mz;
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
Sin embargo, requieren una matemática vectorial básica en 3D. Consulte aquí las ecuaciones:
- Comprender las matrices de transformada homogénea 4x4