program SBE

! SBE é o acrônimo p/ Simulador de Balística de Espaçomodelos.
! O objetivo deste programa é efetuar cálculos de balística de
! espaçomodelos e prever, teoricamente, o desempenho de projetos
! do Grupo de Aerodinâmica de Veículos Auto Propulsados - GREAVE
! do curso de Engenharia Mecânica da Universidade Positivo.

! Implementação da versão inicial (1.0)
!                  versão atual: 2.0
!      Programadores: Jeremy Paraná Blavier e Antonio Carlos Foltran
!      Local:         GREAVE/Eng.Mec/UP, Curitiba, PR
!      Período:       07 a 19 de setembro de 2015
!                     10 a 21 de janeiro de 2016
!      Última alteração: 13 de abril de 2018: Inclusão da abertura de
! paraquedas e novas funções p/ calcular as propriedades da atmosfera.
! Alterada lógica do Sonda IV p/ simular tempo de ação dos estágios.

! 1) O referencial utilizado neste programa para localizar a posição
! e medir a velocidade e aceleração do espaço-modelo é localizado na
! base da plataforma de lançamento e o eixo coordenado onde a altitude
! é medida aumenta na direção da atmosfera. Por isso a altitude é
! considerada nula quando o foguete está sobre a rampa de lançamento
! independentemente da altitude da rampa em relação ao nível do mar
! ou ao centro geométrico do planeta;
! 2) Modelo para descrever a massa específica do ar pode ser a equação
! de estado ou a Atmosfera Padrão Internacional, sendo que considerou-
! se nula a massa específica em altitudes maiores que 47 km;
! 3) Modelo para a aceleração gravitacional variável (ver subrotina
! "propriedades").

use portlib	! Biblioteca básica do FORTRAN.
implicit none

! Declaração das variáveis:
character*50 :: sim_ID  ! Identificador da simulação;
character*50 :: caso_ID ! Identificador específico da simulação;
character*50 :: tab_E   ! Nome arquivo de dados de empuxo em função do tempo;
character*50 :: tab_Cd  ! Nome arquivo de dados de Cd em função do tempo;
integer :: ID_E	 ! Seletor do modelo para descrever a força de empuxo;
integer :: ID_Cd ! Seletor do modelo para descrever o coeficiente de arrasto;
integer :: ID_ro ! Seletor do modelo para descrever a massa específica do ar;
integer :: ID_g  ! Seletor do modelo para descrever a aceleração da gravidade;
integer :: ID_pq ! Seletor da variável de abertura do paraquedas;
integer :: i     ! Índice identificador de um ponto da trajetória;
integer :: j     ! Índice de um ponto qualquer da curva empuxo x tempo;
integer :: NpE   ! Número total de pontos da curva empuxo x tempo;
integer :: NpCd  ! Número total de pontos da curva Cd x V;
integer :: ite   ! Contador;
integer :: itmax ! Número máximo de iterações para processo iterativo;
integer :: ver	 ! Variável de visualização dos resultados;
integer :: w	 ! Frequência de escrita dos resultados no arquivo de saída;
integer :: apogeu ! Assume valor zero e muda p/ 1 quando o apogeu é detectado;
integer :: fim_de_queima !Assume valor zero e muda p/ 1 quando o final da queima é detectado;
integer :: abertura_paraquedas ! Assume valor 0 e muda p/ 1 qdo o paraquedas inicia a abertura;
real*8  :: Mp    ! Massa de propelente [kg];
real*8  :: Mf    ! Massa final [kg];
real*8  :: dr    ! Diâmetro de referência do foguete [m];
real*8  :: dr_pq ! Diâmetro efetivo do paraquedas [m];
real*8  :: pi    ! Constante matemática;
real*8  :: Ar    ! Área de referência [m2];
real*8  :: Ar_f  ! Área de referência do foguete[m2];
real*8  :: Ar_pq ! Área de referência do paraquedas[m2];
real*8  :: Ec    ! Empuxo constante do motor [N];
real*8  :: tp    ! Tempo de queima do propelente [s];
real*8  :: It    ! Impulso total produzido pelo propulsor [Ns];
real*8  :: Cd    ! Coeficiente de arrasto [adim.];
real*8  :: Cd_pq ! Coeficiente de arrasto do paraquedas aberto[adim.];
real*8  :: g0    ! Aceleração gravitacional local na superfície [m/s2];
real*8  ::  g    ! Aceleração gravitacional em função da altitude [m/s2];
real*8  :: ro    ! Massa específica do ar [kg/m3];
real*8  ::p_atm  ! Pressão atmosférica local [Pa];
real*8  ::T_atm  ! Temperatura atmosférica local [°C];
real*8  ::R_ar   ! Constante dos gases perfeitos do ar [J/kg.K];
real*8	:: dt    ! Intervalo de tempo utilizado na discretização [s];
real*8  :: Vms   ! Velocidade máxima na subida [m/s];
real*8  :: Vmd   ! Velocidade máxima na descida [m/s];
real*8  :: Vi    ! Velocidade do impacto [m/s];
real*8	::   t   ! Instante de tempo [s];
real*8	::   th  ! Instante de tempo em que o apogeu é atingido [s];
real*8  :: t_pq  ! Tempo após o lançamento no qual o paraquedas é aberto [s];
real*8  :: ta_pq ! Tempo que o paraquedas leva para abrir [s];
real*8  ::   ti  ! Instante de tempo em que ocorre o impacto [s];
real*8	::   h   ! Altitude do espaçomodelo [m];
real*8  ::  hmax ! Apogeu alcançado pelo espaçomodelo [m];
real*8  ::  h_pq ! Altitude após o apogeu na qual o paraquedas abre [m];
real*8	::   V   ! Velocidade do espaçomodelo [m/s];
real*8	::   a   ! Aceleração do espaçomodelo [m/s2];
real*8	::   M   ! Massa instantânea do espaçomodelo [kg];
real*8	::   P   ! Força peso do espaçomodelo [N];
real*8	::   E   ! Força de empuxo do propulsor [N];
real*8	::   D   ! Força de arrasto [N];
real*8	::  Cv   ! Cd [adim.];
real*8	::  t0   ! Instante de tempo na iteração anterior [s];
real*8	:: t00   ! Instante de tempo duas iterações anteriores [s];
real*8	::  a0   ! Aceleração do espaçomodelo na iteração anterior [m/s2];
real*8	::  V0   ! Velocidade do espaçomodelo na iteração anterior [m/s];
real*8	::  h0   ! Altitude do espaçomodelo na iteração anterior [m];
real*8	::  h00  ! Altitude do espaçomodelo duas iterações anteriores [m];
real*8	:: L1_a, L1_V, L1_h ! Norma L1 da aceleração, velocidade e altitude;
real*8	:: a_aux,V_aux,h_aux ! Variáveis auxiliares p/ cálculo da norma;
real*8	:: tol   ! Tolerância para encerrar o processo iterativo;
real*8,dimension(:),allocatable :: Et ! empuxo da tabela empuxo x tempo [N];
real*8,dimension(:),allocatable	:: tt ! Instante de tempo da tabela empuxo x tempo [s];
real*8,dimension(:),allocatable :: Ct ! Cd da tabela Cd x V [adim.];
real*8,dimension(:),allocatable :: Vt ! V da tabela Cd x V [adim.];
! Variáveis auxiliares
integer :: bandeira, bandeira1, io
real*8  :: t1, e1, C1, V1
real*8  :: tcpu ! tempo de processamento

! Leitura do arquivo inicial do programa SBE:
call ler_dados_de_entrada

! Inicialização das variáveis:
pi = dacos(-1.0d0)
tol = 1.0d-15
apogeu = 0
fim_de_queima = 0
abertura_paraquedas = 0
Ar_f  = 2.5d-1*pi*dr*dr
Ar_pq = 1.25d-1*pi*dr_pq*dr_pq
if (ID_E  == 0) It=Ec*tp
R_ar = 286.9d0 ! J/(kgK)
if (ID_ro == 0) ro = p_atm / (R_ar*(T_atm+273.15d0))
if (ID_pq == 0) then
    t_pq = 0.0d0
    h_pq = 0.0d0
    ta_pq = 0.0d0
end if
if (ID_E == 1) then
	write(*,*) 'O arquivo de texto da CURVA DE EMPUXO foi aberto.'
	write(*,*) 'Se a curva de empuxo mostrada precisa ser alterada, efetue as alteracoes e salve o arquivo'
	write(*,*) 'Primeira coluna: Instante de tempo [s]. Segunda coluna: Empuxo [N]'
	write(*,*)
	ver = system('notepad '// tab_E)
	open(1,file = tab_E)
	bandeira = 0
	NpE = 0
	do while (bandeira == 0)
		read(1,*,iostat=io) t1, E1
		if (io < 0) then
			bandeira = 1
		else
			NpE = NpE + 1
		end if
	end do
	close (1)
	open(1,file = tab_E)
	allocate(tt(NpE), Et(NpE))
	It = 0.0d0
	do j = 1, NpE
		read(1,*) tt(j), Et(j)
		if (j /= 1) then
			It = It + 0.5d0*(Et(j)+Et(j-1))*(tt(j)-tt(j-1))
		end if
	end do
	close(1)
	tp = tt(NpE)
end if

if (ID_Cd == 1) then
	write(*,*) 'O arquivo de texto da CURVA DE Cd foi aberto.'
	write(*,*) 'Se a curva de Cd mostrada precisa ser alterada, efetue as alteracoes e salve o arquivo'
	write(*,*) 'Primeira coluna: Velocidade do foguete [m/s]. Segunda coluna: Cd [adm]'
	write(*,*)
	ver = system('notepad '// tab_Cd)
	open(1,file = tab_Cd)
	bandeira = 0
	NpCd = 0
	do while (bandeira == 0)
		read(1,*,iostat=io) V1, C1
		if (io < 0) then
			bandeira = 1
		else
			NpCd = NpCd + 1
		end if
	end do
	close (1)
	open(1,file = tab_Cd)
	allocate(Vt(NpCd), Ct(NpCd))
	do j = 1, NpCd
		read(1,*) Vt(j), Ct(j)
	end do
	close(1)
end if

! Abertura do arquivo de saída de dados e registro de alguns dados de saída e da simulação:
open(2,file='Saida_de_dados_SBE.txt')
call gravar_dados_de_saida_1

! Desenvolvimento dos cálculos principais
tcpu = timef()
bandeira = 0
bandeira1 = 0
VMS = 0.0d0
VMD = 0.0d0
itmax = 100000000
t0 = 0.0d0
a0 = 0.0d0
V0 = 0.0d0
h0 = 0.0d0

write(2,18)
18 format('  >>> Fase Propulsada <<<', /,  &
      5x,'  i', t13,'t [s]', t29,'h [m]',  t45,'V [km/h]', t61,'a [g]', &
      t77,'M [g]', t93,'P [N]', t109,'D [N]',   t125,'E [N]', t141,'Cd [adim.]')
t = t0
V = V0
do i = 1, itmax

	t = t0 + dt

	if (t < tp) then
		M = Mf + Mp*(1.0d0 - t/tp)
 	    !if (t < 40.0d0) then !Executa, de forma bastante improvisada, o cálculo para foguetes de dois estágios. No caso, o Sonda IV.
		!	M = Mf + Mp * ( 1 - t/40.0d0 )
		!else
		!    if (t < 61.0d0) then
		!        M = 2702.0d0
		!        dr = 1.007d+0
		!        Ec = 0.0d0
		!    else
		!        if (t < 83.0d0) then
		!	        Mf = 403.0d0
		!	        Mp = 879.0d0
		!	        M = Mf + Mp * ( 1 - (t-61.0d0)/22.0d0 )
		!	        Ec = 95.9d3
        !           dr = 555.0d-3
        !        else
        !            M = 403.0d0
        !            Ec = 0.0d0
        !            dr = 555.0d-3
        !        end if
        !    end if
		!end if
	else
		if (fim_de_queima == 1) then
		else
			fim_de_queima = 1
			bandeira = 1 ! usada aqui p/ escrever linha de resultados no arquivo de saída dados do fim da queima.
			t = tp
		end if
		
		M = Mf
	end if

	if ( ID_E == 0 ) then
		if (t > tp) then
			E = 0.0d0
		else
			E = Ec
		end if
	else
		if (t > tt(NpE)) then
			E = 0.0d0
		else
			do j = 1, NpE-1
				if ( t >= tt(j) .and. t < tt(j+1) ) then
					E = Et(j) + (Et(j+1)-Et(j))*(t-tt(j))/(tt(j+1)-tt(j))
					exit
				end if
			end do
		end if
	end if

    if (ID_pq == 1) then
        if (t <= t_pq) then
            Ar = Ar_f
        else
            if (t > (t_pq + ta_pq)) then
                Ar = Ar_pq
            else
                Ar = Ar_f + ((t-t_pq)/ta_pq )  *  (Ar_pq - Ar_f)
                if (abertura_paraquedas == 0) then
                    h_pq = h
                    abertura_paraquedas = 1
                end if
            end if
        end if
    else
        Ar = Ar_f
    end if

	a_aux = 0.0d0
	V_aux = 0.0d0
	h_aux = 0.0d0
	do ite = 1, 100

		h = h0 + V*dt
		
		call propriedades

		P = M*g

        if (abertura_paraquedas == 0) then
            D = Cv*0.5d0*ro*V*V*Ar
        else
            D = Cd_pq*0.5d0*ro*V*V*Ar
        end if
		a = -g + (E - sign(D,V))/M

		if (a<0.0d0.and.h<=0.0d0.and.t<tp) a=0.0d0 ! Condição de contato com a plataforma de lançamento.

		V = V0 + dt * a

		! calcula a norma L1 adimensionalizada de a, V e h.
		L1_a = dabs(a - a_aux)
		L1_V = dabs(V - V_aux)
		L1_h = dabs(h - h_aux)

		if (L1_a<tol .and. L1_V<tol .and. L1_h<tol) exit

		a_aux = a
		V_aux = V
		h_aux = h

	end do
	
	h = h0 + V*dt
	
	if (V > VMS) VMS = V
	if (V < VMD .and. V < -1.0d-15) VMD = V

	if (h < h0 .and. apogeu == 0) then ! Apogeu detectado.
		apogeu = 1
		t = t0 - (t-t0)*V0/(V-V0)
		th = t
		h = h00 + (h0-h00)*(t-t00)/(t0-t00)
		hmax = h
		a = -g
		D = 0.0d0
		V = 0.0d0
		if ( ID_Cd == 1 ) then
		   Cv = Ct(1)
		else
		   Cv = Cd
		end if
		write(2,19)
			19 format('  >>> Fase Balística Descendente <<<', /,  &
			  5x,'  i', t13,'t [s]', t29,'h [m]',  t45,'V [km/h]', t61,'a [g]', &
			  t77,'M [g]', t93,'P [N]', t109,'D [N]',   t125,'E [N]', t141,'Cd [adim.]')
		write(2,12) i, t, h, V*3.6d+0, a/g0, M*1000, P, D, E, Cv
	end if
	
	if(ID_pq == 2) then
	    if (apogeu == 1 .and. h <= h_pq) then
	        if (abertura_paraquedas == 0) then
	            t_pq = t
	            abertura_paraquedas = 1
	        end if
	        if (t >= t_pq + ta_pq) then
    	        Ar = Ar_pq
	        else
                Ar = ((t-t_pq)/ta_pq )  *  (Ar_f + Ar_pq)
	        end if
	    else
	        Ar = Ar_f
	    end if
	end if
	
	if (h<0.0d0) exit ! Impacto com o solo detectado.

	t00 = t0
	t0 = t
	a0 = a
	V0 = V
	h00 = h0
	h0 = h

	if (mod(i,w)==0 .or. i==1 .or. bandeira == 1) then ! Escreve resultados no arquivo de saída.
		write(2,12) i, t, h, V*3.6d+0, a/g0, M*1000, P, D, E, Cv
		if (bandeira == 1) then
			write(2,21)
			21 format('  >>> Fase Balística Ascendente <<<', /,  &
			  5x,'  i', t13,'t [s]', t29,'h [m]',  t45,'V [km/h]', t61,'a [g]', &
			  t77,'M [g]', t93,'P [N]', t109,'D [N]',   t125,'E [N]', t141,'Cd [adim.]')
			bandeira = 0
		end if
	end if

end do
if (i > itmax) then
	write(*,*) 'Numero maximo de intervalos de tempo discretizados.'
	write(*,*) 'Aumentar valor da variavel dt e simular novamente.'
end if

12 format (i8, 9 (1x, 1pe15.6))
17 format (i8, 3 (1x, 1pe12.6))

t = t0 - h0*(t-t0)/(h-h0) ! Encontra o momento do impacto por interpolação.
ti = t
a = a0 - h0*(a-a0)/(h-h0)
V = V0 - h0*(V-V0)/(h-h0)
h = h0 - h0*(h-h0)/(h-h0)
Vi = V
if (V < VMD) VMD = V

write(2,12) i, t, h, V*3.6d+0, a/g0, M*1000, P, D, E, Cv

tcpu = timef() - tcpu

call gravar_dados_de_saida_2

close(2)

ver = system('Saida_de_dados_SBE.txt')

if (ID_E == 1) then
	deallocate(tt, Et)
end if
if (ID_Cd == 1) then
	deallocate(Vt, Ct)
end if
!======================================

contains

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

  subroutine ler_dados_de_entrada
  ver = system('Modelo_matematico.txt')
    open(1,file='Modelo_matematico.txt')
    read(1,*)
    read(1,*)
	read(1,*)   ID_E ! Identificador do modelo para a força de empuxo (0 -> E constante; 1 -> E tabelado);
	read(1,*)  ID_Cd ! Identificador do modelo para o coeficiente de arrasto (0 -> Cd constante; 1 -> Cd tabelado);
	read(1,*)   ID_g ! Modelo para a aceleração da gravidade (0 -> g cte; 1-> função da altitude);    
	read(1,*)  ID_ro ! Modelo para a massa específica do ar (0 -> gás perfeito; 1 -> atmosfera padrão internacional);
	read(1,*)   ID_pq ! Variável de comando de acionamento do paraquedas;
    close(1)
    
    open(2,file="Entrada_de_dados_SBE.txt")
	write(2,*)
    write(2,*)' ARQUIVO DE ENTRADA DE DADOS DO PROGRAMA SBE, versao 2.0'
	write(2,*)' Entre com os dados no início de cada linha abaixo:'
	write(2,*) 'SondinhaII   ...: sim_ID  = Identificador da simulação'
    write(2,*) 'voo_geral    ...: caso_ID = Descrição específica da simulação (estudo de caso)'
    write(2,8) 10.0d-3, ' ...: Mp      = Massa de propelente [kg]'
    write(2,8) 20.0E-3, ' ...: Mf      = Massa final do espaçomodelo [kg]'
    write(2,8) 20.0E-3, ' ...: dr      = Diâmetro de referência [m]'
    if (ID_E == 0) then
    	write(2,8) 3.0E+0, ' ...: Ec      = Empuxo constante do motor [N]'
		write(2,8) 0.7E+0, ' ...: tp      = Tempo de queima do propelente [s]'
    else
        write(2,12) 'Curva_empuxo.txt', ' ....: tab_E    = Arquivo de dados de empuxo em função do tempo'
    end if
	if (ID_Cd == 0) then
        write(2,8)   0.9E+0, ' ...: Cd      = Coeficiente de arrasto constante [adim.]'
    else
	    write(2,12)  'Curva_Cd.txt', ' ....: tab_Cd   = Arquivo de dados de Cd em função do tempo'
    end if
    if (ID_g == 0) then
        write(2,9)  9.806563E+0, ' .......: g0      = Aceleração gravitacional local no nível do solo [m/s2]'
    end if
    if(ID_ro == 0) then
        write(2,10)101.325E+3, ' ......:p_atm    = Pressão atmosférica local [Pa]'
        write(2,8)  25.0E+0, ' ...:T_atm    = Temperatura atmosférica local [°C]'
	else
    end if
    if (ID_pq == 1) then
        write(2,8) 5.0E+0, ' ...:t_pq    = Tempo após o lançamento no qual o paraquedas se abre [s]'
        write(2,8) 1.0E+0, ' ...:ta_pq   = Tempo de abertura do paraquedas [s]'
        write(2,8) 2.0E-1, ' ...:dr_pq   = Diâmetro efetivo do paraquedas quando aberto [m]'
        write(2,8) 1.0E+0, ' ...:Cd_pq   = Cd_pq= Coeficiente de arrasto do paraquedas [adim]'
    end if
    if (ID_pq == 2) then
        write(2,8) 1.0E+2, ' ...:h_pq    = Altitude após o apogeu na qual o paraquedas se abre [s]'
        write(2,8) 1.0E+0, ' ...:ta_pq   = Tempo de abertura do paraquedas [s]'
        write(2,8) 2.0E-1, ' ...:dr_pq   = Diâmetro efetivo do paraquedas quando aberto [m]'
        write(2,8) 1.0E+0, ' ...:Cd_pq   = Cd_pq= Coeficiente de arrasto do paraquedas [adim]'
    end if
    write(2,8) 1.0E-4, ' ...: dt      = Intervalo de tempo utilizado na discretização [s]'
    write(2,11) 100, ' ...: w	   = Frequência de escrita dos resultados no arquivo de saída'
    write(2,*)
    write(2,16) '>>> Simulador de Balística de Espaçomodelos - Versão: 2.0, Data: 13/Abr/2018 <<<'
    close(2)
    8 format  ((1pe13.6),A)
	9 format  (F9.6,A)
	10 format (F10.2,A)
	11 format (i13,A)
	12 format (A16,A)
	16 format (A)
  
  ver = system('Entrada_de_dados_SBE.txt')
	open(3,file="Entrada_de_dados_SBE.txt")
	read(3,*)		     ! Linha em branco no arquivo de entrada de dados;
	read(3,*)		     ! Linha em branco no arquivo de entrada de dados;
	read(3,*)		     ! Linha em branco no arquivo de entrada de dados;
	read(3,*) sim_ID     ! Identificador da simulação;
	read(3,*)caso_ID     ! Identificador específico da simulação (estudo de caso);
	read(3,*)     Mp     ! Massa de propelente [kg];
	read(3,*)     Mf     ! Massa final do espaçomodelo [kg];
	read(3,*)     dr     ! Diâmetro de referência [m];
    if (ID_E == 0) then  ! Empuxo constante;
    	read(3,*)     Ec ! Empuxo constante do motor [N];
		read(3,*)     tp ! Tempo de queima do propelente [s];
    else                 ! Empuxo dado pela tabela de empuxo;
        read(3,*) tab_E  ! Nome da tabela onde estão os dados do empuxo em função do tempo;
    end if
	if (ID_Cd == 0) then ! Cd constante;
        read(3,*)     Cd ! Coeficiente de arrasto [adim.];
    else                 ! Cd variável dado pela tabela de Cd função da velocidade;
	    read(3,*) tab_Cd ! Nome da tabela onde estão os dados do Cd em função da velocidade;
    end if
    if (ID_g == 0) then  ! Aceleração da gravidade constante;
        read(3,*)     g0 ! Aceleração gravitacional local no nível do solo [m/s2];
    else                 ! Aceleração da gravidade variável, função da altitude;
		g0 = 9.806563d0	 ! Aceleração gravitacional local [m/s2];
    end if
    if(ID_ro == 0) then  ! Massa espescífica do ar dado pelo modelo de gás perfeito;
        read(3,*)  p_atm ! Pressão atmosférica local [Pa];
        read(3,*)  T_atm ! Temperatura atmosférica local [°C];
    else                 ! Massa específica do ar dado pela atmosfera padrão internacional;
                         ! ver subrotina "propriedades";
    end if
    if (ID_pq == 1) then
        read(3,*) t_pq   ! Tempo após o lançamento no qual o paraquedas se abre [s]'
        read(3,*) ta_pq  ! Tempo de abertura do paraquedas [s]'
        read(3,*) dr_pq  ! Diâmetro efetivo do paraquedas quando aberto [m]'
        read(3,*) Cd_pq  ! Coeficiente de arrasto do paraquedas [adim]'
    end if
    if (ID_pq == 2) then
        read(3,*) h_pq   ! Altitude após o apogeu na qual o paraquedas se abre [s]'
        read(3,*) ta_pq  ! Tempo de abertura do paraquedas [s]'
        read(3,*) dr_pq  ! Diâmetro efetivo do paraquedas quando aberto [m]'
        read(3,*) Cd_pq  ! Coeficiente de arrasto do paraquedas [adim]'
    end if
	read(3,*)    dt      ! Intervalo de tempo utilizado na discretização [s];
    read(3,*)     w 	 ! Frequência de escrita dos resultados no arquivo de saída.
	close(3)
  end subroutine ler_dados_de_entrada

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

  subroutine gravar_dados_de_saida_1
	write(2,*)'======================================='
	write(2,*)'       >>>  SBE   versão 2.0  <<<'
	write(2,*)'Arquivo de saída de dados da simulação'
	write(2,*)'======================================='
	write(2,*)
	write(2,*) trim(adjustl(sim_ID)) , ' = Identificador da simulação'
	write(2,*) trim(adjustl(caso_ID)), ' = Estudo de caso (identificador específico da simulação)'
	write(2,1) Mf, Mp, Mp+Mf, dr, tp
	1 format(1pe13.6, ' =  Mf  = Massa final [kg]',/, &
			 1pe13.6, ' =  Mp  = Massa de propelente [kg]',/, &
			 1pe13.6, ' =  Mo  = Massa inicial [kg] = Mf + Mp',/, &
			 1pe13.6, ' =  dr  = Diâmetro de referência [m]',/, &
			 1pe13.6, ' =  tp  = Tempo de queima do propelente [s]')
	if (ID_E == 0) then
		write(2,*) 'Modelo p/ a força de empuxo: empuxo constante'
		write(2,5) Ec, It
		5 format(1pe13.6, '  = Ec   = Empuxo constante do motor [N]',/, &
				 1pe13.6, '  = It   = Impulso total [Ns]')
	else
		write(2,*) 'Modelo p/ a força de empuxo: empuxo variável conforme tabela abaixo'
	end if
	if (ID_Cd == 0) then
		write(2,*) 'Modelo p/ o coeficiente de arrasto: Cd constante'
		write(2,6) Cd
		6 format(1pe13.6, '  = Cd  = Coeficiente de arrasto [adim.]')
	else
		write(2,*) 'Modelo p/ Cd: Cd função da velocidade conforme tabela abaixo'
	end if
	if (ID_g == 0) then
		write(2,*) 'Modelo p/ a aceleração gravitacional: g constante'
		write(2,2) g0
		2 format(1pe13.6, ' =  g0  = Aceleração gravitacional local [m/s2]')
	else
		write(2,3)
		3 format(' Modelo p/ a aceleração gravitacional: g função da altitude')
		write(2,2) g0
	end if
	if (ID_ro == 0) then
		write(2,*) 'Modelo p/ a massa específica do ar: equação de gás perfeito'
		write(2,4) p_atm, T_atm, R_ar, ro
		4 format(1pe13.6, ' = p_atm= Pressão atmosférica local [Pa]',/, &
				 1pe13.6, ' = T_atm= Temperatura atmosférica local [°C]',/, &
				 1pe13.6, ' = R_ar = Constante dos gases perfeitos do ar [J/kg.K]',/, &
				 1pe13.6, ' = ro  = Massa específica do ar [kg/m3]')
	else
		write(2,*)'Modelo p/ a massa específica do ar: atmosfera internacional padrão'
	end if
    if (ID_pq == 0) write(2,*) 'Modelo p/ comando de abertura do paraquedas: sem paraquedas.'
    if (ID_pq == 1) then
        write(2,*) 'Modelo p/ comando de abertura do paraquedas: tempo após lançamento'
        write(2,22) t_pq,   ' = t_pq =  Tempo após o lançamentono no qual o paraquedas se abre [s]'
        write(2,22) ta_pq,  ' = ta_pq= Tempo de abertura do paraquedas [s]'
        write(2,22) dr_pq,  ' = dr_pq= Diâmetro efetivo do paraquedas quando aberto [m]'
        write(2,22) Cd_pq,  ' = Cd_pq= Coeficiente de arrasto do paraquedas [adim]'
    end if
    if (ID_pq == 2) then
        write(2,*) 'Modelo p/ comando de abertura do paraquedas: altitude após o apogeu'
        write(2,22) h_pq,   ' = h_pq = Altitude após o apogeu na qual o paraquedas se abre [s]'
        write(2,22) ta_pq,  ' = ta_pq= Tempo de abertura do paraquedas [s]'
        write(2,22) dr_pq,  ' = dr_pq= Diâmetro efetivo do paraquedas quando aberto [m]'
        write(2,22) Cd_pq,  ' = Cd_pq= Coeficiente de arrasto do paraquedas [adim]'
    end if
    22 format (1pe13.6,A)
	write(2,7) dt, w
	7 format(1pe13.6, ' =  dt = Intervalo de tempo utilizado na discretização [s]',/, &
             i13,1x, '=  w   = Freqüência de escrita dos resultados de campo',/)
	if (ID_E == 1) then
		write(2,*) '======================================='
		write(2,*) 'Tabela de empuxo utilizada na simulação:'
		write(2,13) 't [s]','E [N]'
		do j = 1, NpE
			write(2,14) tt(j), Et(j)
		end do
		write(2,*)'Impulso total (regra retangulo):', It, '[Ns]'
		write(2,*)
	end if
	13 format (4x, A5, 20x, A5)
	14 format (2(1pe25.15))
	if (ID_Cd == 1) then
		write(2,*) '======================================='
		write(2,*) 'Tabela de coeficiente de arrasto utilizada na simulação:'
		write(2,15) 'V [m/s]','Cd [adm.]'
		do j = 1, NpCd
			write(2,14) Vt(j), Ct(j)
		end do
		write(2,*)
	end if
	15 format (4x, A7, 18x, A9)
	write(2,*) '======================================='

  end subroutine gravar_dados_de_saida_1

  !-------------------------------
  subroutine gravar_dados_de_saida_2

	write(2,*) '======================================='
	write(2,*)
	write(2,*) 'Dados globais do voo:'
	write(2,*)
	write(2,20) 'Tempo de voo até o apogeu:', th, ' s   '
	if (ID_pq /= 0) then
	    write(2,20) 'Início da abertura do paraquedas:', t_pq, ' s   '
	    write(2,20) 'Altitude na abertura do paraquedas:', h_pq, ' m   '
	end if
	write(2,20) 'Altitude no apogeu:', hmax, ' m   '
	write(2,20) 'Velocidade máxima na subida:', VMS*3.6d0, ' km/h'
	write(2,20) 'Velocidade máxima na descida:', VMD*3.6d0, ' km/h'
	write(2,20) 'Velocidade de impacto:', Vi*3.6d0, ' km/h'
	write(2,20) 'Tempo de voo até o impacto:', ti, ' s   '
	write(2,20) 'Tempo de processamento:', tcpu, ' s   '
	write(2,*)
	write(2,*)'============================================='
	write(2,*)'          >>>  SBE   versão 2.0  <<<'
	write(2,*)'  Suportado pela Coordenação de Pesquisa e'
	write(2,*)' pela Coordenação de Engenharia Mecânica da'
	write(2,*)'  Universidade Positivo - UP, Curitiba/PR'
	write(2,*)'Programa de Iniciação Científica - PIC1434/15'
	write(2,*)'   Programadores: Jeremy Blavier Paraná'
	write(2,*)'                  Antonio Carlos Foltran'
	write(2,*)'============================================='
	20 format(A35, 1pe25.15, A5)

  end subroutine gravar_dados_de_saida_2

  !-------------------------------
  subroutine propriedades
    real*8 :: Temp ! temperatura [K];
	real*8 :: pres ! pressão [Pa];
    real*8 :: aux
	if ( ID_Cd == 1 ) then
		do j = 1, NpCd-1
			if ( dabs(V) >= Vt(j) .and. dabs(V) < Vt(j+1) ) then
				Cv = Ct(j) + (Ct(j+1)-Ct(j))*(dabs(V)-Vt(j))/(Vt(j+1)-Vt(j))
				exit
			end if
		end do   
		if ( dabs(V) > Vt(NpCd) .and. bandeira1 == 0) then
			write(*,*)'************************************************'
			write(*,*)'ATENCAO: Cd requerido alem dos valores da tabela'
			write(*,*)'************************************************'
			bandeira1 = 1
		end if
	else
		Cv = Cd
	end if

	if (ID_g == 1) then
		g = 9.806563d0*(6378137.0d0 / (6378137.0d0 + h))**2.0d0
		!g = 3.989371d14 / ((6378137d0 + h)**2.0d0)
	else
		g = g0
	end if

	if (ID_ro == 1) then
!		if (h < 11.0d3) then ! troposfera   ! Fórmulas extraídas de CAVCAR, Mustafa. The International Standard Atmosphere (ISA). Eskisehir, TR: Anadolu University. Documento não publicado.
!			Temp = 288.15d0 - 6.5d0*h/1.0d3
!			pres = 101325.0d0*(1.0d0 - 6.5d-3*h/288.15d0)**5.2561d0
!			ro = pres/(287.04d0*Temp)
!		else
!			if (h >= 11.0d3 .and. h < 20.0d3) then ! tropopausa
!				Temp = 216.65d0
!				pres = 22630.61892d0*dexp(-9.80665d0*(h - 11.0d3)/(287.04d0*216.65d0))
!				ro = pres/(287.04d0*Temp)
!			else ! acima da tropopausa
!				ro = 0.0d0
!			end if
!		end if
		if (h < 47.0d+3) then ! Fórmulas extraídas do capítulo 1 de DRELA, Mark. Flight Vehicle Aerodynamics, Cambridge: The MIT Press, 2014.
            aux = 1.0d-3*h
			Temp = 216.65d0 + 2.0d0*dlog(1.0d0 + dexp(35.75d0 - 3.25d0*aux) + dexp(-3.0d0 + 3.0d-4*aux**3.0d0))
			pres = 1.0132d5*dexp(-0.118d0*aux - 0.0015d0*aux**2.0d0/(1.0d0-0.018d0*aux+0.0011d0*aux**2.0d0))
            ro = pres/(287.04d0*Temp)
        else
            ro = 0.0d0
        end if
	end if

  end subroutine propriedades
end program SBE