program lcwcyl
c
c	Les Wood's method for free streamline
c	solution toflow about cylinder
c
	real *8 g(0:100),q(0:100),c(0:100),b(0:100,0:100),om(0:100)
	real *8 p(0:100),phi(0:100),xi(0:100),th(0:100)
	real *8 phis,dg,pi,angle
	integer iterations,i
	complex *8 ey,z(0:10,-50:200),a1,a2,a3,omz(-50:200)
	pi=4*atan(1.)
	iterations=20
	write(6,1)
1	format('Input separation angle in degrees')
	read(5,*)angle
	xis=angle*pi/180.
	n=40
	dg=pi/n
c
c	set up initial guess
c
	do i=0,n
		g(i)=i*dg
		p(i)=(sin(g(i)/2))**2
		q(i)=1.
		do j=0,n-1
			gj=(j+0.5)*dg
			arg1=(gj-g(i))/4
			arg2=(gj+g(i))/4
			fac=sin(arg1)*sin(arg2)/(cos(arg1)*cos(arg2))
			b(i,j)=-log(abs(fac))*dg/pi
		enddo
	enddo
	phis=xis
	oms=-.4077
	tau=pi
	ak=0.
c
c     loop a fixed number of times, convergence usually occurs in
c	only a few iterations
c
	do loop=1,iterations
		do i=0,n-1
			gm=(g(i)+g(i+1))/2
			qm=(q(i)+q(i+1))/2
			c(i)=0.5*phis*sin(gm)/qm
		enddo
		do i=1,n
			
			f1=tau*log(abs(tan(g(i)/4)))/pi
			om(i)=oms/(1.+ak*cos(g(i)/2))-f1
			do j=0,n-1
				om(i)=om(i)-b(i,j)*c(j)
			enddo
			
			q(i)=exp(-om(i))
		enddo
		q(0)=0
		s1=0
		do i=0,n-1
			s1=s1+c(i)
		enddo
c		tau=pi-2*xis+2*pi*s1/n
		s2=0
		do i=0,n-1
			gm=(i+0.5)*dg
			qm=(q(i)+q(i+1))/2
			s2=s2+dg*sin(gm)/qm
		enddo
		phis=2*xis/s2
		do i=0,n
			phi(i)=phis*p(i)
		enddo
		xi(0)=0
		th(0)=pi/2
		do i=1,n
			qm=0.5*(q(i-1)+q(i))
			xi(i)=xi(i-1)+(phi(i)-phi(i-1))/qm
			th(i)=pi/2-xi(i)
		enddo
		s3=0.
		do i=0,n-1
			thm=0.5*(th(i)+th(i+1))
			gm=0.5*(g(i)+g(i+1))
			s3=s3+thm*sin(gm/2)
		enddo
		s3=s3*dg/pi
		ak=-oms/s3
		write(6,*)phis,tau,ak

	enddo
	write(6,5)(i,180*th(i)/pi,i=0,n)	
5	format(i5,f15.6)	
	write(7,10)(180*xi(i)/pi,1-q(i)*q(i),i=0,n)
10	format(2f15.6)
	stop
	end