c
c Simulation of velocity vs. time for a bicyclist
c        - Euler or 2nd order Runge-Kutta
c Program to accompany "Computational Physics" by N. Giordano and H. Nakanishi
c Fortran version written by H. Nakanishi, writes to files
c
	program bicycle
c
c       Declare the arrays we will need
c
	dimension t(5003), velocity(5003)
	real mass
	character ans,yes
	character*12 ofile
	yes='y'
c
c       Use subroutines to do the work
c
20	continue
	call initialize(t,velocity,dt,power,mass,nmax,c,area,density,
     c		method,ofile)   
	call calculate(t,velocity,dt,power,mass,nmax,c,area,density,
     c		method)
	call record(t,velocity,nmax,ofile)
	print *,'another try?'
	read(5,10) ans
10	format(a1)
	if(ans.eq.yes) go to 20
	stop
	end
c
	subroutine initialize(t,velocity,dt,power,mass,nmax,
     c		c,area,density,method,ofile)
c
c       Initialize variables
c
	dimension t(1),velocity(1)
	real mass
	character*12 ofile
	print *,'Euler (1) or Runge-Kutta (2)? -> '
	read(5,*) method
	if(method.ne.1.and.method.ne.2) then
		print *,'must select 1 or 2 ..'
		stop
	endif
	print *,'initial velocity -> '
	read(5,*) velocity(1)
	t(1) = 0
	print *,'time step -> '
	read(5,*) dt 
	print *,'max time -> '
	read(5,*) tmax
	nmax=min(int(tmax/dt),5000)
	print *,'constant power -> '
	read(5,*) power
	print *,'c -> '
	read(5,*) c
	mass = 70
	area = 0.33
	density = 1.29
	print *,'mass = 70 kg, area = 0.33 m^2, c =',c
	print *,'output file name?'
	read(5,14) ofile
14	format(a12)
	return
	end
c
	subroutine calculate(t,v,dt,power,mass,nmax,c,area,density,
     c		method)
	real mass,f
c
c       Now use the Euler method or the Runge-Kutta (2nd order)
c
	dimension t(1),v(1)
	if(method.eq.1) then
		do 10 i = 1,nmax-1
		v(i+1)=v(i)+dt
     c			*f(v(i),power,mass,c,density,area)
		t(i+1) = t(i) + dt
10		continue
	else
		do 30 i = 1,nmax-1
		v1=v(i)+0.5*dt*f(v(i),power,mass,c,density,area)
	        v(i+1)=v(i)+dt*f(v1,power,mass,c,density,area)
		t(i+1) = t(i) + dt
30		continue
	endif
	return
	end
c
	function f(v,p,x,c,rho,a)
	f=(p/v-c*rho*a*v**2)/x
	return
	end
c
	subroutine record(t,velocity,n,ofile)
	character*12 ofile
	dimension t(1),velocity(1)
	open(1,file=ofile)
	do 10 i=1,n
	write(1,15) t(i),velocity(i)
10	continue
15	format(1x,2e12.5)
	close(1)
	return
	end
