-
Notifications
You must be signed in to change notification settings - Fork 106
/
NormalGravity.f95
65 lines (48 loc) · 2.13 KB
/
NormalGravity.f95
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
function NormalGravity(geocentric_lat, GM, omega, a, b)
!------------------------------------------------------------------------------
!
! Compute the total predicted gravity normal to the ellipsoid at a given
! GEOCENTRIC latitude (input in DEGREES), using Somigliana's formula.
! This is taken from Physical Geodesy (Hofmann-Wellenhof and Moritz,
! sec. ed.), sections 2.7 and 2.8.
!
! Copyright (c) 2005-2019, SHTOOLS
! All rights reserved.
!
!------------------------------------------------------------------------------
use ftypes
implicit none
real(dp) :: NormalGravity
real(dp), intent(in) :: geocentric_lat, gm, omega, a, b
real(dp) :: geodetic_lat, pi, ga, gb, m, ep, bigE, q0, q0p
if (a < b) then
print*, "Warning --- NormalGravity"
print*, "The semimajor axis A should be greater than the " // &
"semiminor axis B."
end if
if (a == b .and. omega == 0.0_dp) then
NormalGravity = GM / a**2
else if (a == b .and. omega /= 0.0_dp) then
print*, "Warning --- NormalGravity"
print*, "A can not be equal to B when OMEGA is non zero."
print*, "Setting OMEGA equal to zero."
NormalGravity = GM / a**2
else
pi = acos(-1.0_dp)
m = omega**2 * a**2 * b / GM
bigE = sqrt(a**2 - b**2) ! linear eccentricity
ep = bigE / b ! second eccentricity
q0 = 0.50_dp * ((1.0_dp + 3.0_dp * (b / bigE)**2) * atan(bigE / b) &
- 3.0_dp * b / bigE)
q0p = 3.0_dp * (1.0_dp + (b / bigE)**2) &
* (1.0_dp - b / bigE * atan(bigE / b) ) - 1.0_dp
! gravity on equator
ga = GM / (a*b) * (1.0_dp - m - m * ep * q0p / 6.0_dp / q0)
gb = GM / a**2 * (1.0_dp + m * ep * q0p / 3.0_dp / q0) ! gravity at poles
geodetic_lat = atan((a / b)**2 * tan(geocentric_lat * pi / 180.0_dp))
NormalGravity = a * ga * cos(geodetic_lat)**2 + b * gb &
* sin(geodetic_lat)**2
NormalGravity = NormalGravity / sqrt(a**2 * cos(geodetic_lat)**2 &
+ b**2 * sin(geodetic_lat)**2 )
end if
end function NormalGravity