	program lcwcyl
c
c	lcw cylinder solution
c
	real *8 g(0:200),q(0:400),c(0:200),b(0:200,0:200),om(0:200)
	real *8 p(0:200),phi(0:200),xi(0:400),th(0:400),dth(0:200)
	real *8 phis,dg,pi
	complex *8 ey,z(-400:400,0:10),a1,a2,a3,omz,w,zeta,omp
	pi=3.14159
	ey=cmplx(0.,1.)
	write(6,1)
1	format('Input separation angle')
	read(5,*)angle
	write(6,2)
2	format('Input Cp behind body')
	read(5,*)cpw
	qwake=sqrt(1.-cpw)
	oms=-log(qwake)
	write(6,33)oms
33	format('Omega_0 = ',f15.6)
	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
	tau=pi
	ak=0.

	do loop=1,20
		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
c		write(6,*)phis,ak

	enddo
	do i=0,n-1
		dth(i)=th(i+1)-th(i)
	enddo
c
c	output angle theta to fort.9
c
	write(9,39)(-i*dg,-180*th(i)/pi,i=n,0,-1)
	write(9,39)(i*dg,180*th(i)/pi,i=0,n)
39	format('data theta woods',/,(2f15.6))
	write(10,41)(i*dg,q(i),i=0,n)
	write(11,41)(180*(pi/2-th(i))/pi,q(i),i=0,n)
41	format('data speed woods',/,(2f15.6))
c
c	calculate extent of singularity in curvature
c
	gam=0.
		do i=0,n-1
			gm2=(i+0.5)*dg/2
			qm=(q(i)+q(i+1))/2
			gam=gam+dg*sin(gm2)/qm
		enddo
	
	gam=-(0.5-phis*gam/pi)
	critk=sqrt(2*s3*gam)
	critcp=1-exp(2*critk)
	dqs=(q(n)-q(n-1))/dg
	write(6,3)gam,dqs,critk,critcp
3	format(' A = ',f15.6, ' slope= ' ,f15.6,' Kc=',f15.6,' Cpc=',f15.6)
	write(6,4)gam/sqrt(phis)
4	format(' Singularity strength = ',f15.7)
	gam2=gam+ak*oms/2
	write(6,29)gam, ak,oms,ak*oms/2,gam2
29	format('gamma2 = ',5f15.6)
c	write(6,28)(q(i),(q(i)-q(i-1))/dth(i),i=1,n-1)
28	format(2f15.6)
		
c	write(6,5)(i,180*xi(i)/pi,180*th(i)/pi,i=0,n)	
5	format(i5,2f15.6)	
	write(7,10)(180*xi(i)/pi,1-q(i)*q(i),i=0,n)
10	format(2f15.6)
	asum=0.
	do i=0,n-1
		gm2=(i+0.5)*dg/2
		thm=(th(i)+th(i+1))/2
		asum=asum+dg*thm*sin(gm2)
	enddo
	asum=2*sqrt(phis)*asum/pi
	write(6,6)asum
6	format('wake growth = ',f15.6)
c
c	calculate free streamlines
c
	write(8,20)(-sin(i*pi/100),cos(i*pi/100),i=0,200)
20	format('cylinder data',/,(2f12.5))
	omp=cmplx(0.,th(n))
	do i=0,n
		z(i,0)=cmplx(-sin(th(i)),cos(th(i)))
	enddo
c
c	have to calculate starting points as psi varies for phi=0
c
	write(6,25)
25	format('calculating free streamlines')
	dphi=0.1
	n2=n+160
		phiz=phis-0.5*dphi
		do i=n+1,n2
			phiz=phiz+dphi
			w=cmplx(phiz,0.)
			zeta=pi*ey+2*log(sqrt(phiz/phis)-sqrt(phiz/phis-1))
			a1=(1-exp(zeta/2))/(1+exp(zeta/2))
			a1=log(a1)-pi*ey/2
			a2=exp(zeta/2)
			a2=(a2+1./a2)/2
			omz=omp+oms/(1.+ak*a2)-a1
			do j=0,n-1
			  gm=(g(j)+g(j+1))/2
			  a3=sin((ey*zeta+gm)/4)*sin((ey*zeta-gm)/4)
			  a3=a3/(cos((-ey*zeta+gm)/4)*cos((-ey*zeta-gm)/4))
			  a3=log(a3)
			  omz=omz-a3*dth(j)/pi
			enddo
			omz=omz+pi*ey
			z(i,0)=z(i-1,0)+exp(omz)*dphi
c
c	calculate q, cp and theta on free streamline
c
		eta=2*log(sqrt(phiz/phis)+sqrt(phiz/phis-1))
		sinhetahalf=0.5*(exp(eta/2)-exp(-eta/2))
		q(i)=exp(-oms/(1.+ak*ak*sinhetahalf*sinhetahalf))
		xi(i)=pi-atan2(aimag(z(i,0)),real(z(i,0)))
		enddo
		write(7,26)(180*xi(i)/pi,1-q(i)*q(i),i=n+1,n2)
26	format(2f15.6)
		phiz=0.5*dphi
		do i=-1,-n,-1
			phiz=phiz-dphi
			w=cmplx(phiz,0.)
			zeta=2*log(ey*sqrt(w/phis)+sqrt(1.-w/phis))
			a1=exp(zeta/2)
			a1=(1-a1)/(1+a1)
			a1=log(a1)-pi*ey/2
			a2=exp(zeta/2)
			a2=(a2+1./a2)/2
			omz=omp+oms/(1.+ak*a2)-a1
			do j=0,n-1
			  gm=(g(j)+g(j+1))/2
			  a3=sin((gm+ey*zeta)/4)*sin((ey*zeta-gm)/4)
			  a3=a3/(cos((gm+ey*zeta)/4)*cos((gm-ey*zeta)/4))
			  a3=log(a3)
			  omz=omz-a3*dth(j)/pi
			enddo
			z(i,0)=z(i+1,0)+exp(omz)*dphi
		enddo
		write(8,30)(z(i,0),i=-n,n2)
		write(8,30)(real(z(i,0)),-aimag(z(i,0)),i=-n,n2)

30	format('data',/,(2f12.5))
	dpsi=0.2
		k=0
c		write(6,*)k,z(0,k)
		do k=1,10
		w=cmplx(0.,(k-0.5)*dpsi)
			zeta=2*log(ey*sqrt(w/phis)+sqrt(1.-w/phis))
			a1=(1-exp(zeta/2))/(1+exp(zeta/2))
			a1=log(a1)-pi*ey/2
			a2=exp(zeta/2)
			a2=(a2+1./a2)/2
			omz=omp+oms/(1.+ak*a2)-a1
			do j=0,n-1
			  gm=(g(j)+g(j+1))/2
			  a3=sin((ey*zeta+gm)/4)*sin((ey*zeta-gm)/4)
			  a3=a3/(cos((-ey*zeta+gm)/4)*cos((-ey*zeta-gm)/4))
			  a3=log(a3)
			  omz=omz-a3*dth(j)/pi
			enddo
			omz=omz-pi*ey/2
			z(0,k)=z(0,k-1)+exp(omz)*dpsi
		write(6,*)k,z(0,k)
		phiz=-0.5*dphi
		do i=1,n2
			phiz=phiz+dphi
			w=cmplx(phiz,k*dpsi)
			zeta=2*log(ey*sqrt(w/phis)+sqrt(1.-w/phis))
			a1=(1-exp(zeta/2))/(1+exp(zeta/2))
			a1=log(a1)-pi*ey/2
			a2=exp(zeta/2)
			a2=(a2+1./a2)/2
			omz=omp+oms/(1.+ak*a2)-a1
			do j=0,n-1
			  gm=(g(j)+g(j+1))/2
			  a3=sin((ey*zeta+gm)/4)*sin((ey*zeta-gm)/4)
			  a3=a3/(cos((-ey*zeta+gm)/4)*cos((-ey*zeta-gm)/4))
			  a3=log(a3)
			  omz=omz-a3*dth(j)/pi
			enddo
			omz=omz+pi*ey
			z(i,k)=z(i-1,k)+exp(omz)*dphi
		enddo
		phiz=0.5*dphi
		do i=-1,-n,-1
			phiz=phiz-dphi
			w=cmplx(phiz,k*dpsi)
			zeta=2*log(ey*sqrt(w/phis)+sqrt(1.-w/phis))
			a1=exp(zeta/2)
			a1=(1-a1)/(1+a1)
			a1=log(a1)-pi*ey/2
			a2=exp(zeta/2)
			a2=(a2+1./a2)/2
			omz=omp+oms/(1.+ak*a2)-a1
			do j=0,n-1
			  gm=(g(j)+g(j+1))/2
			  a3=sin((gm+ey*zeta)/4)*sin((ey*zeta-gm)/4)
			  a3=a3/(cos((gm+ey*zeta)/4)*cos((gm-ey*zeta)/4))
			  a3=log(a3)
			  omz=omz-a3*dth(j)/pi
			enddo
			z(i,k)=z(i+1,k)+exp(omz)*dphi
		enddo
c		write(8,30)(z(i,k),i=-n,n2)
c		write(8,30)(real(z(i,k)),-aimag(z(i,k)),i=-n,n2)

	         enddo
	y1=aimag(z(n2,0))
	x1=real(z(n2,0))
	y2=aimag(z(n2-1,0))
	x2=real(z(n2-1,0))
	slope=(y1-y2)/(sqrt(x1)-sqrt(x2))
	write(6,9)slope,x1,y1
9	format('calculated expansion rate ' ,3f15.6)

		do k=0,10,2
		write(8,30)(z(i,k),i=-n,n2)
		write(8,30)(real(z(i,k)),-aimag(z(i,k)),i=-n,n2)
		enddo



	stop
	end
		
			


