program lcw4
c
c	program to calculate flow separated from upper wall in channel
c	modified for separation-reattachment
c	modified for variable theta at reattachment
c
c	p(i)=psi(x(i))
c	dp(i)=d psi/dx(x(i))
c	speed(i)=d psi/ds(x(i))
c	th(i)=theta(x(i))
c	b(i)=dx/ds(x(i))

	real*8 p(-1000:10000),dp(-1000:10000)
	real*8 pr(-1000:10000),dpr(-1000:10000)
	real *8 dx, f(-1000:10000),df(-1000:10000),xs,s(-1000:10000)
	real*8 th(-1000:10000),speed(-1000:10000)
	real*8 b(-1000:10000)
	real *8 pi,ominf
	complex *8 ey,om(-1000:10000),arg,a(-1000:10000),z(-1000:10000)
	complex *8 zs(-1000:10000,0:11),w,arg1
	pi=3.14159
	ey=(0.,1.)
	
c---------------------------------------------------------------------------
c	separation and reattachment can be specified arbitrarily
c	but only correct values will give correct downstream conditions.
c----------------------------------------------------------------------------
	write(6,1)
1	format(' Woods Method: Separation and reattachment',/,
     +         'Input separation point (0<x<1)')
	read(5,*)xs
	write(6,2)
2	format('Input xr')
	read(5,*)xr
c----------------------------------------------------------------------------
c	estimate phi0 from q=1 over xs->xr
c----------------------------------------------------------------------------
	np=100
	dx=xs/np
	phi0=xr-xs
c	nrs = start of reattachement
c	nre = end of no zero theta
	nre=4./dx-np
	nrs=xr/dx-np
	ys=-0.5*(cos(pi*xs)-1)
	na=-400
	nb=nre+50
	write(6,4)na,np,nrs,nre,nb
4	format(' na = ',i5,' np = ',i5,' nrs = ',i5,' nre = ',i5,' nb = ',i5)
	write(6,3)na*dx,nb*dx
3	format('Upstream limit   x= ',f8.2,/,'Downstream limit x= ',f8.2)
c----------------------------------------------------------------------------
c	set up wall and intial guess for potential p(i)
c----------------------------------------------------------------------------
	do i=na,nb
		x=i*dx
		p(i)=x
		dp(i)=1.
		if (i.lt.-np)then
			f(i)=1
			df(i)=0.
			th(i)=0.
		else if(x.le.4-xs) then
			f(i)=1-0.5*(cos(0.5*pi*(x+xs))-1)
			df(i)=0.5*0.5*pi*sin(0.5*pi*(x+xs))
			th(i)=atan(df(i))
		else
			f(i)=1
			df(i)=0
			th(i)=0
		endif
		z(i)=cmplx(x,f(i))
	enddo
c	guess theta =0 on free streamline to start
	do i=1,nrs-1
		th(i)=0.
	enddo
c
c	output boundary shape
c
	write(10,30)(z(i),i=na,nb)
c---------------------------------------------------------------------------------
c	iterate to calculate phi(x) for x<0, particularly on wall with theta<>0
c---------------------------------------------------------------------------------
	do loop=1,10
c
c		calculate ominf from integral ---------------------ominf
c
		ominf=0.
		do j=-np,-1
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			os=sqrt((1-exp(pi*ps))/(1-exp(pi*ps0)))
			os=0.5*log((1+os)/(1-os))
			ominf=ominf+2*os*dth/pi
		enddo
		do j=nrs,nre
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1-exp(pi*(ps-phi0)))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			os=real(arg)
			ominf=ominf+2*os*dth/pi
		enddo
		vel=exp(-ominf)
c
c		calculate integral on wall and free streamline
c
		do i=na,-1
			b(i)=1./sqrt(1+df(i)*df(i))
			s(i)=0.
			do j=-np,-1
				dth=th(j+1)-th(j)
				ps=0.5*(p(j)+p(j+1))
				ps0=ps  -phi0
				pi0=p(i)-phi0
				arg=(1-exp(pi*ps))/(1-exp(pi*p(i)))
				arg=arg*(1-exp(pi*pi0))/(1-exp(pi*ps0))
				arg=sqrt(arg)
				arg=0.5*log((1+arg)/(1-arg))
				arg=2*arg*dth/pi
				s(i)=s(i)+real(arg)
			enddo
			do j=nrs,nre
				dth=th(j+1)-th(j)
				ps=0.5*(p(j)+p(j+1))
				ps0=ps  -phi0
				pi0=p(i)-phi0
				arg=(1-exp(pi*ps))/(1-exp(pi*p(i)))
				arg=arg*(1-exp(pi*pi0))/(1-exp(pi*ps0))
				arg=sqrt(arg)
				arg=0.5*log((1+arg)/(1-arg))
				arg=2*arg*dth/pi
				s(i)=s(i)+real(arg)
			enddo
			s(i)=ominf-s(i)
			speed(i)=exp(-s(i))
			dp(i)=speed(i)/b(i)
			
		enddo
		b(0)=1./sqrt(1+df(0)*df(0))
		speed(0)=exp(-ominf)
		dp(0)=speed(0)/b(0)
		p(0)=0.
		do i=1,nrs-1
			w=cmplx(p(i),1.)
			om(i)=cmplx(ominf,th(nrs))
			do j=-np,-1
			   dth=th(j+1)-th(j)
			   ps=0.5*(p(j)+p(j+1))
			   ps0=ps-phi0
			   arg=(1-exp(pi*ps))/(1+exp(pi*w))
			   arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			   arg=sqrt(arg)
			   arg=0.5*log((1+arg)/(1-arg))
			   arg=2*arg*dth/pi
			   om(i)=om(i)-arg
			enddo
			do j=nrs,nre
			   dth=th(j+1)-th(j)
			   ps=0.5*(p(j)+p(j+1))
			   ps0=ps-phi0
			   arg=(1-exp(pi*ps))/(1+exp(pi*w))
			   arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			   arg=sqrt(arg)
			   arg=0.5*log((1+arg)/(1-arg))
			   arg=2*arg*dth/pi
			   om(i)=om(i)-arg
			enddo
			th(i)=imag(om(i))
			speed(i)=exp(-real(om(i)))
			dp(i)=speed(i)/cos(th(i))	
		enddo
		do i=nrs,nb-1
			b(i)=1./sqrt(1+df(i)*df(i))
			s(i)=0.
			do j=-np,-1
				dth=th(j+1)-th(j)
				ps=0.5*(p(j)+p(j+1))
				ps0=ps  -phi0
				pi0=p(i)-phi0
				arg=(1-exp(pi*ps))/(1-exp(pi*p(i)))
				arg=arg*(1-exp(pi*pi0))/(1-exp(pi*ps0))
				arg=sqrt(arg)
				arg=0.5*log((1+arg)/(1-arg))
				arg=2*arg*dth/pi
				s(i)=s(i)+real(arg)
			enddo
			do j=nrs,nre
				dth=th(j+1)-th(j)
				ps=0.5*(p(j)+p(j+1))
				ps0=ps  -phi0
				pi0=p(i)-phi0
				arg=(1-exp(pi*ps))/(1-exp(pi*p(i)))
				arg=arg*(1-exp(pi*pi0))/(1-exp(pi*ps0))
				arg=sqrt(arg)
				arg=0.5*log((1+arg)/(1-arg))
				arg=2*arg*dth/pi
				s(i)=s(i)+real(arg)
			enddo
			s(i)=ominf-s(i)
			speed(i)=exp(-s(i))
			dp(i)=speed(i)/b(i)
			
		enddo
c------------------------------------------------------------------------------
c	calculate potential p(i) from dp(i)
c------------------------------------------------------------------------------
		p(0)=0.
		do i=-1,na,-1
			p(i)=p(i+1)-0.5*(dp(i+1)+dp(i))*dx
		enddo
		do i=1,nb
			p(i)=p(i-1)+0.5*(dp(i)+dp(i-1))*dx
		enddo
c	update guess for phi0
		phi0=p(nrs)
		write(6,10)exp(-ominf),phi0
10		format('V = ',2f12.5)
	enddo
c-----------------------------------------------------------------------------
c	end of iteration to calculate potential on wall and free streamline
c-----------------------------------------------------------------------------
	write(12,20)(i*dx,th(i),i=na,nb)
	write(7,20)(i*dx,p(i),i=na,nb)
	write(13,20)(i*dx,dp(i),i=na,nb)
20	format(2f15.6)
	write(8,20)(i*dx,speed(i),i=na,nb)
c------------------------------------------------------------------------------
c	base values p(i), ominf now calculated
c------------------------------------------------------------------------------
c	calculate singularity strength
		sing=0.
		do j=-np,-1
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			os=(exp(pi*phi0)-exp(pi*ps))
			os=os/(1-exp(pi*ps))
			os=os/(exp(pi*phi0)-1)
			os=sqrt(os)
			sing=sing+2*os*dth/sqrt(pi)
		enddo
		do j=nrs,nre
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			os=(exp(pi*phi0)-exp(pi*ps))
			os=os/(1-exp(pi*ps))
			os=os/(exp(pi*phi0)-1)
			os=sqrt(os)
			sing=sing+2*os*dth/sqrt(pi)
		enddo
		singp=vel*vel*sqrt(vel)*sing/2
		write(6,24)singp
24	format('Pressure singularity strength:',f15.6)

c-----------------------------------------------------------------------------	
c	now calculate omega (complex)for positive values of phi
c	up to now, index i has refered to i*dx, now it will refer
c	to i*dphi, index j still refers to p(j) at points x(j)
c-----------------------------------------------------------------------------
	dphi=0.01
	phi=0
	do i=1,nb
		phi=(i-0.5)*dphi
		w=cmplx(phi,1.)
		om(i)=cmplx(ominf,th(nrs))
		do j=-np,-1
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
		enddo
		do j=nrs,nre
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
			speed(i)=exp(-real(om(i)))
		enddo
	write(15,15)i,om(i),exp(om(i))
15	format(i5,4f12.3)
		z(i)=z(i-1)+exp(om(i))*dphi
	enddo
	write(9,20)(z(i),i=na,nb)
	write(11,20)(real(z(i)),speed(i),i=na,nb)
c
c	calculate streamlines
c
	do i=na,nb
		zs(i,10)=z(i)
	enddo
	dpsi=0.1
	do k=9,1,-1
	  psi=dpsi*(k+0.5)
	  w=cmplx(0.,psi)
	  arg1=cmplx(ominf,th(nrs))
	  do j=-np,-1
		dth=th(j+1)-th(j)
		ps=0.5*(p(j)+p(j+1))
		ps0=ps-phi0
		arg=(1-exp(pi*ps))/(1+exp(pi*w))
		arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
		arg=sqrt(arg)
		arg=0.5*log((1+arg)/(1-arg))
		arg=2*arg*dth/pi
		arg1=arg1-arg
	  enddo
	  do j=nrs,nre
		dth=th(j+1)-th(j)
		ps=0.5*(p(j)+p(j+1))
		ps0=ps-phi0
		arg=(1-exp(pi*ps))/(1+exp(pi*w))
		arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
		arg=sqrt(arg)
		arg=0.5*log((1+arg)/(1-arg))
		arg=2*arg*dth/pi
		arg1=arg1-arg
	  enddo
		zs(0,k)=zs(0,k+1)-exp(arg1)*dpsi*ey
	write(6,25)k,zs(0,k)
25	format(i5,2f15.6)
	  do i=1,nb
		phi=(i-0.5)*dphi
		w=cmplx(phi,psi)
		om(i)=cmplx(ominf,th(nrs))
		do j=-np,-1
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
		enddo
		do j=nrs,nre
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
		enddo
		zs(i,k)=zs(i-1,k)+exp(om(i))*dphi
	  enddo
	  do i=-1,na,-1
		phi=(i+0.5)*dphi
		w=cmplx(phi,psi)
		om(i)=cmplx(ominf,th(nrs))
		do j=-np,-1
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
		enddo
		do j=nrs,nre
			dth=th(j+1)-th(j)
			ps=0.5*(p(j)+p(j+1))
			ps0=ps-phi0
			arg=(1-exp(pi*ps))/(1+exp(pi*w))
			arg=arg*(1+exp(pi*(w-phi0)))/(1-exp(pi*ps0))
			arg=sqrt(arg)
			arg=0.5*log((1+arg)/(1-arg))
			arg=2*arg*dth/pi
			om(i)=om(i)-arg
		enddo
		zs(i,k)=zs(i+1,k)-exp(om(i))*dphi
	  enddo
	enddo
	do k=1,10
	write(10,30)(zs(i,k),i=na,nb)
30	format('data',/,(2f15.7))
	enddo
	q=-log(ominf)
	write(6,40)z(nb),q
40	format('Limiting streamline:',2f15.6,/,'Speed',f15.6)
	stop
	end