
Program curva_empuxo_2p0

! Analisador de curvas de empuxo experimentais

! Carlos H. Marchi
!   DEMEC/UFPR, Curitiba, PR

! versão 1.0: 11 Set 2006
! versão 1.1: 28 Nov 2006
! versão 1.2:  9 Nov 2008

! versão atual 2.0   :  7 Jun 2012
! 
!   subversão  2.0.2 : 10 Jul 2012
!   subversão  2.0.3 : 18 Jul 2012
!   subversão  2.0.4 : 19 Jul 2012
!   subversão  2.0.5 : 25 Jul 2012
!   subversão  2.0.6 : 30 Out 2012

!              2.0.6.2 : 19 Nov 2012

! última alteração     : 19 Nov 2012

! Precisão dupla
! Linguagem Fortran
! Aplicativo usado: Compaq Fortran 6.6
! Tipo de projeto: Console Application

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

! *** DEFINIÇÃO DAS VARIÁVEIS DO PROGRAMA ***

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

implicit none

!integer, parameter :: nmax = 10000      ! número máximo de pontos para cada vetor

real*8,  parameter :: go   =  9.80665d0 ! aceleração gravitacional ao nível do mar (m/s2)

! Vetores com dimensão NMAX
real*8, allocatable, dimension(:) :: t    ! instante de tempo da curva processada em cada F (s)
real*8, allocatable, dimension(:) :: F    ! valor do empuxo da curva processada em cada t (N)
real*8, allocatable, dimension(:) :: Fp   ! F de projeto (N)
real*8, allocatable, dimension(:) :: Ii   ! impulso total até t (Ns)
real*8, allocatable, dimension(:) :: Mi   ! massa total de propelente consumida até t (g)
real*8, allocatable, dimension(:) :: fmi  ! fluxo de massa em cada t (g/s)
real*8, allocatable, dimension(:) :: p    ! pressão de combustão em cada t (Pa)
real*8, allocatable, dimension(:) :: Cf   ! coeficiente de empuxo real (.)
real*8, allocatable, dimension(:) :: Ca   ! velocidade característica (m/s)
real*8, allocatable, dimension(:) :: Rg   ! constante de gás perfeito dos gases do propelente (J/kg.K)
real*8, allocatable, dimension(:) :: ga   ! razão entre os calores específicos dos gases do propelente (.)
real*8, allocatable, dimension(:) :: Tg   ! temperatura dos gases do propelente na câmara (K)
real*8, allocatable, dimension(:) :: pe   ! pressão na saída da tubeira (Pa)
real*8, allocatable, dimension(:) :: Ve   ! velocidade de exaustão média real (m/s)
real*8, allocatable, dimension(:) :: Vep  ! velocidade de exaustão média de projeto (m/s)
real*8, allocatable, dimension(:) :: ci   ! velocidade de exaustão efetiva média (m/s)
real*8, allocatable, dimension(:) :: Me   ! número de Mach médio real na saída da tubeira (.)
real*8, allocatable, dimension(:) :: Mep  ! número de Mach médio de projeto na saída da tubeira (.)
real*8, allocatable, dimension(:) :: Te   ! temperatura média dos gases real na saída da tubeira (K)
real*8, allocatable, dimension(:) :: Tep  ! temperatura média dos gases de projeto na saída da tubeira (K)
real*8, allocatable, dimension(:) :: RAep ! razão entre a área de projeto na saída da tubeira e a área da garganta (.)
real*8, allocatable, dimension(:) :: Dep  ! diâmetro de projeto na saída da tubeira (m)
real*8, allocatable, dimension(:) :: at   ! velocidade do som na garganta (m/s)
integer,allocatable, dimension(:) :: som  ! =0, saída subsônica;     =1, saída supersônica

! Vetores com dimensão NTAB
real*8, allocatable, dimension(:) :: ptab  ! pressão de combustão da tabela de propriedades (Pa)
real*8, allocatable, dimension(:) :: Tgtab ! temperatura de combustão da tabela de propriedades (K)
real*8, allocatable, dimension(:) :: gatab ! ga da tabela de propriedades (.)
real*8, allocatable, dimension(:) :: Rgtab ! Rg da tabela de propriedades (J/kg.K)

integer :: w     ! freqüência da aquisição de dados (Hz)
integer :: nc    ! número de pontos da curva processada
integer :: i     ! número de cada ponto da curva processada
integer :: j     ! número da iteração no cálculo de p, pe
integer :: k     ! número do ponto da tabela de dados de To, gama, Rg x po
integer :: nmax  ! número máximo de pontos para a curva de empuxo
integer :: itmax ! número máximo de iterações para obter po e pe supersônicos
integer :: ntab  ! número de pontos da tabela de dados de To, gama, Rg x po
integer :: ver   ! variável para ativar o comando system
integer :: flag  ! variável de controle da leitura da curva original
integer :: ctor  ! =1: usa tempo original; =0: usa tempo corrigido
integer :: cgraf ! =1: mostra gráficos;    =0: não
integer :: clis  ! =1: mostra listagens auxiliares;    =0: não

real*8  :: Fcor  ! valor de corte do empuxo da curva de empuxo
real*8  :: Fo    ! valor do empuxo da curva original 
real*8  :: ti    ! tempo original para o início da curva de empuxo (s) 
real*8  :: tf    ! tempo original para o fim da curva de empuxo (s) 
real*8  :: tor   ! tempo original de cada ponto da curva de empuxo (s) 
real*8  :: dt    ! intervalo de tempo entre duas medidas (s) 
real*8  :: tmin  ! tempo original de Fmin (s) 
real*8  :: tmax  ! tempo original de Fmax (s)
real*8  :: taux  ! tempo auxiliar para ctor 
real*8  :: M     ! massa de propelente (kg)
real*8  :: tq    ! tempo de queima (s)
real*8  :: It    ! impulso total real (Ns)
real*8  :: Itp   ! impulso total de projeto (Ns)
real*8  :: c     ! velocidade de exaustão efetiva média (m/s)
real*8  :: Is    ! impulso específico médio (s)
real*8  :: fm    ! fluxo de massa médio (g/s)
real*8  :: Ag    ! área da garganta da tubeira (m2)
real*8  :: Dg    ! diâmetro da garganta da tubeira (m)
real*8  :: Ae    ! área da saída da tubeira (m2)
real*8  :: De    ! diâmetro da saída da tubeira (m)
real*8  :: pa    ! pressão atmosférica local (Pa)
real*8  :: aux   ! variável auxiliar
real*8  :: pi    ! constante matemática
real*8  :: tolp  ! tolerância para obter po e pe supersônicos (Pa)
real*8  :: tolM  ! tolerância para obter Me (.)

real*8  :: Fmin,    Fmed,    Fmax    ! valores mínimo, médio e máximo de F (N)

character*4  :: classe   ! nome da classe do motor
character*50 :: original ! nome do arquivo de dados originais de empuxo versus tempo
character*50 :: tabela   ! nome do arquivo de dados de To, gama, Rg x po
character*90 :: titulo   ! título de identificação da análise
character*90 :: dados    ! nome do arquivo de dados gerais
character*90 :: caso     ! rótulo (prefixo) do nome dos arquivos de saída

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

pi = dacos(-1.0d0)

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

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

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

! Lê os dados numéricos

open(6,file=trim(adjustl(dados)))

read(6,*) original
read(6,*) ti
read(6,*) tf
read(6,*) Fcor
read(6,*) M
read(6,*) w
read(6,*) nmax
read(6,*) ctor
read(6,*) cgraf
read(6,*) clis
read(6,*) caso
read(6,*) titulo

close(6)

! abertura do arquivo H para registro do avanço da execução do programa

open(80,file=trim(adjustl(caso))//"_saida_H.txt")

write(80,11) trim(adjustl(caso)), trim(adjustl(titulo))

11 format('#', /, &
          '#', t4, 'Saída do programa CURVA EMPUXO 2.0.6.2, versão de 19 Nov 2012', /, &
          '#', /, &
          '#', t4, 'Caso   = ', a, /, &
          '#', /, &
          '#', t4, 'Título = ', a, /, &
          '#' )

write(80,*) 'Leu os dados principais'

write(80,*) 'Vai alocar memória para os vetores'

allocate ( t(nmax),  F(nmax),   Fp(nmax),   Ii(nmax),  Mi(nmax), fmi(nmax), &
           p(nmax),  Cf(nmax),  Ca(nmax),   Rg(nmax),  ga(nmax), Tg(nmax),  &
           pe(nmax), Ve(nmax),  Vep(nmax),  ci(nmax),  Me(nmax), Mep(nmax), &
           Te(nmax), Tep(nmax), RAep(nmax), Dep(nmax), at(nmax), som(nmax)  )

write(80,*) 'vai começar o cálculo de tq, fm, dt, Ag e Ae'

tq = tf - ti
fm = 1000 * M / tq
dt = 1.0d0 / w

! *** geração da curva de empuxo corrigida ***

nc   = 0
flag = 0
Fmin = nmax
Fmax = 0.0d0

write(80,*) 'Vai iniciar a leitura de t e F do arquivo original da célula de carga'
write(80,*) ' ponto    t (s)      F (N)'

open(8,file=trim(adjustl(original)))

do while ( flag == 0 )

   read(8,*) tor, Fo

   if ( tor > ti-dt/2 .and. tor < tf+dt/2 ) then
      nc = nc + 1
      t(nc) = (nc - 1) * dt ! = tor - ti
      F(nc)  = Fo
      if ( Fo < Fcor ) F(nc) = Fcor
         write(80,501) nc, t(nc), F(nc)
         501 format (i5, 3x, f12.7, 1pe15.7 )
      if ( F(nc) < Fmin ) then
         Fmin = F(nc)
         tmin = tor
      end if
      if ( F(nc) > Fmax ) then
         Fmax = F(nc)
         tmax = tor
      end if
   end if

   if ( tor > tf+dt/2 ) exit

end do

close(8)

write(80,*) 'Terminou a leitura de t e F do arquivo original da célula de carga'

write(80,502) fm, Fmin, tmin-ti, tmin, Fmax, tmax-ti, tmax, nmax, nc
502 format("fm  = fluxo de massa médio = M/tq (g/s) = ", 1pe15.6, &
         /,"Fmin = empuxo mínimo (N) = ", 1pe15.6, &
         /,"tmin-ti = tempo corrigido do empuxo mínimo (s) = ", 1pe15.7, &
         /,"tmin    = tempo original do empuxo mínimo (s)  = ", 1pe15.7, &
         /,"Fmax = empuxo máximo (N) = ", 1pe15.6, &
         /,"tmax-ti = tempo corrigido do empuxo máximo (s) = ", 1pe15.7, &
         /,"tmax    = tempo original do empuxo máximo (s)  = ", 1pe15.7, &
         /,"nmax = número máximo de pontos para a curva de empuxo = ", i8, &
         /,"nc = número de pontos da curva processada = ", i6 )

! cálculo de Ii e It

write(80,*) 'vai começar o cálculo de Ii e It'

Ii(1) = 0.0d0
do i = 2, nc
   Ii(i) = Ii(i-1) + dt * ( F(i) + F(i-1) ) / 2
end do
It = Ii(nc)

! cálculo de Mi (estimativa inicial) e Fmed

write(80,*) 'vai começar o cálculo de Mi (estimativa inicial) e Fmed'

Mi = 1000 * M * Ii / It

!tq = tf - ti

Fmed  = It / tq

! cálculo de fmi (estimativa inicial)

write(80,*) 'vai começar o cálculo de fmi (estimativa inicial)'

fmi(1) = 0.0d0
do i = 2, nc
   fmi(i) = ( Mi(i) - Mi(i-1) ) / dt
end do

! cálculo de ci (estimativa inicial)

write(80,*) 'vai começar o cálculo de ci (estimativa inicial)'

ci(1) = 0.0d0
do i = 2, nc
   ci(i) =  1000 * ( Ii(i) - Ii(i-1) ) / ( Mi(i) - Mi(i-1) )
end do

! cálculo de c e Is

write(80,*) 'vai começar o cálculo de c e Is'

c = It / M

Is = c / go

! cálculo da classe

write(80,*) 'vai começar a determinação da classe do motor'

if (     0.0000d0 < It .and. It <=     0.3125d0 ) classe = '1/8A'
if (     0.3125d0 < It .and. It <=     0.6250d0 ) classe = '1/4A'
if (     0.6250d0 < It .and. It <=     1.2500d0 ) classe = '1/2A'
if (     1.2500d0 < It .and. It <=     2.5000d0 ) classe = 'A'
if (     2.5000d0 < It .and. It <=     5.0000d0 ) classe = 'B'
if (     5.0000d0 < It .and. It <=    10.0000d0 ) classe = 'C'
if (    10.0000d0 < It .and. It <=    20.0000d0 ) classe = 'D'
if (    20.0000d0 < It .and. It <=    40.0000d0 ) classe = 'E'
if (    40.0000d0 < It .and. It <=    80.0000d0 ) classe = 'F'
if (    80.0000d0 < It .and. It <=   160.0000d0 ) classe = 'G'
if (   160.0000d0 < It .and. It <=   320.0000d0 ) classe = 'H'
if (   320.0000d0 < It .and. It <=   640.0000d0 ) classe = 'I'
if (   640.0000d0 < It .and. It <=  1280.0000d0 ) classe = 'J'
if (  1280.0000d0 < It .and. It <=  2560.0000d0 ) classe = 'K'
if (  2560.0000d0 < It .and. It <=  5120.0000d0 ) classe = 'L'
if (  5120.0000d0 < It .and. It <= 10240.0000d0 ) classe = 'M'
if ( 10240.0000d0 < It .and. It <= 20480.0000d0 ) classe = 'N'


! *** ESCRITA DOS RESULTADOS EM ARQUIVOS ***

! ----- saída A -----

write(80,*) 'vai começar a escrita da saída A'

open(10,file=trim(adjustl(caso))//"_saida_A.txt")

write(10,11) trim(adjustl(caso)), trim(adjustl(titulo))

write(10,12) trim(adjustl(original)), ti, tf, Fcor, M*1000, w, dt, go, nmax, nc
12 format(//, " ---------- DADOS ----------", /, &
          //, "   original = nome do arquivo de dados originais de empuxo versus tempo = ", a, &
          //, "   ti = tempo original inicial (s) = ", 1pe15.7, &
          //, "   tf = tempo original final (s)   = ", 1pe15.7, &
          //, "   Fcor = valor de corte do empuxo da curva de empuxo (N) = ", 1pe15.6, &
          //, "   M = massa de propelente (g) = ", 1pe15.6, &
          //, "   w = freqüência da aquisição de dados (Hz) = ", i6, &
          //, "   dt = intervalo de tempo entre duas medidas (s) = ", 1pe15.7, &
          //, "   go = aceleração gravitacional ao nível do mar (m/s2) = ", 1pe15.6, &
          //, "   nmax = número máximo de pontos para a curva de empuxo = ", i8, &
           /, "   nc = número de pontos da curva processada = ", i6 )

write(10,13) It, Fmed, tq, trim(adjustl(classe)), idnint(Fmed)
13 format(//, " ---------- RESULTADOS PRINCIPAIS ----------", /, &
         //,  "   It  = impulso total real (Ns) = ", 1pe15.6, &
         //,  "   Fmed   = empuxo médio (N) = ", 1pe15.6, &
         //,  "   tq     = tempo de queima (s) = ", 1pe15.7, &
         //,  "   classe = classificação do motor =   ", a, i2 )

write(10,14) fm, c, Is
14 format(//, " ---------- OUTROS RESULTADOS ----------", /, &
         //,  "   fm  = fluxo de massa médio = M/tq (g/s) = ", 1pe15.6, &
         //,  "   c   = velocidade de exaustão efetiva média (m/s) = ", 1pe15.6, &
         //,  "   Is  = impulso específico médio (s) = ", 1pe15.6 )

write(10,15) Fmed,        Fmax,        Fmin,        Fmax/Fmed,       Fmin/Fmed
        
15 format(//,  1x, '--- variável ---',  3x, '--- média -----',      &
               3x, '--- máximo ----',   3x, '--- mínimo ----',      &
               3x, '-- max / med --',   3x, '-- min / med --',  //, &
               t2, 'F real [N] .....', t18,  5 ( 3x, 1pe15.6 ) )

close(10)

! ----- saída B -----

write(80,*) 'vai começar a escrita da saída B'

open(20,file=trim(adjustl(caso))//"_saida_B.txt")

write(20,11) trim(adjustl(caso)), trim(adjustl(titulo))

write(20,22)
22 format('#', t4, 't [s]', t16,'F [N]', t31,'I [Ns]', t46,'fm [g/s]', t61,'M [g]' )

if ( ctor == 1 ) t = t + ti
do i = 1, nc
   write(20,23) t(i), F(i), Ii(i), fmi(i), Mi(i)
   23 format ( f12.7, 8 ( 1pe15.6 ) )
end do

close(20)

! --------  GRÁFICOS  ---------------------------------------------------------------------

if ( cgraf == 1 ) then

write(80,*) 'vai começar a gerar o gráfico de F-t.gnu'

! *** gera o gráfico da curva de empuxo ***

open(18,file="curva_empuxo_2p0_curva_F-t.gnu")

do i = 1, 6
   read(18,*)
end do 

write(18,211) trim(adjustl(titulo))
211 format("set title '", a, "'")

!write(18,212) trim(adjustl(caso))//"_saida_C.txt"
write(18,212) trim(adjustl(caso))//"_saida_B.txt"
212 format("plot '", a, "' using 1:2 notitle")
!212 format("plot '", a, "' using 1:4 title 'F real'")

close(18)

ver = system("Wgnuplot curva_empuxo_2p0_curva_F-t.gnu")

end if

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

! *** MOSTRA O CONTEÚDO DOS ARQUIVOS DE SAÍDA ***

write(80,*) 'vai começar a abrir os arquivos de saída'

close(80)

!ver = system("Notepad " // trim(adjustl(caso)) // "_saida_H.txt" )

ver = system("Notepad " // trim(adjustl(caso)) // "_saida_A.txt" )

if ( clis == 1 ) then

ver = system("Notepad " // trim(adjustl(caso)) // "_saida_B.txt" )

end if

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

end program curva_empuxo_2p0
