
Program curva_empuxo_2p1

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

! versão atual 2.1: 23 Set 2014

! última alteração: 29 Set 2014

! Precisão dupla
! Linguagem Fortran
! Aplicativo usado: Microsoft Visual Studio 2008 versão 9.0.30729.1 SP
! Compilador usado: Intel Fortran 11.1.065
! Tipo de projeto: Console Application

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

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

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

implicit none

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)

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 :: nmax  ! número máximo de pontos para a curva de empuxo
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 mínimo de corte do empuxo na 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  :: 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  :: aux   ! variável auxiliar

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

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

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

! Lê o nome do arquivo com os dados numéricos
open(5,file="curva_empuxo_2p1_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)

11 format('#', /, &
          '#', t4, 'Saída do programa CURVA EMPUXO 2.1, versão de 29 Set 2014', /, &
          '#', /, &
          '#', t4, 'Caso   = ', a, /, &
          '#', /, &
          '#', t4, 'Título = ', a, /, &
          '#' )

allocate ( t(nmax),  F(nmax) )

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

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

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

Fmed  = It / tq

c = It / M

Is = c / go

! cálculo da classe

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

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 mínimo de corte do empuxo na 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)), 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, 1pe6.0 )

if ( ctor == 0 ) tmax = tmax - ti

write(10,14) fm, c, Is, Fmax, tmax
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, &
         //,  "   Fmax = empuxo máximo (N) = ", 1pe15.6, &
          /,  "   tmax = tempo de Fmax (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 -----

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

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

close(20)

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

if ( cgraf == 1 ) then

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

open(18,file="curva_empuxo_2p1_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_B.txt"
212 format("plot '", a, "' using 1:2 notitle")

close(18)

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

end if

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

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

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_2p1
