
      PROGRAM MACH1

!***********************************************************************

!    Escoamentos MACH e RANKINE
 
!    Solução analítica do escoamento quase-unidimensional
!    permanente de fluido compressível invíscido, 
!    sem troca de calor, com ou sem onda de choque normal.
 
!    CARLOS HENRIQUE MARCHI - DEMEC/UFPR
!    CURITIBA, MAIO/96                                
!                                           última alteração: 18 Dez 97
!***********************************************************************
 
!*** NOMENCLATURA
 
!   A = AREA
!   AI = AREA NA ENTRADA
!   AE = AREA NA SAIDA
!   AT = AREA NA GARGANTA
!   EMPUXO = EMPUXO
!   FMI = FLUXO DE MASSA NA ENTRADA
!   FME = FLUXO DE MASSA NA SAIDA
!   FMT = FLUXO DE MASSA NA GARGANTA
!   GAMA = RAZAO DOS CALORES ESPECIFICOS
!   IC = NUMERO DE PONTOS + 1 ATE O CHOQUE
!   ISHOCK = CONSTANTE QUE DEFINE SE HA (1) OU NAO (0) CHOQUE NORMAL
!   L = COMPRIMENTO TOTAL DO DOMINIO DE CALCULO
!   M = NUMERO DE MACH
!   MACH1DEX.D = DISPOSITIVO 5 DE I/O, ARQUIVO DE DADOS
!   MACH1DEX.S = DISPOSITIVO 6 DE I/O, ARQUIVO DE RESULTADOS
!   MC1 = NUMERO DE MACH ANTES DA ONDA DE CHOQUE
!   MC2 = NUMERO DE MACH APOS A ONDA DE CHOQUE
!   NG = NUMERO DE PONTOS + 1 ATE A GARGANTA
!   NI = NUMERO DE PONTOS TOTAL + 1
!   NITM = NUMERO MAXIMO DE ITERACOES NO METODO DE NEWTON-RAPHSON
!   PO1 = PRESSAO TOTAL NA ENTRADA
!   PO2 = PRESSAO TOTAL NA SAIDA
!   P = PRESSAO ESTATICA
!   PC1 = PRESSAO ESTATICA ANTES DA ONDA DE CHOQUE
!   PC2 = PRESSAO ESTATICA APOS A ONDA DE CHOQUE
!   PI = CONSTANTE MATEMATICA
!   PO1 = PRESSAO TOTAL NA ENTRADA
!   PO2 = PRESSAO TOTAL NA SAIDA
!   R = CONSTANTE DO GAS
!   RA = RAZAO DE AREAS ENTRE A AREA LOCAL E A AREA DA GARGANTA
!   RO = MASSA ESPECIFICA ESTATICA
!   ROC1 = MASSA ESPECIFICA ESTATICA ANTES DA ONDA DE CHOQUE
!   ROC2 = MASSA ESPECIFICA ESTATICA APOS A ONDA DE CHOQUE
!   ROI = MASSA ESPECIFICA ESTATICA NA ENTRADA
!   ROE = MASSA ESPECIFICA ESTATICA NA SAIDA
!   ROT = MASSA ESPECIFICA ESTATICA NA GARGANTA
!   RO1 = MASSA ESPECIFICA TOTAL NA ENTRADA
!   RO2 = MASSA ESPECIFICA TOTAL NA SAIDA
!   T = TEMPERATURA ESTATICA
!   TC1 = TEMPERATURA ESTATICA ANTES DA ONDA DE CHOQUE
!   TC2 = TEMPERATURA ESTATICA APOS A ONDA DE CHOQUE
!   TO1 = TEMPERATURA TOTAL NA ENTRADA
!   TO2 = TEMPERATURA TOTAL NA SAIDA
!   TOL = TOLERANCIA NO METODO DE NEWTON-RAPHSON
!   U = VELOCIDADE LOCAL
!   UI = VELOCIDADE NA ENTRADA
!   UE = VELOCIDADE NA SAIDA
!   UT = VELOCIDADE NA GARGANTA
!   X = COORDENADA ESPACIAL NA DIRECAO DO ESCOAMENTO
!   Y = RAIO
 
!   OBS: EMPREGA-SE O SISTEMA INTERNACIONAL DE UNIDADES (SI)
    
!***********************************************************************
      use msflib
!***********************************************************************
 
      PARAMETER (   NV = 1004  )

!***********************************************************************
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION X(NV),T(NV),P(NV),Y(NV),RO(NV),U(NV),cpa(nv),cva(nv), &
                ve(nv),aa(nv)
      REAL*8 M(NV),L,MC1,MC2,mm,LC,LT
      logical chamada
      integer*4 int4
 
      chamada = systemqq('notepad mach1.d')
 
      OPEN(5,FILE='mach1.d',FORM='FORMATTED')
      OPEN(6,FILE='mach1.s',FORM='FORMATTED')
 
      ru = 8.314510d+3  ! constante universal dos gases (J/kmol.K)
      g  = 9.80665d+0	! aceleração gravitacional ao nível do mar (m/s)

      CALL DADOS(RU,TO1,PO1,GAMA,NI,R,RO1,MM,PA,ISHOCK,IC,IS,LC,LT,DT,RCOS)

      NIP1 = NI + 1
      PI = DACOS(-1.0D0)
      CALL GRID(X,Y,NIP1,PI,NG,L,LC,LT,DT,RCOS)
 
      IF(ISHOCK.EQ.1) THEN
        NI = NI + 1
        NIP1 = NI + 1
        DO I = NIP1,IC+2,-1
          X(I) = X(I-1) 
          Y(I) = Y(I-1)
    	END DO
        X(IC+1) = X(IC)
        Y(IC+1) = Y(IC)
      END IF
 
      AG = PI * (Y(NG) ** 2)
      AI = PI * (Y(2) ** 2)
      AE = PI * (Y(NI) ** 2)
 
      write(6,64)
   64 format(///,5x,'Solução com Newton-Raphson',/)

 
      DO 10 I = 2, NI
 
        IF(I.EQ.IC+1.AND.ISHOCK.EQ.1) GO TO 5
 
        A = PI * (Y(I) ** 2)
        RA = A / AG
        IF(I.GT.IC.AND.ISHOCK.EQ.1) RA = A / AG2
 
        IF(ISHOCK.EQ.1) THEN
          IF(I.LT.NG.OR. I.GT.IC) IND = 0
          IF(I.GT.NG.AND.I.LE.IC) IND = 1
        ELSE
          IF(I.LT.NG) IND = 0
          IF(I.GT.NG.AND.IS==0) IND = 1
          IF(I.GT.NG.AND.IS==1) IND = 0
        END IF
        IF(I.NE.NG) CALL NEWTON(GAMA,RA,IND,M(I),I)
        IF(I.EQ.NG) M(NG) = 1.0D0
 
        IF(ISHOCK.EQ.1) THEN
          IF(I.LE.IC) CALL PROPRI(M(I),TO1,PO1,T(I),P(I),GAMA)
          IF(I.GT.IC) CALL PROPRI(M(I),TO2,PO2,T(I),P(I),GAMA)
          IF(I.EQ.IC) CALL CHOQUE(M(I),T(I),P(I),R,GAMA,MC1,TC1,PC1,ROC1, &
                                  MC2,TC2,PC2,ROC2,AG2,A)
        ELSE 
          CALL PROPRI(M(I),TO1,PO1,T(I),P(I),GAMA)
        END IF
 
        GO TO 10
    5   M(I) = MC2
        T(I) = TC2
        P(I) = PC2
        CALL ESTAG(MC2,TC2,PC2,GAMA,TO2,PO2)
        RO2 = PO2 / (R * TO2)
 
   10 CONTINUE

 
      ROI = P(2) / (R * T(2))
      ROT = P(NG) / (R * T(NG))
      ROE = P(NI) / (R * T(NI))
 
      UI = M(2) * DSQRT(GAMA * R * T(2))
      UT = M(NG) * DSQRT(GAMA * R * T(NG))
      UE = M(NI) * DSQRT(GAMA * R * T(NI))
      ao = DSQRT(GAMA * R * TO1)
 
      FMI = ROI * AI * UI
      FMT = ROT * AG * UT
      FME = ROE * AE * UE
      fmo = ro1 * ag * ao

      fo = po1*ag

      cest = fo / fme
 
      yt = y(ng)

      DO I = 2, NI
        RO(I) = P(I)/(R*T(I)*RO1)
        aa(i) = PI * (Y(I)**2)
        U(I) = FMI / (RO(I)*RO1*aa(i)*ao)
        ve(I) = FMI / (RO(I)*RO1*aa(i))
        if(m(i) >= 1.0d0) then
          pe = p(i)
        else
          pe = pa
        end if
        cpa(i) = ve(i) + (pe-pa)*aa(i)/fmi
        cva(i) = ve(i) + p(i)*aa(i)/fmi
        aa(i) = aa(i)/ag
        X(I) = X(I) / L
        Y(I) = Y(I) / yt
        T(I) = T(I) / TO1
        P(I) = P(I) / PO1
      END DO


      WRITE(6,90) FMI,FMT,FME,FMI/FMo,ag,4.0d0*rcos+dt,   &
                  cest,fo,ui,ao,cva(ni)/ve(ni)
   90 FORMAT(///,5X,'FLUXOS DE MASSA (kg/s)',//, &
                 1pe20.9,' = entrada',/,	  &
                 1pe20.9,' = garganta',/,	  &
                 1pe20.9,' = saída',//,  &
                 1pe20.9,' = fluxo de massa (adimensional)',//,  &
                 1pe20.9,' = área da garganta (m2)',/,  &
                 1pe20.9,' = diâmetro da câmara de combustão (m)',//,  &
                 1pe20.9,' = velocidade característica (m/s)',//,  &
                 1pe20.9,' = empuxo de referência (N)',//, &
                 1pe20.9,' = velocidade na entrada (m/s)',/,  &
                 1pe20.9,' = velocidade de referência (m/s)',//,  &
                 1pe20.9,' = Is no vácuo / Is ao nível do mar (adim.)')

      write(6,91) x(2),x(ng),x(ni),     &
                  aa(2),aa(ng),aa(ni),  &
                  y(2),y(ng),y(ni),  &
                  m(2),m(ng),m(ni),  &
                  u(2),u(ng),u(ni),  &
                  p(2),p(ng),p(ni),  &
                  t(2),t(ng),t(ni),  &
                  ro(2),ro(ng),ro(ni),  &
                  ve(2),ve(ng),ve(ni),  &
                        cpa(ng),cpa(ni),  &
                        cva(ng),cva(ni),  &
                  ve(2)*fmi/fo,ve(ng)*fmi/fo,ve(ni)*fmi/fo,  &
                        cpa(ng)*fmi/fo,cpa(ni)*fmi/fo,  &
                        cva(ng)*fmi/fo,cva(ni)*fmi/fo,  &
                  ve(2)*fmi,ve(ng)*fmi,ve(ni)*fmi,  &
                        cpa(ng)*fmi,cpa(ni)*fmi,  &
                        cva(ng)*fmi,cva(ni)*fmi,  &
                  ve(2)/g,ve(ng)/g,ve(ni)/g,  &
                        cpa(ng)/g,cpa(ni)/g,  &
                        cva(ng)/g,cva(ni)/g
   91 format(///,5x,'RESULTADOS PRINCIPAIS',//,     &
             t3,'variável',t20,'entrada',t37,'garganta',t54,'saída',/,  &
             t3,'x / L',t18,3(1pe17.9),/,  &
             t3,'A / At',t18,3(1pe17.9),/,  &
             t3,'y / Dt',t18,3(1pe17.9),/,  &
             t3,'Mach',t18,3(1pe17.9),/,  &
             t3,'u / ao',t18,3(1pe17.9),/,  &
             t3,'p / poi',t18,3(1pe17.9),/,  &
             t3,'T / Toi',t18,3(1pe17.9),/,  &
             t3,'ro / roi',t18,3(1pe17.9),/, &
             t3,'ve (m/s)',t18,3(1pe17.9),/,  &
             t3,'c <pa> (m/s)',t35,2(1pe17.9),/,  &
             t3,'c <vácuo> (m/s)',t35,2(1pe17.9),/,  &
             t3,'Cf <ve>',t18,3(1pe17.9),/,  &
             t3,'Cf <pa>',t35,2(1pe17.9),/,  &
             t3,'Cf <vácuo>',t35,2(1pe17.9),/,  &
             t3,'F <ve> (N)',t18,3(1pe17.9),/,  &
             t3,'F <pa> (N)',t35,2(1pe17.9),/,  &
             t3,'F <vácuo) (N)',t35,2(1pe17.9),/,  &
             t3,'Is <ve> (s)',t18,3(1pe17.9),/,  &
             t3,'Is <pa> (s)',t35,2(1pe17.9),/,  &
             t3,'Is <vácuo> (s)',t35,2(1pe17.9))

      OPEN(7,FILE='saida1.dat',FORM='FORMATTED')
      write(6,110)
  110 format(///,t6,'x / L',t26,'A / At',t46,'y / Dt',/)
      do i=2,ni
        write(6,111) x(i),aa(i),y(i)
        write(7,111) x(i),aa(i),y(i)
      end do
  111 format(6(1pe20.9))
      close(7)

      write(6,112)
  112 format(///,t6,'x / L',t26,'Mach',t46,'u / ao',/)
      do i=2,ni
        write(6,111) x(i),m(i),u(i)
      end do

      OPEN(7,FILE='saida3.dat',FORM='FORMATTED')
      write(6,113)
  113 format(///,t6,'x / L',t26,'p / poi',t46,'T / Toi',t66,'ro / roi',/)
      do i=2,ni
        write(6,111) x(i),p(i),t(i),ro(i)
        write(7,111) x(i),p(i),t(i),ro(i),m(i),u(i)
      end do
      close(7)

      OPEN(7,FILE='saida4.dat',FORM='FORMATTED')
      write(6,114)
  114 format(///,t6,'x / L',t26,'ve',t46,'c (pa)',t66,'c (vácuo)',/)
      do i=2,ni
        write(6,111) x(i),ve(i),cpa(i),cva(i)
        if(i >= ng) write(7,111) x(i),ve(i),cpa(i),cva(i)
      end do
      close(7)

      OPEN(7,FILE='saida5.dat',FORM='FORMATTED')
      write(6,115)
  115 format(///,t6,'x / L',t26,'Cf (ve)',t46,'Cf (pa)',t66,'Cf (vácuo)',/)
      do i=2,ni
        write(6,111) x(i),ve(i)*fmi/fo,cpa(i)*fmi/fo,cva(i)*fmi/fo
        if(i >= ng) write(7,111) x(i),ve(i)*fmi/fo,cpa(i)*fmi/fo,cva(i)*fmi/fo
      end do
      close(7)

      OPEN(7,FILE='saida6.dat',FORM='FORMATTED')
      write(6,116)
  116 format(///,t6,'x / L',t26,'empuxo (ve) [N]',t46,'empuxo (pa) [N]',   &
                t66,'empuxo (vácuo) [N]',/)
      do i=2,ni
        write(6,111) x(i),ve(i)*fmi,cpa(i)*fmi,cva(i)*fmi
        if(i >= ng) write(7,111) x(i),ve(i)*fmi,cpa(i)*fmi,cva(i)*fmi
      end do
      close(7)

      OPEN(7,FILE='saida7.dat',FORM='FORMATTED')
      write(6,117)
  117 format(///,t6,'x / L',t26,'Is (ve)',t46,'Is (pa)',t66,'Is (vácuo)',/)
      do i=2,ni
        write(6,111) x(i),ve(i)/g,cpa(i)/g,cva(i)/g
        if(i >= ng) write(7,111) x(i),ve(i)/g,cpa(i)/g,cva(i)/g
      end do
      close(7)
 
      IF(ISHOCK.EQ.1) THEN
        WRITE(6,77) TO1/TO1,PO1/PO1,RO1/RO1
        WRITE(6,79) TO2/TO1,PO2/PO1,RO2/RO1
        WRITE(6,80)
        WRITE(6,82) X(IC)
        WRITE(6,84) MC1,MC2,MC2-MC1
        WRITE(6,86) TC1/TO1,TC2/TO1,(TC2-TC1)/TO1
        WRITE(6,88) PC1/PO1,PC2/PO1,(PC2-PC1)/PO1
        WRITE(6,89) ROC1/RO1,ROC2/RO1,(ROC2-ROC1)/RO1
      END IF
 
   77 FORMAT(///,5X,'CONDIÇÕES DE ESTAGNAÇÃO ANTES DO CHOQUE',//, &
            1pe20.9,' = temperatura: T/To1',/,			  &
            1pe20.9,' = pressão: p/po1',/,			  &
            1pe20.9,' = massa específica: Ro/Ro1')
   79 FORMAT(//,5X,'CONDIÇÕES DE ESTAGNAÇÃO DEPOIS DO CHOQUE',//, &
            1pe20.9,' = temperatura: T/To1',/,			   &
            1pe20.9,' = pressão: p/po1',/,			   &
            1pe20.9,' = massa específica: Ro/Ro1')
   80 FORMAT(///,5X,'INFORMAÇÕES SOBRE A ONDA DE CHOQUE',/)
   82 FORMAT(1pe20.9,' = posição: x/L',/)
   84 FORMAT(1pe20.9,' = Mach antes do choque',/, &
             1pe20.9,' = Mach após o choque',/, &
             1pe20.9,' = variação',/)
   86 FORMAT(1pe20.9,' = temperatura antes do choque: T/To1',/, &
             1pe20.9,' = temperatura após o choque: T/To1',/, &
             1pe20.9,' = variação',/)
   88 FORMAT(1pe20.9,' = pressão antes do choque: p/po1',/, &
             1pe20.9,' = pressão após o choque: p/po1',/, &
             1pe20.9,' = variação',/)
   89 FORMAT(1pe20.9,' = massa específica antes do choque: Ro/Ro1',/, &
             1pe20.9,' = massa específica após o choque: Ro/Ro1',/, &
             1pe20.9,' = variação',/)
 
      CLOSE(5)
      CLOSE(6)
 
      chamada = systemqq('notepad mach1.s')

      chamada = systemqq('wgnuplot saida1.gnu')
      chamada = systemqq('wgnuplot saida3.gnu')
      chamada = systemqq('wgnuplot saida4.gnu')
      chamada = systemqq('wgnuplot saida5.gnu')
      chamada = systemqq('wgnuplot saida6.gnu')
      chamada = systemqq('wgnuplot saida7.gnu')

      int4 = setexitqq(qwin$exitnopersist)
 
      END
 
!***********************************************************************
! *** DADOS ***
!              LE OS DADOS DO PROBLEMA
!***********************************************************************
      SUBROUTINE DADOS(RU,TO1,PO1,GAMA,NI,R,RO1,MM,PA,ISHOCK,IC,IS,  &
                       LC,LT,DT,RCOS)
      IMPLICIT REAL*8(A-H,O-Z)
	  real*8 mm,LC,LT
 
      READ(5,*) TO1
      READ(5,*) PO1
      READ(5,*) PA
      READ(5,*) MM
      READ(5,*) GAMA
      READ(5,*) NI
      READ(5,*) LC
      READ(5,*) LT
      READ(5,*) DT
      READ(5,*) RCOS
      READ(5,*) ISHOCK
      READ(5,*) IC
      READ(5,*) IS

      R = RU / MM 
      RO1 = PO1/(R*TO1)
 
      WRITE(6,35)
      WRITE(6,40) TO1,PO1,RO1,RU,MM,R,GAMA,PA,LC,LT,DT,RCOS
      WRITE(6,50) NI,ISHOCK,IC,IS
 
   35 FORMAT(///,5X,'DADOS DE ENTRADA',//)
   40 FORMAT(3X,1PE17.9,' = Temperatura total na entrada (K)',/, &
             3X,1PE17.9,' = Pressão total na entrada (Pa)',/, &
             3X,1pe17.9,' = Massa específica total na entrada (kg/m3)',/, &
             3X,1pe17.9,' = Constante universal dos gases (J/kmol.K)',/, &
             3X,1pe17.9,' = Massa molecular do gás (kg/kmol)',/, &
             3X,1pe17.9,' = Constante do gás (J/kg.K)',/, &
             3X,1pe17.9,' = Razão entre os calores específicos (adim.)',/, &
             3X,1pe17.9,' = Pressão ambiente (Pa)',/,  &
             3X,1pe17.9,' = Comprimento da câmara de combustão (m)',/, &
             3X,1pe17.9,' = Comprimento da tubeira (m)',/, &
             3X,1pe17.9,' = Diâmetro da garganta (m)',/, &
             3X,1pe17.9,' = Amplitude da cossenóide (m)',/)
   50 FORMAT(14x,i6,' = Número total de pontos + 1',/, &
             14x,i6,' = ishock (0=sem choque, 1=com)',/,	&
             14x,i6,' = Número de pontos + 1 até o choque',/, &
             14x,i6,' = is (0=saída supersônica, 1=subsônica)')
 
      END

!***********************************************************************
! *** PROPRIEDADES ***
!                     CALCULA A PRESSAO E TEMPERATURA ESTATICAS EM FUNCAO     
!                     DAS CONDICOES DE ESTAGNACAO E DO NUMERO DE MACH LOCAL
!***********************************************************************
      SUBROUTINE PROPRI(M,TO,PO,T,P,GAMA)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 M
 
      EXP = GAMA / (GAMA - 1.0D0)
      CTE = 1.0D0 + (GAMA - 1.0D0) * (M ** 2) / 2.0D0
      T = TO / CTE
      P = PO / (CTE ** EXP)
 
      END

!***********************************************************************
! *** NEWTON ***
!               DETERMINA O NUMERO DE MACH EM FUNCAO DA RAZAO DE AREAS
!               USANDO O METODO DE NEWTON-RAPHSON
!***********************************************************************
      SUBROUTINE NEWTON(GAMA,A,IND,M,INP)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 M
 
      FDEN(M,E,C) =  M * ( 1.0D0 + E * M * M) ** C
      F(DEN,D,A) = D / DEN - A
      G(E,M) = 1.0D0 + E * M * M
 
      NITM = 500
      TOL = 1.0D-14
 
      IF(IND.EQ.1)M = 1.2D0
      IF(IND.EQ.0)M = 0.4D0

      E = (GAMA - 1.0D0) / 2.0D0
      C = (GAMA + 1.0D0) / ( 2.0D0 - 2.0D0 * GAMA)
      D = ((GAMA + 1.0D0) / 2.0D0) ** C

      DO I = 1,NITM
        FAT1 = M * ( C * G(E,M) ** (C - 1.0D0) * 2.0D0 * E * M)
        DEN = FDEN(M,E,C)
        FLINHA = - D * (G(E,M) ** C + FAT1) / (DEN * DEN)
        DELTAM = - F(DEN,D,A) / FLINHA
        M = M + DELTAM
        IF(M.LT.0.0D0) M = 0.01D0
        IF(DABS(DELTAM).LT.TOL) EXIT
      END DO

      WRITE(6,5) INP,I,DELTAM
    5 FORMAT(6X,i6,' = ponto',5X,i6,' = iterações', &
                              5X,1pe10.2,' = variação')
 
      END

!***********************************************************************
! *** CHOQUE ***
!               DETERMINA AS PROPRIEDADES ESTATICAS APOS A ONDA 
!               DE CHOQUE NORMAL
!***********************************************************************
      SUBROUTINE CHOQUE(M,T,P,R,GAMA,MC1,TC1,PC1,ROC1,MC2,TC2,PC2,  &
                        ROC2,AG,A)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 M,MC1,MC2
 
      MC1 = M
      TC1 = T
      PC1 = P
      ROC1 = P/(R*T)
 
      C1 = 2.D0 + (GAMA - 1.D0) * (MC1 ** 2)
      C2 = 2.D0 * GAMA * (MC1 ** 2) - (GAMA - 1.D0)
      C3 = ((GAMA + 1.D0) ** 2) * (MC1 ** 2)
      C4 = (GAMA + 1.D0) * (MC1 ** 2)
 
      MC2 = DSQRT(C1/C2)
      TC2 = TC1 * C1 * C2 / C3
      PC2 = PC1 * C2 / (GAMA + 1.D0)
      ROC2 = ROC1 * C4 / C1
 
      EXP = 0.5D0 * (GAMA + 1.0D0)/ (GAMA - 1.0D0)
      AUX = 1.D0 + 0.5D0 * (GAMA - 1.D0) * (MC2**2)
      AG = A/((((2.D0/(GAMA+1.D0))*AUX)**EXP)/MC2)
 
      END

!***********************************************************************
! *** ESTAGNACAO ***
!                   CALCULA AS PROPRIEDADES DE ESTAGNACAO APOS A ONDA
!                   DE CHOQUE NORMAL
!***********************************************************************
      SUBROUTINE ESTAG(M1,T1,P1,GAMA,TO,PO)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 M1
 
      EXP = GAMA / (GAMA - 1.0D0)
      CTE = 1.0D0 + (GAMA - 1.0D0) * (M1 ** 2)/ 2.0D0
      TO = T1 * CTE
      PO = P1 * (CTE ** EXP)
 
      END

!***********************************************************************
! *** GRID ***
!              TUBEIRA COSSENOIDAL
!***********************************************************************
      SUBROUTINE GRID(X,Y,NIP1,PI,NIG,L,LC,LT,DT,RCOS)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(NIP1),X(NIP1)
      REAL*8 LC,LT,L
 
      NI = NIP1 - 1
      NIG = 3*(NI-2)/5+2
      NIC = (NI-2)/5+2
 
      RT = 0.5D0*DT
      L  = LC + LT
      RCA = 2.0D0 * RCOS + RT
 
      RC = (LT**2)/(4.0D0 * RCOS * (PI**2))
 
      DX = L/(NI-2)
      X(1) = - DX
      DO NP = 2, NI
        NPW = NP - 1
        X(NP) = X(NPW) + DX
        IF(NP.LE.NIC) THEN
          Y(NP) = RCA
        ELSE
          A1 = 2.0D0 * PI * (X(NP) - LC) / LT
          Y(NP) = RT + RCOS * (1.0D0 + DCOS(A1))
        END IF
      END DO
 
      Y(1) = Y(2)
      X(NIP1) = X(NI) + DX
      Y(NIP1) = Y(NI)
 
      END

