program pohlhausen
c
c	implement pohlhausen 1928 solution
c
	implicit real *8(a-h,o-z)
	real *8 cp(0:200),ue(0:200),due(0:200),th(0:200)
	real *8 ddue(0:200),d(0:200),z(0:200)
	real *8 pi
	pi=3.14159
	radius=9.75/2
c	radius in cm
c	velocity scale 19 cm/sec
	uscale=19.
	u1=7.151*radius/uscale
	u3=-0.04497*radius**3/uscale
	u5=-0.000330*radisu**5/uscale
c	u_e(alpha )=u1*alpha+u3*alpha^3+u5*alpha^5	
c
c	Pohlhausens result: separation at x=6.94cm
c
	acrit=6.94*180/(pi*radius)
	write(6,1)acrit
1	format('Pohlhasuen gets separation at alpha = ',f10.4,' degrees')
	amax=120*pi/180
	da=amax/100
	do i=0,100
		alpha=i*da
		a2=alpha*alpha
		th(i)=180*alpha/pi
		ue(i)=alpha*(u1+a2*(u3+a2*u5))
		due(i)=u1+a2*(3*u3+a2*5*u5)
		ddue(i)=alpha*(6*u3+20*u5*a2)
		cp(i)=1-ue(i)**2
	enddo
	write(7,10)(th(i),ue(i),i=0,100)
	write(8,10)(th(i),cp(i),i=0,100)
	write(9,10)(th(i),due(i),i=0,100)
10	format(f8.2,f15.6)
	c0=-2
	c1=116./315
	c2=-1./7560
	c3=-1./4536
	f0=-37./315
	f1=1./315
	f2=5./9027
	p0=c0
	p1=c1*due(0)
	p2=c2*79*due(0)**2
	p3=c3*due(0)**3
	write(6,5)p3/p3,p2/p3,p1/p3,p0/p3
5	format('Input solution of:',f15.5,'x^3+',f15.5,
     +        'x^2+',f15.5,'x+',f15.5)
	read(5,*)z(0)
	d(0)=sqrt(z(0))
	z(1)=z(0)
	d(1)=d(0)
	alpha=0.
	do i=1,100
	  alpha=alpha+da
	  u=(ue(i)+ue(i-1))/2
	  du=(due(i)+due(i-1))/2
	  ddu=(ddue(i)+ddue(i-1))/2
	  z(i)=z(i-1)+(2./u)*g1(alpha,d(i-1),u,du,ddu)*da
	  d(i)=sqrt(z(i))
	write(6,2)i,z(i-1),d(i-1),u,du,ddu,z(i),d(i)
2	format(i5,7f12.3)
	enddo
	write(10,10)(th(i),d(i),i=0,100)
	write(11,10)(th(i),due(i)*z(i),i=0,100)
	stop
	end
	real *8 function g1(a,d,u,du,ddu)
	implicit real*8 (a-h,o-z)
	gn=-2+d*d*(116*du/315-d*d*((79*du*du+8*u*ddu)/7560
     +     +d*d*du*(du*du+u*ddu)/4536))
	gd=-37./315+d*d*(du/315+5*d*d*du*du/9072)
	g1=gn/gd
	return
	end