// ***************************
// **   Pendulum animation  **
// ** Computing 2 Lecture 3 **
// **  Ver 1.0 02/02/1998   **
// ***************************

#pragma windows 
#include <stdio.h>
#include <math.h>
#include <windows.h>
#include <DBOS\graphics.h>

#define g 9.8

float d_omega_dt(float); 
int update_func();

float theta[2000];                     

int main(void)
       {     
         float omega[2000],time=0.0,time_step=0.2;
         int t;
          
         theta[0]=M_PI/4.0;             // Starting angle 
         omega[0]=0.00;                 // Starting velocity is 0 rad/s
                                   
         for (t=0;t<1999;t++)
           {                  
             omega[t+1]=omega[t]+d_omega_dt(theta[t])*time_step; 
             theta[t+1]=theta[t]+(time_step*omega[t+1]);         
             
             // NB Euler-Cromer method used to calculate theta[t+1].
             // That is, the NEW value of omega (omega[t+1]) is used to
             // calculate theta[t+1]. Euler's method uses the previous
             // value of omega (i.e. omega[t]).
  
             // Modify the program so that Euler's method is used - you'll
             // see a dramatic difference in the behaviour of the
             // pendulum.See Giordano, Appendix 1 and A. Cromer, "Stable Solutions using the Euler
             // Approximation", Am. J. Phys. vol. 49 p. 455 (1981) for details of the
             // Euler-Cromer method.
                                                                                                                                                                
             time+=time_step;
            }
         
         // Graphics commands - see lab. session notes
         
         winio("%ca[Simple pendulum: Euler-Cromer method]%gr[black]&",500,500);
         winio("%dl",0.1,update_func); 
    }
         
// Function to calculate dv/dt using
// formula dv/dt = -(g/l)*sin(theta)

float d_omega_dt(float angle)
        {         
      
          float l=9.8; 
          float dveldt;       

          dveldt=-(g/l)*sin(angle);                    

          return(dveldt);
         } 
         
// ***********************************
// Function to update graphics display
// ***********************************

int update_func()
{
   long l;
   short errcode;
   static int i=0;
   static int xdist=196*sin(theta[0]),ydist=196*cos(theta[0]);
   
   // i, xdist and ydist are declared as static variables
   // NB..Declaring a variable as static means that the value 
   //     of the variable is "remembered" from function call
   //     to function call. Otherwise the variable is reset
   //     every time the function is called. Try removing
   //     the "static" declaration from before "int xdist="
   //     and observe what happens to the graphics display!
   
   // "blank out" pendulum at old position
   
   draw_line(250,250,250+xdist,250+ydist,0);
   fill_ellipse(250+xdist,250+ydist,20,20,0);
   
   // calculate new position of pendulum from value of theta
   
   i++;
   xdist=196*sin(theta[i]);
   ydist=196*cos(theta[i]);
    
   // draw pendulum at new position
   
   draw_line(250,250,250+xdist,250+ydist,10);
   fill_ellipse(250+xdist,250+ydist,20,20,10);
   
   if (i!=1999)
   return 2;        // If update_function() returns 2 the
   else             // Windows operating system keeps the window open
   return 0;        // If update_function() returns 0 the window is closed
}