program voo
use portlib
implicit none
character*128 :: arquivo, traj, temp 
real*8 :: tP, TBA, tBD, De, A, Pa, T, mo, mej, g
real*8, parameter :: Rg = 286.9d0 ! Constante de gás perfeito do ar (J/kg.K) 
real*8, parameter :: Pi = dacos(-1.d0)
real*8 :: rho, h, ha

real*8,allocatable, dimension(:,:) :: Alt
! Alt(:,1) = tempo (s)
! Alt(:,2) = altitude (m)

integer :: i, j, k, l
integer :: n ! Número de linhas da trajetória
integer :: nP ! Ponto final da fase propulsada
integer :: nBA ! Ponto final da fase Balistica Ascendente
integer :: nBD ! Ponto final da fase Balística Descendente

real*8 :: ThetaP, BP

integer :: na ! Número de parâmetros
integer :: ene ! Número de equações

real*8, allocatable, dimension(:) :: pFP ! parâmetros da fase propulsada
real*8, allocatable, dimension(:) :: pBA
real*8, allocatable, dimension(:) :: pBD
real*8, allocatable, dimension(:) :: D ! Vetor com as diferenças entre as medidas e os valores da função
real*8, allocatable, dimension(:,:) :: Drag
real*8, allocatable, dimension(:) :: Forca, ForcaSemArrasto

real*8 :: sumq, res, erro, pant
real*8 :: dt, tempo, vel, acel, Fd, It, ItFP, c, ItSA
real*8 :: SStot, SSres

real*8 :: t1, t2
integer :: Nm, Nint, Nres
logical :: existe
    
    ! Abre o arquivo que dará a entrada da trajetoria
    open(1,file='.\entrada\voo.in')
        read(1,*) Nint
        read(1,*) Nres
        read(1,*) arquivo  
    close(1)
    
    ! Abre o arquivo da entrada da trajetória
    open(1,file='.\entrada\' // trim(adjustl(arquivo)))
        read(1,*) traj
        read(1,*) tP
        read(1,*) tBA
        read(1,*) tBD        
        read(1,*) De
        read(1,*) Pa
        read(1,*) T
        read(1,*) mo
        read(1,*) mej
        read(1,*) g
        read(1,*) Nm
    close(1)
        
    allocate(Forca(0:Nm), ForcaSemArrasto(0:Nm))
    Forca=0;ForcaSemArrasto=0
    
    ! Confere se existe a pasta de saída e as pastas auxiliares, senão cria-as
    INQUIRE (DIRECTORY='.\saida', EXIST=EXISTE)
    If(.not. Existe) existe = system("md .\saida")
    
    INQUIRE (DIRECTORY='.\saida\FBA', EXIST=EXISTE)
    If(.not. Existe) existe = system("md .\saida\FBA")
    
    INQUIRE (DIRECTORY='.\saida\FBD', EXIST=EXISTE)
    If(.not. Existe) existe = system("md .\saida\FBD")
    
    INQUIRE (DIRECTORY='.\saida\FP', EXIST=EXISTE)
    If(.not. Existe) existe = system("md .\saida\FP")    

    ! Calculo da massa específica do ar (kg/m³)
    rho = Pa/Rg/T
    
    ! Cálculo da área da seção do foguete
    A = Pi/4*De*De
    
    ! Início da leitura dos dados da trajetória
    ! Conta o número de linhas presentes no arquivo
    open(1,file='.\entrada\' // trim(adjustl(traj)) )
        do i=1, 1000000         
            read(1,*,iostat=n) 
            if(n /= 0) exit
        enddo   
    ! Calcula-se o número de linhas        
    n = i-1
    allocate( Alt(n,2), D(n) )
    Alt = 0; D=0
    rewind(1)    
    ! Leitura da trajetória
        do i=1, n
            read(1,*) Alt(i,1), Alt(i,2)
        enddo
    close(1)
    
    ! Cria gráfico do voo completo    
    open(1,file='.\saida\voo.txt')
        do i=1, N
            write(1,'(10(1PE25.13))') Alt(i,1), Alt(i,2)
        enddo        
    close(1)
    
    open(1,file='.\saida\voo.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Altitude (m)'"
        write(1,'(A)') "set out '.\saida\voo.png'"
        write(1,'(A)') "plot '.\saida\voo.txt' using 1:2 title '' lw 2 "    
    close(1)
    
    i = system("Wgnuplot ./saida/voo.gnu")
    
    t1 = timef()
    ! Separação do voo nas 3 fases
    ! Fase Propulsada (FP)
    ! Definição: enquanto o motor ainda estiver ligado
    ! na fase propulsada: E-D-P=m.a, enquanto E>0
    ! ou seja, enquanto m.(a+g)+D>0
    
    ! Fase Balística Ascendente (FBA)
    ! Definição: enquanto a velocidade Vertical não for Nula
    ! -D-mf*g=mf*a
    ! mf*(a+g)+D=0
    
    ! Fase Balística Descendente (FBD)
    ! Definição : enquanto a velocidade não for nula após a FBA 
    
    do i=1, N
        if(Alt(i,1) <= tP) nP = i
        if(Alt(i,1) <= tBA) nBA = i
        if(Alt(i,1) <= tBD) nBD = i        
    enddo   
    
    ! Cria gráfico das fases separadas   
    open(1,file='.\saida\FP\FP.txt')
        do i=1, nP
            write(1,'(10(1PE25.13))') Alt(i,1), Alt(i,2)
        enddo        
    close(1)
    open(1,file='.\saida\FBA\FBA.txt')
        do i=np, nBA
            write(1,'(10(1PE25.13))') Alt(i,1), Alt(i,2)
        enddo        
    close(1)
    open(1,file='.\saida\FBD\FBD.txt')
        do i=nBA, nBD
            write(1,'(10(1PE25.13))') Alt(i,1), Alt(i,2)
        enddo        
    close(1)    
    
    open(1,file='.\saida\fases1.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Altitude (m)'"
        write(1,'(A)') "set out '.\saida\fases.png'"
        write(1,'(A)') "plot '.\saida\FP\FP.txt' using 1:2 title 'Propulsada' lw 2 \"    
        write(1,'(A)') "   , '.\saida\FBA\FBA.txt' using 1:2 title 'Balística Ascendente' lw 2 \"    
        write(1,'(A)') "   , '.\saida\FBD\FBD.txt' using 1:2 title 'Balística Descendente' lw 2 "    
    close(1)
    
    i = system("Wgnuplot ./saida/fases1.gnu")   
    
    ! Ajuste da fase propulsada
    write(*,'(A9, 10(A13))') "It", "R2", "Sum Quad"
    write(*,*) "FP"
    ene = nP
    na  = 2 ! Apenas 2 parâmetros
    allocate( pFp(na) )    
    ! Estimativa dos parâmetros
    pFp(1) = 1*rho*A/2/(mo-mej/2) ! ThetaP, admitindo Cd=1
    pFp(2) = acosh( sqrt(exp(2*Alt(nP,2)*pFp(1)) ) )**2/pFp(1)/tP**2;  ! B = (F-g*m)/m, usando o ponto final da curva para interpolar
    !write(*,*) pFp;pause
    pant = pFp(1)*pFp(2)
    
    !call ajusteFP
    do k=1, nres
        res = 10d0**(-1-k)
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"
        open(1,file=".\saida\FP\"//trim(adjustl(temp)) )
        do j=1, nint            
                call diff_sq_FP
                sumq = sum(D(1:np)**2)
                ha = res*pFp(1)
                pFp(1) = pFp(1) + ha/2
                call diff_sq_FP     
                erro = sumq - sum(D(1:nP)**2)
                h = 1.*sign(ha, erro )
                pFp(1) = pFp(1) + h*3
    
                call diff_sq_FP
                sumq = sum(D(1:nP)**2)       
                ha = res*pFp(2)
                pFp(2) = pFp(2) + ha/2   
                call diff_sq_FP                
                erro = sumq - sum(D(1:nP)**2)
                h = 1.*sign(ha, erro )
                pFp(2) = pFp(2) + h*3  
    
    
                SStot = sum( (Alt(1:nP,2) - sum(Alt(1:nP,2) )/nP )**2)
                SSres = 0
                do l=1, nP
                    SSres = SSres + (Alt(l,2)-hP(Alt(l,1), pFp(1),pFp(2)) )**2
                enddo                                   
            write(1,'(I9,100(1PE25.15))') (k-1)*nint+j, pFp(1)*pFp(2), abs(pFp(1)*pFp(2)-pant), sumq, 1-SSres/SStot
            pant = pFp(1)*pFp(2)
        enddo
        close(1)
        !write(*,'(i5, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    enddo
    write(*,'(i9, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    
    ! Cria gráfico das fases separadas   
    open(1,file='.\saida\FP\FPa.txt')
        dt = (tP-0)/Nm
        do i=0, Nm   
            tempo = 0+dt*i
            write(1,'(10(1PE25.13))') tempo, hP(tempo, pFp(1),pFp(2))
        enddo        
    close(1)
    
    open(1,file='.\saida\fasesP.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Altitude (m)'"
        write(1,'(A)') "set out '.\saida\FP\FP.png'"
        write(1,'(A)') "plot '.\saida\FP\FP.txt' using 1:2 title 'Propulsada' lw 2 \"    
        write(1,'(A)') "   , '.\saida\FP\FPa.txt' using 1:2 title 'Propulsada - ajustada' lw 2 "   
    close(1)
    
    i = system("Wgnuplot ./saida/fasesP.gnu")   
    
    ! Cria gráfico dos resíduos
    open(1,file='.\saida\FP\FPr.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'Iteracao'"
        
        write(1,'(A)') "set out '.\saida\FP\FPr2.png'"
        write(1,'(A)') "set ylabel 'R²'"
        
        !k = 0
        !write(temp,*) k
        !write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"        
        !write(1,'(A,A,A,I1,A)') "plot '.\saida\FP\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 "  
        !
        !write(1,'(A)') "set ylabel 'Soma dos quadrados dos resíduos'"
        !write(temp,*) k
        !write(1,'(A)') "set out '.\saida\FP\FPr_" // trim(adjustl(temp)) //  ".png'"        
        !write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"
        !write(1,'(A,A,A,I2,A)') "plot '.\saida\FP\",trim(adjustl(temp)),"' using 1:4 title '",k,"' lw 2 "        
             
        k=1
        !res = 10d0**(-1-k)
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"        
        write(1,'(A,A,A,I1,A)') "plot '.\saida\FP\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 \"            
        do k=2, nres-1
            !res = 10d0**(-1-k)
            write(temp,*) k+1
            write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "   , '.\saida\FP\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 \"
        enddo 
        k = nres
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"
        write(1,'(A,A,A,I2,A)') "   , '.\saida\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 "         
        
        !write(1,'(A)') "set format y '10**{%L}'"       
        write(1,'(A)') "set ylabel 'Soma dos quadrados dos resíduos'"
        !write(1,'(A)') "set logscale y"
                   
        do k=1, nres
            write(temp,*) k+1
            write(1,'(A)') "set out '.\saida\FP\FPr_" // trim(adjustl(temp)) //  ".png'"        
            write(temp,'(A,A,A)') "FP_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "plot '.\saida\FP\",trim(adjustl(temp)),"' using 1:4 title '",k,"' lw 2 "
        enddo           
    
    close(1)
    
    i = system("Wgnuplot ./saida/FP/FPr.gnu")       
    
    
    
    !pause
    
    
    ! Ajuste da fase balistica ascendente
    write(*,*) "FBA"
    ene = nBA-np
    na  = 4 
    allocate( pBA(na) )    
    ! Estimativa dos parâmetros
    pBA(1) = rho*A/2/(mo-mej) ! Theta, Cd = 1
    pBA(2) = tP !1d-0 ! tq
    pBA(3) = hP(tP, pFp(1),pFp(2)) !1d-0 ! hq
    pBA(4) = vP(tP, pFp(1),pFp(2)) ! vq
    pant = pBA(1)*pBA(2)*pBA(3)*pBA(4)
    !write(*,'(i5, 10(1PE15.3))') j, pBA,sumq ;pause
    !call ajusteFBA
    do k=1, nres
        res = 10d0**(-1-k)
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBA_",trim(adjustl(temp)),".txt"
        open(1,file=".\saida\FBA\"//trim(adjustl(temp)) )        
        do j=1, nint 
            
    !! Cria gráfico das fases separadas   
    !open(1,file='.\saida\FBA\FBAa.txt')
    !    dt = (tBA-tP)/Nm
    !    do i=0, Nm
    !        tempo = tP+dt*i
    !        write(1,'(10(1PE25.13))') tempo, hBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4))
    !    enddo        
    !close(1) 
    !
    !open(1,file='.\saida\fasesBA.gnu')
    !    write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
    !    write(1,'(A)') "set style data lines"
    !    write(1,'(A)') "set key right box top"
    !    write(1,'(A)') "set time"
    !    write(1,'(A)') "set grid"
    !    write(1,'(A)') "set title '' "
    !    write(1,'(A)') "set xlabel 'tempo (s)'"
    !    write(1,'(A)') "set ylabel 'Altitude (m)'"
    !    write(1,'(A)') "set out '.\saida\FBA\FBA.png'"
    !    write(1,'(A)') "plot '.\saida\FBA\FBA.txt' using 1:2 title 'Balística Ascendente' lw 2 \"    
    !    write(1,'(A)') "   , '.\saida\FBA\FBAa.txt' using 1:2 title 'Balística Ascendente - ajustada' lw 2 "   
    !close(1)
    !i = system("Wgnuplot ./saida/fasesBA.gnu")  
    !pause            
            
                call diff_sq_BA
                sumq = sum(D(np:nBA)**2)
                ha = res*pBA(1)
                pBA(1) = pBA(1) + ha/2    
                call diff_sq_BA
                erro = sumq - sum(D(np:nBA)**2)
                h = 1.*sign(ha, erro )
                pBA(1) = pBA(1) + h*3                               
    
    
                call diff_sq_BA
                sumq = sum(D(np:nBA)**2) 
                ha = res*pBA(2)
                pBA(2) = pBA(2) + ha/2    
                call diff_sq_BA
                erro = sumq - sum(D(np:nBA)**2)
                h = 1.*sign(ha, erro )
                pBA(2) = pBA(2) + h*3    
                
                
                call diff_sq_BA
                sumq = sum(D(np:nBA)**2) 
                ha = res*pBA(3)
                pBA(3) = pBA(3) + ha/2    
                call diff_sq_BA
                erro = sumq - sum(D(np:nBA)**2)
                h = 1.*sign(ha, erro )
                pBA(3) = pBA(3) + h*3    
                
                call diff_sq_BA
                sumq = sum(D(np:nBA)**2) 
                ha = res*pBA(4)
                pBA(4) = pBA(4) + ha/2    
                call diff_sq_BA
                erro = sumq - sum(D(np:nBA)**2)
                h = 1.*sign(ha, erro )
                pBA(4) = pBA(4) + h*3  
                
                
                   
                
                !write(*,'(i5, 10(1PE15.3))') j, pBA,sumq ;pause          
            !if(abs(pBA(1)*pBA(2)*pBA(3)*pBA(4)-pant) <= 1e-14) exit 
                SStot = sum( (Alt(nP:nBA,2) - sum(Alt(np:NBA,2) )/nP )**2)
                SSres = 0
                do i=nP, NBA
                    SSres = SSres + (Alt(i,2)- hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) )**2
                enddo
            write(1,'(I9,100(1PE25.15))') (k-1)*nint+j, pBA(1)*pBA(2)*pBA(3)*pBA(4), abs(pBA(1)*pBA(2)*pBA(3)*pBA(4)-pant), sumq, 1-SSres/SStot                
            pant = pBA(1)*pBA(2)*pBA(3)*pBA(4)
        enddo
        close(1)
        !write(*,'(i5, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    enddo
    write(*,'(i9, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    
    ! Cria gráfico das fases separadas   
    open(1,file='.\saida\FBA\FBAa.txt')
        dt = (tBA-tP)/Nm
        do i=0, Nm
            tempo = tP+dt*i
            write(1,'(10(1PE25.13))') tempo, hBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4))
        enddo        
    close(1) 
    
    open(1,file='.\saida\fasesBA.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Altitude (m)'"
        write(1,'(A)') "set out '.\saida\FBA\FBA.png'"
        write(1,'(A)') "plot '.\saida\FBA\FBA.txt' using 1:2 title 'Balística Ascendente' lw 2 \"    
        write(1,'(A)') "   , '.\saida\FBA\FBAa.txt' using 1:2 title 'Balística Ascendente - ajustada' lw 2 "   
    close(1)
    
    i = system("Wgnuplot ./saida/fasesBA.gnu")       
    
    ! Cria gráfico dos resíduos
    open(1,file='.\saida\FBA\FBAr.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'Iteracao'"
        write(1,'(A)') "set ylabel 'R²'"
        write(1,'(A)') "set out '.\saida\FBA\FBAr2.png'"
        k=1
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBA_",trim(adjustl(temp)),".txt"        
        write(1,'(A,A,A,I1,A)') "plot '.\saida\FBA\",trim(adjustl(temp)),"' using 1:5 title '",k+1,"' lw 2 \"            
        do k=2, nres-1
            write(temp,*) k+1
            write(temp,'(A,A,A)') "FBA_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "   , '.\saida\FBA\",trim(adjustl(temp)),"' using 1:5 title '",k+1,"' lw 2 \"
        enddo 
        k = nres
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBA_",trim(adjustl(temp)),".txt"
        write(1,'(A,A,A,I2,A)') "   , '.\saida\FBA\",trim(adjustl(temp)),"' using 1:5 title '",k+1,"' lw 2 "         
        
        !write(1,'(A)') "set format y '10**{%L}'"
        write(1,'(A)') "set ylabel 'Soma dos quadrados dos resíduos'"
        !write(1,'(A)') "set logscale y"                   
        do k=1, nres
            write(temp,*) k+1
            write(1,'(A)') "set out '.\saida\FBA\FBAr_" // trim(adjustl(temp)) //  ".png'" 
            write(temp,'(A,A,A)') "FBA_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "plot '.\saida\FBA\",trim(adjustl(temp)),"' using 1:4 title '",k+1,"' lw 2 "
        enddo   
        
        
    close(1)
    
    i = system("Wgnuplot ./saida/FBA/FBAr.gnu")       
    
    !pause
    ! Ajuste da fase balistica descendente
    write(*,*) "FBD"
    ene = nBD-nBA
    na  = 3
    allocate( pBD(na) )    
    ! Estimativa dos parâmetros
    pBD(1) = 1*rho*A/2/(mo-mej) ! Theta, Cd=1
    pBD(2) = tBA ! tH
    pBD(3) = hBA(g,tBA, pBA(1), pBA(2), pBA(3), pBA(4)) ! HH
    pant = pBD(1)*pBD(2)*pBD(3)
    do k=1, nres
        res = 10d0**(-1-k)
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBD_",trim(adjustl(temp)),".txt"
        open(1,file=".\saida\FBD\"//trim(adjustl(temp)) )          
        do j=1, nint
            
                call diff_sq_BD
                sumq = sum(D(nBA:nBD)**2) 
                ha = res*pBD(1)
                pBD(1) = pBD(1) + ha/2    
                call diff_sq_BD
                erro = sumq - sum(D(nBA:nBD)**2)
                h = 1.*sign(ha, erro )
                pBD(1) = pBD(1) + h*3   
                
                call diff_sq_BD
                sumq = sum(D(nBA:nBD)**2)
                ha = res*pBD(2)
                pBD(2) = pBD(2) + ha/2    
                call diff_sq_BD
                erro = sumq - sum(D(nBA:nBD)**2)
                h = 1.*sign(ha, erro )
                pBD(2) = pBD(2) + h*3  
                
                call diff_sq_BD
                sumq = sum(D(nBA:nBD)**2) 
                ha = res*pBD(3)
                pBD(3) = pBD(3) + ha/2    
                call diff_sq_BD
                erro = sumq - sum(D(nBA:nBD)**2)
                h = 1.*sign(ha, erro )
                pBD(3) = pBD(3) + h*3                  

                
                !write(*,'(i5, 10(1PE15.3))') j, pBD,sumq ;pause
            !if(abs(pBD(1)*pBD(2)*pBD(3)-pant) <= 1e-14) exit 
                SStot = sum( (Alt(nBA:nBD,2) - sum(Alt(nBA:NBD,2) )/nP )**2)
                SSres = 0
                do i=nBA, NBD
                    SSres = SSres + (Alt(i,2)- hBD(g,Alt(i,1), pBD(1), pBD(2), pBD(3)) )**2
                enddo
            write(1,'(I9,100(1PE25.15))') (k-1)*nint+j, pBD(1)*pBD(2)*pBD(3), abs(pBD(1)*pBD(2)*pBD(3)-pant), sumq, 1-SSres/SStot                                
                
            pant = pBD(1)*pBD(2)*pBD(3)
        enddo
        close(1)
        !write(*,'(i5, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    enddo
    write(*,'(i9, 10(1PE13.3))') k*nint, 1-SSres/SStot,sumq
    
    ! Cria gráfico das fases separadas   
    open(1,file='.\saida\FBD\FBDa.txt')
        dt = (tBD-tBA)/Nm
        do i=0, Nm
            tempo = tBA+dt*i
            write(1,'(10(1PE25.13))') tempo, hBD(g,tempo, pBD(1), pBD(2), pBD(3))
        enddo        
    close(1)

    open(1,file='.\saida\fasesBD.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Altitude (m)'"
        write(1,'(A)') "set out '.\saida\FBD\FBD.png'"
        write(1,'(A)') "plot '.\saida\FBD\FBD.txt' using 1:2 title 'Balística Descendente' lw 2 \"    
        write(1,'(A)') "   , '.\saida\FBD\FBDa.txt' using 1:2 title 'Balística Descendente - ajustada' lw 2 "   
    close(1)
    
    i = system("Wgnuplot ./saida/fasesBD.gnu")      
    
    ! Cria gráfico dos resíduos
    open(1,file='.\saida\FBD\FBDr.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'Iteracao'"
        write(1,'(A)') "set ylabel 'R²'"
        write(1,'(A)') "set out '.\saida\FBD\FBDr2.png'"
        k=1
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBD_",trim(adjustl(temp)),".txt"        
        write(1,'(A,A,A,I1,A)') "plot '.\saida\FBD\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 \"            
        do k=2, nres-1
            write(temp,*) k+1
            write(temp,'(A,A,A)') "FBD_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "   , '.\saida\FBD\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 \"
        enddo 
        k = nres
        write(temp,*) k+1
        write(temp,'(A,A,A)') "FBD_",trim(adjustl(temp)),".txt"
        write(1,'(A,A,A,I2,A)') "   , '.\saida\FBD\",trim(adjustl(temp)),"' using 1:5 title '",k,"' lw 2 "         
        
        !write(1,'(A)') "set format y '10**{%L}'"
        write(1,'(A)') "set ylabel 'Soma dos quadrados dos resíduos'"
        !write(1,'(A)') "set logscale y"                   
        do k=1, nres
            write(temp,*) k+1
            write(1,'(A)') "set out '.\saida\FBD\FBDr_" // trim(adjustl(temp)) //  ".png'"             
            write(temp,'(A,A,A)') "FBD_",trim(adjustl(temp)),".txt"
            write(1,'(A,A,A,I2,A)') "plot '.\saida\FBD\",trim(adjustl(temp)),"' using 1:4 title '",k,"' lw 2 "
        enddo    
        
        
    close(1)
    
    i = system("Wgnuplot ./saida/FBD/FBDr.gnu")      

    write(*,*) "Drag"
    allocate(Drag(0:Nm,2))
    Drag = 0
    ! Calcula a curva de força de arrasto na fase balística ascendente
    ! -Drag-P=m*a
    ! Drag = -mf*(g+a)
    dt = (tBA-tP)/Nm
    do j=0, Nm
        tempo = tP+dt*j     
        vel  = vBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4))
        acel = aBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4))
        Drag(j,1) = vel
        Drag(j,2) = -(mo-mej)*(g + acel  )

        !write(*,*) vBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4)), -(mo-mej)*(g+aBA(g,tempo, pBA(1), pBA(2), pBA(3), pBA(4))  )
    enddo
    open(2,file='.\saida\Drag.txt')
        do i=0, Nm
            tempo = tP+dt*i
            write(2,'(10(1PE25.13))') tempo,  Drag(i,1), Drag(i,2), Drag(i,2)/rho/A*2/Drag(i,1)/Drag(i,1)
        enddo        
    close(2)
    open(1,file='.\saida\Drag.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'Velocidade (m/s)'"
        write(1,'(A)') "set ylabel 'Força de Arrasto (N)'"
        write(1,'(A)') "set out '.\saida\Drag.png'"
        write(1,'(A)') "plot '.\saida\Drag.txt' using 2:3 title '' lw 2 "     
        !write(1,'(A)') "set yrange [0:10]"
        write(1,'(A)') "set ylabel 'Coeficiente de arrasto'"
        write(1,'(A)') "set out '.\saida\Cd.png'"
        write(1,'(A)') "plot '.\saida\Drag.txt' using 2:4 title '' lw 2 "           
    close(1)
    
    i = system("Wgnuplot ./saida/Drag.gnu") 
    ! Calcula a força de empuxo com:
    ! E - Drag - P = m*a
    ! E = Drag + m*(a+g)
    write(*,*) "Thrust"    
    open(1,file='.\saida\Forca.txt')    
    dt = (tP-0)/Nm
    do i=0, Nm
        tempo = 0+dt*i
        vel = vP(tempo, pFp(1),pFp(2))
        Fd=0
        do j=1, Nm
            ! Interpolação linear para a força de arrasto dado a velocidade
            if( vel< Drag(j-1,1) .and. vel>=Drag(j,1)) then
                Fd = (Drag(j,2)-Drag(j-1,2))/(Drag(j,1)-Drag(j-1,1))*(vel-Drag(j-1,1)) + Drag(j-1,2)                
                exit
            endif            
        enddo
        !if(mod(i,100)==0) then
        !    write(*,*) tempo, vel, Fd, (mo-mej/2)*(aP(tempo, pFp(1),pFp(2))+g); pause
        !endif

        Forca(i) = Fd +( aP(tempo, pFp(1),pFp(2))+g )*(mo-mej/2)!((mo-mej) + mej*(tP-tempo)) 
        ForcaSemArrasto(i) = ( aP(tempo, pFp(1),pFp(2))+g )*(mo-mej/2)!((mo-mej) + mej*(tP-tempo))
        

        write(1,'(10(1PE25.13))') tempo,  Forca(i), ForcaSemArrasto(i), Fd, vel             
    enddo                
    close(1)
    
    open(1,file='.\saida\Forca.gnu')
        write(1,'(A)') "set term png nocrop enhanced size 1280, 720"
        write(1,'(A)') "set style data lines"
        write(1,'(A)') "set key right box top"
        write(1,'(A)') "set time"
        write(1,'(A)') "set grid"
        write(1,'(A)') "set title '' "
        write(1,'(A)') "set xlabel 'tempo (s)'"
        write(1,'(A)') "set ylabel 'Forças (N)'"
        write(1,'(A)') "set out '.\saida\Forca.png'"
        write(1,'(A)') "plot '.\saida\Forca.txt' using 1:2 title 'Empuxo (com arrasto)' lw 2  \"     
        write(1,'(A)') ", '.\saida\Forca.txt' using 1:3 title 'Empuxo (sem arrasto)' lw 2  \"  
        write(1,'(A)') ", '.\saida\Forca.txt' using 1:4 title 'Arrasto' lw 2  "  
    close(1)
    
    i = system("Wgnuplot ./saida/Forca.gnu")  
    write(*,*) "It"
    ! Calculo do Impulso total
    dt = (tP-0)/Nm
    It = 0
    ItSA=0
    do i=1, Nm
        It = It + dt*(Forca(i) + Forca(i-1))/2
        ItSA = ItSA+dt*(ForcaSemArrasto(i) + ForcaSemArrasto(i-1))/2
    enddo  
    
    t2 = timef()
    open(1,file='.\saida\resultados_' // trim(adjustl(arquivo)))
    
        write(1,'(/A)') " --- Dados de entrada ---"
        write(1,'(A,A)') "'"//trim(adjustl(traj))//"'", " = nome do arquivo de trajetória utilizado"
        write(1,'(1PE23.13,A)') tP, " = tempo até o fim da fase propulsada (s)"
        write(1,'(1PE23.13,A)') tBA, " = tempo até o fim da fase balística ascendente (s)"
        write(1,'(1PE23.13,A)') tBD, " = tempo até o fim da fase balística descendente (s)"
        write(1,'(1PE23.13,A)') De, " = diâmetro de referência (externo) do minifoguete (m)"
        write(1,'(1PE23.13,A)') Pa, " = pressão atmosférica no momento do lançamento (Pa)"
        write(1,'(1PE23.13,A)') T, " = temperatura atmosférica no momento do lançamento (K)"
        write(1,'(1PE23.13,A)') mo, " = massa inicial do minifoguete (kg)"
        write(1,'(1PE23.13,A)') mej, " = massa ejetada do minifoguete (kg)"
        write(1,'(1PE23.13,A)') g, " = aceleração da gravidade utilizad (m/s**2)"
        write(1,'(I23,A/)') Nm, " = Número de pontos para a criação dos gráficos das fases"
        
    
        write(1,'(/A)') " --- Parâmetros das curvas ajustadas ---"
        write(1,'(/A)') ":::: Fase Propulsada ::::"
        write(1,'( 1PE23.13,A)') pFp(1)," = Theta = Cd*rho*A/2/m [1/m]"
        write(1,'( 1PE23.13,A)') pFp(2)," = B = (Fmed-g*m)/m [m/s**2]"  
        write(1,'( 1PE23.13,A)') (pFp(2)+g)*(mo-mej/2)," = Fmed, considerando a massa média ao longo da fase propulsada [N]"        
        write(1,'( 1PE23.13,A)') pFp(1)*2*(mo-mej/2)/rho/A," = Cd, considerando a massa média ao longo da fase propulsada [adim]"

        
        write(1,'(/A)') ":::: Fase Balística Ascendente ::::"
        write(1,'( 1PE23.13,A)') pBA(1)," = Theta = Cd*rho*A/2/m [1/m]"        
        write(1,'( 1PE23.13,A)') pBA(2)," = tq [s]"
        write(1,'( 1PE23.13,A)') pBA(3)," = hq [m]"
        write(1,'( 1PE23.13,A)') pBA(4)," = vq [m/s]" 
        write(1,'( 1PE23.13,A)') pBA(1)*2*(mo-mej)/rho/A," = Cd [adim]"        
        
        
        write(1,'(/A)') ":::: Fase Balístia Descendente ::::" 
        write(1,'( 1PE23.13,A)') pBD(1)," = Theta = Cd*rho*A/2/m [1/m]" 
        write(1,'( 1PE23.13,A)') pBD(2)," = tH [s]"
        write(1,'( 1PE23.13,A)') pBD(3)," = hH [m]"                
        write(1,'( 1PE23.13,A)') pBD(1)*2*(mo-mej)/rho/A," = Cd [adim]"         
        
        write(1,*)  
        
        write(1,'(/A)') " --- Resultados ---"
        write(1,'(/A)') ":::: Força de empuxo (N) ::::" 
        
        write(1,'( 1PE23.13,A)') ItSA/tP," = Fmed (Sem força de Arrasto) "
        write(1,'( 1PE23.13,A)') It/tP ," = Fmed (Com força de Arrasto) "
        
        write(1,'(/A)') ":::: Impulso total (N.s) ::::" 
        write(1,'( 1PE23.13,A)') ItSA," = It (Sem força de Arrasto) "
        write(1,'( 1PE23.13,A)') It ," = It (Com força de Arrasto) "
        
        write(1,'(/A)') ":::: Velocidade de ejeção (m/s) ::::" 
        
        write(1,'(/1PE23.13,A )') ItSA/mej ," = c (Sem força de Arrasto) (m/s)"        
        write(1,'( 1PE23.13,A )') It/mej ," = c (Com força de Arrasto) (m/s)"

        write(1,'(/A)') ":::: Coeficiente de arrasto ao longo da fase balística ascendente (adim) ::::"         
        write(1,'(/1PE23.13,A )') sum(Drag(0:nm,2)/rho/A*2/Drag(0:nm,1)/Drag(0:nm,1))/nm ," = Cd"
        
        write(1,'(/A)') ":::: Outros ::::"  
        write(1,'(/1PE23.13,A )') t2-t1," = Tempo computacional (s)"        
    close(1)
    
!---------------------------------------    
    contains    
!---------------------------------------
real*8 function hP(t, theta,B)    
    real*8 t, theta, B    
    hP = 0.5d0/theta*log(cosh(sqrt(theta*B)*t)**2)    

end function hP
!---------------------------------------
real*8 function dhPdtheta(t, theta,B)    
    real*8 t, theta, B    
    dhPdtheta = 0.5d0*(sinh(sqrt(Theta*B)*t)*t*B*Theta-&
        &log(cosh(sqrt(Theta*B)*t)**2)*cosh(sqrt(Theta*B)*t)*sqrt(Theta*B))/&
        &cosh(sqrt(Theta*B)*t)/sqrt(Theta*B)/Theta**2

end function dhPdtheta

!---------------------------------------
real*8 function dhPdB(t, theta,B)    
    real*8 t, theta, B    
    dhPdB = 0.5d0*cosh(sqrt(Theta*B)*t)*sinh(sqrt(Theta*B)*t)/sqrt(Theta*B)*t

end function dhPdB
!---------------------------------------
real*8 function vP(t, theta,B)    
    real*8 t, theta, B    
    vP = 1d0/theta/cosh(sqrt(theta*B)*t)*sinh(sqrt(theta*B)*t)*sqrt(theta*B)

end function vP
!---------------------------------------
real*8 function aP(t, theta,B)    
    real*8 t, theta, B    
    aP = -1d0/cosh(sqrt(theta*B)*t)**2*sinh(sqrt(theta*B)*t)**2*B+B

end function aP
!---------------------------------------
real*8 function hBA(g, t, theta, tq, hq, vq)    
    real*8 g, t, theta, hq, tq, vq
    
hBA= 0.5d0/theta*log( exp(2*hq*theta)/g/(sin(tq*sqrt(g*theta)))**2*(sin(sqrt(g*theta)*t-tq*sqrt(g*theta))*&
    (sqrt(g)*cos(tq*sqrt(g*theta))-vq*sqrt(theta)*sin(tq*sqrt(g*theta)))-sqrt(g)*sin(sqrt(g*theta)*t))**2 )
end function hBA
!---------------------------------------
real*8 function hBD(g, t, theta, tH,HH)    
    real*8 g, t, theta, tH, HH
    
hBD=0.5d0*(2*sqrt(g)*sqrt(Theta)*t+2*log(2d0)-log(exp(-2*sqrt(g)*sqrt(Theta)*tH-2*HH*Theta)*(exp(2*sqrt(g)*sqrt(Theta)*t)+exp(2*sqrt(g)*sqrt(Theta)*tH))**2))/Theta
end function hBD
!---------------------------------------
real*8 function aBA(g, t, theta, tq, hq, vq)    
implicit none
    real*8 g, t, theta, hq, tq, vq
    
aBA=(-2*g**2-2*g*Theta*vq**2)/&
    (-Theta*vq**2*cos(2*sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+&
    g*cos(2*sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+2*sqrt(g)*sqrt(Theta)*vq*sin(2*sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+g+Theta*vq**2)
end function aBA
!---------------------------------------
real*8 function vBA(g, t, theta, tq, hq, vq)
implicit none
    real*8 g, t, theta, hq, tq, vq
    
vBA =-(-g*cos(sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+&
    cos(sqrt(g)*sqrt(Theta)*t)*g-sqrt(Theta)*sqrt(g)*vq*sin(sqrt(g)*&
    sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+sin(sqrt(g)*sqrt(Theta)*t)*&
    sqrt(Theta)*sqrt(g)*vq)/(sqrt(Theta)*sqrt(g)*sin(sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)-&
    sin(sqrt(g)*sqrt(Theta)*t)*sqrt(Theta)*sqrt(g)-vq*Theta*cos(sqrt(g)*sqrt(Theta)*t-2*sqrt(g)*sqrt(Theta)*tq)+cos(sqrt(g)*sqrt(Theta)*t)*vq*Theta)
end function vBA
!---------------------------------------
subroutine diff_sq_FP
implicit none
        do i=1, np
            ! Cálculo de D (Diferença entre as variáveis medidas e os valores da função)
            D(i) = Alt(i,2) - hP(Alt(i,1), pFp(1),pFp(2))
        enddo                    

end subroutine diff_sq_FP
!---------------------------------------
subroutine diff_sq_BA
implicit none
        do i=Np, NBA
            ! Cálculo de D (Diferença entre as variáveis medidas e os valores da função)
            D(i) = Alt(i,2) - hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4))
        enddo
end subroutine diff_sq_BA
!---------------------------------------
subroutine diff_sq_BD
implicit none
        do i=NBA, NBD
            ! Cálculo de D (Diferença entre as variáveis medidas e os valores da função)
            D(i) = Alt(i,2) - hBD(g,Alt(i,1), pBD(1), pBD(2), pBD(3))
        enddo
end subroutine diff_sq_BD
!---------------------------------------  
subroutine ajusteFP
implicit none
integer i, j
real*8 :: somaN(2), somaD(2)
pFp(1) = .1
pFp(2) = 100
pant = pFp(1)*pFp(2)
write(temp,'(A,A,A)') "FP_0.txt"
open(1,file=".\saida\FP\"//trim(adjustl(temp)) ) 
do i=1, Nint*10   
    somaN= 0
    somaD= 0
    do j=1, np
        somaN(1) = somaN(1) +                    Alt(j,2)*dhPdtheta(Alt(j,1), pFp(1),pFp(2))
        somaD(1) = somaD(1) + hP(Alt(j,1), pFp(1),pFp(2))*dhPdtheta(Alt(j,1), pFp(1),pFp(2))
        
        somaN(2) = somaN(2) +                    Alt(j,2)*dhPdB(Alt(j,1), pFp(1),pFp(2))
        somaD(2) = somaD(2) + hP(Alt(j,1), pFp(1),pFp(2))*dhPdB(Alt(j,1), pFp(1),pFp(2))                    
    enddo
    
    pFp(1) = pFp(1)*somaN(1)/somaD(1)
    pFp(2) = pFp(2)*somaN(2)/somaD(2)

    SStot = sum( (Alt(1:nP,2) - sum(Alt(1:nP,2) )/nP )**2)
    SSres = 0
    do j=1, nP
        SSres = SSres + (Alt(j,2)-hP(Alt(j,1), pFp(1),pFp(2)) )**2
    enddo 

    call diff_sq_FP
    sumq = sum(D(1:np)**2)    
    
    write(1,'(I9,100(1PE25.15))') i, pFp(1)*pFp(2), abs(pFp(1)*pFp(2)-pant), sumq, 1-SSres/SStot
    if(abs(pFp(1)*pFp(2)-pant)<1E-15) exit
    pant = pFp(1)*pFp(2)
enddo
close(1)
write(*,'(i9, 10(1PE13.3))') i, 1-SSres/SStot,pFp
end subroutine ajusteFP
!---------------------------------------  
subroutine ajusteFBA
implicit none
integer i, j
real*8 :: somaN(4), somaD(4)
real*8 :: re,dt
!pBA(1) = .1
!pBA(2) = 100
!pBA(3) = 100
!pBA(4) = 100
pant = pBA(1)*pBA(2)*pBA(3)*pBA(4)
write(temp,'(A,A,A)') "FBA_0.txt"
open(1,file=".\saida\FBA\"//trim(adjustl(temp)) ) 

re = 1d-10

do i=1, Nint*10   
    somaN= 0
    somaD= 0
    do j=1, np
        dt = re*pBA(1)
        somaN(1) = somaN(1) +                                        Alt(j,2)*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1)-dt, pBA(2), pBA(3), pBA(4)))/dt
        somaD(1) = somaD(1) + hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4))*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1)-dt, pBA(2), pBA(3), pBA(4)))/dt
        
        dt = re*pBA(2)
        somaN(2) = somaN(2) +                                        Alt(j,2)*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2)-dt, pBA(3), pBA(4)))/dt
        somaD(2) = somaD(2) + hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4))*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2)-dt, pBA(3), pBA(4)))/dt
        
        dt = re*pBA(3)
        somaN(3) = somaN(3) +                                        Alt(j,2)*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3)-dt, pBA(4)))/dt
        somaD(3) = somaD(3) + hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4))*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3)-dt, pBA(4)))/dt
        
        dt = re*pBA(4)
        somaN(4) = somaN(4) +                                        Alt(j,2)*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)-dt))/dt
        somaD(4) = somaD(4) + hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4))*(hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)) - hBA(g,Alt(i,1), pBA(1), pBA(2), pBA(3), pBA(4)-dt))/dt
        
    enddo
    write(*,*) SomaN/SomaD
    stop
    pBA(1) = pBA(1)*somaN(1)/somaD(1)
    pBA(2) = pBA(2)*somaN(2)/somaD(2)
    pBA(3) = pBA(3)*somaN(3)/somaD(3)
    pBA(4) = pBA(4)*somaN(4)/somaD(4)

    SStot = sum( (Alt(nP:nBA,2) - sum(Alt(np:NBA,2) )/nP )**2)
    SSres = 0
    do j=nP, NBA
        SSres = SSres + (Alt(j,2)- hBA(g,Alt(j,1), pBA(1), pBA(2), pBA(3), pBA(4)) )**2
    enddo

    call diff_sq_BA
    sumq = sum(D(Np:nBA)**2)     
    
    write(*,'(i9, 10(1PE13.3))') i, 1-SSres/SStot,pBA;pause
    
    write(1,'(I9,100(1PE25.15))') i, pBA(1)*pBA(2)*pBA(3)*pBA(4), abs(pBA(1)*pBA(2)*pBA(3)*pBA(4)-pant), sumq, 1-SSres/SStot                
    if(abs(pBA(1)*pBA(2)*pBA(3)*pBA(4)-pant)<1E-15) exit
    pant = pBA(1)*pBA(2)*pBA(3)*pBA(4)
enddo
close(1)
write(*,'(i9, 10(1PE13.3))') i, 1-SSres/SStot,pBA
end subroutine ajusteFBA
!---------------------------------------  
end program voo
    