//linda Cobb
//winter 98/99
//http://www.timestocome.com
//model classical electron orbits in the helium atom
//in 3d using opengl



//written and compiled under redhat linux
//to compile and link
// link in   -lm -lMesaGl -lMesaGLU -lglut -lMesatk -lMesaaux
//#include <gl/glut>

//to compile under OSX 
//gcc helium3d.c -o helium -framework OpenGL -framework GLUT






#include <stdio.h>

#include <math.h>

#include <GLUT/glut.h>              //open gl libraries for graphics



#define N    12                    //number of o.d.e. to solve

#define loop     2000000L          //number of time steps to take

#define DT_INIT  1.0e-5            //initial time step



static GLfloat theta[] = {10.0, 10.0, 10.0};

static GLint axis = 2;





double *dvector(long nl, long nh)

/* allocate a double vector with subscript range v[nl..nh] */

{

	double *v;



	v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));

	if (!v) printf("allocation failure in dvector()");

	return v-nl+1;

}



void free_dvector(double *v, long nl, long nh)

/* free a double vector allocated with dvector() */

{

	free((char*) (v+nl-1));

}





void rk4(double y[], double  dydx[], int n, double x, double h, double yout[],

	void (*derivs)(double, double [], double []))

{



	int i;

	double xPlush, halfstep, h6, *dym, *dyt, *yt;



	dym = dvector(1,n);

	dyt = dvector(1,n);

	yt  = dvector(1,n);

	halfstep  = h*0.5;

	h6  = h/6.0;

	



	for (i=1;i<=n;i++)                       //for each ode do an euler step

	  yt[i] = y[i] + halfstep*dydx[i];       //dydx given for here



	(*derivs)(x + halfstep, yt, dyt);        //2nd rk4 step  

	for (i=1;i<=n;i++) 

	  yt[i] = y[i] + halfstep*dyt[i];



	(*derivs)(x + halfstep, yt, dym);        //3rd step

	for (i=1;i<=n;i++) {

	  yt[i] = y[i] + h*dym[i];

	  dym[i] += dyt[i];

	}



	(*derivs)(x + h, yt, dyt);                 //4th step

	for (i=1;i<=n;i++)                         //weight total

		yout[i] = y[i] + h6*(dydx[i] + dyt[i] + 2.0*dym[i]);





	free_dvector(yt,1,n);

	free_dvector(dyt,1,n);

	free_dvector(dym,1,n);

}



//calculated derivatives  dydt[N] = f';

void derivs(double t, double y[], double dydt[])

{



  double r1, r2, r3, r123, r12;          //magnitudes of distances

  double r1Cubed, r2Cubed, r12Cubed, r123Cubed;

  double r1x, r2x, r1y, r2y, r1z, r2z;

  double r1x_r2x, r1y_r2y, r1z_r2z;



  r1x = y[1];                           //rename these to make debugging easier

  r2x = y[2];

  r1y = y[5];

  r2y = y[6];

  r1z = y[9];

  r2z = y[10];



  r1x_r2x = r1x - r2x;

  r1y_r2y = r1y - r2y;

  r1z_r2z = r1z - r2z;

 



  r1 = sqrt( r1x*r1x + r1y*r1y + r1z*r1z);

  r2 = sqrt( r2x*r2x + r2y*r2y + r2z*r2z);

  



  r12 = sqrt( (r1x - r2x)*(r1x - r2x) + (r1y - r2y)*(r1y - r2y) );

  r123 = sqrt( (r1x - r2x)*(r1x - r2x) + (r1y - r2y)*(r1y - r2y) 

               + (r1z - r2z)*(r1z - r2z) );

  



  r1Cubed = r1*r1*r1;

  r2Cubed = r2*r2*r2;

  r12Cubed = r12*r12*r12;

  r123Cubed = r123*r123*r123;

  

  

  dydt[1] = y[3];                                 //e1x

  dydt[2] = y[4];                                 //e2x

  dydt[3] = -2.0 * r1x/r1Cubed + r1x_r2x/r123Cubed;//e1vx

  dydt[4] = -2.0 * r2x/r2Cubed - r1x_r2x/r123Cubed;//e2vx

  dydt[5] = y[7];                                 //e1y

  dydt[6] = y[8];                                 //e2y

  dydt[7] = -2.0 * r1y/r1Cubed + r1y_r2y/r123Cubed;//e1vy

  dydt[8] = -2.0 * r2y/r2Cubed - r1y_r2y/r123Cubed;//e2vy

  dydt[9] = y[11];                                 //e1z

  dydt[10] = y[12];                                //e2z

  dydt[11] = -2.0 * r1z/r1Cubed + r1z_r2z/r123Cubed;//e1Vz

  dydt[12] = -2.0 * r2z/r2Cubed - r1z_r2z/r123Cubed;//e2Vz

  



}







void drawAxis(){



 

  typedef GLfloat point[3];          //define a point array

  long int k, j;                     //loop control

  point p = {0.0, 0.0, 0.0};         //initialize



 glColor3f(0.0, 0.0, 1.0);

  for(k=-250; k<250; k++){            // blue y axis

    p[0] = 0.0;

    p[1] = k;

    p[2] = 0.0;

    glBegin(GL_POINTS);

    glVertex3fv(p);

    glEnd();

  }



  glColor3f(0.0, 1.0, 0.0);           // green x axis

  for(k=-250; k<250; k++){

    p[0] = k;

    p[1] = 0.0;

    p[2] = 0.0;

    glBegin(GL_POINTS);

    glVertex3fv(p);

    glEnd();

  }





  glColor3f(1.0, 0.0, 0.0);           // red z axis

    for(k=-250; k<250; k++){

    p[0] = 0.0;

    p[1] = 0.0;

    p[2] = k;

    glBegin(GL_POINTS);

    glVertex3fv(p);

    glEnd();

  }

  



}







//plot data in an opengl window previously stored in array 

//while calculating data

void displayOrbit (){ 



    int    i, nsteps, j, k;         //loop controls

	double t, dt, *y, *dydx, *yout;

    typedef GLfloat point3[3];      //define a point array

    point3 p1 = {0.0, 0.0, 0.0};    //e1 x, y, z

    point3 p2 = {0.0, 0.0, 0.0};    //e2 x, y, z

    FILE *fp;                       //file pointer

    double e1velocity, e2velocity, distance;

 

    y = dvector(1,N);               //data before rk4 call

	dydx = dvector(1,N);            //derivative info

	yout = dvector(1,N);            //data after rk4 call



    j = 0; k = 0;

	nsteps = 0;                     //initialize max time steps to take

	t = 0.0;                        //initial time

	dt = DT_INIT;                   //time step size

	

    fp = fopen("data", "w");        //create file, allow writing

//###############################################################

//set initial conditions here.....



	y[1] = 0.2;                     //initial x position e1

	y[2] = -1.0;                     //initial x position e2

	y[3] = 0.75;                     //initial x velocity e1

	y[4] = 0.9;                     //initial x velocity e2

	y[5] = 0.0;                     //initial y position e1

    y[6] = -1.0;                     //initial y position e2

    y[7] = -0.9;                     //initial y velocity e1

    y[8] = -0.75;                    //initial y velocity e2

    y[9] = 1.0;                     //initial z position e1

    y[10] = 0.0;                   //initial z position e2

    y[11] = 0.9;                    //initial z velocity e1

    y[12] = -0.1;                    //initial z velocity e2

    

      

	derivs(t, y, dydx);             //initial value of derivatives

    fprintf(fp, "\ne1 initial cond. pos( %.3lf, %.3lf, %.3lf), vel (%.3lf,%.3lf,%.3lf)",

            y[1],y[5],y[9],y[3],y[7],y[11]);

    fprintf(fp, "\ne2 initial cond. pos( %.3lf, %.3lf, %.3lf), vel (%.3lf,%.3lf,%.3lf)",

            y[2],y[6],y[10],y[4],y[8],y[12]);



    

  

  //set up display window

  glClear(GL_COLOR_BUFFER_BIT);    //clear window

  glLoadIdentity();

 

  //rotate picture

  glRotatef(theta[0], 1.0, 0.0, 0.0);

  glRotatef(theta[1], 0.0, 1.0, 0.0);

  glRotatef(theta[2], 0.0, 0.0, 1.0);

  

  drawAxis();



  //calculate and plot at the same time

    do{	

	

	  t += dt;                     //increase time by timestep

	  rk4( y, dydx, N, t, dt, yout, derivs);

  

           

  	  derivs(t, y, dydx);          //update all values

      for (i=1; i<=N; i++){  y[i] = yout[i];  }       

      

      j++;

     

      e1velocity = sqrt(y[3]*y[3] + y[7]*y[7] + y[11]*y[11]);

      e2velocity = sqrt(y[4]*y[4] + y[8]*y[8] + y[12]*y[12]);

      distance = sqrt( fabs( (y[1]-y[2])*(y[1]-y[2])

                       + (y[5]-y[6])*(y[5]-y[6])

                       + (y[9]-y[10])*(y[9]-y[10]) ) );

      

      

      

      //scale timestep according to velocity

      if( (e1velocity > 9.5)||(e2velocity > 9.5) )  dt = 10e-7;                  

      else if( (e1velocity < 1.0 )&&(e2velocity < 1.0) ) dt = 10e-5;

      else dt = 10e-6;

      



    

      if(j>=1000){





        k++;

        if(k >=100){

          k = 0;

          

          fprintf(fp, "\n time %lf, distance %lf", t, distance);

          

          fprintf(fp, "\n e1 x %lf, y %lf, z %lf, vel %lf ", 

                 yout[1], yout[5], yout[9], e1velocity);   

          fprintf(fp, "\n e2 x %lf, y %lf, z %lf, vel %lf ", 

                  yout[2], yout[6], yout[10], e2velocity);   

        }

        



        j=0;

        

        p1[0] = yout[1]*50;          //e1x position (scale data so can see better...)

        p2[0] = yout[2]*50;          //e2x position

        p1[1] = yout[5]*50;          //e1y position

        p2[1] = yout[6]*50;          //e2y position

        p1[2] = yout[9]*50;          //e1z

        p2[2] = yout[10]*50;         //e2z

        

        glBegin(GL_POINTS);          //plot points

          glColor3f(1.0, 0.0, 1.0);  //purple for e1

          glVertex3fv(p1);

          glColor3f(0.0, 1.0, 1.0);  //turquoise for e2

          glVertex3fv(p2);

        glEnd();

        

      }

      



      nsteps++;                    //increment loop counter

	



    }while(nsteps < loop);



    free_dvector(yout,1,N);        //cleanup

	free_dvector(dydx,1,N);

	free_dvector(y,1,N);		



    fclose(fp);                    //close-save file

    

    glutSwapBuffers();

    



}





//initialize opengl window for plotting

void myinit(void){



  glClearColor(0.0, 0.0, 0.0, 0.0);    //backround color  (r,g,b, alpha)

  glColor3f(1.0, 0.0, 0.0);            //draw color (r, g, b)

  glMatrixMode(GL_PROJECTION);         //projection stack

  glOrtho(-250.0, 250.0, -250.0, 250.0, -250.0, 250.0);  //3d  @(0,0,0)

  glMatrixMode(GL_MODELVIEW);          //model view stack 



}



void reshape(int w, int h){

  

  GLfloat aspect = (float) w / (float) h;

  glViewport(0, 0, w, h);

  glLoadIdentity();

  glFrustum( -aspect, aspect, -100.0, 100.0, 0.0, 250.0);

  glMatrixMode(GL_MODELVIEW);



}



void Key( unsigned char key, int x, int y){

  

  if(key==27)

    exit(0);

}



//change view point....

//using arrow and home keys...

void SpecialKey(int key, int x, int y){





  switch(key){

  case GLUT_KEY_UP:

    theta[0] += 10.0;

    break;

  case GLUT_KEY_DOWN:

    theta[1] += 10.0;

    break;

  case GLUT_KEY_LEFT:

    theta[2] += 10.0;

    break;

   case GLUT_KEY_HOME:

    theta[0] = 0.0;

    theta[1] = 0.0;

    theta[2] = 0.0;

    break;

    

  }

  

  glutPostRedisplay();

}







void plotData(int argc, char **argv){



    glutInit(&argc, argv);         



	  glutInitWindowSize(500,500);

	  glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);

	  (void)glutCreateWindow("3d Helium Atom");

	  glutDisplayFunc(displayOrbit);

	  myinit();

      glutReshapeFunc(reshape);

      glutKeyboardFunc(Key);

      glutSpecialFunc(SpecialKey);

	

    glutMainLoop();



}







void main(int argc, char **argv){

				

	system("clear");               //start w/ clean user interface



	plotData(argc, argv);          //plot data in opengl windows

	

}