/* Pendulum motion, y''+sen(y)=0. 
compile with "gcc -o main main.c -lm"
*/

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

FILE *gv, *tk, *cf ; 
 
# define DIM        2

# define omega  1.

# define  sqr(x)  ((x)*(x))
# define min(x,y,z) ((x)<(y)? ((x)<(z)? (x):(z)) : ((y)<(z)? (y):(z)))

# include "../lib/pendulum_eqns1.c" 
# include "../lib/runge_kutta_step.c"

/*********************************************************************
 ******************            main        ****************************
 **********************************************************************/

main()   
     
{float y[2],t, dt=.05 ,v; 
char command;

 void start_geomview(void);
 void update_geomview(float);
 void quit(void);

 /*  make pipes to Geomview and the  Tk interface */

 gv=popen("geomview  -c - -wpos 500,500@285,70", "w");
 
 tk=popen("tk","r");
 
 start_geomview();
 
 /* set  initial position and velocity */
 y[0] = 1.; y[1] = 0.;

 /* tell geomview to draw ang vel vector and some other stuff */
 update_geomview(y[0]);

 /* tell geomview to adjust size of picture  */
 fprintf(gv,"(look world)"); fflush(gv);

 while(1) {
   switch (fgetc(tk)) {
     
   case 'v': /* change velocity of integration */
     fscanf(tk, "%f", &v);
     dt=v/20.;
     break;

   case 'p': /* change initial position */
     fscanf(tk, "%f",y);
     update_geomview(y[0]);
     break;
     
   case 'V': /* change initial velocity */
     fscanf(tk, "%f",y+1);
     break;
     
     
   case 'q': 
     quit();
     
   case 'g': /* start integrfation loop */
     command='g';
     for(t=0.;command=='g'; t+=dt)
       {runge_kutta_step(t, y, dt);
       update_geomview(y[0]);
       cf=fopen("/tmp/command_file", "r");
       fscanf(cf,"%c",&command) ;
       fclose(cf);
       if (command=='q') quit();
       if (command=='v') 
	 {fscanf(tk, "%c %f", &command, &v); 
	 dt=v/20.;
	 command='g';
	 };
       }; 
     break;
   };
 };
}

/*********************************************************************
 ****************** start_geomview ***********************************
 **********************************************************************/

void start_geomview()
     
{fprintf(gv,"
(progn
(geometry base {<pendulum_base})
(geometry pendulum
{INST transform:pendulum_transform 
geom <pendulum_geometry})
(backcolor allcams .5 .5 .7)
(normalization allgeoms none)
(bbox-draw world no)
)
");

fflush(gv);
}



/*********************************************************************
****************** update_geomview  *****************************
**********************************************************************/

void update_geomview(float angle)

{fprintf(gv, "
(read transform 
{define pendulum_transform {
%f          %f          0 0", 
cos(angle), sin(angle));

fprintf(gv, "
%f           %f           0 0 \n", 
-sin(angle), cos(angle));

fprintf(gv, " 
0  0 1 0
0  0 0 1 
}})
");

fflush(gv);
}

/*********************************************************************
******************            quit        ****************************
**********************************************************************/
void quit()

{fprintf(gv,"(exit)\n");
pclose(tk); 
fclose(cf);
system("rm -f /tmp/command_file");
printf("BYE\n"); 
exit(0);}