	program brodgen
c
c	calculate brodetskys solution for arbitrary constants
c
c	z is the body & free streamline flow path
c	circ is the actual body
c	om is omega = log(zeta)
c	w is the complex potential
c	
c	capA are coefficents in Brodetsky's expansion, calculated using
c	matlab routine solvebrod.m
c
c	q is speed along z streamline
c	th is angle relative to -x axis of point z
c	nb = number of points on body
c	nfs = number of points on free streamline

	complex *8 z(0:5000),ey,arg1,arg2,arg3,circ(0:200)
	complex *8 om(0:5000),zeta(0:1000),w(0:5000)
	real *8 capA(100),a(100),q(0:5000),th(0:5000)
	pi=3.1415927
	nb=500
	nfs=500
	dt=pi/(2*nb)
	z(0)=(0.,0.)
	ey=(0.,1.)
	read(5,*)k,(capA(i),i=1,k)
	sum=1.
	kk=1
	do i=1,k
		sum=sum+kk*capA(i)
		kk=-kk
	enddo
c------------------------------------------------------------------------
c	generate circle
c------------------------------------------------------------------------
	c=0
	do i=1,k
		c=c+(-1)**i*capA(i)*(2*i-1)
	enddo
	r=2/c
	do i=0,200
		circ(i)=r*exp(ey*pi*i/200.)
	end do
	write(7,10)(circ(i)/r,i=0,200)
	z(0)=cmplx(-r,0.)
	pgamma=-sum*sqrt(r)/2
	write(6,1)pgamma,r
1	format('Brodetsky pressure singularity= ',f15.6,' radius ',f15.6)
	write(6,21)(2+capA(1))/sqrt(r)
21	format('Wake growth ',f15.6)
	
c------------------------------------------------------------------------
c	calculate body before separation
c------------------------------------------------------------------------
	do i=1,nb
		t=(i-0.5)*dt
		
		arg1=exp(ey*t)
		arg2=(0.,0.)
		do j=1,k
			arg2=arg2+capA(j)*exp((2*j-1)*ey*t)/(2*j-1)
		enddo
		om(i)=log((1+arg1)/(1-arg1))+arg2
		zeta(i)=exp(om(i))
		arg2=exp(arg2)
		arg3=ey*(1+arg1)*(1+arg1)*(1+arg1*arg1)/(2*arg1**2)
		z(i)=z(i-1)+arg3*arg2*dt
		q(i)=exp(-real(om(i)))
		th(i)=pi-atan2(aimag(z(i)),real(z(i)))
	end do
	dq=(q(nb)-q(nb-1))/dt
	ddth=(th(nb)-2*th(nb-1)+th(nb-2))/(dt*dt)
c	write(6,2)sum,dq
2	format('gamma =',f15.6,' dq= ',f15.6)
c	write(6,3)2/r,ddth
3	format('2/r   =',f15.6,' ddth=',f15.6)
c------------------------------------------------------------------------
c	calculate free streamline after separation
c------------------------------------------------------------------------
	dt=0.9/nfs
	do i=nb+1,nb+nfs
		t=1.-dt*(i-nb-0.5)
		arg2=(0.,0.)
		do j=1,k
			n=2*j-1
			arg2=arg2+ey*capA(j)*(-1)**(j-1)*t**n/n
		enddo
		om(i)=log((1+t*ey)/(1-t*ey))+arg2
		zeta(i)=exp(om(i))
		arg2=exp(arg2)
		arg3=(1+ey*t)**2*(1-t*t)/t**3
		z(i)=z(i-1)+arg3*arg2*dt/2
		q(i)=exp(-real(om(i)))
		th(i)=pi-atan2(aimag(z(i)),real(z(i)))
	enddo
	y1=aimag(z(nb+nfs))
	x1=real(z(nb+nfs))
	y2=aimag(z(nb+nfs-1))
	x2=real(z(nb+nfs-1))
	slope=(y1-y2)/(sqrt(x1)-sqrt(x2))
	write(6,9)slope,x1,y1
9	format('calculated expansion rate ' ,3f15.6)


	write(7,10)(z(i)/r,i=0,nb)
10	format('data',/,(2f15.6))
	write(7,10)(z(i)/r,i=nb+1,nb+nfs)

	write(8,11)(i,abs(z(i)-r),i=0,nb)
11	format('data',/,(i5,f15.8))

	write(9,10)(om(i),i=1,nb)
	write(9,10)(om(i),i=nb+1,nb+nfs)

	write(10,10)(zeta(i),i=1,nb)
	write(10,10)(zeta(i),i=nb+1,nb+nfs)

	write(11,10)(180*th(i)/pi,1.-q(i)**2,i=1,nb)
	write(11,10)(180*th(i)/pi,1.-q(i)**2,i=nb+1,nb+nfs)



	stop
	end

