program leigh
c
c	implementation of Leigh's 1955 paper
c
	implicit real *8 (a-h,o-z)
	parameter (ix=400,iy=80)
	real *8 uy(0:ix),w(0:iy),u(0:ix,0:iy),dp(0:ix),ue(0:ix)
	real *8 v(0:ix,0:iy)
	real *8 a(iy,iy),z(iy)
	real *8 h,dx
	iterations=10
	read(5,*)xs,h
	dx=xs/ix
	ix_print=ix/5
	pi=3.14159
	h2=1./(h*h)
	do i=0,ix
		x=i*dx
		ue(i)=1.-x/8
c		ue(i)=2*sin(pi*x)
c		ue(i)=1-x*x/8
		dp(i)=ue(i)/8
c		dp(i)=-ue(i)*2*pi*cos(pi*x)
c		dp(i)=ue(i)*x/4
		do j=0,iy
			y=j*h
			if(y.lt.3.1519/2)then
				u(i,j)=ue(i)*sin(y)
			else
				u(i,j)=ue(i)
			endif
c			u(i,j)=ue(i)*(1.-exp(-j*h))
		enddo
	
	enddo
c
c	starting from approximate flow - remember
c
	do i=1,ix
	write(6,5)i
5	format('i=',i5)
		pg=dp(i)+dp(i-1)
		do j=0,iy
		     w(j)=2*u(i-1,j)
		enddo
		do k=1,iterations
	           do i1=1,iy
	             z(i1)=0.
		     do j=1,iy
			a(i1,j)=0.
		     enddo
                   enddo
		   s=0.
		   i1=1
		   dw=(w(i1+1)-w(i1-1))/2
		   a(i1,i1)=-2*dx+h*h*dw/2-w(i1)*h*h
		   a(i1,i1+1)=dx
                   z(i1)=h*h*dx*pg-2*u(i-1,i1)*w(i1)*h*h+
     +                                 (s+u(i-1,i1))*dw*h*h
		   s=s+2*u(i-1,i1)
		   
		   do i1=2,iy-1
		      dw=(w(i1+1)-w(i1-1))/2
		      do j=1,i1-1
                        a(i1,j)=dw*h*h
		      enddo
                      a(i1,i1-1)=a(i1,i1-1)+dx
		      a(i1,i1)=-2*dx + h*h*dw/2 - w(i1)*h*h
		      a(i1,i1+1)=dx
                      z(i1)=h*h*dx*pg-2*u(i-1,i1)*w(i1)*h*h+
     +                                 (s+u(i-1,i1))*dw*h*h
		      s=s+2*u(i-1,i1)
		   enddo
		   i1=iy
		   a(i1,i1)=1.
		   z(i1)=ue(i)+ue(i-1)

		   call invert(a,iy,iy,z)

		   diff=0.

		   do j=1,iy
		      diff=diff+abs(z(j)-w(j))
		      w(j)=z(j)
		   enddo
		   diff=diff/(iy*(ue(i)+ue(i-1)))
c		   write(6,10)i,diff
c10	           format('looping:',i5,f15.6)
		   
		enddo
		do j=1,iy
		  u(i,j)=w(j)-u(i-1,j)
		enddo
	enddo
	do i=0,ix,ix_print
	write(7,20)(j*h,u(i,j),j=0,iy)
20	format('data',/,(f10.4,f15.6))
	
	enddo
c
c	output wall shear
c
	write(8,20)(i*dx,(u(i,1)-(1-i*dx/8)*h*h/16)/h,i=0,ix)
c
c	calculte v velocity at midpoints
c
	do i=1,ix
		v(i,0)=0.
		do j=1,iy
		  du0=(u(i,j-1)-u(i-1,j-1))/dx
		  du1=(u(i,j)-u(i-1,j))/dx
		  v(i,j)=v(i,j-1)-0.5*h*(du0+du1)
		enddo
	enddo
	do j=0,iy,iy/5
	write(9,20)((i-0.5)*dx,v(i,j),i=1,ix)
	enddo
	stop
	end