Cuánto pesa la Tierra, o cómo medir la masa de nuestro planeta usando R

1 Objetivos

El principal objetivo de esta práctica es obtener el valor de la masa de la Tierra usando datos reales del movimiento de satélites terrestres. Así mismo se estimará la incertidumbre de esta medida para poder concluir si el valor experimental obtenido concuerda con el valor aceptado de la masa de la Tierra. Para el análisis de los datos experimentales se usará el software estadístico R.¹

2 Marco teórico

Comenzaremos la actividad haciendo un repaso de la tercera ley de Kepler y analizando cómo se puede aplicar al movimiento de satélites terrestres.

2.1 ¿Cómo se mide la masa de la Tierra?

Todos nos podemos hacer una idea de cuánto pesan media docena de naranjas. En caso de duda nos acercamos a la frutería más cercana, y si la balanza marca alrededor de un kilo nos quedamos muy satisfechos con el resultado. Sin embargo, sabemos que estamos siendo engañados –aunque este engaño no nos preocupe lo más mínimo–: en realidad la balanza no nos está dando el peso de las naranjas, sino su masa.

Usar una balanza para medir la masa de un objeto es algo muy habitual. Pero si lo que queremos medir es la masa de la Tierra este método resulta, por motivos evidentes, totalmente inadecuado. Entonces, ¿qué podemos hacer? La solución está en recurrir a la influencia gravitatoria que ejerce la Tierra sobre otros cuerpos. En concreto, mediremos la masa de la Tierra a partir del efecto que produce nuestro planeta sobre el movimiento de los satélites que orbitan a su alrededor.

2.2 La tercera ley de Kepler

Utilizando los exactos y completos datos recabados durante años por el astrónomo Tycho Brahe y basándose en conceptos geométricos básicos, Johannes Kepler estableció un modelo sencillo y preciso para describir el movimiento de los planetas alrededor del Sol. En el capítulo 3 de su obra Harmonices mundi. Libri V (La armonía de los mundos. Libro V), publicada en Linz en el año 1619, Kepler enuncia el resultado que hoy conocemos como tercera ley de Kepler.

Johannes Kepler: Harmonices mundi. Libri V (La armonía de los mundos. Libro V). Linz, 1619

Escrita en términos modernos, esta ley afirma que la relación entre el cuadrado del período orbital de un planeta, \(T\), y el cubo de su distancia media al Sol, \(r\), es la misma para todos los planetas que orbitan alrededor de nuestra estrella:

\[ \frac {T^2}{r^3} = cte\] Recurriendo a los posteriores descubrimientos de Newton sobre la fuerza gravitatoria, la expresión de la tercera ley de Kepler tal y como la solemos escribir hoy en día es:

\[ \frac {T^2}{r^3} = \frac{4\pi^2}{G M_S} \]

siendo:

• \(G\): constante de gravitación universal
  
• \(M_S\): masa del Sol
  
• \(T\): período orbital del planeta
  
• \(r\): distancia media del planeta al Sol

2.3 Satélites terrestres

Kepler enunció sus tres famosas leyes del movimiento planetario estudiando las trayectorias de los planetas de nuestro Sistema Solar. Debemos a Newton el descubrimiento de la universalidad de la fuerza gravitatoria que provoca este movimiento, lo que nos permite generalizar las leyes de Kepler para describir el movimiento de los cuerpos que orbitan alrededor de cualquier astro, no importa que sean planetas girando alrededor de una estrella o satélites orbitando un planeta.

Pensemos en la Tierra. Nuestro planeta tiene un único satélite natural, la Luna. Pero, a cambio, hay miles de satélites artificiales dando vueltas a su alrededor: satélites de comunicaciones, meteorológicos, de investigación, de posicionamiento (los famosos GPS)… Pero no se trata únicamente de satélites operativos, sino que orbitando la Tierra hay una ingente cantidad de desperdicios que se han convertido en basura espacial.

En esta práctica estudiaremos el movimiento orbital de varios satélites terrestres. No importa que estemos hablando de la Estación Espacial Internacional o de un trozo de basura del tamaño de una manzana; en ambos casos se mantiene la validez de la tercera ley de Kepler. Eso sí, debemos reescribir la expresión anterior de la tercera ley de Kepler para reflejar que los cuerpos están orbitando alrededor de la Tierra: \[\frac {T^2}{r^3} = \frac{4\pi^2}{G M_T}\]

o, lo que es lo mismo:

\[r^3 = \frac{G M_T}{4\pi^2}T^2\] donde \(G\) es la constante de gravitación universal, \(M_T\) la masa de la Tierra, \(T\) el período orbital del satélite considerado y \(r\) la distancia media del satélite a la Tierra.

2.4 Un apunte sobre las trayectorias elípticas

De acuerdo con la primera ley de Kepler, los planetas describen órbitas elípticas alrededor del Sol, que está situado en uno de los focos de la elipse; por tanto, lo mismo sucederá con los satélites que orbitan la Tierra. Según esto, la distancia entre el satélite y la Tierra varía a medida que el satélite se mueve (si la órbita fuese perfectamente circular esto no sucedería y el satélite estaría siempre a la misma distancia de la Tierra… pero la órbita es una elipse).

Hay un punto en el que el satélite está a la menor distancia posible del planeta; es el perigeo. Y hay otro punto, el apogeo, en el que el satélite está en el punto más alejado de la Tierra. En consecuencia, la distancia entre el planeta y la Tierra varía continuamente entre el valor mínimo en el perigeo, \(r_P\), y el valor máximo en el apogeo, \(r_A\). Entonces, ¿cuál es el valor de la distancia \(r\) que aparece en la tercera ley de Kepler?

Primera ley de Kepler: el satélite describe una órbita elíptica con la Tierra situada en el foco. El dibujo, obviamente, no está a escala (Imagen: Freepik / Beatriz Padín)

La tercera ley de Kepler habla de la distancia media. Es decir, \(r\) es la media de todas las distancias a las que se encuentra el satélite con respecto a la Tierra a lo largo de su movimiento. Esta distancia media no podría ser más fácil de calcular: por las propiedades de las elipses se puede demostrar que su valor coincide con el del semieje mayor de la elipse. Y este semieje mayor es la semisuma de las distancias en el apogeo y en el perigeo. Es decir:

\[r = \frac{r_A + r_P}{2}\]

donde \(r_A\) y \(r_P\) son, respectivamente, la distancia del satélite al centro de la Tierra en el apogeo y en el perigeo.

3 Datos experimentales

Una parte fundamental en el trabajo experimental es la toma de datos. En el experimento que nos ocupa, sin embargo, no vamos a poder tomar medidas directamente (por razones evidentes). A cambio recurriremos a fuentes online que nos proporcionarán todos los datos que necesitamos. Además de recoger los datos y prepararlos para su posterior análisis, haremos un breve análisis con R para conocer un poco más sobre los satélites seleccionados.

3.1 Toma de datos

Como ya hemos dicho, para calcular la masa de la Tierra vamos a aplicar la tercera ley de Kepler al movimiento de satélites terrestres. Pero nos falta algo fundamental para poder aplicar esta ley: los datos. En primer lugar, ¿dónde vamos a conseguir los datos? Y, no menos importante, ¿qué datos necesitamos?

Para contestar la primera pregunta, existen diversas fuentes en internet que suministran datos reales sobre los miles de objetos que orbitan la Tierra. Por ejemplo, la aplicación Satellite Map ofrece un mapa en tres dimensiones donde se pueden seleccionar los satélites y obtener información básica sobre ellos. Otra opción es la página de seguimiento de satélites N2YO.com, en la que se puede hacer una búsqueda por tipo de satélite, país u organización responsable, fecha de lanzamiento y otras muchas opciones.

Una vez escogido un satélite debemos decidir qué información se necesita sobre él. Una de las variables que aparece en la tercera ley de Kepler es el período orbital del satélite. La otra es la distancia media del satélite a la Tierra, para la cual se necesitan las distancias a las que se encuentra el satélite en el afelio y en el perihelio. Pero, ¡cuidado! Los valores que ofrecen estas fuentes de afelio y perihelio se refieren a la altura sobre la superficie de la Tierra; por tanto, cuando procesemos los datos debemos tener en cuenta que a cada uno de estos valores hay que sumarle el radio terrestre.

En resumen, estos son los tres datos imprescindibles que debemos recoger para cada satélite: período orbital, afelio y perihelio. Y, por supuesto, no nos podemos olvidar de las unidades.

Datos del satélite del sistema de navegación europeo Galileo 24, en N2YO.com

Datos de Gorizont 21, satélite de comunicaciones geoestacionario de la Comunidad de Estados Independientes, en Satellite Map

Por ejemplo, hemos obtenido los siguientes valores de perigeo, apogeo y período orbital para diecisiete satélites:

Satélite	hA (km)	hP (km)	T (min)
INTELSAT 2-F3	35850	35714	1435,49
INTELSAT 3-F4	39037	38822	1599,25
GSTAR 4	36140	36088	1452,49
COSMOS 2382	19151	19122	675,71
ARIANE 5	35172	247	619,81
GALILEO 24	23233	23225	844,69
TACSAT 4	12300	463	238,87
ATLAS 75E	20493	1046	372,22
TELSTAR 1	5652	950	157,78
RBSP B	28863	219	501,59
NAVSTAR 18	21244	21164	759,5
MOLNIYA 1-90	25309	958	451,82
SATCOM K1	33248	394	585,81
BREEZE-M	38407	36200	1513,91
FREGAT	35853	34297	1399,55
SYNCOM 1	34568	29003	1236,26
STARLINK-2558	556	554	95,6

Estos son los datos que usaremos para hacer los cálculos en esta actividad.

3.2 Jugando con los datos

Lo habitual es usar una hoja de cálculo para recoger los datos en una tabla como la anterior. Guardamos los datos en un archivo .csv al que hemos llamado satelites.csv. Hemos cambiado el nombre a las columnas ya que es conveniente usar nombres de variables que sean válidos en R; así, la columna satelite almacena el nombre del satélite, apogeoy perigeo son las alturas en kilómetros en el apogeo y en el perigeo, respectivamente, y periodo es el período orbital en minutos.

Una vez creado el archivo de datos debemos leerlo en R, ya que el procesamiento de los datos lo haremos en este entorno. Para importar una tabla de datos se usa la función read.table(). En nuestro caso vamos a guardar los datos en un dataframe llamado datos, para lo que hemos usado la instrucción datos <- read.table("satelites.csv", header=TRUE, sep=","). Los datos tienen ahora el siguiente aspecto:

datos

        satelite apogeo perigeo periodo
1  INTELSAT 2-F3  35850   35714 1435.49
2  INTELSAT 3-F4  39037   38822 1599.25
3        GSTAR 4  36140   36088 1452.49
4    COSMOS 2382  19151   19122  675.71
5       ARIANE 5  35172     247  619.81
6     GALILEO 24  23233   23225  844.69
7       TACSAT 4  12300     463  238.87
8      ATLAS 75E  20493    1046  372.22
9      TELSTAR 1   5652     950  157.78
10        RBSP B  28863     219  501.59
11    NAVSTAR 18  21244   21164  759.50
12  MOLNIYA 1-90  25309     958  451.82
13     SATCOM K1  33248     394  585.81
14      BREEZE-M  38407   36200 1513.91
15        FREGAT  35853   34297 1399.55
16      SYNCOM 1  34568   29003 1236.26
17 STARLINK-2558    556     554   95.60

El dataframe anterior está formado por 17 observaciones (las filas), en cada una de las cuales se han medido cuatro variables (las columnas). Cada elemento de la tabla queda definido por dos índices, de manera que datos[i,j] representa el valor que está en la fila iy en la columna j. Si alguno de los índices está en blanco, eso quiere decir que se seleccionarán todos los valores posibles de ese índice (es decir, todas las filas o todas las columnas).

R nos permite seleccionar datos, ordenar las observaciones, añadir nuevas columnas, etc. de cualquier manera que se nos ocurra. Veamos unos cuantos ejemplos:

1. Accedemos a los datos de una observación concreta (es decir, de un satélite), mostrando los valores de todas las variables. Por ejemplo, la siguiente instrucción devuelve todos los datos correspondientes a la observación 13:

datos[13, ]

    satelite apogeo perigeo periodo
13 SATCOM K1  33248     394  585.81

2. Análogamente, mostramos todos los valores que toma la segunda variable, donde están almacenados los valores del apogeo de los satélites:

datos[ , 2]

 [1] 35850 39037 36140 19151 35172 23233 12300 20493  5652 28863 21244 25309
[13] 33248 38407 35853 34568   556

Como esta segunda variable se llama “apogeo” también podemos usar su nombre:

datos[ , "apogeo"]

 [1] 35850 39037 36140 19151 35172 23233 12300 20493  5652 28863 21244 25309
[13] 33248 38407 35853 34568   556

E incluso es muy habitual usar la siguiente notación para conseguir el mismo resultado:

datos$apogeo

 [1] 35850 39037 36140 19151 35172 23233 12300 20493  5652 28863 21244 25309
[13] 33248 38407 35853 34568   556

3. Mostramos el nombre y el período de todos los satélites:

datos[ , c("satelite", "periodo")]

        satelite periodo
1  INTELSAT 2-F3 1435.49
2  INTELSAT 3-F4 1599.25
3        GSTAR 4 1452.49
4    COSMOS 2382  675.71
5       ARIANE 5  619.81
6     GALILEO 24  844.69
7       TACSAT 4  238.87
8      ATLAS 75E  372.22
9      TELSTAR 1  157.78
10        RBSP B  501.59
11    NAVSTAR 18  759.50
12  MOLNIYA 1-90  451.82
13     SATCOM K1  585.81
14      BREEZE-M 1513.91
15        FREGAT 1399.55
16      SYNCOM 1 1236.26
17 STARLINK-2558   95.60

4. Ordenamos los satélites por período, en orden creciente:

datos[order(datos$periodo), ]

        satelite apogeo perigeo periodo
17 STARLINK-2558    556     554   95.60
9      TELSTAR 1   5652     950  157.78
7       TACSAT 4  12300     463  238.87
8      ATLAS 75E  20493    1046  372.22
12  MOLNIYA 1-90  25309     958  451.82
10        RBSP B  28863     219  501.59
13     SATCOM K1  33248     394  585.81
5       ARIANE 5  35172     247  619.81
4    COSMOS 2382  19151   19122  675.71
11    NAVSTAR 18  21244   21164  759.50
6     GALILEO 24  23233   23225  844.69
16      SYNCOM 1  34568   29003 1236.26
15        FREGAT  35853   34297 1399.55
1  INTELSAT 2-F3  35850   35714 1435.49
3        GSTAR 4  36140   36088 1452.49
14      BREEZE-M  38407   36200 1513.91
2  INTELSAT 3-F4  39037   38822 1599.25

5. Mostramos los satélites por tipo. Los satélites de órbita baja (LEO) son aquellos que orbitan a una altura menor de 2000 km:

datos[datos$apogeo <= 2000, ]

        satelite apogeo perigeo periodo
17 STARLINK-2558    556     554    95.6

Los de órbita alta orbitan a más de 35786 km de altura en su apogeo (algunos de estos tienen órbitas con una excentricidad tan grande que llegan a pasar por zonas LEO, pero aún así se consideran de órbita alta):

datos[datos$apogeo > 35786, ]

        satelite apogeo perigeo periodo
1  INTELSAT 2-F3  35850   35714 1435.49
2  INTELSAT 3-F4  39037   38822 1599.25
3        GSTAR 4  36140   36088 1452.49
14      BREEZE-M  38407   36200 1513.91
15        FREGAT  35853   34297 1399.55

Y, finalmente, los de órbita media (MEO) tienen un apogeo entre 2000 km y 35786 km:

datos[datos$apogeo > 2000 & datos$apogeo <= 35786, ]

       satelite apogeo perigeo periodo
4   COSMOS 2382  19151   19122  675.71
5      ARIANE 5  35172     247  619.81
6    GALILEO 24  23233   23225  844.69
7      TACSAT 4  12300     463  238.87
8     ATLAS 75E  20493    1046  372.22
9     TELSTAR 1   5652     950  157.78
10       RBSP B  28863     219  501.59
11   NAVSTAR 18  21244   21164  759.50
12 MOLNIYA 1-90  25309     958  451.82
13    SATCOM K1  33248     394  585.81
16     SYNCOM 1  34568   29003 1236.26

6. Contamos el número de satélites de órbita baja en nuestra tabla:

sum(datos$apogeo < 2000)

[1] 1

7. Añadimos una columna con el período en horas, escrito con 2 cifras decimales, para hacernos una mejor idea de lo rápido que orbitan los satélites:

horas = round(datos$periodo/60, 2)
cbind(datos, horas)

        satelite apogeo perigeo periodo horas
1  INTELSAT 2-F3  35850   35714 1435.49 23.92
2  INTELSAT 3-F4  39037   38822 1599.25 26.65
3        GSTAR 4  36140   36088 1452.49 24.21
4    COSMOS 2382  19151   19122  675.71 11.26
5       ARIANE 5  35172     247  619.81 10.33
6     GALILEO 24  23233   23225  844.69 14.08
7       TACSAT 4  12300     463  238.87  3.98
8      ATLAS 75E  20493    1046  372.22  6.20
9      TELSTAR 1   5652     950  157.78  2.63
10        RBSP B  28863     219  501.59  8.36
11    NAVSTAR 18  21244   21164  759.50 12.66
12  MOLNIYA 1-90  25309     958  451.82  7.53
13     SATCOM K1  33248     394  585.81  9.76
14      BREEZE-M  38407   36200 1513.91 25.23
15        FREGAT  35853   34297 1399.55 23.33
16      SYNCOM 1  34568   29003 1236.26 20.60
17 STARLINK-2558    556     554   95.60  1.59

Y cualquier otra combinación que se nos pueda ocurrir.

3.3 Preprocesado de los datos

Para estar en disposición de aplicar la tercera ley de Kepler debemos primero hacer una serie de transformaciones en los datos:

• en lo que se refiere a las unidades, necesitamos pasar los datos al Sistema Internacional;
  
• los datos de apogeo y perigeo están medidos desde la superficie de la Tierra, por lo que debemos sumarles el valor del radio de la Tierra;
  
• tenemos que calcular la distancia media entre satélite y la Tierra.

En lo que sigue usaremos los siguientes nombre de variables: Tes el período del satélite en segundos; RT representa el radio de la Tierra; rA y rP son la distancia entre el satélite y el centro de la Tierra en el afelio y en el perihelio, respectivamente; y r es la distancia media entre satélite y la Tierra.

Cálculo de T

Para calcular el período en segundos simplemente debemos multiplicar el valor medido, periodo, por 60:

T = datos$periodo*60

Cálculo de rA y rP

Primero declaramos la constante RT, en la que almacenamos el valor del radio de la Tierra, \(R_T=6,371 \cdot 10^6 \, m\). Para la distancia en el apogeo multiplicamos la variable apogeo por 1000 para expresarla en metros y le sumamos el radio de la Tierra, y lo mismo para el perigeo:

RT = 6.371e6
rA = datos$apogeo*1000 + RT
rP = datos$perigeo*1000 + RT

Cálculo de r

La distancia media es la semisuma de las distancias en el afelio y en el perihelio:

r = (rA + rP)/2

Datos

Una vez obtenidos estos tres valores, podemos mostrar la tabla con las nuevas variables:

cbind(datos, T, rA, rP, r)

        satelite apogeo perigeo periodo       T       rA       rP        r
1  INTELSAT 2-F3  35850   35714 1435.49 86129.4 42221000 42085000 42153000
2  INTELSAT 3-F4  39037   38822 1599.25 95955.0 45408000 45193000 45300500
3        GSTAR 4  36140   36088 1452.49 87149.4 42511000 42459000 42485000
4    COSMOS 2382  19151   19122  675.71 40542.6 25522000 25493000 25507500
5       ARIANE 5  35172     247  619.81 37188.6 41543000  6618000 24080500
6     GALILEO 24  23233   23225  844.69 50681.4 29604000 29596000 29600000
7       TACSAT 4  12300     463  238.87 14332.2 18671000  6834000 12752500
8      ATLAS 75E  20493    1046  372.22 22333.2 26864000  7417000 17140500
9      TELSTAR 1   5652     950  157.78  9466.8 12023000  7321000  9672000
10        RBSP B  28863     219  501.59 30095.4 35234000  6590000 20912000
11    NAVSTAR 18  21244   21164  759.50 45570.0 27615000 27535000 27575000
12  MOLNIYA 1-90  25309     958  451.82 27109.2 31680000  7329000 19504500
13     SATCOM K1  33248     394  585.81 35148.6 39619000  6765000 23192000
14      BREEZE-M  38407   36200 1513.91 90834.6 44778000 42571000 43674500
15        FREGAT  35853   34297 1399.55 83973.0 42224000 40668000 41446000
16      SYNCOM 1  34568   29003 1236.26 74175.6 40939000 35374000 38156500
17 STARLINK-2558    556     554   95.60  5736.0  6927000  6925000  6926000

4 Análisis de los datos

Para estar en condiciones de extraer información de un experimento es necesario hacer un análisis de los datos medidos. Eso es lo que haremos en este apartado. En primer lugar comprobaremos que los satélites seleccionados verifican la tercera ley de Kepler para, a continuación, calcular la masa de la Tierra. Haremos también un análisis de las incertidumbres implicadas en la actividad.

4.1 Comprobación experimental de la tercera ley de Kepler

Ya tenemos los datos preparados para empezar a trabajar con ellos. Lo primero que debemos hacer es comprobar si realmente se puede aplicar la tercera ley de Kepler a nuestros datos. Para ello no hay nada mejor que hacer una representación gráfica.

Recordemos que esta ley afirma que la distancia \(r\) al cubo y el período \(T\) al cuadrado son directamente proporcionales: \[r^3 \propto T^2\]

Para comprobar esta relación representamos \(r^3\) frente a \(T^2\), y esperamos que la gráfica obtenida sea una recta. Para hacer esta representación creamos dos nuevas variables, r3 y T2, cuyo significado es obvio:

r3 = r^3
T2 = T^2

Representar una gráfica básica en R no puede ser más sencillo. Para representar la gráfica de dispersión de r3 (variable dependiente) frente a T^2 (variable independiente) solo es necesaria la siguiente instrucción:

plot(T2, r3)

Por supuesto esta gráfica necesita un poco de refinamiento, pero como primera aproximación sirve a las mil maravillas: nos deja perfectamente claro que los datos medidos cumplen la tercera ley de Kepler.

4.2 Recta de ajuste

Ya que, a simple vista, parece que los datos experimentales se ajustan a una recta, en este apartado procederemos a obtener la ecuación de dicha recta.

La manera de obtener la recta que mejor se ajusta a los datos experimentales es utilizando el método de los mínimos cuadrados. En R la función lm() se encarga aplicar dicho método (para los curiosos, lm quiere decir linear model). Usando esta función, la instrucción lm(r3 ~ T2 - 1) calcula la recta de mejor ajuste de r3 frente a T2 (el valor - 1 indica que hemos forzado que la recta pase por el origen). Habitualmente el resultado del ajuste se guarda en una variable para poder utilizarlo posteriormente:

ajuste = lm(r3 ~ T2 - 1)

La variable ajuste contiene toda la información sobre la recta de mejor ajuste de r3 frente a T2. Por ejemplo, podemos usar la variable para representar la recta junto a los puntos experimentales:

plot(T2, r3)
abline(ajuste)

Ahora vemos aún mejor lo bien que se cumple la relación lineal supuesta inicialmente.

La representación gráfica anterior se puede mejorar hasta límites insospechados. Una muestra de una mejora sustancial, aunque un tanto compleja de programar, es la siguiente:

4.3 Pendiente de la recta

El interés por obtener la recta de ajuste (además de servirnos para hacer las representaciones gráficas anteriores) radica en la información que guardan los coeficientes del ajuste: la pendiente de la recta y la ordenada en el origen. Habitualmente estos coeficientes nos van a permitir obtener los valores de ciertas magnitudes relacionadas con el experimento que estemos llevando a cabo; en el caso que nos ocupa, obtendremos la masa de la Tierra a partir de la pendiente de la recta de ajuste.

Para obtener los coeficientes de nuestro ajuste veamos primero un resumen de la información almacenada en la variable ajuste:

summary(ajuste)

Call:
lm(formula = r3 ~ T2 - 1)

Residuals:
       Min         1Q     Median         3Q        Max 
-1.504e+18 -8.323e+16  3.937e+16  5.215e+17  8.479e+17 

Coefficients:
    Estimate Std. Error t value Pr(>|t|)    
T2 1.010e+13  3.045e+07  331593   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.811e+17 on 16 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.1e+11 on 1 and 16 DF,  p-value: < 2.2e-16

De todos los valores que nos muestra el resumen del ajuste nos interesan especialmente dos: los que se llaman Estimate y Std. Error del apartado Coefficients. Para acceder individualmente a cada uno de estos valores debemos fijarnos en que el valor de Estimate es el primer valor dentro de los coeficientes del resumen; por tanto podemos obtenerlo con la instrucción:

summary(ajuste)$coefficients[1]

[1] 1.009666e+13

Análogamente, el valor de Std. Error es el segundo:

summary(ajuste)$coefficients[2]

[1] 30448965

Ahora bien, ¿qué significado tienen estos valores? Estimate es el valor estimado de la pendiente de la recta, y Std. Error es la incertidumbre de este valor. Es decir, si llamamos \(m\) a la pendiente de la recta y \(u(m)\) a su incertidumbre, estos valores son:

\[ \begin{matrix} m=1,010 \cdot 10^{13} \, m^3 \, s^{-2} \\ u(m)=3,045 \cdot 10^{7} \, m^3 \, s^{-2} \end{matrix} \] De esta manera hemos obtenido no sólo la pendiente de la recta de ajuste, sino que tenemos también su incertidumbre. Y esto es precisamente lo que necesitamos para calcular la masa de la Tierra y la incertidumbre en su medida. Como estos valores los vamos a utilizar en cálculos posteriores, es una buena idea almacenarlos en sendas variables, que llamaremos m (la pendiente) y um (la incertidumbre en la pendiente):

m = summary(ajuste)$coefficients[1]
um = summary(ajuste)$coefficients[2]

4.4 Cálculo de la masa de la Tierra

Ahora ya estamos en condiciones de calcular el valor de la masa de la Tierra, que llamaremos MT. Anteriormente escribimos la tercera ley de Kepler como: \[r^3 = \frac{G M_T}{4\pi^2}T^2\] De esta expresión es inmediato deducir que \(m\), la pendiente de la recta de ajuste de \(r^3\) frente a \(T^2\), es:

\[m = \frac{G M_T}{4\pi^2}\] Despejando \(M_T\) en la expresión anterior se obtiene:

\[M_T = \frac{4\pi^2m}{G}\] Es decir, para calcular \(M_T\) necesitamos conocer dos valores: \(m\) (la pendiente de la recta de ajuste de \(r^3\) frente a \(T^2\)) y \(G\) (la constante de gravitación universal).

El valor de \(m\) lo hemos obtenido en el apartado anterior usando el método de los mínimos cuadrados. Nos falta por tanto \(G\). Para conocer su valor recurrimos a la información que ofrece el Comité de Datos (CODATA) del Consejo Científico Internacional (ISC). El valor recomendado para la constante \(G\) es:

Valor	Incertidumbre
\(6,67430 \cdot 10^{-11} \, m^{3} \, kg^{-1} \, s^{-2}\)	\(0,00015 \cdot 10^{-11} \, m^{3} \, kg^{-1} \, s^{-2}\)

De nuevo, almacenamos estos dos valores en las variables G y uG, respectivamente, para poder utilizarlos posteriormente:

G = 6.67430e-11
uG = 0.00015e-11

Ya tenemos todo lo que necesitamos para calcular la masa de la Tierra:

MT = 4*pi^2*m/G
MT

[1] 5.972164e+24

Es decir, el valor de la masa de la Tierra obtenido en este experimento es: \[M_T=5,972164 \cdot 10^{24} \, kg\]

El objetivo de esta práctica era calcular la masa de la Tierra, y acabamos de obtener ese valor. ¿Quiere esto decir que el trabajo ya se ha terminado? La respuesta es no. Falta algo fundamental: el cálculo de la incertidumbre de esta medida.

4.5 Incertidumbre

Siempre que midamos una magnitud, junto con el valor experimental debemos expresar su incertidumbre. Pero, ¿para qué necesitamos conocer esta incertidumbre?

• Por un lado, la incertidumbre cuantifica la duda que tenemos sobre la validez del resultado obtenido; al fin y al cabo, el valor de la masa de la Tierra del apartado anterior no deja de ser una estimación.
  
• Por otro lado, la incertidumbre también va a determinar con cuántos decimales debemos expresar la medida: ¿el valor correcto es \(5,972164 \cdot 10^{24} \, kg\)? ¿O \(5,972 \cdot 10^{24} \, kg\)? ¿O debemos redondear a \(6 \cdot 10^{24} \, kg\)? La respuesta viene determinada por el valor de la incertidumbre.
  
• Además, para poder concluir si nuestra medida es o no consistente con el valor aceptado de la masa de la Tierra, es decir, para ver si hay o no discrepancias significativas entre el valor medido y el valor “real”, debemos primero conocer la incertidumbre en la medida.

En este apartado calcularemos cuánto vale la incertidumbre con que hemos medido la masa de la Tierra en este experimento.

Cuando calculamos la masa de la Tierra vimos que su valor dependía de dos magnitudes: la pendiente \(m\) y la constante de gravitación universal \(G\):

\[M_T = \frac{4\pi^2m}{G}\] Por tanto, la incertidumbre en \(M_T\) vendrá determinada por la incertidumbre en estas dos magnitudes (el valor de \(4\pi^2\) es exacto, y por tanto no está sujeto a incertidumbre).

El primer paso consiste en analizar, por propagación de incertidumbres, cuánto contribuye la incertidumbre de cada una de las magnitudes \(m\) y \(G\) a la incertidumbre en la masa de la Tierra.

Incertidumbre debida a \(m\)

Comenzamos calculando la incertidumbre en la masa de la Tierra debida a la pendiente \(m\), a la que llamaremos \(u_m(M_T)\):

\[u_m(M_T)= \left| \frac{\partial M_T}{\partial m}\right|u(m)=\frac{4\pi^2}{G}u(m)\] Su valor es:

uMT.m = (4*pi^2/G)*um
uMT.m

[1] 1.801053e+19

Incertidumbre debida a \(G\)

Por su parte, la incertidumbre en la masa de la Tierra debida a la \(G\), \(u_G(M_T)\), es:

\[u_G(M_T)= \left| \frac{\partial M_T}{\partial G}\right|u(G)=\frac{4\pi^2 m}{G^2}u(G)\] y vale:

uMT.G = (4*pi^2*m/G^2)*uG
uMT.G

[1] 1.3422e+20

Antes de calcular la incertidumbre combinada, comparemos la contribución de cada magnitud a la incertidumbre total. Estos son los resultados que hemos obtenido:

\[ \begin{matrix} u_m(M_T)= 1,8011 \cdot 10^{19} \, kg \\ u_G(M_T)= 1,3422 \cdot 10^{20} \, kg \end{matrix} \] Viendo estos valores se observa que la incertidumbre en \(M_T\) debida a la constante \(G\) es un orden de magnitud mayor que la debida a la pendiente \(m\). Es decir, la incertidumbre debida a \(G\) es la contribución dominante a la incertidumbre en la medida de la masa de la Tierra. Por tanto, a la vista de los resultados, podemos concluir que \(G\) es el factor limitante a la hora de hacer medidas más precisas de la masa de la Tierra.

Incertidumbre combinada

Las dos contribuciones a la incertidumbre (debidas a la pendiente \(m\) y la constante \(G\)) las podemos considerar independientes y aleatorias; por tanto, la incertidumbre combinada, \(u(M_T)\), es:

\[u(M_T)= \sqrt{[u_m(M_T)]^2 + [u_G(M_T)]^2}\]

Esta expresión nos da el valor de la incertidumbre total de \(M_T\), cuyo valor es:

uMT = sqrt(uMT.m^2 + uMT.G^2)
uMT

[1] 1.35423e+20

En conclusión, la incertidumbre de \(M_T\) es \(u(M_T)= 1,3542 \cdot 10^{20} \, kg\). Como habíamos analizado, esta incertidumbre está determinada en mayor medida por la incertidumbre debida a \(G\).

4.6 Cómo expresar correctamente el resultado de la medida

Hasta ahora hemos realizado los cálculos (en realidad, el ordenador los ha hecho por nosotros) tomando gran cantidad de decimales para evitar la propagación de errores de redondeo. Sin embargo, aunque matemáticamente podamos obtener ocho, diez o incluso más cifras significativas, no tiene sentido dar el resultado final con tal precisión. De hecho, es la incertidumbre la que va a determinar con cuántas cifras significativas debemos expresar el resultado de la medida.

Cifras significativas en la incertidumbre

En la Guía para la expresión de la incertidumbre de medida (GUM)² se recoge que es suficiente expresar las incertidumbres con dos cifras significativas como máximo.

Siguiendo este criterio, la incertidumbre en la masa de la Tierra será, con dos cifras significativas:

\[u(M_T) = 1,4 \cdot 10^{20} \, kg\]

Este es el valor de la incertidumbre con el que trabajaremos de ahora en adelante.

En R, la función signif() permite expresar un número con una determinada cantidad de cifras significativas. Reescribimos entonces la incertidumbre uMT con dos cifras significativas:

uMT = signif(uMT, 2)
uMT

[1] 1.4e+20

Decimales en la medida

La GUM también establece que los valores medidos deben redondearse de acuerdo con sus incertidumbres. Esto quiere decir que la medida debe expresarse con el mismo número de decimales que tenga la incertidumbre.

Dado que el valor medido de la masa de la Tierra es \(M_T=5,972164 \cdot 10^{24} \, kg\), primero expresaremos la incertidumbre con la misma potencia de diez para poder comparar ambos números:

\[u(M_T) = 0,00014 \cdot 10^{24} \, kg\]

De esta manera, aparte de la potencia de diez, la incertidumbre queda expresada con cinco decimales; por tanto, la masa de la Tierra la escribiremos también con cinco decimales:

\[M_T=5,97216 \cdot 10^{24} \, kg\]

Dado que la masa queda expresada con seis cifras significativas, reescribimos MT de manera análoga a lo que hicimos con uMT:

MT = signif(MT, 6)
MT

[1] 5.97216e+24

Resultado de la medida

Para finalizar debemos poner juntos ambos valores –la medida y su incertidumbre–, sin olvidarnos de las unidades:

\[M_T=(5,97216 \pm 0,00014)\cdot 10^{24} \, kg\]

Por fin, este es el resultado de nuestro experimento. Este resultado viene a decirnos que el valor que estimamos para la masa de la Tierra se encuentra entre \(5,97202 \cdot 10^{24} \, kg\) y \(5,97230 \cdot 10^{24} \, kg\). Dicho de otra manera, en el contexto del desarrollo experimental realizado, cualquier valor de la masa de la Tierra entre \(5,97202 \cdot 10^{24} \, kg\) y \(5,97230 \cdot 10^{24} \, kg\) es consistente con los datos medidos.

Si decidimos expresar la incertidumbre únicamente con una cifra significativa, el resultado quedaría como sigue:

\[M_T=(5,9722 \pm 0,0001)\cdot 10^{24} \, kg\]

lo cual indica que el valor de \(M_T\) está entre \(5,9721 \cdot 10^{24} \, kg\) y \(5,9723 \cdot 10^{24} \, kg\).

5 Conclusiones

Para finalizar la actividad sacaremos unas conclusiones a partir de los resultados anteriores. Por ejemplo, la incertidumbre en nuestra medida de la masa de la Tierra, ¿es grande o es pequeña? ¿Hay mucha diferencia entre el valor que hemos obtenido y el habitualmente aceptado para la masa de la Tierra? Plantearse estas preguntas es una parte muy importante del trabajo experimental, ya que nos sirven para analizar de manera crítica los resultados de un experimento.

5.1 Incertidumbre relativa

La incertidumbre con que medimos la masa de la Tierra en nuestro experimento es de \(1,4 \cdot 10^{20} \, kg\). La pregunta que vamos a contestar en este apartado es si esta incertidumbre es grande o es pequeña. Nos gustaría que fuese pequeña, ya que cuanto menor sea la incertidumbre menos dispersión habrá en los posibles valores que puede tomar la magnitud. Sin embargo, en términos absolutos esta incertidumbre parece muy grande… aunque por otro lado, comparando la incertidumbre con la masa de la Tierra, puede que no sea tan grande. Esta comparación es precisamente lo que necesitamos hacer para cuantificar el “tamaño” de la incertidumbre. Calcularemos lo que llamamos incertidumbre relativa dividiendo la incertidumbre entre el valor de la medida:

\[\frac{u(M_T)}{M_T}\]

Si este valor lo multiplicamos por 100 obtenemos la incertidumbre relativa en tanto por ciento.

Veamos qué dicen los números:

uMT/MT*100

[1] 0.00234421

Ahora se ve claramente que la incertidumbre en la medida de la masa de la Tierra es muy pequeña, de tan solo el 0,002 %. Este resultado lo podemos ver como una consecuencia de lo bien que cumplían los datos la tercera ley de Kepler.

5.2 Discrepancia entre el valor experimental y el valor aceptado

Acabamos de comprobar que la incertidumbre en nuestra medida de la masa de la Tierra es pequeña. Pero que la incertidumbre sea pequeña no quiere decir que el resultado esté necesariamente cerca del valor “real” de lo que estamos midiendo. Dicho de otra manera, hemos realizado una medida precisa, sí, pero ¿es también exacta?

En multitud de situaciones no podremos responder esta pregunta, ya que el valor real del mensurando (lo que estamos midiendo) suele ser desconocido. Pero en este caso podemos confiar en las medidas de la masa de la Tierra que se han hecho en mejores condiciones: con instrumental más preciso, utilizando un método experimental más adecuado, dejando la medida en manos de científicos con más experiencia que nosotros, etc. De acuerdo con la International Astronomical Union (IAU)³, el valor recomendado de la masa de la Tierra es \(5,9722 \cdot 10^{24} \, kg\), así que este es el valor que vamos a aceptar como “correcto”. Además, habitualmente podemos considerar que los valores “oficiales” van a tener una incertidumbre despreciable en comparación con la que podamos obtener en nuestras medidas.

Dicho esto, ¿es nuestro valor compatible con el valor de la IAU? La respuesta a esta pregunta pasa por calcular la discrepancia que existe entre el valor medido y el valor real. La discrepancia mide, como su nombre indica, cuánto se separa nuestra medida del valor aceptado. Es decir, es la diferencia entre estos dos valores (en valor absoluto, ya que no importa que estemos midiendo un valor mayor o menor que el aceptado).

\[ discrepancia = | M_{T_{medida}} - M_{T_{real}} |\]

Calculemos este valor:

MT.real = 5.9722e24
discrepancia = abs(MT - MT.real)
discrepancia

[1] 4e+19

La discrepancia entre el valor medido y el valor aceptado es \(4 \cdot 10^{19} \,kg\). Este valor es menor que la incertidumbre de nuestra medida (que era, recordemos, \(1,4 \cdot 10^{20} \,kg\)), es decir, el valor real está dentro del intervalo de nuestro valor medido. En consecuencia podemos concluir que no se trata de una discrepancia significativa. Es decir, podemos estar satisfechos con el resultado de nuestra medida.

6 Anexo I: Incertidumbres

En el siguiente esquema se muestra un resumen de cómo calcular la incertidumbre en la masa de la Tierra.

7 Anexo II: Programa completo

A lo largo de la actividad hemos ejecutado las líneas de código una a una para ir viendo los resultados según los íbamos necesitando. Adjuntamos aquí el programa completo, con unos ligeros cambios para mejorar la presentación de los resultados. Es de resaltar que, en lugar de leer los datos de un archivo, los hemos introducido “a mano”; esta no es lo habitual, pero es necesario hacerlo así porque en este documento online no tenemos manera de acceder a un archivo local de datos.

# ============================ DATOS ============================

# Radio de la Tierra (m)
RT = 6.371e6

# Constante de gravitación universal y su incertidumbre (m3·kg-1·s-2)
G = 6.67430e-11
uG = 0.00015e-11

# Valor aceptado de la masa de la Tierra (kg)
MT.real = 5.9722e24

# Satélites
satelite <- c("INTELSAT 2-F3", "INTELSAT 3-F4", "GSTAR 4",
              "COSMOS 2382", "ARIANE 5", "GALILEO 24", "TACSAT 4",
              "ATLAS 75E", "TELSTAR 1", "RBSP B", "NAVSTAR 18",
              "MOLNIYA 1-90", "SATCOM K1", "BREEZE-M", "FREGAT",
              "SYNCOM 1", "STARLINK-2558")
apogeo <- c(35850, 39037, 36140, 19151, 35172, 23233, 12300, 20493,
            5652, 28863, 21244, 25309, 33248, 38407, 35853, 34568, 556)
perigeo <- c(35714, 38822, 36088, 19122, 247, 23225, 463, 1046,
             950, 219, 21164, 958, 394, 36200, 34297, 29003, 554)
periodo <- c(1435.49, 1599.25, 1452.49, 675.71, 619.81, 844.69, 238.87,
             372.22, 157.78, 501.59, 759.5, 451.82, 585.81, 1513.91,
             1399.55, 1236.26, 95.6)
datos <- data.frame(satelite, apogeo, perigeo, periodo)

# ===============================================================

# Período (s)
T = datos$periodo*60

# Distancia en apogeo y perigeo (m)
rA = datos$apogeo*1000 + RT
rP = datos$perigeo*1000 + RT

# Distancia media (m)
r = (rA + rP)/2

# Tabla de datos con T y r
datos <- cbind(datos, T, r)
datos

# Variables para la tercera ley de Kepler
r3 = r^3
T2 = T^2

# Recta de ajuste
ajuste = lm(r3 ~ T2 - 1)

# Pendiente y su incertidumbre
m = summary(ajuste)$coefficients[1]
um = summary(ajuste)$coefficients[2]
paste("Pendiente:", m, "m3/s2")
paste("Incertidumbre en la pendiente:", um, "m3/s2")

# Masa de la Tierra
MT = 4*pi^2*m/G
paste("Masa de la Tierra:", MT, "kg")

# Incertidumbre en la masa de la Tierra debida a la pendiente
uMT.m = (4*pi^2/G)*um
paste("Incertidumbre debida a la pendiente:", uMT.m, "kg")

# Incertidumbre en la masa de la Tierra debida a G
uMT.G = (4*pi^2*m/G^2)*uG
paste("Incertidumbre debida a G:", uMT.G, "kg")

# Incertidumbre total de la masa de la Tierra
uMT = sqrt(uMT.m^2 + uMT.G^2)
paste("Incertidumbre combinada:", uMT, "kg")

# Masa de la Tierra y su incertidumbre
paste("MT = (", signif(MT, 6), " +- ", signif(uMT, 2), ") kg", sep="")
paste("MT = (", signif(MT, 5), " +- ", signif(uMT, 1), ") kg", sep="")

# Tercera ley de Kepler
plot(T2, r3)
abline(ajuste)

        satelite apogeo perigeo periodo       T        r
1  INTELSAT 2-F3  35850   35714 1435.49 86129.4 42153000
2  INTELSAT 3-F4  39037   38822 1599.25 95955.0 45300500
3        GSTAR 4  36140   36088 1452.49 87149.4 42485000
4    COSMOS 2382  19151   19122  675.71 40542.6 25507500
5       ARIANE 5  35172     247  619.81 37188.6 24080500
6     GALILEO 24  23233   23225  844.69 50681.4 29600000
7       TACSAT 4  12300     463  238.87 14332.2 12752500
8      ATLAS 75E  20493    1046  372.22 22333.2 17140500
9      TELSTAR 1   5652     950  157.78  9466.8  9672000
10        RBSP B  28863     219  501.59 30095.4 20912000
11    NAVSTAR 18  21244   21164  759.50 45570.0 27575000
12  MOLNIYA 1-90  25309     958  451.82 27109.2 19504500
13     SATCOM K1  33248     394  585.81 35148.6 23192000
14      BREEZE-M  38407   36200 1513.91 90834.6 43674500
15        FREGAT  35853   34297 1399.55 83973.0 41446000
16      SYNCOM 1  34568   29003 1236.26 74175.6 38156500
17 STARLINK-2558    556     554   95.60  5736.0  6926000
[1] "Pendiente: 10096659380643.2 m3/s2"
[1] "Incertidumbre en la pendiente: 30448964.7748575 m3/s2"
[1] "Masa de la Tierra: 5.97216390389984e+24 kg"
[1] "Incertidumbre debida a la pendiente: 1.8010532145726e+19 kg"
[1] "Incertidumbre debida a G: 1.34220005930955e+20 kg"
[1] "Incertidumbre combinada: 1.35423001223123e+20 kg"
[1] "MT = (5.97216e+24 +- 1.4e+20) kg"
[1] "MT = (5.9722e+24 +- 1e+20) kg"

---

1. The R Project for Statistical Computing: https://www.r-project.org/.↩︎

2. Evaluation of measurement data. Guide to the expression of uncertainty in measurement, publicación del Bureau International des Poids et Mesures (BIPM). El Centro Español de Metrología ha publicado una traducción de esta guía al castellano: Evaluación de datos de medición. Guía para la expresión de la incertidumbre de medida.↩︎

3. The Astronomical Almanac es una publicación conjunta del United States Naval Observatory (USNO) de Estados Unidos y Her Majesty’s Nautical Almanac Office (HMNAO) del Reino Unido que, entre otros muchos otros datos, recoge los valores de las constantes astronómicas publicadas por la International Astronomical Union (IAU).↩︎

---

Cuánto pesa la Tierra, o cómo medir la masa de nuestro planeta usando R.
Beatriz Padín