253 lines
8.1 KiB
C
253 lines
8.1 KiB
C
/*
|
|
This module contains a set of functions that calculate the gravitational
|
|
potential and its first and second derivatives for the sphere in spherical
|
|
coordinates.
|
|
|
|
The position of the sphere and computation point are in spherical coordinates.
|
|
|
|
The derivatives of the potential are made with respect to the local coordinate
|
|
system x->North, y->East, z->out. So it would be normal for a sphere of positive
|
|
density to have negative gz
|
|
|
|
References
|
|
----------
|
|
|
|
* Grombein, T.; Seitz, K.; Heck, B. (2010): Untersuchungen zur effizienten
|
|
Berechnung topographischer Effekte auf den Gradiententensor am Fallbeispiel der
|
|
Satellitengradiometriemission GOCE.
|
|
KIT Scientific Reports 7547, ISBN 978-3-86644-510-9, KIT Scientific Publishing,
|
|
Karlsruhe, Germany.
|
|
*/
|
|
|
|
|
|
#include <math.h>
|
|
#include "geometry.h"
|
|
#include "constants.h"
|
|
#include "grav_sphere.h"
|
|
|
|
|
|
/* Calculates the potential caused by a sphere */
|
|
double sphere_pot(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*(
|
|
sinlatp*sinlatc + coslatp*coslatc*coslon);
|
|
|
|
return G*mass/sqrt(l_sqr);
|
|
}
|
|
|
|
|
|
/* Calculates the gx component of gravitational attraction caused by a sphere */
|
|
double sphere_gx(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., kphi, coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*(
|
|
sinlatp*sinlatc + coslatp*coslatc*coslon);
|
|
|
|
kphi = coslatp*sinlatc - sinlatp*coslatc*coslon;
|
|
|
|
return G*SI2MGAL*mass*(sphere.rc*kphi)/pow(l_sqr, 1.5);
|
|
}
|
|
|
|
|
|
/* Calculates the gy component of gravitational attraction caused by a sphere */
|
|
double sphere_gy(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., cospsi, coslatc, kern;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatc = cos(d2r*sphere.latc);
|
|
|
|
cospsi = sin(d2r*latp)*sin(d2r*sphere.latc) + cos(d2r*latp)*coslatc*
|
|
cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
kern = (sphere.rc*coslatc*sin(d2r*(sphere.lonc - lonp)))/pow(l_sqr, 1.5);
|
|
|
|
return G*SI2MGAL*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculates the gz component of gravitational attraction caused by a sphere */
|
|
double sphere_gz(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., cospsi;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
cospsi = sin(d2r*latp)*sin(d2r*sphere.latc) + cos(d2r*latp)*
|
|
cos(d2r*sphere.latc)*cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
return G*SI2MGAL*mass*(sphere.rc*cospsi - rp)/pow(l_sqr, 1.5);
|
|
}
|
|
|
|
|
|
/* Calculate the xx component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gxx(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., kphi, coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon, kern;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*(sinlatp*sinlatc +
|
|
coslatp*coslatc*coslon);
|
|
|
|
kphi = coslatp*sinlatc - sinlatp*coslatc*coslon;
|
|
|
|
kern = (3*sphere.rc*kphi*sphere.rc*kphi - l_sqr)/pow(l_sqr, 2.5);
|
|
|
|
return G*SI2EOTVOS*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculate the xy component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gxy(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., kphi, coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon, kern;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*(sinlatp*sinlatc +
|
|
coslatp*coslatc*coslon);
|
|
|
|
kphi = coslatp*sinlatc - sinlatp*coslatc*coslon;
|
|
|
|
kern = (3*sphere.rc*sphere.rc*kphi*coslatp*sin(d2r*(sphere.lonc - lonp)))/
|
|
pow(l_sqr, 2.5);
|
|
|
|
return G*SI2EOTVOS*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculate the xz component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gxz(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., kphi, coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon, kern, cospsi;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
|
|
cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon;
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
kphi = coslatp*sinlatc - sinlatp*coslatc*coslon;
|
|
|
|
kern = 3*sphere.rc*kphi*(sphere.rc*cospsi - rp)/pow(l_sqr, 2.5);
|
|
|
|
return G*SI2EOTVOS*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculate the yy component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gyy(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon, sinlon, kern, cospsi;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
sinlon = sin(d2r*(sphere.lonc - lonp));
|
|
|
|
cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon;
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
kern = (3*sphere.rc*sphere.rc*coslatc*coslatc*sinlon*sinlon - l_sqr)/
|
|
pow(l_sqr, 2.5);
|
|
|
|
return G*SI2EOTVOS*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculate the yz component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gyz(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., coslatp, coslatc, sinlatp, sinlatc,
|
|
coslon, sinlon, kern, cospsi;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
coslatp = cos(d2r*latp);
|
|
coslatc = cos(d2r*sphere.latc);
|
|
sinlatp = sin(d2r*latp);
|
|
sinlatc = sin(d2r*sphere.latc);
|
|
coslon = cos(d2r*(lonp - sphere.lonc));
|
|
sinlon = sin(d2r*(sphere.lonc - lonp));
|
|
|
|
cospsi = sinlatp*sinlatc + coslatp*coslatc*coslon;
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
kern = 3*sphere.rc*coslatc*sinlon*(sphere.rc*cospsi - rp)/pow(l_sqr, 2.5);
|
|
|
|
return G*SI2EOTVOS*mass*kern;
|
|
}
|
|
|
|
|
|
/* Calculate the zz component of gravity gradient tensor cause by a sphere */
|
|
double sphere_gzz(SPHERE sphere, double lonp, double latp, double rp)
|
|
{
|
|
double mass, l_sqr, d2r = PI/180., deltaz, cospsi;
|
|
|
|
mass = (double)(sphere.density*4.*PI*sphere.r*sphere.r*sphere.r)/3.;
|
|
|
|
cospsi = sin(d2r*latp)*sin(d2r*sphere.latc) + cos(d2r*latp)*
|
|
cos(d2r*sphere.latc)*cos(d2r*(lonp - sphere.lonc));
|
|
|
|
l_sqr = rp*rp + sphere.rc*sphere.rc - 2*rp*sphere.rc*cospsi;
|
|
|
|
deltaz = sphere.rc*cospsi - rp;
|
|
|
|
return G*SI2EOTVOS*mass*(3*deltaz*deltaz - l_sqr)/pow(l_sqr, 2.5);
|
|
}
|