//linda Cobb

//http://www.timestocome.com

//winter 98/99



//model planetary orbits



#include <stdio.h>

#include <math.h>

#define NRANSI                    //'numerical recipes' libraries

#include "nr.h"

#include "nrutil.h"

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

#include <gl/glu.h>

#include <gl/gl.h>



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

#define pi   3.14159265358979324  //Pi

#define loop 50000                //number of time steps to take

#define mass 1.0                  //mass of earth in a.u.





static int GM = 4.0*pi*pi;        //gravitational constant in A.U.

double displayData[2][loop];      //hold data to plot w/ opengl routines (position)

double dataArray[3][loop];        //opengl data energy, L





//calculated derivatives

//dydt[N] = f';

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

{



    double r3;                          //radius cubed

     r3 = sqrt(y[1]*y[1] + y[2]*y[2]) * (y[1]*y[1] + y[2]*y[2]);



    // r3 = pow( (y[1]*y[1]+ y[2]*y[2]), 3.5);



    dydt[1] = y[3];                     //x position              

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

    dydt[3] = -GM*(y[1]/r3);            //x velocity       

    dydt[4] = -GM*(y[2]/r3);            //y velocity 

   

}







//calculate orbit, position, etc

//double data stores velocity and position information

//double dataArray stores angular momentum in 1st position 

//                        and energy in 2nd position



void calculateOrbit(double velX, double velY, double posX, double posY, double data[N+1][loop]){



        int    i, k,  nsteps;        //loop controls

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

	double j;                    //print data control



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

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

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



	k = 0; j = 0.0;              //loop control used when storing data

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

	t = 0.0;                     //initial time

	dt = 0.0001;                 //time step size

	

	y[1] = posX;                 //initial x position

	y[2] = posY;                 //initial y position

	y[3] = velX;                 //initial x velocity

	y[4] = velY;                 //initial y velocity	



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



        do{	

	

	  t += dt;                   //increase time by timestep

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

	  

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

       	  y[1] = yout[1];           

	  y[2] = yout[2];

	  y[3] = yout[3]; 

	  y[4] = yout[4];

	  

	  displayData[0][nsteps] = yout[1];   //save x,y position for opengl to plot later

	  displayData[1][nsteps] = yout[2];



	  

	  // if(j>=1){                  //collect data in increments

	    data[1][k] = yout[1];    //x position saved

	    data[2][k] = yout[2];    //y position saved

	    data[3][k] = yout[3];    //x velocity saved

	    data[4][k] = yout[4];    //y velocity saved

	    

	    dataArray[0][k] = yout[3]*yout[2] - yout[4]*yout[1]; //angular momentum

	    

	    dataArray[1][k] = (yout[3]*yout[3] + yout[4]*yout[4])*0.50 

	      - GM/(sqrt(yout[1]*yout[1] + yout[2]*yout[2]));    //energy

	    

	    dataArray[2][k] = k;     //time passed

	    

	    //j=0.0;



	    k++;                     //increment array pointer for saved data

	   

	    // }

	    //j+=dt;



	  nsteps++;                   //increment loop counter

	    	    

	}while(nsteps < loop);

	



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

	free_dvector(dydx,1,N);

	free_dvector(y,1,N);		

}









//print data to screen and save to file on hard drive

void printData(double dataout[N+1][loop]){



  FILE *fp;                           //file pointer

  int i;                              //print loop control

  char *heading = " pos(x)  pos(y)   vel(x)   vel(y)      L      energy";



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

  fprintf(fp, heading);               //print file heading





  for(i = 0; i < loop*0.001; i++){          //write data to file

        fprintf(fp,"\n %.1lf      %.1lf      %.1lf      %.1lf,         %.1lf,         %.1lf   ", 

	       dataout[1][i], dataout[2][i], dataout[3][i], dataout[4][i], 

		dataArray[0][i], dataArray[1][i]);

    }





  fclose(fp);                         //close-save file



}





void drawAxis(){



  typedef GLfloat point2[2];          //define a point array

  long int k, j;                      //loop control

  point2 p = {0.0, 0.0};              //begin here



 glColor3f(0.0, 0.0, 1.0);

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

    p[0] = 0.0;

    p[1] = k;

    glBegin(GL_POINTS);

    glVertex2fv(p);

    glEnd();

  }



  glColor3f(0.0, 0.0, 1.0);           //horizontal axis

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

    p[0] = k;

    p[1] = 0.0;

    glBegin(GL_POINTS);

    glVertex2fv(p);

    glEnd();

  }



  //tick marks

  for(k=-250; k<250; k+=10){         //horizontal tick marks

    for(j=-5; j<5; j++){

      p[0] = k;

      p[1] = j;

      glBegin(GL_POINTS);

      glVertex2fv(p);

      glEnd();

    }

    

  }



  for(k=-250; k<250; k+=10){         //vertical tick marks

    for(j=-5; j<5; j++){

      p[0] = j;

      p[1] = k;

      glBegin(GL_POINTS);

      glVertex2fv(p);

      glEnd();

    }

    

  }



}







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

//while calculating data

void displayOrbit (){



  typedef GLfloat point2[2];          //define a point array

  long int k;                         //loop control

  point2 p = {0.0, 0.0};              //begin here



  glClear(GL_COLOR_BUFFER_BIT);       //clear window



  drawAxis();



  glColor3f(1.0, 0.0, 0.0);           //plot color (red)

  

  for(k=0; k < loop; k++){            //plot data

    p[0] = displayData[0][k]*100;     //x position (scale data so can see better...)

    p[1] = displayData[1][k]*100;     //y position

    glBegin(GL_POINTS);               //plot points

    glVertex2fv(p);

    glEnd();

  }



}







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

//while calculating data

void displayEnergy (){



  typedef GLfloat point2[2];          //define a point array

  long int k;                         //loop control

  point2 p = {0.0, 0.0};              //begin here



  glClear(GL_COLOR_BUFFER_BIT);       //clear window

  

  drawAxis();



  glColor3f(1.0, 0.0, 0.0);           //plot color (red)

  for(k=0; k < loop; k++){            //plot data

    p[0] = dataArray[1][k];           //energy

    p[1] = dataArray[3][k];           //time

    glBegin(GL_POINTS);               //plot points

    glVertex2fv(p);

    glEnd();

  }



}





void displayL (){



  typedef GLfloat point2[2];          //define a point array

  long int k;                         //loop control

  point2 p = {0.0, 0.0};              //begin here



  glClear(GL_COLOR_BUFFER_BIT);       //clear window

  

  drawAxis();



  glColor3f(1.0, 0.0, 0.0);           //draw color



  for(k=0; k < loop; k++){            //plot data

    p[0] = dataArray[2][k];           //L

    p[1] = dataArray[3][k];           //time

    glBegin(GL_POINTS);               //plot points

    glVertex2fv(p);

    glEnd();

  }



}











//initialize opengl window for plotting

void myinit(void){



  glClearColor(1.0, 1.0, 1.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

  gluOrtho2D(-250.0, 250.0, -250.0, 250.0);  //2d coordinate system centered at (0,0)

  glMatrixMode(GL_MODELVIEW);          //model view stack 



}







//open 3 windows - one each for position, energy and angular momentum

plotData(int argc, char **argv){





     glutInit(&argc, argv);         //plot position

	glutInitWindowSize(500,500);

	glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);

	(void)glutCreateWindow("Orbit Trajectory");

	glutDisplayFunc(displayOrbit);

	myinit();

	

       	glutInitWindowSize(500,500);  //plot energy

	glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);

	(void)glutCreateWindow("Energy");

	glutDisplayFunc(displayEnergy);

	myinit();

	        

       	glutInitWindowSize(500,500);  //plot angular momentum

	glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);

	(void)glutCreateWindow("Angular Momentum");

	glutDisplayFunc(displayL);

	myinit();

     glutMainLoop();



}







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



        double velocityX, velocityY;   //initial velocity

	double positionX, positionY;   //initial position

	double dataholding[N+1][loop]; //saved data position and velocity



	//set initial values

	velocityX = 0;

	velocityY = 2.0*pi;

	positionX = 1;

	positionY = 0;

		

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



        calculateOrbit(velocityX, velocityY, positionX, positionY, dataholding);



	printData(dataholding);        //save data to file 'data'

       

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



	

}



#undef NRANSI