TIERRA + Saturn V (Apollo 11) + LUNA + SOL [multistage3.c]

#define AA (5.0*9.81) // aceleración constante soportada por los astronautas 5*g (letal)

#define SPECIFIC_IMPULSE_F1 265.4
#define SPECIFIC_IMPULSE_J2 426.0
#define SPECIFIC_IMPULSE_SPACECRAFT 452.0 // inventado por el momento (no afecta al resultado obtenido)
#define EV_F1 (SPECIFIC_IMPULSE_F1*9.81)
#define EV_J2 (SPECIFIC_IMPULSE_J2*9.81)
#define EV_XX (SPECIFIC_IMPULSE_SPACECRAFT*9.81)

#define THRUST_F1 6.7725E6
#define THRUST_J2 1033.5E3

#define POUND2KG 0.45359

#define MC (6484070.0*POUND2KG) // masa del cohete Saturn V (misión Apollo 11)

#define S_IC_STAGE_EXPENDABLES (4739320.0*POUND2KG)
#define FIRST_STAGE_INERT_WEIGHT ((288750.0+11465.0)*POUND2KG)

#define S_II_STAGE_EXPENDABLES (980510.0*POUND2KG)
#define SECOND_STAGE_INERT_WEIGHT ((79920.0+8080.0)*POUND2KG)

#define S_IVB_STAGE_EXPENDABLES (237155.0*POUND2KG)
#define THIRD_STAGE_INERT_WEIGHT ((25000.0+4305.0)*POUND2KG)

#define SPACECRAFT_EXPENDABLES ((23680.0+40605.0)*POUND2KG)
#define SPACECRAFT_INERT_WEIGHT ((4045.0+9520.0+10555.0+12250.0+8910.0)*POUND2KG)

#define DT 0.001 // incremento de tiempo por iteración de cálculo 0.001 segundos

#define MT 5.97E24 // masa de la Tierra
#define RT 6.379E6 // radio de la Tierra
#define ML 7.35E22 // masa de la Luna
#define RL 1.739E6 // radio de la Luna
#define MS 1.99E30 // masa del Sol
#define RS 695700.0E3 // radio del Sol
#define DTL 356400.0E3 // distancia Tierra-Luna (entre los centros de masas)
#define DTS 146.0E9 // distancia Tierra-Sol (entre los centros de masas)
#define GG 6.67259E-11 // constante de gravitación universal



#include <math.h>
#include <stdio.h>
#include <stdlib.h>



int main(int argc, char **argv) {
  double mm; // masa restante del cohete
  double mg; // masa de gas expulsada en iteración
  double mc; // masa de combustible consumido en la etapa
  double mi; // masa inerte en la etapa
  double ee; // empuje del cohete
  double ff; // fuerza sobre el cohete
  double aa; // aceleración del cohete
  double vvi; // velocidad inicial del cohete
  double vvf; // velocidad final del cohete
  double rri; // posición inicial del cohete
  double rrf; // posición final del cohete
  double tt; // tiempo
  double rr2; // posición donde se anulan las fuerzas
  double vescape; // velocidad de escape
  double EV; // "exhaust velocity" de los gases
  int etapa;

  rr2 = 2.448529223E8;

  tt = 0.0;
  mm = MC;
  mc = 0.0;
  vvi = 0.0;
  rri = RT;
  etapa = 1;
  EV = EV_F1;
  mi = FIRST_STAGE_INERT_WEIGHT + SECOND_STAGE_INERT_WEIGHT + THIRD_STAGE_INERT_WEIGHT + SPACECRAFT_INERT_WEIGHT;
  do {
    mg = (DT*mm*AA)/(EV+DT*AA);

    mm = mm - mg;
    mc = mc + mg;

    ee = mg*EV/DT;
    ff = ee + GG*(mm*ML)/pow(DTL-rri,2.0) + GG*(mm*MS)/pow(DTS-rri,2.0) - GG*(mm*MT)/pow(rri,2.0);
    aa = ff/mm;
    vvf = aa*DT + vvi;
    rrf = aa*pow(DT,2.0)/2.0 + vvi*DT + rri;

    vescape = sqrt(2) * sqrt(rrf * rr2 * (DTL * DTL * DTS * DTS - DTL * DTL * DTS * rr2 - DTL * DTL * DTS * rrf + DTL * DTL * rr2 * rrf - DTL * DTS * DTS * rr2 - DTL * DTS * DTS * rrf + DTL * DTS * rr2 * rr2 + 2 * DTL * DTS * rr2 * rrf + DTL * DTS * rrf * rrf - DTL * rr2 * rr2 * rrf - DTL * rr2 * rrf * rrf + DTS * DTS * rr2 * rrf - DTS * rr2 * rr2 * rrf - DTS * rr2 * rrf * rrf + rr2 * rr2 * rrf * rrf) * GG * (DTL * DTL * DTS * DTS * MT * rr2 - DTL * DTL * DTS * DTS * MT * rrf - DTL * DTL * DTS * MT * rr2 * rr2 + DTL * DTL * DTS * MT * rrf * rrf + DTL * DTL * MS * rr2 * rr2 * rrf - DTL * DTL * MS * rr2 * rrf * rrf + DTL * DTL * MT * rr2 * rr2 * rrf - DTL * DTL * MT * rr2 * rrf * rrf - DTL * DTS * DTS * MT * rr2 * rr2 + DTL * DTS * DTS * MT * rrf * rrf + DTL * DTS * MT * (double) pow((double) rr2, (double) 3) + DTL * DTS * MT * rr2 * rr2 * rrf - DTL * DTS * MT * rr2 * rrf * rrf - DTL * DTS * MT * (double) pow((double) rrf, (double) 3) - DTL * MS * (double) pow((double) rr2, (double) 3) * rrf + DTL * MS * rr2 * (double) pow((double) rrf, (double) 3) - DTL * MT * (double) pow((double) rr2, (double) 3) * rrf + DTL * MT * rr2 * (double) pow((double) rrf, (double) 3) + DTS * DTS * ML * rr2 * rr2 * rrf - DTS * DTS * ML * rr2 * rrf * rrf + DTS * DTS * MT * rr2 * rr2 * rrf - DTS * DTS * MT * rr2 * rrf * rrf - DTS * ML * (double) pow((double) rr2, (double) 3) * rrf + DTS * ML * rr2 * (double) pow((double) rrf, (double) 3) - DTS * MT * (double) pow((double) rr2, (double) 3) * rrf + DTS * MT * rr2 * (double) pow((double) rrf, (double) 3) + ML * (double) pow((double) rr2, (double) 3) * rrf * rrf - ML * rr2 * rr2 * (double) pow((double) rrf, (double) 3) + MS * (double) pow((double) rr2, (double) 3) * rrf * rrf - MS * rr2 * rr2 * (double) pow((double) rrf, (double) 3) + MT * (double) pow((double) rr2, (double) 3) * rrf * rrf - MT * rr2 * rr2 * (double) pow((double) rrf, (double) 3))) / rrf / rr2 / (DTL * DTL * DTS * DTS - DTL * DTL * DTS * rr2 - DTL * DTL * DTS * rrf + DTL * DTL * rr2 * rrf - DTL * DTS * DTS * rr2 - DTL * DTS * DTS * rrf + DTL * DTS * rr2 * rr2 + 2 * DTL * DTS * rr2 * rrf + DTL * DTS * rrf * rrf - DTL * rr2 * rr2 * rrf - DTL * rr2 * rrf * rrf + DTS * DTS * rr2 * rrf - DTS * rr2 * rr2 * rrf - DTS * rr2 * rrf * rrf + rr2 * rr2 * rrf * rrf);

    tt = tt + DT;

    printf("Etapa : %d\n", etapa);
    printf("Tiempo : %g [s]\n", tt);
    printf("Velocidad de salida de los gases respecto al cohete : %g [m/s]\n", EV);
    printf("Masa de gas expulsada por segundo : %g [kg/s]\n", mg/DT);
    printf("Masa restante del cohete : %g [kg]\n", mm);
    printf("Masa restante del cohete (porcentaje) : %g [%%]\n", mm/MC*100.0);
    printf("Masa inerte del cohete en esta etapa: %g [kg]\n", mi);
    if (etapa == 1) {
      printf("Empuje del cohete : %g [N] (%g motores F-1)\n", ee, ee/THRUST_F1);
    } else if ((etapa == 2) || (etapa == 3)) {
      printf("Empuje del cohete : %g [N] (%g motores J-2)\n", ee, ee/THRUST_J2);
    } else {
      printf("Empuje del cohete : %g [N]\n", ee);
    }
    printf("Fuerza total sobre el cohete : %g [N]\n", ff);
    printf("Aceleración sufrida por los astronautas : %g [m/(s^2)]\n", AA);
    printf("Aceleración del cohete : %g [m/(s^2)]\n", aa);
    printf("Velocidad del cohete : %g [m/s]\n", vvf);
    printf("Velocidad del cohete (km/h) : %g [km/h]\n", vvf/1000.0*3600.0);
    printf("Velocidad de escape : %g [m/s]\n", vescape);
    printf("Velocidad de escape (km/h) : %g [km/h]\n", vescape/1000.0*3600.0);
    printf("Posición del cohete : %g [m]\n", rrf);
    printf("Altura del cohete sobre la superficie de la Tierra : %g [m]\n", rrf-RT);
    printf("Altura del cohete sobre la superficie de la Tierra (km) : %g [km]\n", (rrf-RT)/1000.0);
    printf("Distancia del cohete a la superficie de la Luna : %g [m]\n", DTL-rrf-RL);
    printf("Distancia del cohete a la superficie de la Luna (km) : %g [km]\n\n", (DTL-rrf-RL)/1000.0);

    switch (etapa) {
      case 1:
        if (mc >= S_IC_STAGE_EXPENDABLES) {
          printf("AGOTADO COMBUSTIBLE DE LA PRIMERA ETAPA : %g [kg]\n", S_IC_STAGE_EXPENDABLES);
          printf("APAGADOS LOS 5 MOTORES F-1 DE LA PRIMERA ETAPA\n");
          mm = mm - FIRST_STAGE_INERT_WEIGHT;
          mi = SECOND_STAGE_INERT_WEIGHT + THIRD_STAGE_INERT_WEIGHT + SPACECRAFT_INERT_WEIGHT;
          printf("LIBERADO PESO INERTE DE LA PRIMERA ETAPA : %g [kg]\n", FIRST_STAGE_INERT_WEIGHT);
          printf("Masa restante del cohete : %g [kg]\n", mm);
          printf("Masa restante del cohete (porcentaje) : %g [%%]\n", mm/MC*100.0);
          EV = EV_J2;
          printf("ENCENDIDOS LOS 5 MOTORES J-2 DE LA SEGUNDA ETAPA\n\n");
          mc = 0.0;
          etapa = 2;
        }
        break;

      case 2:
        if (mc >= S_II_STAGE_EXPENDABLES) {
          printf("AGOTADO COMBUSTIBLE DE LA SEGUNDA ETAPA : %g [kg]\n", S_II_STAGE_EXPENDABLES);
          printf("APAGADOS LOS 5 MOTORES J-2 DE LA SEGUNDA ETAPA\n");
          mm = mm - SECOND_STAGE_INERT_WEIGHT;
          mi = THIRD_STAGE_INERT_WEIGHT + SPACECRAFT_INERT_WEIGHT;
          printf("LIBERADO PESO INERTE DE LA SEGUNDA ETAPA : %g [kg]\n", SECOND_STAGE_INERT_WEIGHT);
          printf("Masa restante del cohete : %g [kg]\n", mm);
          printf("Masa restante del cohete (porcentaje) : %g [%%]\n", mm/MC*100.0);
          printf("ENCENDIDO EL MOTOR J-2 DE LA TERCERA ETAPA\n\n");
          mc = 0.0;
          etapa = 3;
        }
        break;

      case 3:
        if (mc >= S_IVB_STAGE_EXPENDABLES) {
          printf("AGOTADO COMBUSTIBLE DE LA TERCERA ETAPA : %g [kg]\n", S_IVB_STAGE_EXPENDABLES);
return 1;
          printf("APAGADO EL MOTOR J-2 DE LA TERCERA ETAPA\n");
          mm = mm - THIRD_STAGE_INERT_WEIGHT;
          mi = SPACECRAFT_INERT_WEIGHT;
          printf("LIBERADO PESO INERTE DE LA TERCERA ETAPA : %g [kg]\n", THIRD_STAGE_INERT_WEIGHT);
          printf("Masa restante del cohete : %g [kg]\n", mm);
          printf("Masa restante del cohete (porcentaje) : %g [%%]\n", mm/MC*100.0);
          EV = EV_XX;
          printf("ENCENDIDO EL MOTOR DE LA CUARTA ETAPA (ASTRONAVE)\n");
          printf("¡¡¡USANDO COMBUSTIBLE NECESARIO PARA EL VIAJE DE RETORNO!!!\n\n");
          mc = 0.0;
          etapa = 4;
        }
        break;

      case 4:
        if (mc >= SPACECRAFT_EXPENDABLES) {
          printf("AGOTADO COMBUSTIBLE DE LA CUARTA ETAPA (ASTRONAVE) : %g [kg]\n", SPACECRAFT_EXPENDABLES);
          printf("APAGADO EL MOTOR DE LA CUARTA ETAPA (ASTRONAVE)\n\n");
          etapa = 5;
        }
    } // switch

    vvi = vvf;
    rri = rrf;
  } while ((vvf < vescape) && (etapa != 5));

  if (vvf >= vescape) {
    printf("ALCANZADA VELOCIDAD DE ESCAPE\n\n");
  }

  return 0;
}


Etapa : 3
Tiempo : 251.312 [s]
Velocidad de salida de los gases respecto al cohete : 4179.06 [m/s]
Masa de gas expulsada por segundo : 739.248 [kg/s]
Masa restante del cohete : 62983.9 [kg]
Masa restante del cohete (porcentaje) : 2.1415 [%]
Masa inerte del cohete en esta etapa: 33831 [kg]
Empuje del cohete : 3.08936e+06 [N] (2.98922 motores J-2)
Fuerza total sobre el cohete : 2.65958e+06 [N]
Aceleración sufrida por los astronautas : 49.05 [m/(s^2)]
Aceleración del cohete : 42.2263 [m/(s^2)]
Velocidad del cohete : 10141.9 [m/s]
Velocidad del cohete (km/h) : 36510.8 [km/h]
Velocidad de escape : 10202.3 [m/s]
Velocidad de escape (km/h) : 36728.4 [km/h]
Posición del cohete : 7.63707e+06 [m]
Altura del cohete sobre la superficie de la Tierra : 1.25807e+06 [m]
Altura del cohete sobre la superficie de la Tierra (km) : 1258.07 [km]
Distancia del cohete a la superficie de la Luna : 3.47024e+08 [m]
Distancia del cohete a la superficie de la Luna (km) : 347024 [km]

AGOTADO COMBUSTIBLE DE LA TERCERA ETAPA : 107571 [kg]

VELOCIDAD DE ESCAPE TIERRA-LUNA MAL

LA EXPRESIÓN DE LA VELOCIDAD DE ESCAPE TIERRA-LUNA ESTÁ MAL EN TODO LO ANTERIOR.

TAMPOCO SE PUEDE DESPRECIAR EL EFECTO DEL SOL.

UN CASO MÁS FAVORABLE: PUNTO DE LANZAMIENTO, LUNA Y SOL ALINEADOS EN ESE ORDEN.

La velocidad de escape Tierra-Luna-Sol creo que puede ser:

sqrt(2)*sqrt(rrf*rr2*(DTL^2*DTS^2-DTL^2*DTS*rr2-DTL^2*DTS*rrf+DTL^2*rr2*rrf-DTL*DTS^2*rr2-DTL*DTS^2*rrf+DTL*DTS*rr2^2+2*DTL*DTS*rr2*rrf+DTL*DTS*rrf^2-DTL*rr2^2*rrf-DTL*rr2*rrf^2+DTS^2*rr2*rrf-DTS*rr2^2*rrf-DTS*rr2*rrf^2+rr2^2*rrf^2)*GG*(DTL^2*DTS^2*MT*rr2-DTL^2*DTS^2*MT*rrf-DTL^2*DTS*MT*rr2^2+DTL^2*DTS*MT*rrf^2+DTL^2*MS*rr2^2*rrf-DTL^2*MS*rr2*rrf^2+DTL^2*MT*rr2^2*rrf-DTL^2*MT*rr2*rrf^2-DTL*DTS^2*MT*rr2^2+DTL*DTS^2*MT*rrf^2+DTL*DTS*MT*rr2^3+DTL*DTS*MT*rr2^2*rrf-DTL*DTS*MT*rr2*rrf^2-DTL*DTS*MT*rrf^3-DTL*MS*rr2^3*rrf+DTL*MS*rr2*rrf^3-DTL*MT*rr2^3*rrf+DTL*MT*rr2*rrf^3+DTS^2*ML*rr2^2*rrf-DTS^2*ML*rr2*rrf^2+DTS^2*MT*rr2^2*rrf-DTS^2*MT*rr2*rrf^2-DTS*ML*rr2^3*rrf+DTS*ML*rr2*rrf^3-DTS*MT*rr2^3*rrf+DTS*MT*rr2*rrf^3+ML*rr2^3*rrf^2-ML*rr2^2*rrf^3+MS*rr2^3*rrf^2-MS*rr2^2*rrf^3+MT*rr2^3*rrf^2-MT*rr2^2*rrf^3))/(rrf*rr2*(DTL^2*DTS^2-DTL^2*DTS*rr2-DTL^2*DTS*rrf+DTL^2*rr2*rrf-DTL*DTS^2*rr2-DTL*DTS^2*rrf+DTL*DTS*rr2^2+2*DTL*DTS*rr2*rrf+DTL*DTS*rrf^2-DTL*rr2^2*rrf-DTL*rr2*rrf^2+DTS^2*rr2*rrf-DTS*rr2^2*rrf-DTS*rr2*rrf^2+rr2^2*rrf^2))

siendo rr2 el punto donde se anula el campo gravitatorio Tierra-Luna-Sol entre la Tierra y la Luna:

rr2 = 2.448529223E8 (aprox.)