220 lines
6.1 KiB
C
220 lines
6.1 KiB
C
|
/*
|
|||
|
|
Functions that calculate the gravitational potential and its first and second
|
|||
|
|
derivatives for the rectangular prism in spherical coordinates.
|
|||
|
|
Uses the formulas in Nagy et al. (2000).
|
|||
|
|
|
|||
|
|
References
|
|||
|
|
----------
|
|||
|
|
|
|||
|
|
* Nagy, D., Papp, G., Benedek, J. (2000): The gravitational potential and its
|
|||
|
|
derivatives for the prism. Journal of Geodesy, 74, 552–560.
|
|||
|
|
*/
|
|||
|
|
|
|||
|
|
|
|||
|
|
#include <math.h>
|
|||
|
|
#include "geometry.h"
|
|||
|
|
#include "constants.h"
|
|||
|
|
#include "grav_prism_sph.h"
|
|||
|
|
#include "grav_prism.h"
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Transform spherical coordinates to local Cartesian coordinates of the prism*/
|
|||
|
|
int global2local(double lon, double lat, double r, PRISM prism, double *x,
|
|||
|
|
double *y, double *z)
|
|||
|
|
{
|
|||
|
|
double cosa, cosb, sina, sinb, d2r, X, Y, Z;
|
|||
|
|
|
|||
|
|
/* degrees to radians */
|
|||
|
|
d2r = PI/180.;
|
|||
|
|
|
|||
|
|
X = r*cos(d2r*lat)*cos(d2r*lon) -
|
|||
|
|
prism.r*cos(d2r*prism.lat)*cos(d2r*prism.lon);
|
|||
|
|
Y = r*cos(d2r*lat)*sin(d2r*lon) -
|
|||
|
|
prism.r*cos(d2r*prism.lat)*sin(d2r*prism.lon);
|
|||
|
|
Z = r*sin(d2r*lat) - prism.r*sin(d2r*prism.lat);
|
|||
|
|
|
|||
|
|
cosa = cos(d2r*(90 - prism.lat));
|
|||
|
|
sina = sin(d2r*(90 - prism.lat));
|
|||
|
|
cosb = cos(d2r*(180 - prism.lon));
|
|||
|
|
sinb = sin(d2r*(180 - prism.lon));
|
|||
|
|
|
|||
|
|
*x = X*cosa*cosb - Y*cosa*sinb + Z*sina;
|
|||
|
|
*y = -X*sinb - Y*cosb;
|
|||
|
|
/* -1 because Nagy et al. (2000) use z->down */
|
|||
|
|
*z = -1*(-X*sina*cosb + Y*sina*sinb + Z*cosa);
|
|||
|
|
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Rotate the gravity vector from the prisms coordinate system to the local
|
|||
|
|
system of the computation point. */
|
|||
|
|
int g_prism2point(double *atprism, PRISM prism, double lon, double lat,
|
|||
|
|
double r, double *atpoint)
|
|||
|
|
{
|
|||
|
|
#define POS(x, y, cols) (((x)*(cols))+(y))
|
|||
|
|
|
|||
|
|
register int i, k;
|
|||
|
|
double R[9], d2r, cosbeta, sinbeta, cosphi, sinphi, cosphil, sinphil;
|
|||
|
|
|
|||
|
|
/* degrees to radians */
|
|||
|
|
d2r = PI/180.;
|
|||
|
|
|
|||
|
|
cosbeta = cos(d2r*(prism.lon - lon));
|
|||
|
|
sinbeta = sin(d2r*(prism.lon - lon));
|
|||
|
|
cosphi = cos(d2r*lat);
|
|||
|
|
sinphi = sin(d2r*lat);
|
|||
|
|
cosphil = cos(d2r*prism.lat);
|
|||
|
|
sinphil = sin(d2r*prism.lat);
|
|||
|
|
|
|||
|
|
/* The transformation matrix */
|
|||
|
|
R[0] = cosbeta*sinphi*sinphil + cosphi*cosphil;
|
|||
|
|
R[1] = sinbeta*sinphi;
|
|||
|
|
R[2] = -cosbeta*sinphi*cosphil + cosphi*sinphil;
|
|||
|
|
R[3] = -sinbeta*sinphil;
|
|||
|
|
R[4] = cosbeta;
|
|||
|
|
R[5] = sinbeta*cosphil;
|
|||
|
|
R[6] = -cosbeta*cosphi*sinphil + sinphi*cosphil;
|
|||
|
|
R[7] = -sinbeta*cosphi;
|
|||
|
|
R[8] = cosbeta*cosphi*cosphil + sinphi*sinphil;
|
|||
|
|
|
|||
|
|
/* Matrix-vector multiplication */
|
|||
|
|
for(i = 0; i < 3; i++)
|
|||
|
|
{
|
|||
|
|
atpoint[i] = 0;
|
|||
|
|
for(k = 0; k < 3; k++)
|
|||
|
|
{
|
|||
|
|
atpoint[i] += R[POS(i, k, 3)]*atprism[k];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
#undef POS
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Rotate the gravity tensor from the prisms coordinate system to the local
|
|||
|
|
system of the computation point. */
|
|||
|
|
int ggt_prism2point(double *atprism, PRISM prism, double lon, double lat,
|
|||
|
|
double r, double *atpoint)
|
|||
|
|
{
|
|||
|
|
#define POS(x, y, cols) (((x)*(cols))+(y))
|
|||
|
|
|
|||
|
|
register int i, j, k;
|
|||
|
|
double R[9], tmp[9], d2r, cosbeta, sinbeta, cosphi, sinphi, cosphil, sinphil;
|
|||
|
|
|
|||
|
|
/* degrees to radians */
|
|||
|
|
d2r = PI/180.;
|
|||
|
|
|
|||
|
|
cosbeta = cos(d2r*(prism.lon - lon));
|
|||
|
|
sinbeta = sin(d2r*(prism.lon - lon));
|
|||
|
|
cosphi = cos(d2r*lat);
|
|||
|
|
sinphi = sin(d2r*lat);
|
|||
|
|
cosphil = cos(d2r*prism.lat);
|
|||
|
|
sinphil = sin(d2r*prism.lat);
|
|||
|
|
|
|||
|
|
/* The transformation matrix */
|
|||
|
|
R[0] = cosbeta*sinphi*sinphil + cosphi*cosphil;
|
|||
|
|
R[1] = sinbeta*sinphi;
|
|||
|
|
R[2] = -cosbeta*sinphi*cosphil + cosphi*sinphil;
|
|||
|
|
R[3] = -sinbeta*sinphil;
|
|||
|
|
R[4] = cosbeta;
|
|||
|
|
R[5] = sinbeta*cosphil;
|
|||
|
|
R[6] = -cosbeta*cosphi*sinphil + sinphi*cosphil;
|
|||
|
|
R[7] = -sinbeta*cosphi;
|
|||
|
|
R[8] = cosbeta*cosphi*cosphil + sinphi*sinphil;
|
|||
|
|
|
|||
|
|
/* Multiply tmp = R*Tensor */
|
|||
|
|
for(i = 0; i < 3; i++)
|
|||
|
|
{
|
|||
|
|
for(j = 0; j < 3; j++)
|
|||
|
|
{
|
|||
|
|
tmp[POS(i, j, 3)] = 0;
|
|||
|
|
for(k = 0; k < 3; k++)
|
|||
|
|
{
|
|||
|
|
tmp[POS(i, j, 3)] += R[POS(i, k, 3)]*atprism[POS(k, j, 3)];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Multiply tmp*R^T */
|
|||
|
|
for(i = 0; i < 3; i++)
|
|||
|
|
{
|
|||
|
|
for(j = 0; j < 3; j++)
|
|||
|
|
{
|
|||
|
|
atpoint[POS(i, j, 3)] = 0;
|
|||
|
|
for(k = 0; k < 3; k++)
|
|||
|
|
{
|
|||
|
|
atpoint[POS(i, j, 3)] += tmp[POS(i, k, 3)]*R[POS(j, k, 3)];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
#undef POS
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gravity gradient tensor caused by a prism. */
|
|||
|
|
int prism_ggt_sph(PRISM prism, double lonp, double latp, double rp, double *ggt)
|
|||
|
|
{
|
|||
|
|
double x = 0, y = 0, z = 0, ggtprism[9], ggtpoint[9];
|
|||
|
|
|
|||
|
|
global2local(lonp, latp, rp, prism, &x, &y, &z);
|
|||
|
|
ggtprism[0] = prism_gxx(prism, x, y, z);
|
|||
|
|
ggtprism[1] = prism_gxy(prism, x, y, z);
|
|||
|
|
/* -1 because the prisms z is Down, but transformation assumes z is Up */
|
|||
|
|
/* z -> Up is the system of the tesseroid */
|
|||
|
|
ggtprism[2] = -1*prism_gxz(prism, x, y, z);
|
|||
|
|
ggtprism[3] = ggtprism[1];
|
|||
|
|
ggtprism[4] = prism_gyy(prism, x, y, z);
|
|||
|
|
/* Same as xz */
|
|||
|
|
ggtprism[5] = -1*prism_gyz(prism, x, y, z);
|
|||
|
|
ggtprism[6] = ggtprism[2];
|
|||
|
|
ggtprism[7] = ggtprism[5];
|
|||
|
|
ggtprism[8] = -(ggtprism[0] + ggtprism[4]);
|
|||
|
|
ggt_prism2point(ggtprism, prism, lonp, latp, rp, ggtpoint);
|
|||
|
|
ggt[0] = ggtpoint[0];
|
|||
|
|
ggt[1] = ggtpoint[1];
|
|||
|
|
ggt[2] = ggtpoint[2];
|
|||
|
|
ggt[3] = ggtpoint[4];
|
|||
|
|
ggt[4] = ggtpoint[5];
|
|||
|
|
ggt[5] = ggtpoint[8];
|
|||
|
|
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gravitational attraction caused by a prism. */
|
|||
|
|
int prism_g_sph(PRISM prism, double lonp, double latp, double rp, double *gx,
|
|||
|
|
double *gy, double *gz)
|
|||
|
|
{
|
|||
|
|
double x = 0, y = 0, z = 0, gprism[3], gpoint[3];
|
|||
|
|
|
|||
|
|
global2local(lonp, latp, rp, prism, &x, &y, &z);
|
|||
|
|
gprism[0] = prism_gx(prism, x, y, z);
|
|||
|
|
gprism[1] = prism_gy(prism, x, y, z);
|
|||
|
|
/* Nagy wants z down, but the transformation assumes z up */
|
|||
|
|
gprism[2] = -prism_gz(prism, x, y, z);
|
|||
|
|
g_prism2point(gprism, prism, lonp, latp, rp, gpoint);
|
|||
|
|
*gx = gpoint[0];
|
|||
|
|
*gy = gpoint[1];
|
|||
|
|
/* Put z back down again to maintain the normal convention for gz */
|
|||
|
|
*gz = -gpoint[2];
|
|||
|
|
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Calculates the potential caused by a prism. */
|
|||
|
|
double prism_pot_sph(PRISM prism, double lonp, double latp, double rp)
|
|||
|
|
{
|
|||
|
|
double x = 0, y = 0, z = 0, res;
|
|||
|
|
|
|||
|
|
global2local(lonp, latp, rp, prism, &x, &y, &z);
|
|||
|
|
res = prism_pot(prism, x, y, z);
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|