
Program curva_empuxo

! 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  3.0   :  30 Jun 2015

! versão  3.1   :  11 Nov 2015

! versão  3.2   :  24 Abr 2018
! Atualização e melhoria do método de obtenção automática da curva de empuxo

! Precisão dupla
! Linguagem Fortran
! Aplicativo usado: Intel Fortran 11.1
! 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(:) :: Ii   ! impulso total até t (Ns)
real*8, allocatable, dimension(:) :: Mi   ! massa total de propelente consumida até t (g)

real*8, allocatable, dimension(:) :: fluxo_massa ! Fluxo de massa (kg/s)
real*8, allocatable, dimension(:) :: velocidade_ejecao ! velocidade de ejecao (m/s)
real*8, allocatable, dimension(:) :: Fp   ! Forca devido a pressao (N)
real*8, allocatable, dimension(:) :: FQML ! Forca devido a QML (N)
real*8, allocatable, dimension(:) :: p0 ! Pressao estagnacao (Pa)

real*8, allocatable, dimension(:) :: rq ! Taxa de queima

real*8, allocatable, dimension(:) :: Iinst


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 :: 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  :: M     ! massa de propelente (kg)
real*8  :: tq    ! tempo de queima (s)
real*8  :: It, It0    ! 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  :: At    ! área da garganta da tubeira (m2)
real*8  :: Dg    ! diâmetro da garganta da tubeira (m)
real*8  :: De    ! diâmetro da saída da tubeira (m)
real*8  :: pa    ! pressão atmosférica local (Pa)
real*8  :: pi    ! constante matemática

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

real*8  :: vdt, mindt=0.4d0    ! tempo mínimo da curva de empuxo

real*8  :: ve_med, mdot_med, aux1
integer  :: nlinhas    ! número total de linhas da curva de empuxo original



character*128  :: 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


real*8 :: pin, p0c,hi, h, y, gam, T0, R, Mcri, pg, x2, y2, x1, y1, ve, tol, qui

real*8 :: Ee, Ds

real*8 :: Mach, fu, dfu

real*8 :: relax, itmax

real*8 :: Fant, Fpost
integer :: iant, ipost

integer :: expansao, propelente

integer :: iinicial, ifinal

integer :: mod

real*8 :: tol_rel, tol_tq, tol_f

real*8 :: Lp, fa, fn, Di, roP, Llivre, racum, dP0dt, ro0, Ai, Ab, v0

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

open(5,file="./entrada/param.txt")
read(5,*) tol
read(5,*) relax
read(5,*) itmax
read(5,*) hi
close(5)

!------------------
pi = dacos(-1.0d0)

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

! Lê o nome do arquivo com os dados numéricos
open(5,file="./entrada/curva_empuxo_dados_geral.txt")
read(5,*) dados
read(5,*) mod
read(5,*) tol_rel
read(5,*) tol_tq
read(5,*) tol_f
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="./entrada/" // 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,*) Dg
read(6,*) Pa
read(6,*) expansao
read(6,*) De
read(6,*) propelente
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='./resultados/' // trim(adjustl(caso))//"_saida_H.txt")

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

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

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


! Pega o número de linhas do arquivo de texto da curva de empuxo original
open(8,file='./entrada/' //trim(adjustl(original)))
    do i=1, nmax
        read(8,*,iostat=flag)
        if(flag /= 0) exit
    enddo
close(8)
nlinhas = i-1
! A variavel nlinhas possui agora o número total de linhas da curva de empuxo original

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

allocate ( t(nlinhas),  F(nlinhas),   Fp(nlinhas),   Ii(nlinhas),  Mi(nlinhas), &
           fluxo_massa (nlinhas), velocidade_ejecao (nlinhas)  , p0 (nlinhas), FQML(nlinhas), rq(nlinhas)  )

t = 0
f= 0
fp = 0
Ii = 0
Mi = 0
p0=0
fluxo_massa = 0
velocidade_ejecao = 0
FQML = 0

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

tq = tf - ti
dt = 1.0d0 / w


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

flag = 0
Fmin = 1.0d6
Fmax = 0.0d0
nc = 0
vdt = 0

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="./entrada/" // trim(adjustl(original)))

do i=1, nlinhas

   read(8,*) tor, Fo
   
    !if(Fcor > 0.d0) then
    !    if(Fo >= Fcor) then
    !        ! Inicio da curva de empuxo, denotado por flag=0
    !        if(flag == 0) then
    !            nc = 0
    !            ti = tor
    !        endif
    !        nc = nc+1
    !        t(nc) = tor-ti
    !        F(nc) = Fo
    !        flag = 1
    !        vdt = t(nc)
    !        if ( F(nc) < Fmin ) then
    !            Fmin = F(nc)
    !            tmin = tor
    !        end if
    !        if ( F(nc) > Fmax ) then
    !            Fmax = F(nc)
    !            tmax = tor
    !        end if
    !        
    !        write(80,501) nc, t(nc), F(nc)
            501 format (i5, 3x, f12.7, 1pe15.7 )
    !
    !    else
    !        if(vdt > mindt) then
    !            tf = tor
    !            exit
    !        else
    !            flag = 0
    !        endif
    !    endif
    !else
    !    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(t(nc) <= ti) iinicial = i
          
          if(t(nc) <= tf) ifinal   = i
          
          if ( Fo < 0.d0 ) F(nc) =0d0
    !      if ( F(nc) < Fmin ) then
    !         Fmin = F(nc)
    !         tmin = tor
    !      end if
    !      if ( F(nc) > Fmax ) then
    !         Fmax = F(nc)
    !         tmax = tor
    !      end if
          
          write(80,501) nc, t(nc), F(nc)
        !
        !end if
        !
        !if ( tor > tf+dt/2 ) exit
    !endif
end do
close(8)

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

! Definicao automatica da curva de empuxo
if (mod == 1) then
    It = 0    
    ! Define inicio e fim da curva (chute inicial com base na metade de força máxima)
    
    Fcor = maxval(F)/2
    
    write(80,*) "Começará a definicao automatica da curva de empuxo"
    write(*,*) "Comecara a definicao automatica da curva de empuxo"
    write(*,'(T2,A,T5,A,T15,A,T25,A,T35,A,T45,A)') "j", "It", "Dif(%)", "Fcor" , "tq", "Fmed"
    write(80,'(T2,A,T5,A,T15,A,T25,A,T35,A,T45,A)') "j", "It", "Dif(%)", "Fcor" , "tq", "Fmed"
    do j=1, 10000
        It0 = It
        flag = 1
        do i=1, nc
            if(F(i) < Fcor ) then
                if(flag == 1) iinicial = i
                if(flag == 0) then
					ifinal = i 
                    if(t(ifinal)-t(iinicial) > tol_tq) exit
                    flag = 1
                endif                
            elseif(F(i) >= Fcor) then
                flag = 0                               
            endif
        enddo
        ti = t(iinicial)
        tf = t(ifinal)
		!write(*,*) ti, tf; pause

    
        ! Calcula o impulso total correspondente
        It = 0
        do i = iinicial+1, ifinal
           It = It + dt * ( F(i) + F(i-1) ) / 2
        end do
    
        write(*,'(I2, 20(1PE10.2))') j, It, abs(It-It0)/It, Fcor, tf-ti, It/(tf-ti)
        write(80,'(I2, 20(1PE10.2))') j, It, abs(It-It0)/It, Fcor, tf-ti, It/(tf-ti)
        if(abs(It-It0)/It <= tol_rel .or. Fcor <= tol_F) exit
        Fcor = Fcor/2
        !pause
    enddo
    
endif
!
!! Se houver dentro da curva de empuxo pontos onde a força de empuxo é zero (ou era negativa e entrou agora como zero, no caso de explosão, etc..) Faz interpolação linear entre os pontos zerados
!do j=1, 10000
!    flag = 1
!    ! Encontra os pontos onde a função força de empuxo possui buracos
!    do i=iinicial+1, ifinal
!        if(F(i) == 0 .and. F(i-1) /= 0 .and. flag == 1) then
!            flag = 0
!            Fant = F(i-1)
!            iant = i-1
!        endif
!    
!        if(flag == 0 .and. F(i) /=0 .and. F(i-1) == 0) then
!    
!            Fpost = F(i)
!            ipost = i
!            exit        
!        endif    
!    enddo
!    if(flag == 1) exit    
!    ! Interpola linearmente entre os buracos
!    do i=iant, ipost
!        F(i) = Fant + (i-iant*1d0)/(ipost-iant*1d0)*(Fpost-Fant)
!    enddo
!    
!enddo

Fmax = maxval( F(iinicial:ifinal))

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'

allocate(Iinst(nc))
Iinst = 0

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



open(10,file='./resultados/' // trim(adjustl(caso))//"_saida_C.txt")
do i = iinicial, ifinal
   write(10,*) t(i)-ti, Iinst(i)
end do
close(10)

!! GRAFICO Do Iinst
!open(18, file="./gnu/curva_empuxo_curva_F-t9.gnu")
!
!write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
!write(18,223) "set out './resultados/",trim(adjustl(original)),"_Iinst.png'"
!write(18,'(A)') "set style data lines"
!write(18,'(A)') "set key left box top"
!write(18,'(A)') "set xlabel 'tempo (s)'"
!write(18,'(A)') "set ylabel 'Iinst (N.s)'"
!write(18,'(A)') "set time"
!write(18,'(A)') "set grid"
!write(18,211) trim(adjustl(titulo))
!
!write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_C.txt", "' using 1:2 title ''"
!close(18)
!
!ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t9.gnu")



! cálculo de c e Is

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

c = It / M

Is = c / go

tq = tf - ti

! Comeca o calculo da curva de velocidade de ejecao efetiva
! calculo da curva de fluxo de massa

write(80,*) 'vai começar o cálculo da curva de velocidade de ejecao efetiva, fluxo de massa e pressão de estagnação'
write(80,*) "k, j (passos até convergir a dada força e ponto de força analisado"

Fmed  = It / tq

At = pi*0.25d0*Dg*Dg
! Estimativa da pressão de estagnação:

p0c= Pa+100    
!---------------------------
do k=iinicial, ifinal
    h = hi
    
    y = F(k)
    
    ! Calculo das variáveis termodinâmicas
    call Fk(propelente,p0c, gam)
    call FT0(propelente,p0c, T0)
    call FR(propelente,p0c, R)

    call M_cri(p0c, Pa, gam, Mcri)

    if(Mcri <= 1) then
        
        call fm_M(propelente,Mcri, At, p0c, gam, R, T0, fm)
        call ve_M(propelente,p0c, Mcri,gam, R, T0, ve)
        
        pg = pa

        x2 = P0c
        y2 = fm*ve
    elseif(Mcri > 1) then
        call fm_M(propelente,1.d0, At, p0c, gam, R, T0, fm)
        call ve_M(propelente,p0c, 1.d0,gam, R, T0, ve)
        pg = P0c*(1.d0+(gam-1.d0)/2.d0)**(gam/(1.d0-gam))
        
        x2 = P0c
        y2 = (fm*ve+(pg-pa)*At)
    endif
        
    x1 = p0c    
    do j=1, itmax
        p0c = x1+h
    ! Calculo das variáveis termodinâmicas
        call Fk(propelente,p0c, gam)
        call FT0(propelente,p0c, T0)
        call FR(propelente,p0c, R)
            
    ! Analise do comportamento do escoamento no interior da tubeira
    ! M=?
    !------------------------------------------
        call M_cri(p0c, Pa, gam, Mcri)
        if(Mcri <= 1) then
            call fm_M(propelente,Mcri, At, p0c, gam, R, T0, fm)
            call ve_M(propelente,p0c, Mcri,gam, R, T0, ve)
            
            pg = pa

!            write(4,'(I3,A,9(1PE25.15))') j, ' 1< ',(fm*ve)/y-1, Mcri
            x1 = P0c
            y1 = fm*ve
        elseif(Mcri > 1) then
            call fm_M(propelente,1.d0, At, p0c, gam, R, T0, fm)
            if(expansao == 0) then
                call ve_M(propelente,p0c,1.d0,gam, R, T0, ve)
                pg = P0c*(1.d0+(gam-1.d0)/2.d0)**(gam/(1.d0-gam))
            elseif(expansao == 1) then
                call ve_M(propelente,p0c,Mcri,gam, R, T0, ve)
                pg = pa   
            elseif(expansao == 2) then
                Ee =((1.d0+0.5d0*(gam-1.d0)*Mcri**2)/(0.5d0*(gam+1.d0)))**((gam+1.d0)/(2.d0*(gam-1.d0)))/Mcri
                Ds = dsqrt(Ee*Dg**2)
                if(Ds <= De) then
                    call ve_M(propelente,p0c,Mcri,gam, R, T0, ve)
                    pg = pa 
                elseif(Ds > De) then
                    Ee = (De/Dg)**2
                    ! Chuta um numero de Mach relativo a nova razao de expansao
                    Mach = Mcri*Ee/((Ds/Dg)**2)
                        do i=1, 10000
                            fu = 1.d0/Mach*(2.d0/(gam+1.d0)*(1.d0+1.d0/2.d0*(gam-1.d0)*Mach**2))**(1.d0/2.d0*(gam+1.d0)/(gam-1.d0))-Ee
    
                            dfu = -1.d0/Mach**2*(2.d0/(gam+1.d0)*(1.d0+1.d0/2.d0*(gam-1.d0)*Mach**2))**(1.d0/2.d0*(gam+1.d0)/(gam-1.d0))+ &
                            1.d0/2.d0*(2.d0/(gam+1.d0)*(1.d0+1.d0/2.d0*(gam-1.d0)*Mach**2))**(1.d0/2.d0*(gam+1.d0)/(gam-1.d0))*(gam+1.d0)/(1.d0+1.d0/2.d0*(gam-1.d0)*Mach**2)
                            
                            Mach = Mach -fu/dfu
                            
                            if(abs(fu/dfu) < tol) exit
                        enddo
                    call ve_M(propelente,p0c,Mach,gam, R, T0, ve)
                    pg = P0c*(1.d0+(gam-1.d0)/2.d0*Mach**2)**(gam/(1.d0-gam))
                endif
                
            endif
            

            
!            write(4,'(I3,A,9(1PE25.15))') j, ' >1 ', (fm*ve+(pg-pa)*At)/y-1, Mcri
            x1 = P0c
            y1 = (fm*ve+(pg-pa)*At)
        endif
        
        ! Novo valor da pressao
        h = (x2-x1)/(y2-y1)*(y-y1)*relax
        !write(*,'(I2,10(1PE15.7))') j, h , p0c, y2-y1, h/p0c

        If( abs(h/p0c)< tol) exit
        x2 = x1
        y2 = y1
    enddo
    
    if(j == itmax + 1) then
         write(*,'(/A,/A)') "ERRO 01- Metodo de calculo da pressao de estagnacao nao convergiu!!", &
                            "Entre em contato com o suporte"
         pause " "
        stop
    endif
    write(80,*) k, j
   
    p0(k) = P0c
    Fluxo_massa(k) = fm
    velocidade_ejecao(k) = ( ve+(pg-pa)*At/fm )
    Fp(k) = (pg-pa)*At
    FQML(k) = fm*ve
    
enddo


! Conferência
Ve_med = 0.0d0
mdot_med = 0.0d0
aux1 = 0

do i = iinicial+1, ifinal
    ! AUX1: Impulso total
    aux1 = aux1 + 0.5d0* dt * ( velocidade_ejecao(i)*Fluxo_massa(i) + velocidade_ejecao(i-1)*Fluxo_massa(i-1) )
    Ve_med = Ve_med + 0.5d0* dt * ( velocidade_ejecao(i) + velocidade_ejecao(i-1) )
    mdot_med = mdot_med + 0.5d0*dt * ( Fluxo_massa(i) + Fluxo_massa(i-1) )
enddo
Ve_med = Ve_med / tq
mdot_med = mdot_med / tq



write(80,'(A,10(1PE13.5))') "Valores numericos", mdot_med, Ve_med, aux1
write(80,'(A,10(1PE13.5))') "Valores experimentais", M/tq, It/M, It
write(80,'(A,10(1PE13.5))') "Erros Percentuais", (M/tq-mdot_med)*100/(M/tq), (It/M-Ve_med)*100/(It/M), (It-aux1)*100/It


! cálculo da classe

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

call class (It, Fmed, classe)


! *** ESCRITA DOS RESULTADOS EM ARQUIVOS ***
fm = M / tq
! ----- saída A -----

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

open(10,file='./resultados/' // 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 eMuxo 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(*,'(F7.1)') Fmed;stop

write(10,13) It, Fmed, tq, trim(adjustl(classe))
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 )

write(10,14) fm*1e3, 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 ) )
               
        
100 format(//, " ---------- RESULTADOS ADMITINDO EXPANSÃO ",A," ----------", /, &
         //,  "   fm  = fluxo de massa médio = M/tq (g/s) = ", 1pe15.6, &
         //,  "   c   = velocidade de exaustão efetiva média (m/s) = ", 1pe15.6, &
         //,  "   It  = impulso total (N.s) = ", 1pe15.6, &
         //,  "   Efm = erro do fluxo de massa (%) = ", 1pe15.6, &
         //,  "   Ec  = erro da velocidade de exaustão efetiva média (%) = ", 1pe15.6, &
         //,  "   EIt = erro do impulso total (%) = ", 1pe15.6 , &
         //,  "   pmed = pressão média (Bar) = ", 1pe15.6, &
         //,  "   pmax = pressão máxima (Bar) = ", 1pe15.6)

if(expansao == 0) then
    write(10,100) 'ATÉ A GARGANTA',MDOT_MED*1e3, Ve_med, aux1, (M/tq-mdot_med)*100/(M/tq), &
    (It/M-Ve_med)*100/(It/M), (It-aux1)*100/It, sum(p0(iinicial:ifinal))/(ifinal-iinicial)/1e5, maxval(p0(iinicial:ifinal))/1e5
elseif(expansao == 1) then
    write(10,100) 'ÓTIMA',MDOT_MED*1e3, Ve_med, aux1, (M/tq-mdot_med)*100/(M/tq), &
    (It/M-Ve_med)*100/(It/M), (It-aux1)*100/It, sum(p0(iinicial:ifinal))/(ifinal-iinicial)/1e5, maxval(p0(iinicial:ifinal))/1e5       
elseif(expansao == 2) then
    write(10,100) 'ATÉ UM DIÂMETRO FIXO',MDOT_MED*1e3, Ve_med, aux1, (M/tq-mdot_med)*100/(M/tq), &
    (It/M-Ve_med)*100/(It/M), (It-aux1)*100/It, sum(p0(iinicial:ifinal))/(ifinal-iinicial)/1e5, maxval(p0(iinicial:ifinal))/1e5  
endif
    
call FQui(propelente,sum(p0(1:nc))/nc,Qui)
write(10,'(2/,75("-"),2/,1PE15.6, A)') MDOT_MED*1e3*tq, " = Massa de gases embarcada estimada (g)"
write(10,'(1PE15.6, A)') MDOT_MED*1e3*tq*qui/(1-qui), " = Massa de condensado embarcada estimada (g)"
write(10,'(1PE15.6, A)') MDOT_MED*1e3*tq/(1-qui), " = Massa de propelente embarcada estimada (g)"    

!write(10,'(1PE15.6, A)') racum, " = Comprimento estimado do grão-propelente (m)"
!write(10,'(1PE15.6, A)') Lp, " = Comprimento real do grão-propelente (m)"
close(10)

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

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

open(20,file='./resultados/' // 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,'fm[g/s]', t46,'c_inst[m/s]', t61,'P0[Bar]', t76,'FP[N]', t91,'FQML[N]', &
   t106,'V_QML[m/s]', t121,'V_Pressão[m/s]', t136,'Cf', t151,'Ee', t166,'Ds', t181,'Iinst', t196, 'Gamma', t211, 'T0', T226, 'R')

if ( ctor == 1 ) t = t + ti
do i = iinicial+1, ifinal
    ! Calculo das variáveis termodinâmicas
    call Fk(propelente,p0(i), gam)
    call FT0(propelente,p0(i), T0)
    call FR(propelente,p0(i), R)

    ! Calculo do Mach na saida
    call M_cri(p0(i), Pa, gam, Mcri)
    if(Mcri > 1.d0) then
        Ee =((1.d0+0.5d0*(gam-1.d0)*Mcri**2)/(0.5d0*(gam+1.d0)))**((gam+1.d0)/(2.d0*(gam-1.d0)))/Mcri
        Ds = dsqrt(Ee*Dg**2)
    else
        Ee = 1.d0
        Ds = dsqrt(Ee*Dg**2)
    endif
    
   write(20,23) t(i)-ti, F(i), fluxo_massa(i)*1e3, velocidade_ejecao(i), P0(i)/1e5, Fp(i), FQML(i), FQML(i)/fluxo_massa(i), Fp(i)/fluxo_massa(i), &
   F(i)/At/P0(i), Ee, Ds*1e3, 0.5d0*dt * ( F(i) + F(i-1) ), gam, T0, R
   23 format ( f12.7, 50 ( 1pe15.6 ) )
end do
close(20)




! --------  GRÁFICOS  ---------------------------------------------------------------------
211 format("set title '",A,"' ")
if ( cgraf == 1 ) then

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

! GRAFICO da CURVA DE EMPUXO
open(18, file="./gnu/curva_empuxo_curva_F-t2.gnu")
write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
223 format(A,A<len(trim(adjustl(original)))-4>,A)
write(18,223) "set out './resultados/",trim(adjustl(original)),"_empuxo.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'empuxo (N)'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

if(expansao == 0 .or. expansao == 2) then
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:2 title 'Empuxo total' \"
    write(18,'(A,A,A)') ",'./resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:6 title 'Empuxo de pressão' \"
    write(18,'(A,A,A)') ",'./resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:7 title 'Empuxo de QML' "
elseif (expansao == 1) then
        write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:2 title 'Empuxo total'"
endif
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t2.gnu")

! GRAFICO DO FLUXO DE MASSA
open(18, file="./gnu/curva_empuxo_curva_F-t3.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_fm.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'fluxo de massa (g/s)'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:3 title 'Fluxo de Massa'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t3.gnu")

! GRAFICO DA VELOCIDADE DE EJECAO
open(18, file="./gnu/curva_empuxo_curva_F-t4.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_c.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'Velocidade de ejeção efetiva (m/s)'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

if(expansao == 0 .or. expansao == 2) then
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:4 title 'Velocidade de ejeção efetiva média'\"
    write(18,'(A,A,A)') ",'./resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:8 title 'Velocidade dos gases devido a QML' \"
    write(18,'(A,A,A)') ",'./resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:9 title 'Velocidade dos gases devido a pressão' "
elseif (expansao == 1) then
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:4 title 'Velocidade de ejeção efetiva média'"
endif

close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t4.gnu")

! GRAFICO DA PRESSAO de ESTAGNACAO
open(18, file="./gnu/curva_empuxo_curva_F-t5.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_p0.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'Pressão de estagnação (Bar)'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:5 title 'Pressão de Estagnação' \"
write(18,'(A,1pe25.15,A)') ",",Pa/1e5," title 'Pressão Ambiente'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t5.gnu")

! GRAFICO Do CF
open(18, file="./gnu/curva_empuxo_curva_F-t6.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_Cf.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'Cf'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:10 title 'Cf'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t6.gnu")

! GRAFICO Do Ee
open(18, file="./gnu/curva_empuxo_curva_F-t7.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_Ee.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"

if(expansao == 0 .or. expansao == 2) write(18,'(A)') "set ylabel 'Ee necessário para expansão ótima'"
if(expansao == 1) write(18,'(A)') "set ylabel 'Ee efetivo'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

if(expansao == 2) then
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:11 title 'Ee' \"
    write(18,'(A,1PE15.7,A)') ", ", (De/Dg)**2, " title 'Ee embarcada'"
else
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:11 title 'Ee'"
endif
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t7.gnu")

! GRAFICO Do Ds
open(18, file="./gnu/curva_empuxo_curva_F-t8.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_Ds.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
if(expansao == 0 .or. expansao == 2) write(18,'(A)') "set ylabel 'Ds (mm) necessário para expansão ótima'"
if(expansao == 1) write(18,'(A)') "set ylabel 'Ds (mm) efetivo'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))


if(expansao == 2) then
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:12 title 'Ds (mm)' \"
    write(18,'(A,1PE15.7,A)') ", ", Dg*1e3, " title 'Dg (mm)' \"
    write(18,'(A,1PE15.7,A)') ", ", De*1e3, " title 'De embarcado (mm)' \"
else
    write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:12 title 'Ds (mm)' \"
    write(18,'(A,1PE15.7,A)') ", ", Dg*1e3, " title 'Dg (mm)'"
endif

close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t8.gnu")

! GRAFICO Do GAMMA
open(18, file="./gnu/curva_empuxo_curva_F-t9.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_gamma.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'Gamma'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:14 title 'Gamma'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t9.gnu")

! GRAFICO Do T0
open(18, file="./gnu/curva_empuxo_curva_F-t10.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_T0.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'T0 (K)'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:15 title 'T0 (K)'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t10.gnu")

! GRAFICO Do R
open(18, file="./gnu/curva_empuxo_curva_F-t11.gnu")

write(18,'(A)') "set term png nocrop enhanced size 1280, 720"
write(18,223) "set out './resultados/",trim(adjustl(original)),"_R.png'"
write(18,'(A)') "set style data lines"
write(18,'(A)') "set key left box top"
write(18,'(A)') "set xlabel 'tempo (s)'"
write(18,'(A)') "set ylabel 'R ( J/(Kg.K) )'"
write(18,'(A)') "set time"
write(18,'(A)') "set grid"
write(18,211) trim(adjustl(titulo))

write(18,'(A,A,A)') "plot './resultados/", trim(adjustl(caso))//"_saida_B.txt", "' using 1:16 title 'R ( J/(Kg.K) )'"
close(18)

ver = system(".\gnuplot\GNUPLOT.exe ./gnu/curva_empuxo_curva_F-t11.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

! -----------------------------------------------------------------------------
contains
!------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------
subroutine FQui(prop,P0, Qui)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: Qui
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    


select case (prop)
    ! KNSU 65/35
    case (0)
A0 =-2.04799653852393E+01
A1 = 3.62113343656072E+01
A2 =-2.77291460981005E+01
A3 = 1.21466806224666E+01
A4 =-3.34936868681469E+00
A5 = 6.02653493751515E-01
A6 =-7.07220352937977E-02
A7 = 5.21876383009207E-03
A8 =-2.19830376383906E-04
A9 = 4.03090411144927E-06
    ! Polvora Roberto 150/50/20
    case (1)
A0 = 7.26043568864702E+01
A1 =-1.44059368611776E+02
A2 = 1.26207419659459E+02
A3 =-6.37177913512304E+01
A4 = 2.04215720064625E+01
A5 =-4.30677147952935E+00
A6 = 5.97321868085242E-01
A7 =-5.25030031896640E-02
A8 = 2.65184174254155E-03
A9 =-5.85854750762796E-05

    ! Polvora KNSU 60/40
    case (2)
A0 = 1.62460275168875E+01
A1 =-2.76930932627632E+01
A2 = 2.07639961634161E+01
A3 =-8.85836008568855E+00
A4 = 2.37098159313768E+00
A5 =-4.12691242859834E-01
A6 = 4.67339057507249E-02
A7 =-3.32430078032562E-03
A8 = 1.35010149479413E-04
A9 =-2.38958104337635E-06
 
    ! Polvora KNSU 70/30
    case (3)
A0 = 1.99837170806653E+00
A1 =-4.70322436642678E+00
A2 = 4.66891945154650E+00
A3 =-2.54946842877398E+00
A4 = 8.56376845398039E-01
A5 =-1.84348060636415E-01
A6 = 2.55205817914470E-02
A7 =-2.19476961416862E-03
A8 = 1.06509629421422E-04
A9 =-2.22494962331610E-06

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
Qui = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine Fqui
!------------------------------------------------------------------------------------------
subroutine FR(prop,P0, R)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: R
    
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    

select case (prop)
    ! KNSU 65/35
    case (0)
A0 = 3.93261975545070E+02
A1 =-2.78261444164441E+02
A2 = 2.21731823032611E+02
A3 =-1.00207260435184E+02
A4 = 2.74321375819001E+01
A5 =-4.67641633470184E+00
A6 = 4.86891326745650E-01
A7 =-2.87651336070012E-02
A8 = 8.00061878913493E-04
A9 =-5.32904201135432E-06
    ! Polvora Roberto 150/50/20
    case (1)
A0 =-4.61816631501257E+04
A1 = 9.26089352007268E+04
A2 =-8.12663522815713E+04
A3 = 4.11056870481495E+04
A4 =-1.32027764825664E+04
A5 = 2.79128507736463E+03
A6 =-3.88244760033038E+02
A7 = 3.42396307890670E+01
A8 =-1.73615214889782E+00
A9 = 3.85331201552525E-02
        
    ! Polvora KNSU 60/40
    case (2)
A0 = 6.71732634788817E+03
A1 =-1.13880372797209E+04
A2 = 8.76949555557240E+03
A3 =-3.87374827622120E+03
A4 = 1.08094863188746E+03
A5 =-1.97577175720993E+02
A6 = 2.36492984951256E+01
A7 =-1.78690783631224E+00
A8 = 7.73150298950068E-02
A9 =-1.45940036826555E-03

    ! Polvora KNSU 70/30
    case (3)
A0 = 2.47645890800131E+02
A1 = 1.16173534468038E+02
A2 =-2.08988520875701E+02
A3 = 1.46939752486607E+02
A4 =-5.73358292473525E+01
A5 = 1.36290650072157E+01
A6 =-2.02380940138111E+00
A7 = 1.83269113930935E-01
A8 =-9.24730277758449E-03
A9 = 1.99039233116305E-04

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
R = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine FR
!------------------------------------------------------------------------------------------
subroutine FT0(prop,P0, T0)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: T0
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0
    
select case (prop)
    ! KNSU 65/35
    case (0)
A0 =-1.37888992619268E+04
A1 = 2.63269770368983E+04
A2 =-2.03742329478610E+04
A3 = 9.02448771642549E+03
A4 =-2.50843498438118E+03
A5 = 4.53324034236469E+02
A6 =-5.31925781292568E+01
A7 = 3.90421818333940E+00
A8 =-1.62672493431325E-01
A9 = 2.93507026655967E-03
    ! Polvora Roberto 150/50/20
    case (1)
A0 =-1.81135947499828E+05
A1 = 3.63202971159722E+05
A2 =-3.18593433119871E+05
A3 = 1.61107084078521E+05
A4 =-5.17303063313623E+04
A5 = 1.09329676609606E+04
A6 =-1.52011083262037E+03
A7 = 1.34001521365648E+02
A8 =-6.79117562963686E+00
A9 = 1.50634004319533E-01

    ! Polvora KNSU 60/40
    case (2)
A0 =-3.93833638581385E+04
A1 = 7.15176408651588E+04
A2 =-5.54248331358298E+04
A3 = 2.46941412204016E+04
A4 =-6.96071798591121E+03
A5 = 1.28700072060828E+03
A6 =-1.55984837441572E+02
A7 = 1.19385199904895E+01
A8 =-5.23114870016792E-01
A9 = 9.99299233223467E-03
        
    ! Polvora KNSU 70/30
    case (3)
A0 =-2.47644986635254E+03
A1 = 4.22534690529101E+03
A2 =-1.60315577906195E+03
A3 =-4.80970623690143E+01
A4 = 2.41319545777364E+02
A5 =-8.83748663032046E+01
A6 = 1.61431174239937E+01
A7 =-1.65434932347325E+00
A8 = 9.06303622630090E-02
A9 =-2.06798110867363E-03

end select

logP0 = dlog10(P0)
! 0.01 bar a 11 kbar
T0 = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0

end subroutine FT0
!------------------------------------------------------------------------------------------
subroutine Fk(prop,P0, k)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: P0
real*8,intent(out) :: k
real*8 :: A0, A1, A2, A3, A4, A5, A6, A7, A8, A9

real*8 :: logP0


select case (prop)
    ! KNSU 65/35
    case (0)
! kgas
!A0=-4.15198741081744E+00
!A1= 9.43179477126738E+00
!A2=-7.17478856172397E+00
!A3= 3.12235584269792E+00
!A4=-8.57362675678530E-01
!A5= 1.54068099574775E-01
!A6=-1.81219194658271E-02
!A7= 1.34583650541077E-03
!A8=-5.72884438851859E-05
!A9= 1.06539054876278E-06

! kmix
!A0= 1.89225403506303E+00
!A1=-1.23412794184640E+00
!A2= 9.73970487442510E-01
!A3=-4.41653674918613E-01
!A4= 1.25118018390715E-01
!A5=-2.29286320619050E-02
!A6= 2.71256388130721E-03
!A7=-1.99431987464383E-04
!A8= 8.26545014672395E-06
!A9=-1.47293882113353E-07

! k2ph
A0= 3.09510840624045E+00
A1=-3.47788297095420E+00
A2= 2.68589018179646E+00
A3=-1.19233536051624E+00
A4= 3.33445002556658E-01
A5=-6.08847172457605E-02
A6= 7.25433618821334E-03
A7=-5.43620288790910E-04
A8= 2.32502840539917E-05
A9=-4.32683123306007E-07
    ! Polvora Roberto 150/50/20
    case (1)
A0= 5.78092194482495E+00
A1=-9.31175283694160E+00
A2= 8.13298643510669E+00
A3=-4.10099875062160E+00
A4= 1.31571165602474E+00
A5=-2.78400484688017E-01
A6= 3.88294967051815E-02
A7=-3.43954747961506E-03
A8= 1.75404902262076E-04
A9=-3.91867465686797E-06

    ! Polvora KNSU 60/40
    case (2)
A0=-6.11043165506052E-01
A1= 2.98058454366955E+00
A2=-2.25564722929034E+00
A3= 9.67756473233875E-01
A4=-2.60214357140101E-01
A5= 4.54419230940925E-02
A6=-5.15377674951999E-03
A7= 3.66450558281284E-04
A8=-1.48496162006620E-05
A9= 2.61844906051524E-07
        
    ! Polvora KNSU 70/30
    case (3)
A0= 2.11131140257238E+00
A1=-1.48872601067561E+00
A2= 1.02571757620582E+00
A3=-4.04537302415782E-01
A4= 9.82597170151715E-02
A5=-1.51035569743172E-02
A6= 1.44690561421264E-03
A7=-8.09074405121067E-05
A8= 2.24135022850270E-06
A9=-1.84145422651684E-08

end select

logP0 = dlog10(P0)

! 0.01 bar a 11 kbar
k = &
    A9*logP0**9 + A8*logP0**8 + A7*logP0**7 + &
    A6*logP0**6 + A5*logP0**5 + A4*logP0**4 + &
    A3*logP0**3 + A2*logP0**2 + A1*logP0    + A0



end subroutine Fk
!------------------------------------------------------------------------------------------
subroutine fm_M(prop,Mcri, At,p0,k, R, T0, fm)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: Mcri, At, p0, k, R, T0
real*8,intent(out) :: fm

real*8 :: Qui

real*8 :: KL1, KL2, KL3, KL4, rcag, CdKL


fm = Mcri*At*p0/sqrt(R*T0/k*(1+(k-1)/2*Mcri**2)**((k+1)/(k-1)))

call FQui(prop,P0, Qui)
fm = fm*(1+qui)

!if(Mcri>=1.d0) then
!    rcag = 0.d0
!
!    KL1 = ( k + 1 ) / ( ( 1 + rcag ) ** 2 )
!    KL2 = ( 8 * k - 27 ) / ( 2304 * ( 1 + rcag ) )
!    KL3 = 754 * ( k ** 2 ) - ( 757 * k ) + 3633
!    KL4 = 276480 * ( ( 1 + rcag ) ** 2 )
!
!    CdKL = 1 - KL1 * ( ( 1 / 9.6d+1 ) - KL2 + ( KL3 / KL4 ) )
!
!    fm = fm * CdKL
!endif

end subroutine fm_M
!------------------------------------------------------------------------------------------
subroutine ve_M(prop,p0,Mcri,k, R, T0, ve)
implicit none
integer, intent(in) :: prop
real*8,intent(in) :: p0, Mcri, k, R, T0
real*8,intent(out) :: ve

real*8 :: Qui

real*8 :: KL1, KL2, KL3, KL4, rcag, CdKL

ve = Mcri*sqrt(k*R*T0/(1.d0+0.5d0*(k-1.d0)*Mcri**2))

call FQui(prop,P0, Qui)
ve = ve*(1-qui)

!if(Mcri>=1.d0) then
!    rcag = 0.d0
!
!    KL1 = ( k + 1 ) / ( ( 1 + rcag ) ** 2 )
!    KL2 = ( 8 * k - 27 ) / ( 2304 * ( 1 + rcag ) )
!    KL3 = 754 * ( k ** 2 ) - ( 757 * k ) + 3633
!    KL4 = 276480 * ( ( 1 + rcag ) ** 2 )
!
!    CdKL = 1 - KL1 * ( ( 1 / 9.6d+1 ) - KL2 + ( KL3 / KL4 ) )
!
!    ve = ve * CdKL
!endif


end subroutine ve_M
!------------------------------------------------------------------------------------------
subroutine M_cri(p0, Pamb, k, Mcri)
implicit none
real*8,intent(in) :: p0, pamb, k
real*8,intent(out) :: Mcri

Mcri = sqrt(((p0/Pamb)**((k-1.d0)/k)-1.d0)*2.d0/(k-1.d0))

end subroutine M_cri
!------------------------------------------------------------------------------------------
subroutine class(Itot, Fmed, texto)

real*8, intent(in) :: Itot, Fmed

character*128, intent(out) :: texto

write(texto,'(F15.2)') Fmed

40 format(A, "  --  ", A)
if(Itot<=2.5d0/64) then
    write(texto,40) "1/64 A", trim(adjustl(texto))
elseif(Itot<=2.5d0/32) then
    write(texto,40) "1/32 A", trim(adjustl(texto))
elseif(Itot<=2.5d0/16) then
    write(texto,40) "1/16 A", trim(adjustl(texto))
elseif(Itot<=2.5d0/8) then
    write(texto,40) "1/8 A", trim(adjustl(texto))
elseif(Itot<=2.5d0/4) then
    write(texto,40) "1/4 A", trim(adjustl(texto))
elseif(Itot<=2.5d0/2) then
    write(texto,40) "1/2 A", trim(adjustl(texto))
elseif(Itot<=2.500d0) then
    write(texto,40) " A", trim(adjustl(texto))
elseif(Itot<=5.000d0) then
    write(texto,40) " B", trim(adjustl(texto))
elseif(Itot<=10.00d0) then
    write(texto,40) " C", trim(adjustl(texto))
elseif(Itot<=20.00d0) then
    write(texto,40) " D", trim(adjustl(texto))
elseif(Itot<=40.00d0) then
    write(texto,40) " E", trim(adjustl(texto))
elseif(Itot<=80.00d0) then
    write(texto,40) " F", trim(adjustl(texto))
elseif(Itot<=160.0d0) then
    write(texto,40) " G", trim(adjustl(texto))
elseif(Itot<=320.0d0) then
    write(texto,40) " H", trim(adjustl(texto))
elseif(Itot<=640.0d0) then
    write(texto,40) " I", trim(adjustl(texto))
elseif(Itot<=1280.d0) then
    write(texto,40) " J", trim(adjustl(texto))
elseif(Itot<=2560.d0) then
    write(texto,40) " K", trim(adjustl(texto))
elseif(Itot<=5120.d0) then
    write(texto,40) " L", trim(adjustl(texto))    
elseif(Itot<=10200d0) then
    write(texto,40) " M", trim(adjustl(texto))  
endif
    
    
end subroutine class
!------------------------------------------------------------------------------------------

    end program curva_empuxo

