program thinair
      parameter(ixx=101)
      real x(ixx),axi(ixx)
      real f(ixx),fp(ixx)
      real e(ixx),ep(ixx)
      real g(ixx),gp(ixx),dg(ixx)
      real zu(ixx),zo(ixx)
      real ga(ixx),gpa(ixx)
      real ue(ixx),uf(ixx),we(ixx),wf(ixx)
      real uu(ixx),uo(ixx)
      real wu(ixx),wo(ixx)
      real uru(ixx),wru(ixx)
      real uro(ixx),wro(ixx)
      real cpu(ixx),cpo(ixx)
      real cpru(ixx),cpro(ixx)
      real ve(ixx),xx(ixx),vgrad(ixx),theta(ixx),cfr(ixx)
      open(unit=14,file='thinair.data',form='formatted')
      open(unit=15,file='thinair.xz',form='formatted')
      open(unit=16,file='thinair.xcp',form='formatted')
      open(unit=17,file='thinair.xu',form='formatted')
      open(unit=18,file='thinair.xw',form='formatted')
      open(unit=19,file='thinair.xg',form='formatted')
      open(unit=20,file='thinair.xgp',form='formatted')
      open(unit=21,file='thinair.xcfu',form='formatted')
      open(unit=22,file='thinair.xcfo',form='formatted')
      open(unit=23,file='thinair.mses',form='formatted')
      open(unit=24,file='thinair.xge',form='formatted')
      open(unit=25,file='thinair.xgpe',form='formatted')
      open(unit=26,file='thinair.xcpu',form='formatted')
c*****constants
      eps=1.e-10
      pi=2.*asin(1.)
      degrad=pi/180.
c*****geometric/flow coefficients
      cxm=1.
      read(14,*)ix
      read(14,*)omega
      read(14,*)dm
      read(14,*)em
      read(14,*)ithick
      Uref=1.
      read(14,*)alphd
      alpha=degrad*alphd
      dtet=pi/(ix-1)
c*****output on listing
      write(6,*)'*************************'
      write(6,*)'dimensionless parameters:'
      write(6,*)'                           X=C*x'
      write(6,*)'                           Z=C*z'
      write(6,*)'                      P-Pinf=0.5*Rho*V**2*Cp'
      write(6,*)'                           L=0.5*Rho*V**2*C*Cl'
      write(6,*)'                           D=0.5*Rho*V**2*C*Cd'
      write(6,*)'                         M,o=0.5*Rho*V**2*C**2*Cm,o'
      write(6,*)
      write(6,*)'******************'
      write(6,*)'main profile data:'
      write(6,*)'              mesh system ix=',ix
      write(6,*)'     relaxation factor omega=',omega
      write(6,*)'               maximum chord=',cxm,' (m)'
      write(6,*)'             relative camber=',dm,' <0 for flying wing'
      write(6,*)'          relative thickness=',em
      write(6,*)'      thickness distribution=',ithick,
     &     '     (ellip/semicubic/q-j/naca00xx/selig)'
      write(6,*)'***********************'
      write(6,*)'aerodynamic parameters:'
      write(6,*)'          reference velocity=',Uref,' (m/s)'
      write(6,*)'             angle of attack=',alpha,' (rd) ='
     &     ,alphd,' (deg)'
      write(6,*)' '
c*****initialization
      write(6,*)
      write(6,*)'do you want to print the mesh, geometry'
     &     ,' and flow? y/n=1/0'
      read(5,*)ipr
      if(ipr.eq.1)write(6,*)' i=  ','     x(i)=   ','     zu(i)=   ',
     &'     zo(i)=   ','     fp(i)=   ','     ep(i)='
      teti=.5*dtet
      axii=.5*cxm*(1.-cos(.5*dtet))
      call camberdis(cxm,dm,axii,teti,fim)
      call camberdis(cxm,dm,0.,0.,fi)
      fim=2.*fi-fim
      call thickdis(ithick,cxm,em,axii,teti,eim)
      axiim=-axii
      eim=-eim
      do 1 i=1,ix
         teti=(i-1)*dtet
         xi=.5*cxm*(1.-cos(teti))
         x(i)=xi
         axii=.5*cxm*(1.-cos(teti+.5*dtet))
         axi(i)=axii
         call camberdis(cxm,dm,axii,teti,f(i))
         fp(i)=(f(i)-fim)/(axii-axiim)
         if(i.ne.ix)then
            call thickdis(ithick,cxm,em,axii,teti+.5*dtet,e(i))
         else
            call thickdis(ithick,cxm,em,axii,teti-.5*dtet,e(i))
         endif
         ep(i)=(e(i)-eim)/(axii-axiim)
         eim=e(i)
         axii=.5*cxm*(1.-cos(teti))
         call thickdis(ithick,cxm,em,axii,teti,e(i))
         call camberdis(cxm,dm,axii,teti,fi)
         zu(i)=fi+e(i)+alpha*(cxm-xi)
         zo(i)=fi-e(i)+alpha*(cxm-xi)
         fim=f(i)
         f(i)=fi
         axiim=axi(i)
         if(i.eq.ix-1)then
            fim=2.*fim
            eim=2.*eim
            axiim=2.*axi(i)-cxm
         endif
         gp(i)=0.
         g(i)=0.
         if(ipr.eq.1)write(6,1000)i,x(i),zu(i),zo(i),fp(i),ep(i)
         write(15,*)x(i),zu(i),zo(i)
 1    continue
      axi(ix)=1.
      iter=0
c*****iteration loop
 100  continue
      write(6,*)'itx=?'
      read(5,*)itx
      if(itx.le.0)goto 300
      do 200 it=1,itx
      iter=iter+1
      idx=0
      dgx=0.
      do 3 i=2,ix-1
         sum=0.
         do 2 j=1,ix-1
            coefj=1./(x(i)-axi(j))
            sum=sum+(g(j+1)-g(j))*coefj
            if(j.eq.i-1)coefi=coefj
            if(j.eq.i)coefi=coefi-coefj
 2       continue
         wf(i)=-sum/(2.*pi)
         coefi=coefi/(2.*pi)
         dg(i)=omega*(wf(i)-Uref*(fp(i)-alpha))/coefi
         g(i)=g(i)+dg(i)
         if(abs(dgx).lt.abs(dg(i)))then
            dgx=dg(i)
            idx=i
         endif
 3    continue
      g(1)=0.
      coefi=(x(ix)-x(ix-1))**2/(x(ix)-x(ix-2))**2
      g(ix)=(coefi*g(ix-2)-g(ix-1))/(coefi-1.)
 200  continue
      write(6,*)'iter=',iter,' dgx=',dgx,' idx=',idx
      goto 100
 300  continue
c*****compute velocity due to camber and thickness
      uf(1)=0.
      uf(ix)=0.
      do 5 i=2,ix-1
         uf(i)=.5*(g(i+1)-g(i-1))/(x(i+1)-x(i-1))
         sum=0.
         do 4 j=1,ix-1
            coefj=1./(x(i)-axi(j))
            sum=sum+(e(j+1)-e(j))*coefj
 4       continue
         ue(i)=Uref*sum/pi
         we(i)=Uref*ep(i)
 5    continue
      we(1)=Uref*ep(1)
      we(ix)=Uref*ep(ix)
c*****exact solution for flat plate at 1.0 rd incidence
      Cle=2.*pi*(1.+0.*2.*dm/cxm)
      sum=0.
      dum=0.
      ermax=0.
      do 6 i=1,ix
         if(i.ge.2)sum=sum+(uf(i-1)+uf(i))*(x(i)-x(i-1))
         teti=(i-1)*dtet
         sq1=sqrt(axi(i)/cxm)
         sq2=sqrt(1.-axi(i)/cxm)
         gpa(i)=2.*Uref*(1.0*sq2/sq1+0.*8.*dm*sq1*sq2)
         if(i.ge.2)dum=dum+.5*(gpa(i-1)+gpa(i))*(x(i)-x(i-1))
         sq1=sqrt(x(i)/cxm)
         sq2=sqrt(1.-x(i)/cxm)
         ga(i)=Uref*cxm*(1.0*(teti+2.*sq1*sq2)
     &        +0.*2.*dm*(teti-2.*sq1*sq2*(1.-2.*x(i)/cxm)))
         gei=2.0*dm*(teti-0.5*sin(2.0*teti))
         write(24,*)x(i),gei
         ermax=ermax+abs(g(i)-gei)
 6    continue
      dum=dum-.5*gpa(ix)*(x(ix)-x(ix-1))
      coefi=(g(ix)-sum)/dum
      ermax=ermax/ix
c*****end of iteration-begin post processing
      ue(1)=ue(2)
     &     +(ue(3)-ue(2))*(x(1)-x(2))/(x(3)-x(2))
      ue(ix)=ue(ix-1)
     &     +(ue(ix-2)-ue(ix-1))*(x(ix)-x(ix-1))/(x(ix-2)-x(ix-1))
      wf(1)=Uref*fp(1)
      wf(ix)=Uref*fp(ix)
      erpmax=0.
      do 7 i=1,ix-1
         teti=(i-1)*dtet
         gp(i)=(g(i+1)-g(i))/(x(i+1)-x(i))
         uf(i)=uf(i)+.5*coefi*gpa(i)
         uu(i)=ue(i)+uf(i)
         uo(i)=ue(i)-uf(i)
         gpei=8.0*dm*sin(teti+0.5*dtet)
         erpmax=erpmax+abs(gp(i)-gpei)
         write(25,*)axi(i),gpei
         cpi=-2.0*ue(i)
         cpi=-8.0*em*(3.0-4.0*x(i))/(3.0*sqrt(3.0))
         cpi=-2.0*(sin(0.5*teti-alpha)*(1.+4.0*em/(3.0*sqrt(3.0))
     &    *(cos(alpha)*sin(0.5*teti)/sin(0.5*teti-alpha+eps)
     &        +2.0*cos(teti))
     &    -2.0*dm*(sin(alpha)*sin(0.5*teti)/sin(0.5*teti-alpha+eps)
     &        -2.0*sin(teti))
     &    )/sin(0.5*teti+eps)-cos(alpha))
         write(26,*)x(i),-cpi,ue(i)
 7    continue
      cpi=-2.0*ue(ix)
      cpi=-8.0*em*(3.0-4.0*x(ix))/(3.0*sqrt(3.0))
      cpi=-2.0*(sin(0.5*pi-alpha)*(1.+4.0*em/(3.0*sqrt(3.0))
     &     *(cos(alpha)*sin(0.5*pi)/sin(0.5*pi-alpha+eps)
     &     +2.0*cos(pi))
     &     -2.0*dm*(sin(alpha)*sin(0.5*pi)/sin(0.5*pi-alpha+eps)
     &     -2.0*sin(pi))
     &     )/sin(0.5*pi+eps)-cos(alpha))
      write(26,*)x(ix),-cpi,ue(ix)
      erpmax=erpmax/(ix-1)
      write(6,*)'ermax=',ermax,' erpmax=',erpmax
      uf(ix)=0.
      uu(ix)=ue(ix)+uf(ix)
      uo(ix)=ue(ix)-uf(ix)
      gp(ix)=0.
c*****regularize velocity field due to thickness
      axiim=0.
      do 8 i=1,ix
         cpu(i)=2.0*(ue(i)+uf(i))/Uref
         cpo(i)=2.0*(ue(i)-uf(i))/Uref
         dx=axi(i)-axiim
         dz=ep(i)*dx
         uru(i)=(Uref+ue(i)+uf(i))*dx*dx/(dx*dx+dz*dz)
         wru(i)=uru(i)*dz/dx
         uru(i)=uru(i)-Uref
         dz=-ep(i)*dx
         uro(i)=(Uref+ue(i)-uf(i))*dx*dx/(dx*dx+dz*dz)
         wro(i)=uro(i)*dz/dx
         uro(i)=uro(i)-Uref
         cpru(i)=2.*uru(i)/Uref
         cpro(i)=2.*uro(i)/Uref
         write(16,*)x(i),cpu(i),cpo(i),cpru(i),cpro(i)
         axiim=axi(i)
         if(i.eq.ix-1)then
            axiim=2.*axiim-cxm
         endif
 8    continue
c*****output on listing
      Cl=2.*g(ix)/(Uref*cxm)
      Cmo=-0.5*pi*(alpha+4.*dm/cxm)
      Cd=0.
      coefi=coefi/degrad
      write(6,*)
      write(6,*)'**************inviscid aerodynamic coefficients:'
      write(6,*)'Cl=',Cl,'   Cmo=',Cmo,'   Cd=',Cd
c      write(6,*)'gp(i)=',gp
c      write(6,*)'g(i)=',g
c*****output for nxyplot
      do 9 i=1,ix
         write(17,*)x(i),uu(i),uo(i),uru(i),uro(i)
 9    continue
c*****output for nxyplot
      do 10 i=1,ix
         wu(i)=Uref*(fp(i)+ep(i)-alpha)
         wo(i)=Uref*(fp(i)-ep(i)-alpha)
         write(18,*)x(i),wu(i),wo(i),wru(i),wro(i)
         write(19,*)x(i),g(i)
         write(20,*)axi(i),gp(i)
 10   continue
c*****boundary layer correction
      write(6,*)
      write(6,*)'          ************************'
      write(6,*)'**********BOUDARY LAYER CORRECTION************'
      write(6,*)'          ************************'
 400  continue
      write(6,*)
      write(6,*)'Reynolds Number Re=?'
      read(5,*)Re
      if(Re.lt.eps.or.Re.gt.1./eps)goto 500
      write(6,*)
      write(6,*)'upper surface'
      do 11 i=1,ix
         ve(i)=1.+uru(i)/Uref
 11   continue
      call moran(ix,itrans,isep,Re,x,zu,ve,xx,vgrad,theta,cfr)
      Cd=0.
      do 12 i=2,ix
         Cd=Cd+.5*(cfr(i-1)+cfr(i))*(x(i)-x(i-1))/cxm
         if(i.gt.isep)then
            cfr(i)=0.
         endif
         write(21,*)x(i),cfr(i)
 12   continue
      write(6,*)'upper surface itrans=',itrans,' isep=',isep,' Cd=',Cd
      Cdt=Cd
      write(6,*)' '
      write(6,*)' '
      write(6,*)'lower surface'
      do 13 i=1,ix
         ve(i)=1.+uro(i)/Uref
 13   continue
      call moran(ix,itrans,isep,Re,x,zo,ve,xx,vgrad,theta,cfr)
      Cd=0.
      do 14 i=2,ix
         Cd=Cd+.5*(cfr(i-1)+cfr(i))*(x(i)-x(i-1))/cxm
         if(i.gt.isep)then
            cfr(i)=0.
         endif
         write(22,*)x(i),cfr(i)
 14   continue
      write(6,*)'lower surface itrans=',itrans,' isep=',isep,' Cd=',Cd
      Cdt=Cdt+Cd
      write(6,*)'total Cd=',Cdt
      goto 400
 500  continue
      do 15 i=1,ix
         ic=ix+1-i
         write(23,*)x(ic),zu(ic)
 15   continue
      do 16 i=2,ix
         write(23,*)x(i),zo(i)
 16   continue
c*****files
      write(6,*)'******input files:'
      write(6,*)'thinair.data                   :geom. & aero. para.'
      write(6,*)'******output files:'
      write(6,*)'thinair.xz                     :x, zu(x) & zo(x)'
      write(6,*)'thinair.xcp                    :x, Cpu(x) & Cpo(x), Cpr
     &u(x) & Cpro(x)'
      write(6,*)'thinair.xu                     :x, uu(x) & uo(x), uru(x
     &) & uro(x)'
      write(6,*)'thinair.xw                     :x, wu(x) & wo(x), wru(x
     &) & wro(x)'
      write(6,*)'thinair.xg                     :x, g(x)'
      write(6,*)'thinair.xgp                    :x, gp(x)'
      write(6,*)'thinair.xcfu                   :x, Cfu(x)'
      write(6,*)'thinair.xcfo                   :x, Cfo(x)'
      write(6,*)'thinair.mses                   :x, zu(x) & zo(x)'
      write(6,*)'thinair.xge                    :x, ge(x)'
      write(6,*)'thinair.xgpe                   :x, gpe(x)'
      write(6,*)'thinair.xcpu                   :x,-cp,ue(x) analytic'
 1000 format(1x,i3,4x,e10.4,4x,e10.4,4x,e10.4,4x,e10.4,4x,e10.4)
      end

      subroutine camberdis(cxm,dm,axii,teti,fim)
      fim=4.*dm*axii*(1.-axii/cxm)
      if(dm.gt.0.)return
      fim=-fim*(7.0-8.0*axii)
      return
      end

      subroutine thickdis(ithick,cxm,em,axii,teti,eim)
      data fasc/.4592793/faqj/.3849002/fasl/.3088162/fass/.2414953/
c
c     1=elliptic distribution
c     2=semi-cubic distribution
c     3=quasi-joukowski distribution
c     4=naca00em distribution
c     5=selig distribution
c     6=super-selig distribution
c
      goto(1,2,3,4,5,6)ithick
 1    eim=.5*em*cxm*sin(teti)
      goto 7
 2    eim=fasc*em*cxm*sqrt(1.+cos(teti))*sin(teti)
      goto 7
 3    eim=faqj*em*cxm*(1.+cos(teti))*sin(teti)
      goto 7
 4    eim=5.*em*cxm*(.2969*sqrt(axii/cxm)-.126*axii/cxm-
     &.3537*(axii/cxm)**2+.2843*(axii/cxm)**3-.1015*(axii/cxm)**4)
      goto 7
 5    eim=fasl*em*cxm*sin(teti)*(1.+cos(teti))**1.5
      goto 7
 6    eim=fass*em*cxm*sin(teti)*(1.+cos(teti))**2
 7    continue
      return
      end


      subroutine moran(ix,itrans,isep,Re,x,y,ve,xx,vgrad,theta,cfr)
      real x(ix),y(ix),ve(ix)
      real xx(ix),theta(ix),vgrad(ix),cfr(ix)
      real yy(2)
      itrans=ix
      isep=ix
      pi=3.1415926535
      eps=1.e-12
      xx(1)=0.
      do 100 i=2,ix
         dx=x(i)-x(i-1)
         dy=y(i)-y(i-1)
         xx(i)=xx(i-1)+sqrt(dx*dx+dy*dy)
 100  continue
      v1=ve(3)
      x1=xx(3)
      v2=ve(1)
      x2=xx(1)
      do 110 i=1,ix
         v3=v1
         x3=x1
         v1=v2
         x1=x2
         if(i.lt.ix)then
            v2=ve(i+1)
            x2=xx(i+1)
         else
            v2=ve(ix-1)
            x2=xx(ix-1)
            v3=ve(ix-2)
            x3=xx(ix-2)
         endif
         fact=(x3-x1)/(x2-x1)
         vgrad(i)=((v2-v1)*fact-(v3-v1)/fact)/(x3-x2)
 110  continue
      write(6,1000)
      i=1
 120  continue
      v0=ve(i)
      if(ve(i).lt.0.)then
         theta(i)=0.
         h=0.
         cf=0.
         cfr(i)=cf
         write(6,1010)x(i),y(i),ve(i),vgrad(i),theta(i),h,cf
         i=i+1
         goto 120
      endif
      dv=eps*vgrad(i)
      if(dv.lt.eps**2)dv=eps**2
      theta(i)=sqrt(eps*(6.*v0**5+15.*v0**4*dv+20.*v0**3*dv**2
     &+15.*v0**2*dv**3+6.*v0*dv**4+dv**5)/(6.*(v0+dv)**6))
 200  continue
      alam=theta(i)**2*vgrad(i)*Re
      if(alam.gt..1)alam=.1
      if(alam.lt.-.0899)alam=-.0899
c      if(alam.lt.-.0842)goto 400
      call thwats(alam,h,el)
      if(ve(i).lt.eps)then
         cf=2.*el/Re/theta(i)/eps
      else
         cf=2.*el/Re/theta(i)/ve(i)
      endif
      cfr(i)=cf
      write(6,1010)x(i),y(i),ve(i),vgrad(i),theta(i),h,cf
      i=i+1
      if(i.gt.ix)return
      if(ve(i).lt.0.)goto 120
      dth2ve6=.225*(ve(i)**5+ve(i-1)**5)*(xx(i)-xx(i-1))/Re
      theta(i)=sqrt(((theta(i-1)**2)*(ve(i-1)**6)+dth2ve6)
     &     /(ve(i)**6))
      Rex=Re*xx(i)*ve(i)
      Ret=Re*theta(i)*ve(i)
      Retmax=1.174*(1.+22400./Rex)*Rex**.46
      if(Ret.lt.Retmax.or.i.le.1)goto 200
      itrans=i
      h=1.47
      i=itrans
      yy(2)=h1ofh(h)
      yy(1)=theta(i-1)
 320  continue
      dx=xx(i)-xx(i-1)
      call runge2(ix,i,Re,dx,yy,ve,vgrad)
      theta(i)=yy(1)
      h=hofh1(yy(2))
      rtheta=Re*ve(i)*theta(i)
      cf=cfturb(rtheta,h)
      cfr(i)=cf
      write(6,1020)x(i),y(i),ve(i),vgrad(i),theta(i),h,cf
      if(h.gt.2.4)goto 410
      i=i+1
      if(i.le.ix)goto 320
      return
 400  continue
      isep=i
      write(6,*)'laminar separation'
      return
 410  continue
      isep=i
      write(6,*)'turbulent separation'
      return
 1000 format(//,9x,'x',8x,'y',8x,'ve',7x,'vdot',6x,'theta',9x,'h',
     &     9x,'cf',/)
 1010 format('l',f10.5,f9.5,1x,e10.4,1x,e10.4,1x,e10.4,1x,e10.4,1x,
     &     e10.4)
 1020 format('t',f10.5,f9.5,1x,e10.4,1x,e10.4,1x,e10.4,1x,e10.4,1x,
     &     e10.4)
      end

      subroutine thwats(alam,h,el)
      if(alam.lt.0.)then
         el=.22+1.402*alam+.018*alam/(.107+alam)
         h=2.088+.0731/(.14+alam)
      else
         el=.22+alam*(1.57-1.8*alam)
         h=2.61-alam*(3.75-5.24*alam)
      endif
      return
      end

      function h1ofh(h)
      if(h.gt.1.6)then
         h1ofh=3.3+1.5501*(h-.6778)**(-3.064)
      else
         h1ofh=3.3+.8234*(h-1.1)**(-1.287)
      endif
      return
      end

      function hofh1(h1)
      if(h1.lt.5.3)then
         if(h1.lt.3.3)then
            hofh1=3.0
         else
            hofh1=.6778+1.1536*(h1-3.3)**(-.326)
         endif
      else
         hofh1=1.1+.86*(h1-3.3)**(-.777)
      endif
      return
      end

      function cfturb(rtheta,h)
      cfturb=.246*(10.**(-.678*h))*(rtheta)**(-.268)
      return
      end

      subroutine derivs(ix,i,Re,yp,yt,ve,vgrad)
      real ve(ix),vgrad(ix)
      real yp(2),yt(2)
      h1=h1ofh(hofh1(yt(2)))
      h=hofh1(h1)
      rtheta=Re*ve(i)*yt(1)
      yp(1)=-(h+2.)*yt(1)*vgrad(i)/ve(i)
     &     +.5*cfturb(rtheta,h)
      yp(2)=-h1*(vgrad(i)/ve(i)+yp(1)/yt(1))
     &     +.0306*(h1-3.)**(-.6169)/yt(1)
      return
      end

      subroutine runge2(ix,i,Re,dx,yy,ve,vgrad)
      real ve(ix),vgrad(ix)
      real ys(2),yy(2)
      real yp(2),yt(2)
      do 1 j=1,2
         yt(j)=yy(j)
 1    continue
      call derivs(ix,i-1,Re,yp,yt,ve,vgrad)
      do 2 j=1,2
         yt(j)=yy(j)+dx*yp(j)
         ys(j)=yy(j)+.5*dx*yp(j)
 2    continue
      call derivs(ix,i,Re,yp,yt,ve,vgrad)
      do 3 j=1,2
         yy(j)=ys(j)+.5*dx*yp(j)
 3    continue
      return
      end