Simple Earth - Reference system

import numpy as np

class grs80:
    """ Defined and derived constans associated with GRS80 as described in
        Geodetic Reference System 1980, H. Moritz (2000)
    """
    a = 6378137.000       # [m] semimajor axis
    GM = 3986005.0e8      # [m^3*s^-2] geocentric gravitational constant
    J2 = 108263.0e-8      # dynamic form factor
    omega = 7292115e-11   # [rad*s^-1] angular velocity
    b = 6356752.3141      # [m] semiminor axis
    R3 = 6371000.7900     # [m] Radius of sphere with same volume as Earth

earthVolume = 4.0/3.0 * np.pi * np.power(grs80.R3, 3)
print("Volume of Earth {:.4e} m^3".format(earthVolume))

earthCylinderHeight = 2 * grs80.b  # Insert the correct calculation
print("Height of EarthC {:.4e} m".format(earthCylinderHeight))

earthCylinderRadius = grs80.R3 # Insert the correct calculation
print("Radius of EarthC {:.4e} m\n".format(earthCylinderRadius))

v_ant = 2.65e16             # [m^3] Volume of Antarctica ice sheet
rho_ice = 915.0             # [kg*m^-3] Density of ice
rho_rock = 5513.0           # [kg*m^-3] Density of rock

antarcticaHeight = v_ant / (np.pi * np.power(grs80.R3, 2))  # Insert the correct calculation
print("Height of Antarctica {:.4f} m".format(antarcticaHeight))

antarcticaWeight = rho_ice * v_ant # Insert the correct calculation
print("Weight of Antarctica {:.4e} kg".format(antarcticaWeight))

earthRockWeight = rho_rock * earthVolume   # Insert the correct calculation
print("Weight of rock Earth {:.4e} kg\n".format(earthRockWeight))

antarcticaCMX = 0  # Determine the correct value
antarcticaCMY = 0  # Determine the correct value
antarcticaCMZ = (1/2) * antarcticaHeight  # Determine the correct value
print("Center of mass for Antarctica:")
print(" X: {:.4f} m".format(antarcticaCMX))
print(" Y: {:.4f} m".format(antarcticaCMY))
print(" Z: {:.4f} m\n".format(antarcticaCMZ))

earthRockCMX = 0  # Determine the correct value
earthRockCMY = 0  # Determine the correct value
earthRockCMZ = (1/2) * earthCylinderHeight  # Determine the correct value
print("Center of mass for the Rock Earth:")
print(" X: {:.4f} m".format(earthRockCMX))
print(" Y: {:.4f} m".format(earthRockCMY))
print(" Z: {:.4f} m\n".format(earthRockCMZ))

# It is obvius that X and Y are zero due to the simetrical shape of the cilynder
earthCMX = ((antarcticaCMX * antarcticaWeight)/antarcticaWeight) +  ((earthRockCMX * earthRockWeight)/earthRockWeight) # Insert the correct calculation
earthCMY = ((antarcticaCMY * antarcticaWeight)/antarcticaWeight) +  ((earthRockCMY * earthRockWeight)/earthRockWeight)  # Insert the correct calculation
earthCMZ = ((antarcticaCMZ * antarcticaWeight)/antarcticaWeight) +  ((earthRockCMZ * earthRockWeight)/earthRockWeight)  # Insert the correct calculation
print("Center of mass for the entire Earth:")
print(" X: {:.4f} m".format(earthCMX))
print(" Y: {:.4f} m".format(earthCMY))
print(" Z: {:.4f} m\n".format(earthCMZ))

# North Pole is placed in the X and Y 0
# For Z axis it is necessary to change the coordinates, since the heigh of the cylinder is expressed from the beneath point

nPX = 0  # Insert the correct calculation
nPY = 0  # Insert the correct calculation
nPZ = (earthCylinderHeight + antarcticaHeight) - earthCMZ  # Insert the correct calculation

print("Coordinate of the North Pole:")
print(" X: {:.4f} m".format(nPX))
print(" Y: {:.4f} m".format(nPY))
print(" Z: {:.4f} m\n".format(nPZ))

# After melted ¿Does the origo changes?
# Después de que se funda el hielo el agua hace que cambie el sistema de coordenadas. La gracia del ejercicio es determinar si el centro de coordenadas ha cambiado

ro_ocean = 1000 # kg/m^3

v_ant_pri = (ro_ocean/(rho_ice)) * v_ant # Antartica volume after melting
print("Antarctica Volume: {}".format(v_ant_pri))

v_sys = earthVolume + v_ant_pri



## If antarctica melts then the cylinder will be bigger, and the reference system will change. 
## The fluid will placed uniformly in the shape of the cylinder, that is... the origo of our system will be placed
## just in the middle of the cylinder.


antarcticaHeight_pri = v_ant_pri / (np.pi * np.power(grs80.R3, 2))  # Insert the correct calculation
print(antarcticaHeight_pri - antarcticaHeight)

# X and Y still 0 because it is a simetrical object and the ice melts homogeniously

nPMeltX = 0  # Insert the correct calculation
nPMeltY = 0  # Insert the correct calculation
nPMeltZ = nPZ + antarcticaHeight_pri # Insert the correct calculation

print("Coordinate of the North Pole without Antarctica:")
print(" X: {:.4f} m".format(nPMeltX))
print(" Y: {:.4f} m".format(nPMeltY))
print(" Z: {:.4f} m\n".format(nPMeltZ))

print("Change in Z: {:.4f} m".format(nPZ - nPMeltZ))

rho_ocean = 1000.0          # [kg*m^-3] Density of ocean water

earthArea = 2 * np.pi * earthCylinderRadius * (earthCylinderRadius + earthCylinderHeight)

oceanArea = 0   # Insert the correct calculation
oceanDepth = v_ant_pri / earthArea  # Insert the correct calculation

print("Ocean depth: {:.4f} m".format(oceanDepth))