Note
The following appendices present code written in the R language for calculation of Earth’s radius, effective g and theta. I am happy to provide the code in a text file to anyone who wishes to use it. To obtain the code please e-mail me at: sigma@outlook.co.nz
APPENDIX 1
R code for calculation of Earth’s radius at a given latitude, variation in g by latitude, local circumference and local speed due to Earth’s rotation
Instructions: Create a folder to store the pdf graphs of Appendix 4 and navigate your R working directory to that folder. Copy and paste all of the code below into the R workspace together.
rm(list=ls()) # Clears the R workspace
R1 <- 6378.137 # Established radius in km at the equator
R2 <- 6356.752 # Established radius in km at the poles
phideg <- 40 # I have used 40 degrees here, but enter any value of latitude as an angle in degrees, simply by replacing the 40 in this line of code
phirad <- phideg * pi/180 # Automatic conversion from degrees to radians, as required by R
NUMERATOR <- (R1**2 * cos(phirad))**2 + (R2**2 * sin(phirad))**2
DENOMINATOR <- (R1 * cos(phirad))**2 + (R2 * sin(phirad))**2
RADIUS <- (NUMERATOR / DENOMINATOR)**0.5
RADIUS
LOCALCIRCUMFERENCE = 2 * pi * RADIUS * cos(phirad) # Calculates the circumference in km at the given latitude
LOCALCIRCUMFERENCE
LOCALVELOCITY = LOCALCIRCUMFERENCE / 24 # Speed in km per hour
LOCALVELOCITY
LOCALVELOCITYMS = LOCALVELOCITY*1000/3600 # Speed in m per second
LOCALVELOCITYMS
APPENDIX 2
R code for calculation of the magnitude of radial g and effective g at a given latitude. Also calculates radial and effective weight for a given body mass in kg.
Instructions: Copy and paste all of the code below into the R workspace together.
rm(list=ls()) # Clears R workspace
# RADIUS <- 6378.137 # Reminder of radius at equator.
# RADIUS <- 6356.752 # Reminder of radius at poles.
W <- 7.292e-5 # Angular velocity of Earth in radians per second
RADIUS <- 6369.345 # Radius at 40 latitude. For other latitudes calculate the radius using the code of Appendix 1 and replace the radius value above
RADIUS <- RADIUS * 1000 # Convert radius from km to m. g is measured in m/sec2 so radii must also be in m
phideg <- 40 # Enter any value of latitude as an angle in degrees by replacing the 40 in this line of code
phirad <- phideg * pi/180 # Automatic conversion to radians, as required in R
# Now we calculate g for the given latitude and corresponding Earth radius
M <- 5.972e24 # Mass of Earth
G <- 6.674e-11 # Gravitational constant
g <- G * M / RADIUS **2 # Using Newton's Law as an approximation. In fact, the increase in g due to oblateness is not as large as expected from the inverse squared law, and is approximately 0.33 of the observed value. The observed value is larger because the Earth’s density increases toward the center. Thus, our calculated values for local gravitational acceleration slightly overestimate the true values. In this simple model g is overestimated by approximately 0.3%.
g # Prints our estimate of local g on the R console
GEFFRAD <- g - W**2 * RADIUS * cos(phirad)**2 # Calculate radial component of g – normal component of centrifugal. This is the resultant acceleration we use to calculate weight on the solid earth
GEFFTANG <- W**2 * RADIUS * cos(phirad) * sin(phirad) # Calculate tangential component of centrifugal force. On earth, this component is cancelled out by friction but is useful for oceans etc.
GEFF <- (GEFFRAD**2 + GEFFTANG**2)**0.5 # Calculate vector sum for total effective g
GEFFRAD # Prints the magnitude of radial component of g for the given angle on the R console
GEFFTANG # Prints the magnitude of the tangential component of acceleration for the given angle on the R console
GEFF # Prints the magnitude of effective g for the given latitude on the R console
# Now we calculate the weight of a body at the chosen latitude
MASS <- 70 # Enter a body mass in kg – here it’s 70 kg but you can enter your own value
WEIGHT <- MASS * GEFFRAD # Calculates radial weight using the radial component only. OK for terrestrial measurements because friction cancels out the tangential component
WEIGHT # Prints the weight found from the radial component
WEIGHTEFF <- MASS * GEFF # Calculates ‘weight’ using total effective g (OK for ocean water etc)
WEIGHTEFF # Prints the weight using effective g. Because the tangential component is so small, it is essentially the same as that found from the radial component only.
APPENDIX 3
R code for calculation of theta, the angle between g and effective g at a given latitude
The angle theta between g and effective g is given by the inverse tan of the ratio of the tangential and radial components of geff. R uses the atan() function to calculate the inverse tan. The angle is usually very small but reaches nearly 0.1 degree around latitude 45.
rm(list=ls()) # Clears R workspace
# RADIUS <- 6378.137 # Reminder of radius at equator.
# RADIUS <- 6356.752 # Reminder of radius at poles.
W <- 7.292e-5 # Angular velocity of Earth in radians per second
RADIUS <- 6369.345 # Radius at latitude 40 found using the code of Appendix 1 and used in the code of Appendix 2. Again, for any other latitude use your calculated value from the code of Appendix 1
RADIUS <- RADIUS * 1000 # Converting to m
phideg <- 40 # Enter any value of latitude as an angle in degrees by replacing the 40 in this line of code
phirad <- phideg * pi/180 # Automatic conversion to radians, as required in R
# Now we calculate g for the given latitude and corresponding Earth radius
M <- 5.972e24 # Mass of Earth
G <- 6.674e-11 # Gravitational constant
g <- G * M / RADIUS **2 # Using Newton's Law as an approximation got purely gravitational acceleration. In fact, the increase in g due to oblateness is not as large as expected from the inverse squared law, and is approximately 0.33 of the observed value. The observed value is larger because the Earth’s density increases toward the center. Thus, our calculated values for local gravitational acceleration slightly overestimate the true values. In this simple model g is overestimated by approximately 0.3%.
g # Prints our estimate of local g on the R console
GEFFRAD <- g - W**2 * RADIUS * cos(phirad)**2 # Calculate radial component of g - normal component of centrifugal force. This is the resultant acceleration we use to calculate weight on the solid earth
GEFFTANG <- W**2 * RADIUS * cos(phirad) * sin(phirad) # Calculate tangential component of centrifugal force. On earth, this component is cancelled out by friction but is useful for oceans etc.
thetarad <- atan(GEFFTANG / GEFFRAD)
thetarad # Prints the angle in radians on the R console
theta <- thetarad*180/pi # Convert from radians back to degrees
theta # Prints the angle between g and effective g in degrees on the R console
APPENDIX 4
R code for modelling and plotting Earth’s radius, effective g and angle theta by latitude
Remember to create a folder to store the pdf graphs of Appendix 4 and navigate your R working directory to that folder. Copy and paste all of the code below into the R workspace together.
# Cut and paste all the following g code at one time and enter on the R console.
# First we model the variation of Earth’s radius by latitude due to oblateness.
rm(list=ls()) # Clears the R workspace
# Set up arrays of key variables for calculations
phideg <- array(100) # Latitudes in degrees
phirad <- array(100) # Latitudes in radians
NUMERATOR <- array(100) # Used in various calculations
DENOMINATOR <- array(100) # Used in various calculations
RADIUS <- array(100) # Earth radius across latitudes
g <- array(100) # Gravitational acceleration across latitudes
GEFFRAD <- array(100) # Radial component of geff
GEFFTANG <- array(100) # Tangential component of geff
GEFF <- array(100) # Effective g
THETARAD <- array(100) # Theta in radians
THETA <- array(100) # Theta in degrees
W <- 7.292e-5 # Angular velocity of Earth in radians per second
R1 <- 6378.137 # Established radius in km at the equator
R2 <- 6356.752 # Established radius in km at the poles
# We use a loop to evaluate Earth’s radius across the range of latitudes from one degree to the pole. The equator must be treated separately.
for (k in 1:90) {
phideg[k] <- k
phirad[k] <- phideg[k] * pi/180 # Automatic conversion from degrees to radians, as required by R
NUMERATOR[k] <- (R1**2 * cos(phirad[k]))**2 + (R2**2 * sin(phirad[k]))**2
DENOMINATOR[k] <- (R1 * cos(phirad[k]))**2 + (R2 * sin(phirad[k]))**2
RADIUS[k] <- (NUMERATOR[k] / DENOMINATOR[k])**0.5
}
RADIUS # Prints out the radii for 1 DEGREE to 90 DEGREES.
# WE NOW CALCULATE THE RADIUS FOR THE EQUATOR (CANNOT BE DONE IN THE ABOVE LOOP)
NUMERATOREQUATOR <- (R1**2 * cos(0))**2 + (R2**2 * sin(0))**2
DENOMINATOREQUATOR <- (R1 * cos(0))**2 + (R2 * sin(0))**2
RADIUSEQUATOR <- (NUMERATOREQUATOR / DENOMINATOREQUATOR)**0.5
RADIUSEQUATOR
RADIUS <- c(RADIUSEQUATOR, RADIUS) # NOW INCLUDING THE EQUATORIAL VALUE
RADIUS
# CREATE A PNG GRAPH OF EARTH'S RADIUS WITH LATITUDE
png("RADIUS_lat.png")
plot(1:91, RADIUS, pch=16, xaxt="n", type = "l", lwd=2, col = "red", xlab = expression(paste("Latitude (",phi, ")")), ylab = "Earth's Radius (km)", main= "Variation of Earth's Radius by Latitude", cex.lab=1.3, xlim= c(0,95), ylim=c(6355, 6380))
axis(1, at=c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90),las=1)
dev.off()
# Now we use a second loop to calculate radial g (GEFFRAD), the tangential component of the centrifugal acceleration (GEFFTANG) and effective g (GEFF)I use multiple loops for clarity for other users of the code
# Now we calculate g an d geff for across latitudes and corresponding Earth radii
M <- 5.972e24 # Mass of Earth
G <- 6.674e-11 # Gravitational constant
# WE MUST ELIMINATE THE FIRST ELEMENT FOR THE EQUATOR SO WE CAN USE A LOOP FOR 1 DEGREE TO 90 DEGREES
RADIUS <- RADIUS[-c(1)] # Removes the equatorial radius
for (k in 1:90) {
RADIUS[k] <- RADIUS[k]*1000 # Convert km to m
g[k] <- G * M / RADIUS[k]**2 # Purely gravitational acceleration
GEFFRAD[k] <- g[k] - W**2 * RADIUS[k] * cos(phirad[k])**2 # Calculates radial component of g – normal component of centrifugal
GEFFTANG[k] <- W**2 * RADIUS[k] * cos(phirad[k]) * sin(phirad[k]) # Calculate tangential component of centrifugal force
GEFF[k] <- (GEFFRAD[k]**2 + GEFFTANG[k]**2)**0.5 # Calculate vector sum for total effective g
}
g # Prints the magnitude of g on the R console
GEFFRAD # Prints the magnitude of the radial component of acceleration on the R console
GEFFTANG # Prints the magnitude of the tangential component of acceleration on the R console
GEFF # Prints the magnitude of effective g on the R console
# CALCULATE EFFECTIVE G FOR THE EQUATOR (CANNOT BE DONE WITHIN THE ABOVE LOOP WHERE k RUNS FROM 1 TO 90)
RADIUSEQUATOR <- R1*1000
geq <- G * M / (RADIUSEQUATOR)**2 # Estimate of g at the equator
GEFFRADEQUATOR <- geq - W**2 * RADIUSEQUATOR * cos( 0 )**2 # Calculates radial component of g – normal component of centrifugal (FOR EQUATOR)
GEFFTANGEQUATOR <- W**2 * RADIUSEQUATOR * cos( 0 ) * sin( 0 ) # Calculate tangential component of centrifugal force
GEFFEQUATOR <- (GEFFRADEQUATOR **2 + GEFFTANGEQUATOR **2)**0.5 # Calculate vector sum for total effective g (FOR EQUATOR)
# NOW INCLUDE THE EQUATORIAL VALUES
g <- c(geq, g)
GEFFRAD <- c(GEFFRADEQUATOR, GEFFRAD)
GEFFTANG <- c(RADIUSEQUATOR, GEFFTANG)
GEFF <- c(GEFFEQUATOR, GEFF)
GEFF
# CREATE A PNG GRAPH OF PURELY GRAVITATIONAL ACCELERATION g WITH LATITUDE
png("G_lat.png")
plot(1:91, g, pch=16, xaxt="n", yaxt="n", col = "blue", type = "l", lwd=2, xlab = expression(paste("Latitude (",phi, ")")), ylab= "", main= "Variation of Purely Gravitational Acceleration g by Latitude")
axis(1, at=c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90),las=1)
axis(2, at=c(9.76, 9.77, 9.78, 9.79, 9.80, 9.81, 9.82, 9.83, 9.84, 9.85, 9.86, 9.87),las=1)
title(ylab = expression(paste("Purely Gravitational Acceleration g m/", s^2,")")), line= 2.8, cex.lab=1.0, cex.axis = 0.8)
dev.off()
# CREATE A PNG GRAPH OF RADIAL AND EFFECTIVE g WITH LATITUDE (essentially, effective g is the same as the radial component)
png("GEFF_lat.png")
plot(1:91, GEFF, pch=16, xaxt="n", yaxt="n", col = "blue", type = "l", lwd=2, xlab = expression(paste("Latitude (",phi, ")")), ylab= "", main= "Variation of Radial and Effective g by Latitude")
axis(1, at=c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90),las=1)
axis(2, at=c(9.76, 9.77, 9.78, 9.79, 9.80, 9.81, 9.82, 9.83, 9.84, 9.85, 9.86, 9.87),las=1)
title(ylab = expression(paste("Radial and Effective g (m/", s^2,")")), line= 2.8, cex.lab=1.0, cex.axis = 0.8)
dev.off()
# CALCULATE THE VARIATION OF THETA WITH LATITUDE USING ANOTHER LOOP FOR CLARITY
# Remove the equatorial value of GEFF so we can loop from 1 degree to 90 degrees
GEFF <- GEFF[-c(1)]
for (k in 1:90) {
THETARAD[k] <- atan((W**2 * RADIUS[k]*cos(phirad[k])* sin(phirad[k])) / ( GEFF[k] - W**2 * RADIUS[k]*cos(phirad[k])))
THETARAD # Prints the angle in radians on the R console
THETA[k] <- THETARAD[k]*180/pi # Convert from radians back to degrees
}
THETA # Prints the magnitudes of theta on the R console
# CALCULATE THETA FOR THE EQUATOR (CANNOT BE DONE IN THE ABOVE LOOP). CALCULATION IN FULL, THOUGH IT MUST BE ZERO BECAUSE OF THE SINE FUNCTION TAKING THE VALUE ZERO AT ZERO LATITUDE
THETAEQUATOR <- atan((W**2 * RADIUSEQUATOR * cos( 0 )* sin( 0 )) / (GEFFEQUATOR - W**2 * RADIUSEQUATOR *cos( 0 )))
# THETARAD # Prints the angle in radians on the R console
THETAEQUATOR <- THETAEQUATOR *180/pi # Convert from radians back to degrees
THETA <- c(THETAEQUATOR, THETA)
THETA
# CREATE A PNG GRAPH OF THETA WITH LATITUDE
png("THETA_lat.png")
plot(1:91, THETA, pch=16, xaxt="n", yaxt="n", col = "darkgreen", type = "l", lwd=2, xlab = expression(paste("Theta (",theta, ")")),
ylab= "", main= "Variation of Theta for Effective g by Latitude")
axis(1, at=c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90),las=1)
axis(2, at=c(0.0, 0.02, 0.04, 0.06, 0.08, 0.1),las=1)
title(ylab = expression(paste("Theta for Effective g in degrees")), line= 2.8, cex.lab=1.0, cex.axis = 0.8)
dev.off()