PROGRAM Bicycle_Racer
!------------------------------------------------------------------------!
! *** 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)                       !
!                                                                        !
! Numerical solution method is second order Runge-Kutta Method.          !
! Ahmet Bingul, 25/01/2005                                               !
!------------------------------------------------------------------------!
IMPLICIT NONE
REAL    :: v,P,m,C,r,A
REAL    :: t,dt,tm
REAL    :: k1,k2,f


!-- Function --!
f(v) = P/(m*v) - C*r*A*v**2/m

!-- Typical parameters --!
 P  = 400.00  ! watt
 m  =  70.00  ! kg
 C  =   0.50  ! drag coefficient
 r  =   1.29  ! density of air
 A  =   0.33  ! frontal area

 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 --!

DO WHILE(t<tm)

   k1 = dt*f(v)
   k2 = dt*f(v+dt*k1)
   v  = v + (k1+k2)/2.0
   t  = t + dt

   PRINT *,t,v

END DO

END PROGRAM