!************************************************************
!SEGUNDA VERSÃO DO PROGRAMA PARA CÁLCULO DO CD	(17/03/2016)
!BASEADA NOS RELATÓRIOS DE C.P. HOULT, NISKANEN E CANNON
!AUTOR: TOBIAS PINHEIRO QUELUZ
!************************************************************
PROGRAM Coeficiente_de_arrasto

!EXTERNAL Cd_atrito, Cd_base, Nariz_geometria, Area_laminar, Cd_tuboguia, PORTLIB 
!
USE PORTLIB

!DEFINIÇÃO DAS VARIÁVEIS

INTEGER motorligado, Na, Nt, It, i, U, motorexterno, Ac, Nest 

REAL*8 To, Po, Retc, Rete, M, Pi, L, Swet, Dc, Sw, Sref, Cf, Cfl, Cft, C, Cdac, &
       Le, Cdae, Cfe, Cfle, Cfte, Swete, e, Rc, Tc, S, VER, N, Cda, Cdbc, Cdbe, Cdb, &
	   Dt, St, Lna, Dna, Lc, Swetc, Swetna, Cdtb, Detb, Ditb, Ltb, Cbc, Cdpe, Cdi, &
       etb, Mmax, Cd, Sc, Sa, Cdpn, Soma, Cdm, V, Lm, Dm, Stf, Cdpv1, Sm, Dref, &
       Dest, Lest, Sest, Cdest, Cdpa, Sna, Cdam, Cdpv2
       
REAL*8 To1, Po1

CHARACTER (100) SAIDA, TEXTO, TEXTO1, TEXTO2, DADOS

!CHAMADA DO ARQUIVO DE ENTRADA

VER = SYSTEM('Notepad dados.txt')

OPEN(2, file ='dados.txt')

READ(2,*) DADOS                    !Nome do arquivo de entrada

CLOSE(2)

VER = SYSTEM('Notepad ./Entrada/'//DADOS)

OPEN(1, file ='./Entrada/'//DADOS)

READ(1,*) SAIDA                    !Nome do arquivo de saída

!CONSTANTES FÍSICAS
READ(1,*) U                        !Unidade utilizada nas medidas do foguete (1=milímetros, 2=pés, 3=polegadas)
READ(1,*) To1                      !Temperatura atmosférica (°C)
READ(1,*) Po1                      !Pressão atmosférica (kPa)
READ(1,*) Retc                     !Número de Reynolds de transição para corpos de revolução
READ(1,*) Rete	                   !Número de Reynolds de transição para empenas
READ(1,*) Mmax                     !Número de mach máximo
READ(1,*) It                       !Número de pontos do gráfico
To   = (To1+273.15d0)*9d0/5d0      !Temperatura atmosférica absoluta (deg R)
Po   = Po1*20.885434273d0          !Pressão atmosférica (lb/ft²)
Pi   = dacos(-1d0)                 !Número de Pi

!GEOMETRIA DO NARIZ
READ(1,*) Na                       !Tipo de nariz (1=elpsoide, 2=ogiva, 3=cônico, 4=parabólico)
READ(1,*) Lna                      !Comprimento total do corpo (mm)
READ(1,*) Dna                      !Diâmetro do corpo (mm)

!GEOMETRIA DO TUBO-FOGUETE
READ(1,*) Lc                       !Comprimento do tubo (mm)
READ(1,*) Dc                       !Diâmetro do tubo (mm)

!GEOMETRIA DA EMPENA
READ(1,*) N                        !Número de empenas
READ(1,*) Ac                       !Tipo de seção transversal(0=quadrada, 1=arredondada)
READ(1,*) e                        !Espessura (mm)
READ(1,*) Rc                       !Comprimento da base (mm)
READ(1,*) Tc                       !Comprimento do topo (mm)
READ(1,*) S                        !Altura (mm)
READ(1,*) Sc                       !Comprimento do contorno (mm)

!GEOMETRIA DO MOTOR
READ(1,*) motorexterno	           !Posição do motor (Externo/Interno)
READ(1,*) Lm	                   !Comprimento do motor (mm)
READ(1,*) Dm	                   !Diâmetro do motor (mm)
READ(1,*) motorligado	           !Condição do motor (ligado/desligado)
READ(1,*) Dt                       !Diâmetro da tubeira (mm)

!GEOMETRIA DO TUBO-GUIA
READ(1,*) Nt                       !Presença de tubo-guia (1=sim, 0=não)
READ(1,*) Detb                     !Diâmetro externo do tubo-guia (mm)
READ(1,*) etb                      !Espessura do tubo-guia (mm)
READ(1,*) Ltb                      !Comprimento do tubo-guia (mm)

!GEOMETRIA DO ESTOPIM
READ(1,*) Nest                     !Presença de estopim externo (1=sim, 0=não)
READ(1,*) Dest                     !Diâmetro externo do estopim (mm)
READ(1,*) Lest                     !Comprimento do estopim (mm)

CLOSE(1)

!CONVERSÃO DE UNIDADES
SELECT CASE (U)
    
CASE(1) !Milímetros
    
CASE(2) !Pés para milímetros
    CALL pes(Lna, Dna, Lc, Dc, e, Rc, Tc, S, Sc, Dt, Detb, etb, Ltb, Lm, Dm)
    
CASE(3) !Polegadas para milímetros
    CALL polegadas(Lna, Dna, Lc, Dc, e, Rc, Tc, S, Sc, Dt, Detb, etb, Ltb, Lm, Dm)
    
END SELECT   

!GEOMETRIA CALCULADA

!Nariz
CALL Narizgeo (Pi, Na, Lna, Dna, Swetna)    !Área do nariz em contato com o ar (mm²)
Sna = pi*Dna**2d0/4d0                       !Área de referência do nariz (mm²)

!Tubo-foguete
Swetc = Pi*Dc*Lc                            !Área do tubo em contato com o ar (mm²)
Stf = Pi*Dc**2d0/4d0                        !Área de referência do tubo-foguete (mm²)

!Corpo
L   = Lna + Lc                              !Comprimento total do corpo (mm) 
Swet = Swetna + Swetc                       !Área em contato com o ar (mm²)
IF (motorexterno==0) THEN
    Dref = max(Dc,Dna)                      !Diâmetro de referência
    ELSE
    Dref = max(Dc,Dna,Dm)                   !Diâmetro de referência  
END IF                     
Sref = Pi*Dref**2d0/4d0                     !Área de referência (igual à seção tranversal do corpo) (mm²)

!Empena
Sa = datan(Sc/S)                            !Ângulo do contorno (rad)
Swete = (Rc+Tc)*S                           !Área em contato com o ar (mm²)
Le = (Rc+Tc)/2d0                            !Comprimento médio da empena (mm)
    
!Motor
St = Pi*Dt**2d0/4d0	                        !Área de saída dos gases (mm²)
Sm = Pi*Dm**2d0/4d0                         !Área de referência do motor (mm²)
    
!tubo-guia
Ditb = Detb - 2d0*etb                       !Diâmetro interno do tubo-guia (mm)

!Estopim externo
Sest = Pi*Dest**2d0/4d0                     !Área de referência do estopim (mm²)

!GRAVAÇÃO DOS COEFICIENTES DE ARRASTO
TEXTO1 = 'Curva_'//trim(adjustl(SAIDA))//'.txt'
TEXTO2 = 'Curva_trajetória_'//trim(adjustl(SAIDA))//'.txt'
OPEN(7, file='./Saida/'//TEXTO1)
OPEN(9, file='./Saida/'//TEXTO2)
WRITE(7,*) " # Número de Mach            V(km/h)           Coeficiente de arrasto"
TEXTO1 = "Notepad "//'./Saida/'//TEXTO1
TEXTO2 = "Notepad "//'./Saida/'//TEXTO2
Soma = 0d0

DO i= 1, It, 1
    
    M = (Mmax/(It*1d0))*(i*1d0)
    V = M*1224d0
    
    !COEFICIENTE DE ARRASTO DE ATRITO
    CALL CdAtritoc (To,Po,Retc,M,Cf,Swet,Cfl,Cft,C,L,Sref, Cdac, Pi, Lna, Dna, Na, Swetna, Dc)
    CALL CdAtritoe (To,Po,Rete,M,Cfe,Swete,Cfle,Cfte,C,Le,Sref, Cdae, N)
    SELECT CASE (motorexterno)
        CASE(1)
        CALL CdAtritom(To, Po, M, Lm, Dm, Cdam, Sref)
        CASE(0)
        Cdam = 0d0    
    END SELECT      
    Cda = Cdac + Cdae + Cdam   !Coeficiente de arrasto de atrito total
    

    !COEFICIENTE DE ARRASTO DE BASE
    CALL Cdbasec (Sref, Pi, Dc, Cdbc, St, M, motorligado, Cbc, Cda, Cf)
    CALL Cdbasee (N, e, S, Sref, Cfe, Cdbe)
    
    Cdb = Cdbc          !Coeficiente de arrasto de atrito total
    !Cdb = Cdbc +Cdbe   !Coeficiente de arrasto de atrito total

    !COEFICIENTE DE ARRASTO PARASITA
    CALL Cdtubo (Pi, Detb, Ditb, Ltb, Sref, Cdtb, Nt, M, Swetc)
    CALL Cdestopim (Pi, Dest, Lest, Sref, Cdest, Nest, M)
    
    Cdpa = Cdtb + Cdest

   !COEFICIENTE DE ARRASTO DE PRESSÃO
    CALL Cdpressaoe(M, Tc, S, Cbc, N, e, Sref, Cdpe, Sa, Ac, Rc, Sc, Pi)
    CALL Cdpressaon(M, Na, Dna, Lna, Sref, Sna, Cdpn)
    CALL Cdpressaov2(M, Dna, Dc, Sna, Stf, Sref, Cdpv2)
    
    SELECT CASE (motorexterno)
        CASE(1)
        CALL Cdpressaov1(M, Dm, Dc, Sm, Stf, Sref, Cdpv1)
        CASE(0)
        Cdpv = 0d0    
    END SELECT      
      
    Cdp = Cdpe + Cdpn + Cdpv1 + Cdpv2
    !write(*,*) Dref
    !Cdp = Cdpn 
    
    !COEFICIENTE DE ARRASTO DE INTERFERÊNCIA
    CALL Cdinterferencia(Cfe, e, Le, Tc, Dc, N, Sref, Cdi) 


    !COEFICIENTE DE ARRASTO TOTAL
    Cd = (Cda+Cdb+Cdpa+Cdp+Cdi)*1.05
   
    !CÁLCULO DO CD MÉDIO
    Soma = Cd + Soma
    Cdm = Soma/i
    
    WRITE(7,*) M, V, Cd
    WRITE(9,*) V, Cd
    
END DO

CLOSE(7)
CLOSE(9)

VER = SYSTEM(TEXTO1)
!VER = SYSTEM(TEXTO2)

!SAÍDA DOS RESULTADOS

!Vizualização do gráfico
CALL graficos(SAIDA)
VER = SYSTEM('.\gnuplot\wgnuplot Comandos_Gnu.gnu')

OPEN(8,file='./Saida/'//trim(adjustl(SAIDA))//".txt")

WRITE(8,*)"COEFICIENTE DE ARRASTO PARA MACHI =", Mmax
WRITE(8,*)""
WRITE(8,*)"Coeficiente de arrasto de atrito:", Cda
WRITE(8,*)""
WRITE(8,*)"Coeficiente de arrasto de base:", Cdb
WRITE(8,*)""
WRITE(8,*)"Coeficiente de arrasto do tubo-guia:", Cdpa
WRITE(8,*)""
WRITE(8,*)"Coeficiente de arrasto de pressão:", Cdp
WRITE(8,*)""
WRITE(8,*)"Coeficiente de arrasto de interferência:", Cdi
WRITE(8,*)""
WRITE(8,*)""
WRITE(8,*)"COEFICIENTE DE ARRASTO TOTAL:", Cd
WRITE(8,*)""
WRITE(8,*)"COEFICIENTE DE ARRASTO MÉDIO:", Cdm

CLOSE(8)

TEXTO = "Notepad "//'./Saida/'//trim(adjustl(SAIDA))//".txt"
VER = System(TEXTO)

END PROGRAM Coeficiente_de_arrasto
