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;
}