
Program TRAJETORIA_1p0

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

! Carlos H. Marchi
! DEMEC/UFPR, Curitiba, PR
! versão original 1.0:  5 Nov 08
! última alteração:     8 Nov 08

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

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

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

implicit none

integer,parameter :: N = 1000000 

! Dados

character*50 :: caso   ! nome do arquivo de saída
character*50 :: nome   ! nome do arquivo de dados
character*90 :: titulo ! título de identificação da análise

integer :: Iq  ! Número de iterações na fase propulsada
integer :: Ib  ! Número de iterações na fase balística
integer :: w   ! Freqüência de escrita dos resultados de campo

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)

! 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

real*8,dimension(:),allocatable :: &
   t, & ! tempo (s)
   h, & ! altura (m)
   V, & ! velocidade (m/s)
   a, & ! aceleração (m/s2)
   M, & ! massa (kg)
   P, & ! peso (N)
   E, & ! empuxo (N)
   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)

integer :: i   ! Número de um ponto qualquer da trajetória
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

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

! *** LÊ OS DADOS ***

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

! Lê o nome do arquivo com os dados numéricos
open(8,file="Trajetoria_1p0_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,*) Ec
read(8,*) tq
read(8,*) Cd
read(8,*) g
read(8,*) par
read(8,*) Tar
read(8,*) Rar
read(8,*) dtq
read(8,*) dtb
read(8,*) Iq
read(8,*) Ib
read(8,*) w
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,*) titulo

close(8)

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

! *** INICIALIZA VARIÁVEIS E CONSTANTES

!dtq = tq / Nq	

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) )

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

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

open(10,file=caso)

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

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

write(10,1) trim(adjustl(caso)), trim(adjustl(titulo)), Mf, Mp, Mf+Mp, dr, &
            Ec, tq, Ec*tq, Cd, g, par, Tar, Rar, ro, dtq, dtb, Iq, Ib, w

1 format(/,1x, 'Saída do programa TRAJETORIA 1.0, versão de 8 Nov 2008',//, &
           1x, 'Caso   = ', a, //, &
           1x, 'Título = ', a, //, &
        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)',/, &
        1pe18.6, ' = E   = Empuxo constante do motor (N)',/, &
        1pe18.6, ' = tq  = Tempo de queima do propelente (s)',/, &
        1pe18.6, ' = It  = Impulso total (Ns)',/, &
        1pe18.6, ' = Cd  = Coeficiente de arrasto (adim.)',/, &
        1pe18.6, ' = g   = Aceleração gravitacional local (m/s2)',/, &
        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, ' = 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')

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

! *** 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
E(i) = 0.0d0
D(i) = 0.0d0

do i = 1, Nq

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

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

   P(i) = M(i) * g

   E(i) = Ec

   V(i) = V(i-1)

   do it = 1, Iq

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

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

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

   end do 

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

end do   

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

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

do i = Nq+1, N

   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

      D(i) = Cd * 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

end do   

!Nba = i
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

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

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

!do i = Nba+1, N
do i = NH+1, N

   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

      D(i) = Cd * 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) < 0.0d0 ) exit

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

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

! *** ESCREVE RESULTADOS ***

open(11,file='Trajetoria_1p0_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)',/)

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)')

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)
   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)

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)
   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)
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)
   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)

12 format( i8, 8(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)',/)

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)')

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)
   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)

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)
   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)
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)
   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)

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)',/)

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)')

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)
   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)

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)
   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)
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)
   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)

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), &
             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, ' = hq/H (adim.)',/, &
        1pe18.6, ' = hq  = altura no fim da queima (m)',/, &
        1pe18.6, ' = H   = apogeu (m)')

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

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

! gráfico de h
if ( g_h /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_h.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_h.gnu')
end if

! gráfico de V
if ( g_V /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_V.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_V.gnu')
end if

! gráfico de a
if ( g_a /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_a.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_a.gnu')
end if

! gráfico de M
if ( g_M /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_M.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_M.gnu')
end if

! gráfico de forças
if ( g_F /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_F.gnu')
    do ver = 1, 9
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_F.gnu')
end if

! gráfico de E
if ( g_E /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_E.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_E.gnu')
end if

! gráfico de P
if ( g_P /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_P.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_P.gnu')
end if

! gráfico de D
if ( g_D /= 0 ) then
    open(12,file='Trajetoria_1p0_grafico_D.gnu')
    do ver = 1, 6
       read(12,*)
    end do
    write(12,*) "set title '", trim(adjustl(titulo)), "'"
    write(12,*) "replot"
    close(12)
    ver = system('wgnuplot Trajetoria_1p0_grafico_D.gnu')
end if

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

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

close (10)

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

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

ver = system('notepad '//caso)

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

end program TRAJETORIA_1p0
