470 lines
12 KiB
C
470 lines
12 KiB
C
|
/*
|
|||
|
|
Functions that calculate the gravitational potential and its first and second
|
|||
|
|
derivatives for the rectangular prism using the formulas in Nagy et al. (2000).
|
|||
|
|
|
|||
|
|
The coordinate system used is that of the article, ie:
|
|||
|
|
|
|||
|
|
x -> North y -> East z -> Down
|
|||
|
|
|
|||
|
|
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 <stdlib.h>
|
|||
|
|
#include "geometry.h"
|
|||
|
|
#include "constants.h"
|
|||
|
|
#include "grav_prism.h"
|
|||
|
|
|
|||
|
|
double safe_atan2(double y, double x)
|
|||
|
|
{
|
|||
|
|
if(y == 0)
|
|||
|
|
{
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
if((y > 0) && (x < 0))
|
|||
|
|
{
|
|||
|
|
return atan2(y, x) - PI;
|
|||
|
|
}
|
|||
|
|
if((y < 0) && (x < 0))
|
|||
|
|
{
|
|||
|
|
return atan2(y, x) + PI;
|
|||
|
|
}
|
|||
|
|
return atan2(y, x);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
double safe_log(double x)
|
|||
|
|
{
|
|||
|
|
if(x == 0)
|
|||
|
|
{
|
|||
|
|
return 0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return log(x);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Calculates the potential cause by a prism. */
|
|||
|
|
double prism_pot(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
kernel = (x[i]*y[j]*safe_log(z[k] + r)
|
|||
|
|
+ y[j]*z[k]*safe_log(x[i] + r)
|
|||
|
|
+ x[i]*z[k]*safe_log(y[j] + r)
|
|||
|
|
- 0.5*x[i]*x[i]*safe_atan2(z[k]*y[j], x[i]*r)
|
|||
|
|
- 0.5*y[j]*y[j]*safe_atan2(z[k]*x[i], y[j]*r)
|
|||
|
|
- 0.5*z[k]*z[k]*safe_atan2(x[i]*y[j], z[k]*r));
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density */
|
|||
|
|
res *= G*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the x component of gravitational attraction cause by a prism. */
|
|||
|
|
double prism_gx(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
|
|||
|
|
kernel = -(y[j]*safe_log(z[k] + r) + z[k]*safe_log(y[j] + r)
|
|||
|
|
- x[i]*safe_atan2(z[k]*y[j], x[i]*r));
|
|||
|
|
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to mGal units */
|
|||
|
|
res *= G*SI2MGAL*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the y component of gravitational attraction cause by a prism. */
|
|||
|
|
double prism_gy(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
|
|||
|
|
kernel = -(z[k]*safe_log(x[i] + r) + x[i]*safe_log(z[k] + r)
|
|||
|
|
- y[j]*safe_atan2(z[k]*x[i], y[j]*r));
|
|||
|
|
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to mGal units */
|
|||
|
|
res *= G*SI2MGAL*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the z component of gravitational attraction cause by a prism. */
|
|||
|
|
double prism_gz(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
|
|||
|
|
kernel = -(x[i]*safe_log(y[j] + r) + y[j]*safe_log(x[i] + r)
|
|||
|
|
- z[k]*safe_atan2(x[i]*y[j], z[k]*r));
|
|||
|
|
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to mGal units */
|
|||
|
|
res *= G*SI2MGAL*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gxx gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gxx(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
kernel = -safe_atan2(z[k]*y[j], x[i]*r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Calculates the gxy gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gxy(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r, xtmp, ytmp;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
if(x[i] == 0 && y[j] == 0 && z[k] < 0)
|
|||
|
|
{
|
|||
|
|
xtmp = 0.0001*(prism.x2 - prism.x1);
|
|||
|
|
ytmp = 0.0001*(prism.y2 - prism.y1);
|
|||
|
|
r = sqrt(xtmp*xtmp + ytmp*ytmp + z[k]*z[k]);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
}
|
|||
|
|
kernel = safe_log(z[k] + r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gxz gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gxz(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r, xtmp, ztmp;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
if(x[i] == 0 && z[k] == 0 && y[j] < 0)
|
|||
|
|
{
|
|||
|
|
xtmp = 0.0001*(prism.x2 - prism.x1);
|
|||
|
|
ztmp = 0.0001*(prism.z2 - prism.z1);
|
|||
|
|
r = sqrt(xtmp*xtmp + ztmp*ztmp + y[j]*y[j]);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
}
|
|||
|
|
kernel = safe_log(y[j] + r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gyy gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gyy(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
kernel = -safe_atan2(z[k]*x[i], y[j]*r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gyz gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gyz(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r, ytmp, ztmp;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
if(z[k] == 0 && y[j] == 0 && x[i] < 0)
|
|||
|
|
{
|
|||
|
|
ytmp = 0.0001*(prism.y2 - prism.y1);
|
|||
|
|
ztmp = 0.0001*(prism.z2 - prism.z1);
|
|||
|
|
r = sqrt(ztmp*ztmp + ytmp*ytmp + x[i]*x[i]);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
}
|
|||
|
|
kernel = safe_log(x[i] + r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
/* Calculates the gzz gravity gradient tensor component cause by a prism. */
|
|||
|
|
double prism_gzz(PRISM prism, double xp, double yp, double zp)
|
|||
|
|
{
|
|||
|
|
double x[2], y[2], z[2], kernel, res, r;
|
|||
|
|
register int i, j, k;
|
|||
|
|
|
|||
|
|
/* First thing to do is make P the origin of the coordinate system */
|
|||
|
|
x[0] = prism.x2 - xp;
|
|||
|
|
x[1] = prism.x1 - xp;
|
|||
|
|
y[0] = prism.y2 - yp;
|
|||
|
|
y[1] = prism.y1 - yp;
|
|||
|
|
z[0] = prism.z2 - zp;
|
|||
|
|
z[1] = prism.z1 - zp;
|
|||
|
|
|
|||
|
|
res = 0;
|
|||
|
|
|
|||
|
|
/* Evaluate the integration limits */
|
|||
|
|
for(k=0; k<=1; k++)
|
|||
|
|
{
|
|||
|
|
for(j=0; j<=1; j++)
|
|||
|
|
{
|
|||
|
|
for(i=0; i<=1; i++)
|
|||
|
|
{
|
|||
|
|
r = sqrt(x[i]*x[i] + y[j]*y[j] + z[k]*z[k]);
|
|||
|
|
kernel = -safe_atan2(x[i]*y[j], z[k]*r);
|
|||
|
|
res += pow(-1, i + j + k)*kernel;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
/* Now all that is left is to multiply res by the gravitational constant and
|
|||
|
|
density and convert it to Eotvos units */
|
|||
|
|
res *= G*SI2EOTVOS*prism.density;
|
|||
|
|
|
|||
|
|
return res;
|
|||
|
|
}
|