#include <stdio.h>

/* ---------------------------------------------------------------------
   *** Simulation of the motion of a bicycle racer in realistic case ***

   The velocity of the racer is governed by the first order
   differential equation:

   Force acting on the bicycle = effort from racer -  air friction

                       m dv/dt =  P/v - CrAv^2
   or
                        dv/dt  = P/mv - CrAv^2/m

   where P is the power of the racer          (watt)
         m is the mass of the racer+bicycle,  (kg)
         C is the drag coefficient,           (unitless)
         r is the density of air and          (kg/m^3)
         A is the frontal area of the racer.  (m^2)

   Typical values are given in the program.

   Numerical solution method is second order Runge-Kutta Method.
   Ahmet Bingul, 25/01/2005
   --------------------------------------------------------------------- */

float  v,P,m,C,r,A;
float  t,dt,tm;
float  k1,k2;

/* Function */
float f(float v)
{
  return (P/(m*v) - C*r*A*v*v/m);
}

void main(void)
{
   /* Typical parameters  */
   P  = 400.00;  /* watt  */
   m  =  70.00;  /* kg    */
   C  =   0.50;  /* -     */
   r  =   1.29;  /* kg/m3 */
   A  =   0.33;  /* m2    */

   v  =   4.0;   /* initial velocity */
   t  =   0.0;   /* initial time     */
   dt =   1.0;   /* time step, s     */
   tm = 200.0;   /* maximum time, s  */


   /* solution apply second order Runge Kutta Method */
   while(t<tm)
   {
     k1 = dt*f(v);
     k2 = dt*f(v+dt*k1);
     v  = v + (k1+k2)/2.0;
     t  = t + dt;
     printf("%f\t%f\n",t,v);
   }

} /* main */