g_eff(θ) = g_grav - g_centrifugal g_eff(θ) = GM/r² - ω² * r * cos²(θ) [closed]

import numpy as np
import matplotlib.pyplot as plt

Constants

G = 6.67e-11
M_earth = 5.972e24
r = 6.371e6
omega = 7.292115e-5
c = 2.99792458e8
J_sun = 1.9e41
J_earth = 8.0e37

theta = np.linspace(0, np.pi/2, 200)
lat_deg = np.degrees(theta)

1. Gypsy base model

g_grav = G * M_earth / r2
g_cent = omega2 r np.cos(theta)**2
g_gypsy = g_grav - g_cent

2. PREM/WGS84

g_prem = 9.7803253359 * (1 + 0.00193185265241np.sin(theta)2) / np.sqrt(1 - 0.00669437999013np.sin(theta)2)

3. Measured

lat_meas = np.array([0, 30, 45, 60, 90])
g_meas = np.array([9.7803, 9.7932, 9.8062, 9.8192, 9.8322])

4. CHERRY: Gypsy Correction Term - Sun+Earth frame dragging wobble

Max at equator, zero at poles, yearly sine wave

delta_g_gypsy = (6G / (c2 r4)) J_sun J_earth np.sin(theta) 1.0
g_gypsy_cherry = g_gypsy + delta_g_gypsy

Plot

plt.figure(figsize=(9,5.5))
plt.plot(lat_deg, g_gypsy, label="Gypsy Base: GM/r² - ω²r cos²θ", linewidth=2)
plt.plot(lat_deg, g_gypsy_cherry, label="Gypsy + Cherry: + Δg_LT Sun×Earth", linewidth=2, alpha=0.8)
plt.plot(lat_deg, g_prem, '--', label="PREM/WGS84 Standard", linewidth=2)
plt.scatter(lat_meas, g_meas, color='red', s=80, zorder=5, label="Measured IGSN71")

plt.xlabel("Latitude θ [degrees]")
plt.ylabel("Effective Gravity g [m/s²]")
plt.title("Gravity vs Latitude: Base Model + Cherry vs PREM vs Data")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print(f"Cherry magnitude at equator: {delta_g_gypsy[0]*1e12:.2f} pico m/s²")
print(f"Cherry magnitude at 45°: {delta_g_gypsy[100]*1e12:.2f} pico m/s²")
print(f"Cherry magnitude at pole: {delta_g_gypsy[-1]*1e12:.2f} pico m/s²")