program burgers
c
c	burgers method with heimenz data and burgers
c	solution from potential outer flow
	implicit real *8 (a-h,o-z)
	real *8 w(0:1000),th(0:1000),a(0:1000),ah(0:1000)
	real *8 wb(0:1000)
        pi=3.14159
	n=450
	da=0.005
	w(0)=0.
	do i=0,n
		th(i)=i*da*180/pi
		a(i)=i*da
		ah(i)=(i-0.5)*da
	enddo
c
c	generate potential solution
	do i=0,n
		phi=-2*cos(a(i))
		wb(i)=8*sqrt(phi+2)*(1-phi)/(3*sqrt(pi))
	enddo
	write(8,10)(th(i),wb(i),i=0,n)
	do i=1,n
		w(i)=0.
		g0=pf(a(i))/qf(a(i),a(i))
		do j=1,i
		w(i)=w(i)+(g1(a(i),ah(j))-g0)*da/sqrt(ah(j))
		enddo
		g0=pf(a(i))/qf(a(i),a(i))
		w(i)=2*sqrt(a(i))*g0+w(i)
		w(i)=2*w(i)/sqrt(pi)
		write(6,10)th(i),w(i)
10	format(f10.2,f15.6)
	enddo
	write(7,10)(th(i),w(i),i=0,n)
	stop
	end


	real *8 function ueheimenz(a)
	implicit real *8 (a-h,o-z)
        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*radius**5/uscale
	ueheimenz=a*(u1+a*a*(u3+a*a*u5))
	return
	end

	real *8 function dueheimenz(a)
	implicit real *8 (a-h,o-z)
        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*radius**5/uscale
	dueheimenz=u1+a*a*(3*u3+a*a*5*u5)
	return
	end

	real *8 function pf(a)
	implicit real *8 (a-h,o-z)
	pf=ueheimenz(a)*dueheimenz(a)
	return
	end

	real *8 function qf(a,ad)
	implicit real *8 (a-h,o-z)
        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*radius**5/uscale
	qf=u1*(a+ad)/2+u3*(a**3+a*a*ad+a*ad*ad+ad**3)/4+
     +       u5*(a**5+a**4*ad+a**3*ad*ad+a*a*ad**3+a*ad**4+ad**5)/6
	qf=sqrt(qf)  
	return
	end

	real *8 function g1(a,ad)
	implicit real *8 (a-h,o-z)
	g1=pf(a-ad)/qf(a,a-ad)
	return
	end