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