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, need to be compiled with "-lpepl".
c
	program bicycle2
c
c       Declare the arrays we will need
c
	dimension t(5003), velocity(5003)
	real mass
c
c       Use subroutines to do the work
c
	call initialize(t,velocity,dt,power,mass,nmax,c,area,density,
     c		lsym,nsym,method)   
	call calculate(t,velocity,dt,power,mass,nmax,c,area,density,
     c		lsym,nsym,method)   
	call display(t,velocity,dt,power,mass,nmax,c,area,density,
     c		lsym,nsym,method)
	stop
	end
c
	subroutine initialize(t,velocity,dt,power,mass,nmax,
     c		c,area,density,lsym,nsym,method)
c
c       Initialize variables
c
	dimension t(1),velocity(1)
	real mass
	character ans,yes
	yes='y'
	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
	mass = 70
	c = 0.5
	area = 0.33
	density = 1.29
	print *,'set line, symbol?'
	read(5,14) ans
14	format(a1)
	if(ans.eq.yes) then
		print *,'line and symbol numbers -> '
		read(5,*) lsym,nsym
	else
		lsym=-1
		nsym=1
	endif
	return
	end
c
	subroutine calculate(t,v,dt,power,mass,nmax,c,area,density,
     c		lsym,nsym,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 display(t,velocity,dt,power,mass,nmax,c,area,density,
     c		lsym,nsym,method)
c
c       First set up title and label axes for graph. Plotting is a lot
c       of work in fortran.
c       This version displays output as well as writes to a file "graph.out".
c
	dimension t(1),velocity(1)
	call usrmon(.true.,.false.,-0.5,9.,-0.5,12.)
	call pltlun(19,.true.,.false.)
	call pltlfn('graph.out')
	call plots
	call plot(0.,3.5,-3)
	if(method.eq.1) then
	call text(0.2,5.8,0.2,'Bicyclist in air: Euler',0.,24,0)
	else
	call text(0.2,5.8,0.2,'Bicyclist in air: Runge-Kutta',0.,29,0)
	endif
	call text(0.2,5.5,0.2,'Velocity versus Time',0.,20,0)
	call scalex(t,5.,nmax,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,-1)
	call axisx(0.,0.,'Time(s)',-7,5.,0.,
     c	t(nmax+1),t(nmax+2),t(nmax+3),4)
	call axisx(0.,5.,' ',1,5.,0.,
     c	t(nmax+1),t(nmax+2),t(nmax+3),20)
	call scalex(velocity,5.,nmax,1)
	call axctl(0.15,0.04,0.15,0.2,0.2,-1)
	call axisx(0.,0.,'Velocity (m/s)',16,5.,90.,
     c	velocity(nmax+1),velocity(nmax+2),velocity(nmax+3),-5)
	call axisx(5.,0.,' ',-1,5.,90.,
     c	velocity(nmax+1),velocity(nmax+2),velocity(nmax+3),-21)
	call line(t,velocity,nmax,1,lsym,nsym)
	call number(2.,3.0,0.2,power,0.,'''power = '',f6.1')
	call number(2.,3.5,0.2,dt,0.,'''dt = '',f5.2')
	call plot(0.,0.,999)
	return
	end