LAW OF ATTRACTION OF THE ENERGY

where:
+ F: attraction force between energy 1 and energy 2 [N]
+ G2: attraction constant of the energy [N*(m^2)/(J^2)]
+ E1: energy 1 [J]
+ E2: energy 2 [J]
+ r: distance between energy 1 and energy 2 [m]
+ G: gravitational constant [N*(m^2)/(kg^2)]
+ c: speed of light [m/s]

Trip scheme of Apollo 11 mission

"Un cohete Saturno 5 coloca en órbita terrestre la nave Apolo (16 julio 1969). Después de dar varias vueltas alrededor de la Tierra, la nave es transferida a una órbita lunar por la tercera fase del cohete."

Source: National Geographic DVD "Destino, el espacio".

Presión en baja estratosfera (11-20 km)

Fuente: Sebastián Franchini y Óscar López García. "Introducción a la Ingeniería Aeroespacial". Garceta, 2012.

Presión en troposfera (0-11 km)

Fuente: Sebastián Franchini y Óscar López García. "Introducción a la Ingeniería Aeroespacial". Garceta, 2012.

multistage2.c

//#define AA 46.5975 // aceleración soportada por los astronautas 4.75*g (letal en pocos segundos)
#define AA 44.145 // aceleración soportada por los astronautas 4.5*g (letal en pocos segundos)

#define SPECIFIC_IMPULSE_F1 265.4
#define SPECIFIC_IMPULSE_J2 426.0
#define SPECIFIC_IMPULSE_SPACECRAFT 452.0 // inventado por el momento (se ha buscado el valor más alto posible)
#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 DTL 356400.0E3 // distancia Tierra-Luna (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 = (DTL*(MT-sqrt(MT*ML)))/(MT-ML);

  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*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.0*GG*(MT*(1.0/rrf - 1.0/rr2) + ML*(1.0/(DTL-rrf) - 1.0/(DTL-rr2)))); // MAL

    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 : 279.229 [s]
Velocidad de salida de los gases respecto al cohete : 4179.06 [m/s]
Masa de gas expulsada por segundo : 665.373 [kg/s]
Masa restante del cohete : 62988.6 [kg]
Masa restante del cohete (porcentaje) : 2.14166 [%]
Masa inerte del cohete en esta etapa: 33831 [kg]
Empuje del cohete : 2.78063e+006 [N] (2.6905 motores J-2)
Fuerza total sobre el cohete : 2.36204e+006 [N]
Aceleración sufrida por los astronautas : 44.145 [m/(s^2)]
Aceleración del cohete : 37.4994 [m/(s^2)]
Velocidad del cohete : 9918.18 [m/s]
Velocidad del cohete (km/h) : 35705.4 [km/h]
Velocidad de escape : 10008.6 [m/s] MAL
Velocidad de escape (km/h) : 36031.1 [km/h]
Posición del cohete : 7.74225e+006 [m]
Altura del cohete sobre la superficie de la Tierra : 1.36325e+006 [m]
Altura del cohete sobre la superficie de la Tierra (km) : 1363.25 [km]
Distancia del cohete a la superficie de la Luna : 3.46919e+008 [m]
Distancia del cohete a la superficie de la Luna (km) : 346919 [km]

AGOTADO COMBUSTIBLE DE LA TERCERA ETAPA : 107571 [kg]

LA EXPRESIÓN DE LA VELOCIDAD DE ESCAPE ESTÁ MAL.
TAMPOCO SE PUEDE DESPRECIAR EL EFECTO DEL SOL.

VELOCIDAD DE ESCAPE TIERRA-LUNA

Potenciales electrodinámicos

// Carga del electrón
#define qe -1.602176487e-19
// Masa del electrón
#define me 9.10938215e-31
// Permitividad del vacío
#define ep0 8.841941281e-12
// Permeabilidad del vacío
#define mu0 1.256637062e-6
// Constante Pi
#define Pi 3.141592654
// Número de cargas
#define N 5
// Incremento de posición para las derivadas parciales
#define hh 1.0e-13
// Tiempo entre marcos
#define dd 4.0e-19



#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <GL/glut.h>



typedef struct {
  double xx;
  double yy;
  double zz;
} vector;



int ancho; // ancho de la ventana gráfica
int alto; // alto de la ventana gráfica

double arista = 0.5e-9; // arista del cubo donde están las cargas
double pvcubo = 0.5e-9; // distancia del punto de vista al cubo
vector vcentro; // posición central del cubo
double qq[N] = { qe, -qe, qe, -qe, qe }; // cargas de las N cargas
double mq[N] = { me, me, me, me, me }; // masas de las N cargas
vector rq[N][4]; // posiciones de las N cargas en instantes de tiempo tt, (tt-dd), (tt-2*dd) y (tt-3*dd)

double tt = 0.0; // tiempo



void manejarTeclado(unsigned char tecla, int xx, int yy);
void reformar(int width, int height);
void exhibir(void);

int pintarEscena();
int actualizarPosiciones();
void printString(char *ltexto);

vector fVector(double xx, double yy, double zz);
vector vsum(vector v1, vector v2);
vector vres(vector v1, vector v2);
vector vepv(double ee, vector vv);
double vmod(vector vv);
vector vprv(vector v1, vector v2);
void imprimeV(char *etiqueta, vector vv);

vector AA(vector rr, double qq, vector rq, vector vq);
double Fi(vector rr, double qq, vector rq);
vector fAA(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]);
double fFi(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]);
vector fAA2(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]);



int main(int argc, char **argv) {
  int i;
  int j;
  double rx;
  double ry;
  double rz;
   
  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
  glutInitWindowSize(300, 300);
  glutInitWindowPosition(100, 100);
 
  glutCreateWindow(argv[0]);
  glClearColor(0.0, 0.0, 0.0, 0.0);
  glEnable(GL_COLOR_MATERIAL);
  glShadeModel(GL_SMOOTH);
  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);
  glEnable(GL_CULL_FACE);
  glCullFace(GL_BACK);
  glEnable(GL_DEPTH_TEST);

  vcentro= fVector(0.0, 0.0, 0.0);
 
  srand( (unsigned)time( NULL ) );
 
  // Posiciones iniciales aleatorias de las cargas en cubo central de lado arista/2
  for(i=0; i<N; i++) {
    rx = arista/2.0*((double)rand()/(double)RAND_MAX)-arista/4.0;
    ry = arista/2.0*((double)rand()/(double)RAND_MAX)-arista/4.0;
    rz = arista/2.0*((double)rand()/(double)RAND_MAX)-arista/4.0;
    for(j=0; j<4; j++) {
      rq[i][j] = fVector(rx, ry, rz);
    }
  }

  glutKeyboardFunc(manejarTeclado); // registro función atiende evento teclado
  glutReshapeFunc(reformar); // registro función atiende evento volver a formar
  glutDisplayFunc(exhibir); // registro función atiende evento frame

  glutMainLoop();

  return 0;
}



void manejarTeclado(unsigned char tecla, int xx, int yy) {

  switch(tecla) {
    case 27:
      exit(0);
      break;
    case '+':
      pvcubo = pvcubo + arista / 100.0;
      break;
    case '-':
      pvcubo = pvcubo - arista / 100.0;
  }

  reformar(ancho, alto);
}



void reformar(int width, int height) {
  double beta; // campo de visión

  ancho = width;
  alto = height;

  glViewport(0, 0, width, height); // dimensiones área rectangular visualización

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();

  beta = (360.0 / Pi) * atan(arista / (2.0 * pvcubo));

  // defino perspectiva de la cámara
  gluPerspective(beta, (GLfloat)width/(GLfloat)height, pvcubo, pvcubo + arista);

  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();

  // defino posición y orientación de la cámara respecto al centro del cubo
  // gluLookAt(punto de vista; punto de referencia; vector arriba)
  gluLookAt(vcentro.xx, vcentro.yy, vcentro.zz + pvcubo + arista / 2.0, vcentro.xx, vcentro.yy, vcentro.zz, 0.0, 1.0, 0.0);
}



void exhibir(void) {

  // borrar buffer de imagen
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

  if (!pintarEscena()) {
      fprintf(stderr, "Error en pintarEscena\n");
      exit(1);
  }

  // intercambiar buffers de imagen
  glutSwapBuffers();


  if (!actualizarPosiciones()) {
    fprintf(stderr, "Error en actualizarPosiciones\n");
    exit(1);
  }

  glutPostRedisplay(); // provoco evento frame
}



int pintarEscena() {
  int i;
  char etiqueta[256];
  char lineatexto[256];

  // Pinto fronteras

  glPushMatrix();
  glTranslatef(vcentro.xx, vcentro.yy, vcentro.zz);
  glColor4f(1.0, 1.0, 1.0, 1.0);
  glutWireCube(arista - arista / 1000.0);
  glPopMatrix();

  // Muestro cargas

  fprintf(stdout, "t = %g\n", tt);
  for(i=0; i<N; i++) {
    sprintf(etiqueta, "vrr_q%d =", i);
    imprimeV(etiqueta, rq[i][0]);

    glPushMatrix();
    glTranslatef(rq[i][0].xx, rq[i][0].yy, rq[i][0].zz);
    glColor4f(1.0, 1.0, 1.0, 1.0);
    glutWireSphere(arista/100.0, 10, 10);
    glPopMatrix();

    glBegin(GL_LINES);
    glVertex3d(+arista/2.0, rq[i][0].yy, rq[i][0].zz);
    glVertex3d(rq[i][0].xx, rq[i][0].yy, rq[i][0].zz);
    glVertex3d(rq[i][0].xx, +arista/2.0, rq[i][0].zz);
    glVertex3d(rq[i][0].xx, rq[i][0].yy, rq[i][0].zz);
    glVertex3d(rq[i][0].xx, rq[i][0].yy, +arista/2.0);
    glVertex3d(rq[i][0].xx, rq[i][0].yy, rq[i][0].zz);
    glEnd();
  }
  fprintf(stdout, "\n");

  // Pinto línea de texto

  sprintf(lineatexto, "t = %g", tt);
  printString(lineatexto);

  return 1;
}



int actualizarPosiciones() { 
  int i;
  int j;
  double rx;
  double ry;
  double rz;
  double BBx;
  double BBy;
  double BBz;
  double EEx;
  double EEy;
  double EEz;
  vector BBi; // campo B
  vector EEi; // campo E
  vector FF; // fuerza
  vector vqi;
  vector aa; // aceleración

  // ACTUALIZO LAS POSICIONES DE LAS CARGAS

  // Desplazo posiciones en el array

  for(i=0; i<N; i++) {
    rq[i][3] = rq[i][2];
    rq[i][2] = rq[i][1];
    rq[i][1] = rq[i][0];
  }

  // Calculo últimas posiciones rq[i][0]

  for(i=0; i<N; i++) { // recorre cargas

    // Calculo B y E sobre carga i

    BBx = ((fAA(rq[i][1].xx   , rq[i][1].yy+hh, rq[i][1].zz   , i, qq, rq).zz-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).zz)/hh)-((fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz+hh, i, qq, rq).yy-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).yy)/hh);
    BBy = ((fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz+hh, i, qq, rq).xx-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).xx)/hh)-((fAA(rq[i][1].xx+hh, rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).zz-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).zz)/hh);
    BBz = ((fAA(rq[i][1].xx+hh, rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).yy-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).yy)/hh)-((fAA(rq[i][1].xx   , rq[i][1].yy+hh, rq[i][1].zz   , i, qq, rq).xx-fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).xx)/hh);

    EEx = ((fFi(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq)-fFi(rq[i][1].xx+hh, rq[i][1].yy   , rq[i][1].zz   , i, qq, rq))/hh)-((fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).xx-fAA2(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).xx)/dd);
    EEy = ((fFi(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq)-fFi(rq[i][1].xx   , rq[i][1].yy+hh, rq[i][1].zz   , i, qq, rq))/hh)-((fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).yy-fAA2(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).yy)/dd);
    EEz = ((fFi(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq)-fFi(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz+hh, i, qq, rq))/hh)-((fAA(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).zz-fAA2(rq[i][1].xx   , rq[i][1].yy   , rq[i][1].zz   , i, qq, rq).zz)/dd);

    BBi = fVector(BBx, BBy, BBz);
    EEi = fVector(EEx, EEy, EEz );

    // Calculo F sobre carga i

    vqi = vepv(1.0/dd, vres(rq[i][1], rq[i][2]));
    FF = vepv(qq[i], vsum(EEi, vprv(vqi, BBi))); // Fuerza de Lorentz

    // Calculo aceleración y nueva posición de la carga i

    aa = vepv(1.0/mq[i] , FF); // a = F/m
    rq[i][0] = vsum(vsum(vepv(pow(dd, 2.0)/2.0, aa), vepv(dd, vqi)), rq[i][1]); // rf = a*t^2/2 + vi*t + ri

  } // for

  tt = tt + dd;

  return 1;
}



void printString(char *ltexto) {
  int longitud = strlen(ltexto);
  int i;
  GLboolean lit = glIsEnabled(GL_LIGHTING);

  glDisable(GL_LIGHTING);
  glColor4f(1.0, 1.0 , 1.0, 1.0);

  glRasterPos3f(vcentro.xx - arista / 2.0 + arista / 100.0, vcentro.yy - arista / 2.0 + arista / 100.0, vcentro.zz + arista / 2.0 - arista / 100.0);

  i = 0;
  while (i < longitud) {
    glutBitmapCharacter(GLUT_BITMAP_8_BY_13, ltexto[i]);
    i++;
  }

  if (lit) {
    glEnable(GL_LIGHTING);
  }
}



// Crea vector a partir de sus componentes
vector fVector(double xx, double yy, double zz) {
  vector vv;

  vv.xx = xx;
  vv.yy = yy;
  vv.zz = zz;

  return vv;
}



// Suma vectores; v1+v2
vector vsum(vector v1, vector v2) {
  vector res;

  res.xx = v1.xx + v2.xx;
  res.yy = v1.yy + v2.yy;
  res.zz = v1.zz + v2.zz;

  return res;
}



// Resta vectores; v1-v2
vector vres(vector v1, vector v2) {
  vector res;

  res.xx = v1.xx - v2.xx;
  res.yy = v1.yy - v2.yy;
  res.zz = v1.zz - v2.zz;

  return res;
}



// Multiplica escalar por vector
vector vepv(double ee, vector vv) {
  vector res;

  res.xx = ee * vv.xx;
  res.yy = ee * vv.yy;
  res.zz = ee * vv.zz;

  return res;
}



// Módulo de vector
double vmod(vector vv) {
  return sqrt(pow(vv.xx, 2.0) + pow(vv.yy, 2.0) + pow(vv.zz, 2.0));
}



// Producto vectorial v1 x v2
vector vprv(vector v1, vector v2) {
  vector res;

  res.xx = v1.yy * v2.zz - v1.zz * v2.yy;
  res.yy = v1.zz * v2.xx - v1.xx * v2.zz;
  res.zz = v1.xx * v2.yy - v1.yy * v2.xx;

  return res;
}



// Imprime vector
void imprimeV(char *etiqueta, vector vv) {
  fprintf(stdout, "%s <%g,%g,%g>\n", etiqueta, vv.xx, vv.yy, vv.zz);
}



// Potencial vector magnético electrodinámico A
vector AA(vector rr, double qq, vector rq, vector vq) {
  return vepv(mu0/(4.0*Pi)*qq/vmod(vres(rr,rq)),vq);
}



// Potencial escalar electrodinámico Fi
double Fi(vector rr, double qq, vector rq) {
  return (1.0/(4.0*Pi*ep0)*qq/vmod(vres(rr,rq)));
}



// Potencial vector magnético electrodinámico A total sobre carga i en último instante
vector fAA(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]) {
  vector rr;
  vector AAi;
  vector vqj;
  int j;

  rr = fVector(xx, yy, zz);
  AAi = fVector(0.0, 0.0, 0.0);

  // Acumulo alfa
  for(j=0; j<N; j++) { // recorre cargas
    if (i != j) {
      vqj = vepv(1.0/dd,vres(rq[j][1],rq[j][2]));
      AAi = vsum(AAi, AA(rr, qq[j], rq[j][1], vqj));
    }
  }

  return AAi;
}



// Potencial escalar electrodinámico Fi total sobre carga i en último instante
double fFi(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]) {
  vector rr;
  double Fii;
  int j;

  rr = fVector(xx, yy, zz);
  Fii = 0.0;

  // Acumulo fi
  for(j=0; j<N; j++) { // recorre cargas
    if (i != j) {
      Fii = Fii + Fi(rr, qq[j], rq[j][1]);
    }
  }

  return Fii;
}



// Potencial vector magnético electrodinámico A total sobre carga i en instante anterior al último
vector fAA2(double xx, double yy, double zz, int i, double qq[N], vector rq[N][4]) {
  vector rr;
  vector AAi;
  vector vqj;
  int j;

  rr = fVector(xx, yy, zz);
  AAi = fVector(0.0, 0.0, 0.0);

  // Acumulo alfa
  for(j=0; j<N; j++) { // recorre cargas
    if (i != j) {
      vqj = vepv(1.0/dd,vres(rq[j][2],rq[j][3]));
      AAi = vsum(AAi, AA(rr, qq[j], rq[j][2], vqj));
    }
  }

  return AAi;
}