program terrill
c
c	implementation of Terrills 1960 paper
c
	implicit real *8 (a-h,o-z)
	parameter (ix=800,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),yr(0:ix,0:iy)
	real *8 phi(0:ix),b(0:ix),x(0:ix)
	real *8 h,xs,dphi,phis,h2
	integer method,plate,cylinder,heimenz
	omega=0.01
	plate=1
	cylinder=2
	heimenz=3
		u1=1.8348
		u3=-0.2742
		u5=-0.0478

	iterations0=500
	write(6,4)
4	format('Implementation of Terrills method',/,
     +         'Input separation distance and eta step')
	read(5,*)xs,h
	write(6,3)
3	format('Input method:')
	read(5,*)method
	if(method.eq.plate)then
	   phis=xs*(1-xs/16)
	else if(method.eq.cylinder)then
	   phis=2*(1-cos(xs))
	else
	   phis=xs*xs*(u1/2+xs*xs*(u3/4+xs*xs*u5/6))
	endif
	dphi=phis/ix
	ix_print=ix/5
	if(ix_print.eq.0)ix_print=1
	pi=3.14159
	h2=1./(h*h)
	if(method.eq.plate)then
	   istart=0
	   do i=istart,ix
		g=i*dphi
		phi(i)=2*i-1.
		x(i)=8-sqrt(64-16*g)
		ue(i)=sqrt(1.-g/4)
		b(i)=-g/(4*(1-g/4))
		do j=0,iy
			y=j*h
			u(i,j)=1.-exp(-y)
			yr(i,j)=y*sqrt(2.*i*dphi)/ue(i)
	        enddo
	  enddo
	else if(method.eq.cylinder)then
           istart=0
	   do i=0,ix
		g=i*dphi
		phi(i)=2*i-1.
		x(i)=acos(1.-g/2)
		ue(i)=2*sin(x(i))
		b(i)=cos(x(i))/(cos(x(i)/2)**2)
		do j=0,iy
			y=j*h
			u(i,j)=1.-exp(-y)
			yr(i,j)=y*sqrt(2*i*dphi+0.00000002)/(ue(i)+0.0001)
	        enddo
	  enddo

	else
c
c	heimenz
c
		istart=0
		u1=1.8348
		u3=-0.2742
		u5=-0.0478
	   do i=0,ix
		g=i*dphi
		phi(i)=2*i-1.
		call cubic(u1,u3,u5,g,root)
		x(i)=root
		ue(i)=x(i)*(u1+x(i)**2*(u3++u5*x(i)**2))
		due=u1+x(i)**2*(3*u3+5*u5*x(i)**2)
		b(i)=2*g*due/ue(i)**2
		do j=0,iy
			y=j*h
			u(i,j)=1.-exp(-y)
			yr(i,j)=y*sqrt(2*i*dphi+0.00000002)/(ue(i)+0.0001)
	        enddo
	  enddo
	  b(0)=1.0
	  

	endif
	phi(0)=0.
	write(6,99)((yr(i,j),j=0,iy,iy/5),b(i),ue(i),i=0,ix,ix/10)
99	format(8f15.3)
	read(5,*)omega
c
c	starting from approximate flow - remember
c
	do i=istart,ix
	write(6,5)i
5	format('i=',i5)
		if(i.gt.0)then
		   do j=0,iy
		       w(j)=2*u(i-1,j)
		   enddo
	        else
		   do j=0,iy
                      w(j)=2*u(0,j)
		   enddo
                endif
		if(i.le.istart)then
			iterations=5*iterations0
		else
			iterations=iterations0
		endif
		do k=1,iterations
	           do i1=1,iy
	             z(i1)=0.
		     do j=1,iy
			a(i1,j)=0.
		     enddo
                   enddo
c------------------------------------------------------------------------------		   
	i1=1
		   if(i.gt.0)then
			s=0.5*u(i-1,i1)
			q1=u(i-1,i1)
			q2=w(i1)-q1
			pg=h*h*(b(i-1)*(1.-q1**2)+b(i)*(1.-q2**2))
		   else
			s=0.
			pg=h*h*b(0)*(1-w(i1)*w(i1))
			q1=0.
		   endif
		   dw=(w(i1+1)-w(i1-1))/2
		   a(i1,i1)=-2.+(0.5+phi(i))*h*h*dw/2
     +                        -phi(i)*w(i1)*h*h
		   a(i1,i1+1)=1.
                   z(i1)=2*phi(i)*(s*dw-w(i1)*q1)*h*h-pg
c     +                         +phi(i)*w(i1)*h*h

	do i1=2,iy-1
		      if(i.gt.0)then
			s=s+0.5*(u(i-1,i1)+u(i-1,i1-1))
			q1=u(i-1,i1)
			q2=w(i1)-q1
			pg=h*h*(b(i-1)*(1.-q1**2)+b(i)*(1.-q2**2))
			
		      else
			s=0.
			pg=h*h*b(0)*(1-w(i1)*w(i1))
			q1=0.
		      endif
		      dw=(w(i1+1)-w(i1-1))/2
		      do j=1,i1-1
                            a(i1,j)=(0.5+phi(i))*dw*h*h
		      enddo
                      a(i1,i1-1)=a(i1,i1-1)+1.
		      a(i1,i1)=-2.+(0.5+phi(i))*h*h*dw/2
     +                        -phi(i)*w(i1)*h*h
		      a(i1,i1+1)=1.
                      z(i1)=2*phi(i)*(s*dw-w(i1)*q1)*h*h-pg
c     +                       +phi(i)*w(i1)*h*h  
	enddo
	i1=iy
		   a(i1,i1)=1.
		   z(i1)=2.
c----------------------------------------------------------------------------

		   call invert(a,iy,iy,z)

		   diff=0.
		   do j=1,iy
		      diff=diff+abs(z(j)-w(j))
		      w(j)=(1-omega)*w(j)+omega*z(j)
		   enddo
		   diff=diff/iy
		   write(6,10)i,diff
10	           format('+looping:',i5,f15.6)
		enddo
		if(i.gt.0)then
		  do j=1,iy
		    u(i,j)=w(j)-u(i-1,j)
		  enddo
		else
		  do j=1,iy
		    u(i,j)=w(j)/2
		  enddo
		endif
	enddo
	do i=ix_print,ix,ix_print
	   write(7,20)(yr(i,j),ue(i)*u(i,j),j=0,iy)
20	   format('data',/,(f10.4,f15.6))	
	enddo
c
c	output wall shear
c
	write(8,20)(x(i),ue(i)*u(i,1)/yr(i,1),i=1,ix)
c
c	calculte v velocity at midpoints
c
        do i=1,ix
                v(i,0)=0.
		ph=sqrt((i-0.5)*dphi)
		um=0.5*(ue(i)+ue(i-1))
                do j=1,iy
                  du0=um*(u(i,j-1)-u(i-1,j-1))/dphi
                  du1=um*(u(i,j)-u(i-1,j))/dphi
                  v(i,j)=v(i,j-1)-0.5*h*ph*(du0+du1)
                enddo
		
        enddo
        do j=0,iy,iy/5
        write(9,20)((x(i-1)+x(i))/2,v(i,j),i=1,ix)
        enddo

	
	stop
	end
c----------------------------------------------------------------------------------
c  invert
c----------------------------------------------------------------------------------
	subroutine invert(a,n,iy1,z)
	implicit real *8 (a-h,o-z)
	real *8 up(200),low(200,200),c(200)
	real *8 a(iy1,iy1),z(iy1)

	up(1)=a(1,1)
	low(1,1)=1.
	do j=2,n
		low(j,1)=a(j,1)/up(1)
	enddo
	do i=2,n-1
		up(i)=a(i,i)-low(i,i-1)
		low(i,i)=1.
		do j=i+1,n
			low(j,i)=(a(j,i)-low(j,i-1))/up(i)
		enddo
	enddo
	up(n)=a(n,n)-low(n,n-1)

	c(1)=z(1)
	do i=2,n
		c(i)=z(i)
		do j=1,i-1
			c(i)=c(i)-low(i,j)*c(j)
		enddo
	enddo
	z(n)=c(n)/up(n)
	do i=n-1,1,-1
		z(i)=(c(i)-z(i+1))/up(i)
	enddo
	return
	end
	subroutine cubic(a,b,c,g,root)
	implicit real *8 (a-h,o-z)
	r=0.
	tol=0.000001
	do i=1,1000
		s=2*(g-b*r*r/4-c*r*r*r/6)/a
		er=abs(s-r)
		r=s
		if(er.lt.tol)go to 1
	enddo
1	root=sqrt(r)
	
	return
	end