Program TRAJETORIA_2p0

! Cálculo de trajetória de espaçomodelos

! Carlos H. Marchi, Antonio C. Foltran, Diego Fernando Moro
! DEMEC/UFPR, Curitiba, PR
! versão original 1.0:  5 Nov 08
! versão atual    2.0:  2 Ago 13
! última alteração:    13 Abr 16: Inserida linha 191 e modificado formato na
                       !subrotina plot (linhas 983 e 984). A. Foltran e D. Moro
! -----------------------------------------------------------------------------

! *** DEFINE AS VARIÁVEIS DO PROGRAMA ***

use portlib ! inserção da biblioteca básica do Fortran

implicit none

integer,parameter :: N = 10000000 ! número de pontos da trajetória 
integer,parameter :: Q =    10000 ! número de pontos de tabelas 

! Dados

character*90 :: caso    ! nome do arquivo de saída
character*90 :: nome    ! nome do arquivo de dados
character*90 :: tabela  ! nome do arquivo com tabela empuxo x tempo
character*90 :: tabela2 ! nome do arquivo com tabela Cd x V
character*90 :: titulo  ! título de identificação da análise

character*90 :: PASTA, ARQUIVO, DADOS ! Variáveis para uso dos gráficos

integer :: Iq  ! Número de iterações na fase propulsada
integer :: Ib  ! Número de iterações na fase balística
integer :: ic  ! =0, Cd constante; =1, Cd variável
integer :: ie  ! =0, E constante; =1, E variável
integer :: ig  ! =0, g constante; =1, g=3.9887d14 / ((6.3776d6 + altura)**2.0d0)
integer :: iro ! =0, ro=par/(Rar*(Tar+273.15)); =1, ro=1.225d0-1.175d-4*h+4.2851d-9*(h**2.0d0)-6.504d-14*(h**3.0d0)
integer :: w   ! Freqüência de escrita dos resultados de campo
integer :: solver ! =0, Euler; =1, Runge-Kuta 4°ordem

real*8  :: Mf  ! Massa final (kg)
real*8  :: Mp  ! Massa de propelente (kg)
real*8  :: dr  ! Diâmetro de referência (m)
real*8  :: Ec  ! Empuxo constante do motor (N)
real*8  :: tq  ! Tempo de queima do propelente (s)
real*8  :: Cd  ! Coeficiente de arrasto (adim.)
real*8  :: g   ! Aceleração gravitacional local (m/s2)
real*8  :: par ! Pressão atmosférica local (Pa)
real*8  :: Tar ! Temperatura atmosférica local (oC)
real*8  :: Rar ! Constante dos gases perfeitos do ar (J/kg.K)
real*8  :: dtq ! Intervalo de tempo na fase propulsada (s)
real*8  :: dtb ! Intervalo de tempo na fase balística (s)
real*8  :: Vms ! velocidade máxima na subida (m/s)
real*8  :: Vmd ! velocidade máxima na descida (m/s)
real*8  :: k1h ! coeficiente k1 do método Runge-Kutta usado para cálculo da altitude
real*8  :: k2h ! coeficiente k2 do método Runge-Kutta usado para cálculo da altitude
real*8  :: k3h ! coeficiente k3 do método Runge-Kutta usado para cálculo da altitude
real*8  :: k4h ! coeficiente k4 do método Runge-Kutta usado para cálculo da altitude
real*8  :: k1v ! coeficiente k1 do método Runge-Kutta usado para cálculo da velocidade
real*8  :: k2v ! coeficiente k2 do método Runge-Kutta usado para cálculo da velocidade
real*8  :: k3v ! coeficiente k3 do método Runge-Kutta usado para cálculo da velocidade
real*8  :: k4v ! coeficiente k4 do método Runge-Kutta usado para cálculo da velocidade

! Resultados

integer :: Nq  ! Número do ponto da trajetória no fim da queima
integer :: NH  ! Número do ponto da trajetória no apogeu
integer :: NI  ! Número do ponto da trajetória no impacto
integer :: is  ! Número do ponto de Vms
integer :: id  ! Número do ponto de Vmd

real*8,dimension(:),allocatable :: &
   t,  & ! tempo (s)
   tt, & ! tempo da tabela empuxo x tempo (s)
   h,  & ! altura (m)
   V,  & ! velocidade (m/s)
   Vt, & ! V da tabela Cd x V (adim.)
   a,  & ! aceleração (m/s2)
   M,  & ! massa (kg)
   P,  & ! peso (N)
   E,  & ! empuxo (N)
   Et, & ! empuxo da tabela empuxo x tempo (N)
   Ct, & ! Cd da tabela Cd x V (adim.)
   Cv, & ! Cd (adim.)
   D     ! arrasto (N)

! Variáveis auxiliares

real*8  :: pi  ! Constante matemática
real*8  :: Ar  ! Área de referência (m2)
real*8  :: ro  ! Massa específica do ar (kg/m3)
real*8  :: t1, e1, C1, V1

integer :: i   ! Número de um ponto qualquer da trajetória
integer :: j   ! Número de um ponto qualquer da curva empuxo x tempo
integer :: Nt  ! Número total de pontos da curva empuxo x tempo
integer :: N2  ! Número total de pontos da curva Cd x V
integer :: it  ! Número da iteração
integer :: ver
integer :: g_h, g_V, g_a, g_M, g_F, g_E, g_P, g_D, g_Ev, g_Cv, g_Cd
integer :: flag, io

real*8  :: tcpu ! tempo de processamento

! -----------------------------------------------------------------------------

! PASTA DOS ARQUIVOS AUXILIARES DO GNUPLOT
PASTA = "GNU"

! *** LÊ OS DADOS ***

! Mostra o conteúdo do arquivo de dados do programa
ver = system("Notepad Trajetoria_2p0_dados.txt")

! Lê o nome do arquivo com os dados numéricos
open(8,file="Trajetoria_2p0_dados.txt")
read(8,*) nome
close(8)

! Mostra o conteúdo do arquivo com os dados numéricos
ver = system('notepad '//nome)

! Lê os dados numéricos

open(8,file=nome)

read(8,*) caso
read(8,*) Mf
read(8,*) Mp
read(8,*) dr
read(8,*) ie
read(8,*) Ec
read(8,*) tq
read(8,*) ic
read(8,*) Cd
read(8,*) ig
read(8,*) g
read(8,*) iro
read(8,*) par
read(8,*) Tar
read(8,*) Rar
read(8,*) dtq
read(8,*) dtb
read(8,*) solver
read(8,*) Iq
read(8,*) Ib
read(8,*) w
read(8,*) g_Ev
read(8,*) g_Cv
read(8,*) g_h
read(8,*) g_V
read(8,*) g_a
read(8,*) g_M
read(8,*) g_F
read(8,*) g_E
read(8,*) g_P
read(8,*) g_D
read(8,*) g_Cd
read(8,*) titulo
read(8,*) tabela
read(8,*) tabela2

close(8)

! -----------------------------------------------

! *** INICIALIZA VARIÁVEIS E CONSTANTES

j = 0

flag = 0

Vms = 0
Vmd = 0

pi = dacos(-1.0d0)

Ar = pi * (dr**2) / 4

ro = par / ( (Tar+273.15d0) * Rar )

Nq = tq / dtq	
if ( Nq*dtq < tq ) Nq = Nq+1

allocate ( a(0:N), D(0:N), V(0:N), h(0:N), t(0:N), P(0:N), M(0:N), E(0:N), Cv(0:N) )
allocate ( tt(0:Q), Et(0:Q), Ct(0:Q), Vt(0:Q) )

! -----------------------------------------------------------------------------

! *** CRIA O ARQUIVO DE SAÍDA ***
caso = trim(adjustl(caso)) // '.txt'
open(10,file=caso)

! -----------------------------------------------------------------------------

! *** ESCREVE OS DADOS NO ARQUIVO DE SAÍDA ***

write(10,1) trim(adjustl(caso)), trim(adjustl(titulo))
1 format(/,1x, 'Saída do programa TRAJETORIA 2.0, versão de 02 Ago 2013',//, &
           1x, 'Caso   = ', a, //, &
           1x, 'Título = ', a  )

if ( ie == 1 ) write(10,2) trim(adjustl(tabela))
2 format(/,1x, 'Tabela  E x t = ', a )

if ( ic == 1 ) write(10,3) trim(adjustl(tabela2))
3 format(/,1x, 'Tabela  Cd x V = ', a )

write(10,4) Mf, Mp, Mf+Mp, dr, ie, Ec, tq, Ec*tq, ic, Cd, ig, &
            g, iro, par, Tar, Rar, ro, dtq, dtb, solver, Iq, Ib, w

4 format(/,1pe18.6, ' = Mf  = Massa final (kg)',/, &
           1pe18.6, ' = Mp  = Massa de propelente (kg)',/, &
           1pe18.6, ' = Mo  = Massa inicial (kg) = Mf + Mp',/, &
           1pe18.6, ' = dr  = Diâmetro de referência (m)',/, &
            i7,11x, ' = ie  : =0, E constante; =1, E variável',/, &
           1pe18.6, ' = E   = Empuxo constante do motor (N)',/, &
           1pe18.6, ' = tq  = Tempo de queima do propelente (s)',/, &
           1pe18.6, ' = It  = Impulso total (Ns)',/, &
            i7,11x, ' = ic  : =0, Cd constante; =1, Cd variável',/, &
           1pe18.6, ' = Cd  = Coeficiente de arrasto (adim.)',/, &
            i7,11x, ' = ig  : =0, g constante; =1, g=3.9887d14 / ((6.3776d6 + altura)**2.0d0)',/,&
           1pe18.6, ' = g   = Aceleração gravitacional local (m/s2)',/, &
            i7,11x, ' = iro : =0, ro=par/(Rar*(Tar+273.15)); =1, ro=0.1225-1.175d-4*h+4.2851d-9*(h**2.0d0)-6.504d-14*(h**3.0d0)',/,&
           1pe18.6, ' = par = Pressão atmosférica local (Pa)',/, &
           1pe18.6, ' = Tar = Temperatura atmosférica local (oC)',/, &
           1pe18.6, ' = Rar = Constante dos gases perfeitos do ar (J/kg.K)',/, &
           1pe18.6, ' = ro  = Massa específica do ar (kg/m3)',/, &
           1pe18.6, ' = dtq = Intervalo de tempo na fase propulsada (s)',/, &
           1pe18.6, ' = dtb = Intervalo de tempo na fase balística (s)',/, &
            i7,11x, ' = solver : =0, Euler; =1, Runge-Kuta 4°ordem',/,&
            i7,11x, ' = Iq  = Número de iterações na fase propulsada',/, &
            i7,11x, ' = Ib  = Número de iterações na fase balística',/, &
            i7,11x, ' = w   = Freqüência de escrita dos resultados de campo')

! -----------------------------------------------------------------------------

if ( ie == 0 ) goto 99

! Mostra o conteúdo do arquivo de dados da curva empuxo x tempo
ver = system('notepad '//tabela)

! Lê os dados da curva empuxo x tempo

open(8,file=tabela)

do while ( flag == 0 )

   read(8,*,iostat=io) t1, e1

   tt(j) = t1
   Et(j) = e1

   if ( io < 0 ) flag = 1

   j = j + 1

end do

close(8)

Nt = j - 2

write(10,103) Nt+1
103 format(/,"Número de pontos na tabela de empuxo =", i5)

write(10,101) 
101 format (/, 3x, "Curva de empuxo x tempo", //, &
      3x,'ponto', t12,'t (s)', t27,'E (N)')
do j = 0, Nt
   write(10,102) j, tt(j), Et(j)
   102 format( i8, 2(1pe15.6) )
end do

! gráfico de E x t (dados)
ARQUIVO = "Trajetoria_2p0_grafico_Ev.gnu"
DADOS = trim(adjustl(tabela))
if ( g_Ev /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 5)


99 continue

! -----------------------------------------------------------------------------

if ( ic == 0 ) goto 98

j = 0
flag = 0

! Mostra o conteúdo do arquivo de dados da curva Cd x V
ver = system('notepad '//tabela2)

! Lê os dados da curva Cd x V

open(8,file=tabela2)

do while ( flag == 0 )

   read(8,*,iostat=io) V1, C1

   Vt(j) = V1 / 3.6d0
   Ct(j) = C1

   if ( io < 0 ) flag = 1

   j = j + 1

end do

close(8)

N2 = j - 2

write(10,104) N2+1
104 format(/,"Número de pontos na tabela de Cd =", i5)

write(10,105) 
105 format (/, 3x, "Curva de Cd x V", //, &
      3x,'ponto', t12,'V (km/h)', t27,'Cd (adim.)')
do j = 0, N2
   write(10,102) j, Vt(j)*3.6d0, Ct(j)
end do

! gráfico de Cd x V (dados)
ARQUIVO = "Trajetoria_2p0_grafico_Cv.gnu"
DADOS = trim(adjustl(tabela2))
if ( g_Cv /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 5)

98 continue

! -----------------------------------------------------------------------------

tcpu = timef()

! *** FASE PROPULSADA ***

i    = 0
t(i) = 0.0d0
V(i) = 0.0d0
a(i) = 0.0d0
h(i) = 0.0d0
M(i) = Mf+Mp
P(i) = M(i)*g
if ( ie == 0 ) then
   E(i) = 0
else
   E(i) = Et(0)
end if
D(i) = 0.0d0
if ( ic == 0 ) then
   Cv(i) = Cd
else
   Cv(i) = Ct(0)
end if

do i = 1, Nq
    
    if ( solver == 0 ) then !Método de Euler 
       
       t(i) = t(i-1) + dtq

       M(i) = Mf + Mp * ( 1 - t(i)/tq )

       if ( ie == 0 ) then
          E(i) = Ec
       else
          do j = 0, Nt
             if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
             end if
          end do   
          if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
       end if

       V(i) = V(i-1)

       do it = 1, Iq

          h(i) = h(i-1) + dtq * V(i)
          
          call propriedades
          
          P(i) = M(i) * g

          D(i) = Cv(i) * ro * ( V(i)**2 ) * Ar / 2

          a(i) = -g + (E(i) - D(i)) / M(i)

          if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0   ! nova condição de contato com a plataforma de lançamento
          
          V(i) = V(i-1) + dtq * a(i)   
          
          !If ( V(i) < 0 ) V(i) = 0.0d0   ! antiga condição de contato com a plataforma de lançamento

       end do 

       if ( V(i) > Vms ) then
          Vms = V(i)
          is  = i
       end if

       h(i) = h(i-1) + dtq * V(i)
    
    else    ! Método Runge-Kutta 4°ordem
        
        t(i) = t(i-1)
        V(i) = V(i-1)
        k1h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k1h
        if ( ie == 0 ) then
            E(i) = Ec
        else
            do j = 0, Nt
            if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
            end if
            end do   
            if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
        end if
        M(i) = Mf + Mp*(1 - t(i)/tq)
            
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0
        k1v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k1v
        k2h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k2h
        t(i) = t(i-1) + 5.0d-1*dtq
        if ( ie == 0 ) then
            E(i) = Ec
        else
            do j = 0, Nt
            if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
            end if
            end do   
            if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
        end if
        M(i) = Mf + Mp*(1 - t(i)/tq)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0
        k2v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k2v
        t(i) = t(i-1) + 5.0d-1*dtq
        k3h = dtq*V(i)
        h(i) = h(i-1) + k3h
        if ( ie == 0 ) then
            E(i) = Ec
        else
            do j = 0, Nt
            if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
            end if
            end do   
            if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
        end if
        M(i) = Mf + Mp*(1 - t(i)/tq)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0
        k3v = dtq*a(i)
        V(i) = V(i-1) + k3v
        t(i) = t(i-1) + dtq
        k4h = dtq*V(i)
        h(i) = h(i-1) + (1.0d0/6.0d0)*(k1h + 2.0d0*(k2h + k3h) + k4h)
        if ( ie == 0 ) then
            E(i) = Ec
        else
            do j = 0, Nt
            if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
            end if
            end do   
            if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
        end if
        M(i) = Mf + Mp*(1 - t(i)/tq)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0
        k4v = dtq*a(i)
        V(i) = V(i-1) + (1.0d0/6.0d0)*(k1v + 2.0d0*(k2v + k3v) + k4v)
        if ( ie == 0 ) then
            E(i) = Ec
        else
            do j = 0, Nt
            if ( t(i) >= tt(j) .and. t(i) < tt(j+1) ) then
                E(i) = Et(j) + (Et(j+1)-Et(j))*(t(i)-tt(j))/(tt(j+1)-tt(j))
                exit
            end if
            end do   
            if ( t(i) > tt(Nt) ) E(i) = Et(Nt)
        end if
        M(i) = Mf + Mp*(1 - t(i)/tq)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        P(i) = M(i)*g
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if (a(i)<0.0d0 .and. h(i)<=0.0d0) a(i) = 0.0d0
        if ( V(i) > Vms ) then
            Vms = V(i)
            is  = i
        end if
    
    end if

end do   

! -----------------------------------------------------------------------------

! *** FASE BALÍSTICA ASCENDENTE ***

do i = Nq+1, N
    
    if ( solver == 0 ) then !Método de Euler
       
        t(i) = t(i-1) + dtb

        M(i) = Mf

        P(i) = M(i) * g

        E(i) = 0.0d0

        V(i) = V(i-1)

        do it = 1, Ib

          h(i) = h(i-1) + dtb * V(i)

          call propriedades
                
          D(i) = Cv(i) * ro * ( V(i)**2 ) * Ar / 2

          a(i) = - g - D(i) / M(i)

          V(i) = V(i-1) + dtb * a(i)   

        end do 

        h(i) = h(i-1) + dtb * V(i)

        if ( h(i) < h(i-1) ) exit
    
    else    ! Método de Runge-Kutta de 4°ordem
    
        t(i) = t(i-1)
        V(i) = V(i-1)
        k1h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k1h
        E(i) = 0.0d0
        M(i) = Mf
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        k1v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k1v
        k2h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k2h
        t(i) = t(i-1) + 5.0d-1*dtq
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        k2v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k2v
        t(i) = t(i-1) + 5.0d-1*dtq
        k3h = dtq*V(i)
        h(i) = h(i-1) + k3h
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        k3v = dtq*a(i)
        V(i) = V(i-1) + k3v
        t(i) = t(i-1) + dtq
        k4h = dtq*V(i)
        h(i) = h(i-1) + (1.0d0/6.0d0)*(k1h + 2.0d0*(k2h + k3h) + k4h)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) - D(i)) / M(i))
        k4v = dtq*a(i)
        V(i) = V(i-1) + (1.0d0/6.0d0)*(k1v + 2.0d0*(k2v + k3v) + k4v)
        
        call propriedades
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        P(i) = M(i)*g
        a(i) = -g + ((E(i) - D(i)) / M(i))
        if ( h(i) < h(i-1) ) exit
    
    end if

end do   

NH = i

t(i) = t(i-1) - (t(i)-t(i-1))*V(i-1)/(V(i)-V(i-1))

h(i) = h(i-2) + (h(i-1)-h(i-2))*(t(i)-t(i-2))/(t(i-1)-t(i-2))

a(i) = -g

D(i) = 0.0d0

V(i) = 0.0d0

if ( ic == 1 ) then
   Cv(i) = Ct(0)
else
   Cv(i) = Cd
end if

! -----------------------------------------------------------------------------

! *** FASE BALÍSTICA DESCENDENTE ***

!Cd = 1.5d0              ! utilizado caso haja sistema de recuperação (pára quedas) Neste caso simula-se sua abertura exatamente no apogeu.
!Ar = 1.9634954085d-1    ! utilizado caso haja sistema de recuperação (pára quedas) Neste caso simula-se sua abertura exatamente no apogeu.

do i = NH+1, N
    
    if ( solver == 0 ) then ! Método de Euler

        t(i) = t(i-1) + dtb

        M(i) = Mf

        E(i) = 0.0d0

        P(i) = M(i) * g

        V(i) = V(i-1)

        do it = 1, Ib

          h(i) = h(i-1) + dtb * V(i)

          v(i) = -V(i)  ! muda-se o sinal para que se consiga calcular o Cd dentro da subrotina propriedades
          call propriedades
          v(i) = -V(i)
          
          D(i) = Cv(i) * ro * ( V(i)**2 ) * Ar / 2

          a(i) = - g + D(i) / M(i)

          V(i) = V(i-1) + dtb * a(i)   

        end do 

        if ( V(i) < Vmd ) then
          Vmd = V(i)
          id  = i
        end if

        h(i) = h(i-1) + dtb * V(i)

        if ( h(i) < 0.0d0 ) exit
    
    else    ! Método de Runge-Kutta de 4°ordem
    
        t(i) = t(i-1)
        V(i) = V(i-1)
        k1h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k1h
        E(i) = 0.0d0
        M(i) = Mf
        
        v(i) = -V(i)  ! muda-se o sinal para que se consiga calcular o Cd dentro da subrotina propriedades
        call propriedades
        v(i) = -V(i)
                  
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) + D(i)) / M(i))
        k1v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k1v
        k2h = dtq*V(i)
        h(i) = h(i-1) + 5.0d-1*k2h
        t(i) = t(i-1) + 5.0d-1*dtq
        
        v(i) = -V(i)
        call propriedades
        v(i) = -V(i)
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) + D(i)) / M(i))
        k2v = dtq*a(i)
        V(i) = V(i-1) + 5.0d-1*k2v
        t(i) = t(i-1) + 5.0d-1*dtq
        k3h = dtq*V(i)
        h(i) = h(i-1) + k3h
        
        v(i) = -V(i)
        call propriedades
        v(i) = -V(i)
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) + D(i)) / M(i))
        k3v = dtq*a(i)
        V(i) = V(i-1) + k3v
        t(i) = t(i-1) + dtq
        k4h = dtq*V(i)
        h(i) = h(i-1) + (1.0d0/6.0d0)*(k1h + 2.0d0*(k2h + k3h) + k4h)
        
        v(i) = -V(i)
        call propriedades
        v(i) = -V(i)
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        a(i) = -g + ((E(i) + D(i)) / M(i))
        k4v = dtq*a(i)
        V(i) = V(i-1) + (1.0d0/6.0d0)*(k1v + 2.0d0*(k2v + k3v) + k4v)
        
        v(i) = -V(i)
        call propriedades
        v(i) = -V(i)
        
        D(i) = 5.0d-1*Cv(i)*ro*V(i)*V(i)*Ar
        P(i) = M(i)*g
        a(i) = -g + ((E(i) + D(i)) / M(i))
        if ( V(i) < Vmd ) then
            Vmd = V(i)
            id  = i
        end if
        
        if ( h(i) < 0.0d0 ) exit
    
    end if

end do   

NI = i

t(i) = t(i-1) - (t(i)-t(i-1))*h(i-1)/(h(i)-h(i-1))

V(i) = V(i-1) - (V(i)-V(i-1))*h(i-1)/(h(i)-h(i-1))

a(i) = a(i-1) - (a(i)-a(i-1))*h(i-1)/(h(i)-h(i-1))

D(i) = D(i-1) - (D(i)-D(i-1))*h(i-1)/(h(i)-h(i-1))

h(i) = 0.0d0

if ( ic == 1 ) then
   do j = 0, N2
      if ( dabs(V(i)) >= Vt(j) .and. dabs(V(i)) < Vt(j+1) ) then
         Cv(i) = Ct(j) + (Ct(j+1)-Ct(j))*(dabs(V(i))-Vt(j))/(Vt(j+1)-Vt(j))
         exit
      end if
   end do   
   if ( dabs(V(i)) > Vt(N2) ) Cv(i) = Ct(N2)
else
    Cv(i) = Cd
end if
Cv(i) = Cd

tcpu = timef()

! -----------------------------------------------------------------------------

! *** ESCREVE RESULTADOS ***

open(11,file='Trajetoria_2p0_graficos.dat')

write(10,11)
11 format(/,1x,'*** FASE PROPULSADA ***',//,  &
      7x,'i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
              t72,'M (g)', t87,'P (N)', t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)',/)

write(11,31)
31 format('# *** FASE PROPULSADA ***', /,  &
      5x,'# i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
                t72,'M (g)', t87,'P (N)', t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)')

i = 0
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)

do i = 1, Nq-1, w
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
end do

i = Nq
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)

12 format( i8, 9(1pe15.6) )

write(10,13)
13 format(/,1x,'*** FASE BALÍSTICA ASCENDENTE ***',//,  &
      7x,'i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
              t72,'M (g)', t87,'P (N)' ,t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)',/)

write(11,33)
33 format('# *** FASE BALÍSTICA ASCENDENTE ***', /,  &
      5x,'# i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
                t72,'M (g)', t87,'P (N)', t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)')

i = Nq+1
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)

do i = Nq+2, NH-1, w
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
end do

i = NH
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), -D(i), E(i), Cv(i)

write(10,14)
14 format(/,1x,'*** FASE BALÍSTICA DESCENDENTE ***',//,  &
      7x,'i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
              t72,'M (g)', t87,'P (N)' ,t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)',/)

write(11,34)
34 format('# *** FASE BALÍSTICA DESCENDENTE ***', /,  &
      5x,'# i', t12,'t (s)', t27,'h (m)',  t42,'V (km/h)', t57,'a (g)', &
                t72,'M (g)', t87,'P (N)', t102,'D (N)',   t117,'E (N)', t132,'Cd (adim.)')

i = NH+1
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)

do i = NH+2, NI-1, w
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)
end do

i = NI
   write(10,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)
   write(11,12) i, t(i), h(i), V(i)*3.6d+0, a(i)/g, M(i)*1000, -P(i), D(i), E(i), Cv(i)

close (11)

! -----------------------------------------------------------------------------

! *** ESCREVE PARÂMETROS PRINCIPAIS ***

write(10,21) t(NH)-t(Nq), t(NI)-t(NH), t(NH), t(NI), t(NH)/t(NI), &
             Vms*3.6d+0, t(is), Vmd*3.6d+0, t(id), h(Nq)/h(NH), h(Nq), h(NH)

21 format(//,1x, '*** PARÂMETROS PRINCIPAIS ***',//, &
        1pe18.6, ' = tba = Tempo entre fim da queima e apogeu (s)',/, &
        1pe18.6, ' = tbd = Tempo entre apogeu e impacto (s)',/, &
        1pe18.6, ' = tH  = Tempo até o apogeu (s)',/, &
        1pe18.6, ' = tI  = Tempo até o impacto (s)',/, &
        1pe18.6, ' = tH/tI (adim.)',/, &
        1pe18.6, ' = Vms = velocidade máxima na subida (km/h)',/, &
        1pe18.6, ' = tms = Tempo de vôo em Vms (s)',/, &
        1pe18.6, ' = Vmd = velocidade máxima na descida (km/h)',/, &
        1pe18.6, ' = tmd = Tempo de vôo em Vmd (s)',/, &
        1pe18.6, ' = hq/H (adim.)',/, &
        1pe18.6, ' = hq  = altura no fim da queima (m)',/, &
        1pe18.6, ' = H   = apogeu (m)')

write(10,39) tcpu
39 format(//,1pe14.3,' = tempo de CPU (s)')

! -----------------------------------------------

! *** GERA GRÁFICOS DA SOLUÇÃO ***

! gráfico de h
ARQUIVO = "Trajetoria_2p0_grafico_h.gnu"
DADOS = ' '
if ( g_h /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6) 

! gráfico de V
ARQUIVO = "Trajetoria_2p0_grafico_V.gnu"
if ( g_V /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de a
ARQUIVO = "Trajetoria_2p0_grafico_a.gnu"
if ( g_a /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de M
ARQUIVO = "Trajetoria_2p0_grafico_M.gnu"
if ( g_M /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de forças
ARQUIVO = "Trajetoria_2p0_grafico_F.gnu"
if ( g_F /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 9)

! gráfico de E
ARQUIVO = "Trajetoria_2p0_grafico_E.gnu"
if ( g_E /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de P
ARQUIVO = "Trajetoria_2p0_grafico_P.gnu"
if ( g_P /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de D
ARQUIVO = "Trajetoria_2p0_grafico_D.gnu"
if ( g_D /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! gráfico de Cd x t
ARQUIVO = "Trajetoria_2p0_grafico_Cd.gnu"
if ( g_Cd /= 0 ) call plot(PASTA, ARQUIVO, DADOS, titulo, 6)

! -----------------------------------------------

! *** FECHA O ARQUIVO DE SAÍDA ***

close (10)

! -----------------------------------------------

! *** MOSTRA O CONTEÚDO DO ARQUIVO DE SAÍDA ***

ver = system('notepad '//caso)

! -----------------------------------------------

contains
!----------------------------------------------------------------
    subroutine propriedades
    
        if ( ic == 1 ) then
            do j = 0, N2
                if ( V(i) >= Vt(j) .and. V(i) < Vt(j+1) ) then
                    Cv(i) = Ct(j) + (Ct(j+1)-Ct(j))*(V(i)-Vt(j))/(Vt(j+1)-Vt(j))
                    exit
                end if
            end do   
            if ( V(i) > Vt(N2) ) Cv(i) = Ct(N2)
        else
            Cv(i) = Cd
        end if
        
        if (iro == 1) then
            ro = 1.225d0 - 1.175d-4*h(i) + 4.2851d-9*(h(i)**2.0d0) - 6.504d-14*(h(i)**3.0d0)
        else
            ro = par / ( (Tar+273.15d0) * Rar )
        end if
            
        if (ig == 1) then
            g = 9.806563d0*(6378137 / (6378137 + h(i)))**2.0d0
            !g = 3.989371d14 / ((6378137d0 + h(i))**2.0d0)
        else
        end if
        
    end subroutine propriedades
!----------------------------------------------------------------
    subroutine plot(PASTA, ARQUIVO, DADOS, titulo, pulos)

    character*90, intent(in) :: PASTA, ARQUIVO, DADOS, titulo
    integer, intent (in) :: pulos
    integer :: pulo
    
    open(79, file = "./" // trim(adjustl(PASTA)) // "/" // trim(adjustl(ARQUIVO)) )
        do pulo = 1, pulos
            read(79,*)
        enddo
        if(DADOS /= " ") write(79,*) "plot '"// "./" // trim(adjustl(PASTA)) // "/", trim(adjustl(DADOS)), "'"
        write(79,'(A,A,A)') "set title '", trim(adjustl(titulo)), "'"
        write(79,'(A)') "replot"
    close(79)
    
    ver = system( '.\wgnuplot\wgnuplot ' // "./" // trim(adjustl(PASTA)) // "/" // trim(adjustl(ARQUIVO)) )

    end subroutine plot
!----------------------------------------------------------------
end program TRAJETORIA_2p0