program bpgen
c
c	program to integrate (10.17) to obtain
c       functions zeta (z here) gamma (g here)
c       and kappa (kappa here) from Borgas and Pedley's
c       variant of Falkner-Skan equation
c
c       NOTE typo in equation (10.19) page 295
c	factor of z should multiply gamma
c


	implicit real *8 (a-h,o-z)
	real *8 eta(0:2000),f(0:2000),df(0:2000),ddf(0:2000)
	real *8 m,deta,kappa(500),z(500),g(500)
	real *8 a2(500),dfl(500)

c
c	Falkner-Skan parameter b=0.5
c
	b=0.5
c
c	input range of tau for which solution calculated
c
	write(6,3)
3	format('Input calculation range:tau_start, tau_end')
	read(5,*)a2start,a2end

	na2=101
	da2=(a2end-a2start)/(na2-1)
	do k=1,na2
		a2(k)=a2start+(k-1)*da2
		deta=0.01
		n=1000
		call rkped(deta,n,eta,f,df,ddf,b,a2(k))
c
c		use limiting values
c		zeta  ~ f''(infinity)
c		gamma ~ f'-z*zeta
c		kappa ~ equation (10.24)
c
		z(k)=ddf(n)
		g(k)=df(n)-eta(n)*z(k)
		kappa(k)= - (3./8)**(1./3)*g(k)*z(k)**(-2./3.)
		write(6,10)a2(k),z(k),g(k),kappa(k)
10	        format(4f15.6)
c
c	output to files fort.7-fort.9
c
	        write(7,10)a2(k),z(k)
	        write(8,10)a2(k),g(k)
	        write(9,10)a2(k),kappa(k)
	enddo
c
c	have now calculated Borgas and Pedleys functions
c

	stop
	end
	subroutine rkped(deta,n,eta,f,df,ddf,b,a2)
	implicit real *8 (a-h,o-z)
	real *8 eta(0:2000),f(0:2000),df(0:2000),ddf(0:2000),a2
c
c	Runge-Kutta solver
c
		eta(0)=0.
		f(0)=0
		df(0)=0
		ddf(0)=a2
		do i=1,n
		  fi=f(i-1)
		  dfi=df(i-1)
	  	  ddfi=ddf(i-1)

		  fa=deta*dfi
		  dfa=deta*ddfi
		  ddfa=deta*(b*(dfi*dfi+1.)-fi*ddfi)

		  fi=f(i-1)+fa/2
		  dfi=df(i-1)+dfa/2
		  ddfi=ddf(i-1)+ddfa/2

		  fb=deta*dfi
		  dfb=deta*ddfi
		  ddfb=deta*(b*(dfi*dfi+1.)-fi*ddfi)

		  fi=f(i-1)+fb/2
		  dfi=df(i-1)+dfb/2
		  ddfi=ddf(i-1)+ddfb/2

		  fc=deta*dfi
		  dfc=deta*ddfi
		  ddfc=deta*(b*(dfi*dfi+1)-fi*ddfi)

		  fi=f(i-1)+fc
		  dfi=df(i-1)+dfc
		  ddfi=ddf(i-1)+ddfc

		  fd=deta*dfi
		  dfd=deta*ddfi
		  ddfd=deta*(b*(dfi*dfi+1)-fi*ddfi)
		  eta(i)=eta(i-1)+deta
		  f(i)  =  f(i-1)+( fa+2* fb+2* fc+ fd)/6
		  df(i) = df(i-1)+(dfa+2*dfb+2*dfc+dfd)/6
		  ddf(i)=ddf(i-1)+(ddfa+2*ddfb+2*ddfc+ddfd)/6
		enddo
	return
	end