program trajeto_analitico
! Programador: Diego Moro 
! Versão 1p0: em 03/10/2016    
! Versão 1p1: em 04/10/2016 - Obtenção de tH, hH, tI e vI de forma numérica
! Versão 1p2: em 06/10/2016 - Opção de fase propulsada considerando força peso, arrasto, 
!    no entanto é considerado massa constante (média ao longo da queima)
! Resolve as 3 fases do voo de um minifoguete por relações analíticas, compreendendo:
! Fase propulsada: arrasto desprezado, apenas força do motor e força peso
! Fase balística ascendente: arrasto e força peso considerados, solução continua é obtida
! Fase balística descendente: arrasto e força peso considerados, solução continua é obtida
    
implicit none

real*8 :: pi
real*8 :: Mf, Mp, dr, tq, Emed, g
real*8 :: Cd, par, Tar, Rar
real*8 :: rho, A, theta
real*8 :: R, hq, Vq
real*8 :: tH, Hh
real*8 :: dtp, dtba, dtbd
real*8 :: t, tI, vI

real*8 :: mref, thetaref


integer :: Ip, efp
integer :: i, j
character*128 :: arquivo, caso, titulo

pi = dacos(-1d0)

! Dados Constantes


call solucao

contains
!--------------------------------------------------------------
subroutine solucao
use portlib

! Abre arquivo que direciona ao arquivo de dados
i = system("Notepad ./entrada/trajeto.in")
open(1,file='.\entrada\trajeto.in')
    read(1,*) arquivo
close(1)

i = system("Notepad ./entrada/" // trim(adjustl(arquivo)) )

! Abre arquivo de dados
open(1,file=".\entrada\"//arquivo)
    read(1,*) caso
    read(1,*) Mf ! em kg
    read(1,*) Mp ! em kg
    read(1,*) dr ! em m
    read(1,*) Emed ! em N
    read(1,*) tq ! em s
    read(1,*) Cd
	read(1,*) g
    read(1,*) par ! em Pa
    read(1,*) Tar ! Em °C
	read(1,*) Rar
    read(1,*) Ip
    read(1,*) efp
close(1)

! Adequando as variaveis
Tar = Tar + 273.15d0  ! em K

! Calculo da massa específica do ar
rho = par/Tar/Rar

! Calcula da area do MF
A = 0.25d0*pi*dr**2

! Calcula da variavel theta
theta = 0.5d0*Cd*rho*A/Mf

! Calculo das variaveis relacionadas a fase propulsada

open(1,file='.\resultados\'//trim(adjustl(caso)) // '_trj.out')
write(1,'("#",5(A25))') "t(s)", "h(m)", "v(m/s)", "a(m/s²)"
    select case (efp)
        case(0) ! Considerando gravidade e massa variavel
            ! Variavel auxiliar
            R = Mp/(Mf+Mp)
            ! Altura no final da queima do motor
            hq = -0.5d0*tq*tq/mp/mp*(2*Emed*mf*log((tq*(mf+mp))/(mf*tq))-2*Emed*mp+g*mp*mp)     
        
            ! Velocidade no final da queima  do motor
            vq = -tq/mp*(g*mp+Emed*(log((mf*tq)/(tq*(mf+mp)) ) )) 
        
        
            dtp = (tq-0)/Ip
            do j=0, Ip
                t = j*dtp
                write(1,'(20(1PE25.15))') t, hp(Emed, tq, mf, mp, t), vp(Emed, tq, mf, mp, t), ap(Emed, tq, mf, mp, t)
            enddo         
    
        case(1) ! Considerando gravidade e força de arrasto mas com massa constante (média entre o início e fim da queima)
            ! Massa de referencia para o calculo da fase propulsada (desconsiderando a massa variavel)
            mref = mf + mp/2
            thetaref = 0.5d0*Cd*rho*A/Mref
        
            ! Altura no final da queima do motor
            hq = hp2(Emed, mref, Thetaref, tq)        
        
            ! Velocidade no final da queima  do motor
            vq = vp2(Emed, mref, Thetaref, tq) 
        
            dtp = (tq-0)/Ip
            do j=0, Ip
                t = j*dtp
                write(1,'(20(1PE25.15))') t, hp2(Emed, mref, Thetaref, t),vp2(Emed, mref, Thetaref, t),ap2(Emed, mref, Thetaref, t)
            enddo            
    end select
close(1)

!--------------------------------------
! Início da fase balística ascendente
!--------------------------------------
   
    
!! Tempo até o apogeu
!tH = -atan((-cos(sqrt(g)*sqrt(Theta)*tq)**2*sqrt(Theta*exp(2*hq*Theta)*g)*sqrt(g)+cos(sqrt(g)*sqrt(Theta)*tq)*sqrt(Theta*exp(2*hq*Theta)*g)*vq*sqrt(Theta)*sin(sqrt(g)*sqrt(Theta)*tq)+sqrt(Theta*cos(sqrt(g)*sqrt(Theta)*tq)**4*exp(2*hq*Theta)*g-2*Theta**(3/2)*cos(sqrt(g)*sqrt(Theta)*tq)**3*exp(2*hq*Theta)*sqrt(g)*vq*sin(sqrt(g)*sqrt(Theta)*tq)+Theta**2*cos(sqrt(g)*sqrt(Theta)*tq)**2*exp(2*hq*Theta)*vq**2*sin(sqrt(g)*sqrt(Theta)*tq)**2+Theta*exp(2*hq*Theta)*g-Theta*exp(2*hq*Theta)*g*cos(sqrt(g)*sqrt(Theta)*tq)**2+2*Theta**(3/2)*exp(2*hq*Theta)*sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)*vq*sin(sqrt(g)*sqrt(Theta)*tq)-Theta**2*exp(2*hq*Theta)*vq**2*sin(sqrt(g)*sqrt(Theta)*tq)**2+Theta*exp(2*hq*Theta)*sin(sqrt(g)*sqrt(Theta)*tq)**2*g*cos(sqrt(g)*sqrt(Theta)*tq)**2-2*Theta**(3/2)*exp(2*hq*Theta)*sin(sqrt(g)*sqrt(Theta)*tq)**3*sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)*vq+Theta**2*exp(2*hq*Theta)*sin(sqrt(g)*sqrt(Theta)*tq)**4*vq**2)*sqrt(g))/(Theta*exp(2*hq*Theta)*g)**(1/2)/sin(sqrt(g)*sqrt(Theta)*tq)/(-sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)+vq*sqrt(Theta)*sin(sqrt(g)*sqrt(Theta)*tq)))/sqrt(g*Theta)
!
!! Altura do apogeu
!hH = hBA(hq, theta, g, tq, vq, tH)

open(1,file='.\resultados\'//trim(adjustl(caso)) // '_trj.out', position='append')
do j=0, 1000000
    t = tq + j*dtp
    if(vBA(hq, theta, g, tq, vq, t) < 0) then
        tH = (t - dtp) + (dtp)/(vBA(hq, theta, g, tq, vq, t)-vBA(hq, theta, g, tq, vq, t-dtp))*(0-vBA(hq, theta, g, tq, vq, t-dtp))
        hH = hBA(hq, theta, g, tq, vq, t - dtp) + (hBA(hq, theta, g, tq, vq, t)-hBA(hq, theta, g, tq, vq, t-dtp))/(vBA(hq, theta, g, tq, vq, t)-vBA(hq, theta, g, tq, vq, t-dtp))*(0-vBA(hq, theta, g, tq, vq, t-dtp))
        write(1,'(20(1PE25.15))') tH, hH, vBA(hq, theta, g, tq, vq, tH), aBA(hq, theta, g, tq, vq, tH)
        exit
    endif
    write(1,'(20(1PE25.15))') t, hBA(hq, theta, g, tq, vq, t), vBA(hq, theta, g, tq, vq, t), aBA(hq, theta, g, tq, vq, t)
enddo
close(1)


!--------------------------------------
! Início da fase balística descendente
!--------------------------------------
! Tempo até o impacto
!tI = (sqrt(g*Theta)*tH+log(exp(hH*Theta)+sqrt(exp(2*hH*Theta)-1)))/sqrt(g*Theta)
!
!vI = -(exp(2*hH*Theta)+exp(hH*Theta)*sqrt(exp(2*hH*Theta)-1)-1)*sqrt(g)/sqrt(Theta)*exp(-hH*Theta)/(exp(hH*Theta)+sqrt(exp(2*hH*Theta)-1))

!dtbd = (tI-tH)/Ibd
open(1,file='.\resultados\'//trim(adjustl(caso)) // '_trj.out', position='append')
do j=0, 1000000
    t = tH + j*dtp
    if(hBD( theta, g, tH, hH, t) < 0) then
        tI = (t - dtp) + (dtp)/(hBD( theta, g, tH, hH, t)-hBD( theta, g, tH, hH, t-dtp))*(0-hBD( theta, g, tH, hH, t-dtp))
        vI = vBD( theta, g, tH, hH, t - dtp) + (vBD( theta, g, tH, hH, t)-vBD( theta, g, tH, hH, t-dtp))/(hBD( theta, g, tH, hH, t)-hBD( theta, g, tH, hH, t-dtp))*(0-hBD( theta, g, tH, hH, t-dtp))
        write(1,'(20(1PE25.15))') t, 0d0, vI, aBD( theta, g, tH, hH, t)
        exit
    endif    
    write(1,'(20(1PE25.15))') t, hBD( theta, g, tH, hH, t), vBD( theta, g, tH, hH, t), aBD( theta, g, tH, hH, t)
enddo
close(1)

!--------------------------------------
! Cria grafico da trajeto

open(1,file='.\resultados\'//'graf.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 'Altura (m)'"
write(1,'(A)') "set out '" // '.\resultados\'//trim(adjustl(caso)) // "_H.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:2 title '' lw 2"
write(1,'(A)') "set ylabel 'Velocidade (km/h)'"
write(1,'(A)') "set out '" //'.\resultados\'// trim(adjustl(caso)) // "_v_kmh.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($3*3.6) title '' lw 2"
write(1,'(A)') "set ylabel 'Aceleração (gee)'"
write(1,'(A)') "set out '" //'.\resultados\'// trim(adjustl(caso)) // "_a_gee.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($4/9.8065) title '' lw 2"

write(1,'(A)') "set ytics nomirror"
write(1,'(A)') "set y2tics nomirror"
write(1,'(A)') "set ylabel 'Altitude (m)'"
write(1,'(A)') "set y2label 'Velocidade (km/h)'"
!write(1,'(A)') "set xrange[0:14] "
!write(1,'(A)') "set yrange[0:130] "
!write(1,'(A)') "set y2range[-150:450] "
write(1,'(A)') "set out '" //'.\resultados\'// trim(adjustl(caso)) // "_HV.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:2 title 'Altitude (m)' axes x1y1 lw 3 \"
write(1,'(A)') ", '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($3*3.6) title 'Velocidade (km/h)' axes x1y2 lw 2 "

write(1,'(A)') "set ylabel 'Altitude (m)'"
write(1,'(A)') "set y2label 'Aceleração (gee)'"
!write(1,'(A)') "set xrange[0:14] "
!write(1,'(A)') "set yrange[0:130] "
!write(1,'(A)') "set y2range[-30:30] "
write(1,'(A)') "set out '" //'.\resultados\'// trim(adjustl(caso)) // "_HA.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:2 title 'Altitude (m)' axes x1y1 lw 3 \"
write(1,'(A)') ", '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($4/9.8065) title 'Aceleração (gee)' axes x1y2 lw 2"

write(1,'(A)') "set ylabel 'Velocidade (km/h)'"
write(1,'(A)') "set y2label 'Aceleração (gee)'"
!write(1,'(A)') "set xrange[0:14] "
!write(1,'(A)') "set yrange[-150:450] "
!write(1,'(A)') "set y2range[-30:30] "
write(1,'(A)') "set out '" //'.\resultados\'// trim(adjustl(caso)) // "_VA.png'"
write(1,'(A)') "plot '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($3*3.6) title 'Velocidade (km/h)' axes x1y1 lw 3 \"
write(1,'(A)') ", '"//'.\resultados\'//trim(adjustl(caso)) // "_trj.out' using 1:($4/9.8065) title 'Aceleração (gee)' axes x1y2 lw 2"
close(1)
i = system(".\wgnuplot\Wgnuplot ./resultados/graf.gnu")

! Cria arquivo de resultados
open(1,file='.\resultados\'//trim(adjustl(caso)) // '_saida.out')
    write(1,'(A,A)') "'" // trim(adjustl(caso)), "' caso  = nome do arquivo de saída"
    write(1,'(1PE25.15,A)') Mf, "  ..... Mf  = Massa final (kg)"
    write(1,'(1PE25.15,A)') Mp, "  ..... Mp  = Massa de propelente (kg)"
    write(1,'(1PE25.15,A)') dr, "  ..... dr  = Diâmetro de referência (m)" 
    write(1,'(1PE25.15,A)') Emed, "  ..... Ec  = Empuxo constante do motor (N)"
    write(1,'(1PE25.15,A)') tq, "  ..... tq  = Tempo de queima do propelente (s) "
    write(1,'(1PE25.15,A)') Cd, "  ..... Cd  = Coeficiente de arrasto (adim.)"
    write(1,'(1PE25.15,A)') par, "  ..... par = Pressão atmosférica local (Pa) "
    write(1,'(1PE25.15,A)') Tar-273.15d0, "  ..... Tar = Temperatura atmosférica local (oC)"
    write(1,'(I25,A)') Ip, "  ..... Ip  = número de pontos para a fase propulsada"
    write(1,'(I25,A)') efp, "  ..... efp = escolha da hipotese da fase propulsada (0=gravidade+massa variavel), (1=gravidade+arrasto+massa fixa)"
    
    write(1,'(/A)') "------------------------------------"
    write(1,'( A )') "------------RESULTADOS--------------"
    write(1,'(A/)') "------------------------------------"

    write(1,'(1PE25.15,A)') g, "  ...... g = Aceleração da gravidade (m/s²)" 
    write(1,'(1PE25.15,A)') rho, "  ..... rho = Massa específica do ar (kg/m³)"    

    write(1,'(/A )') "----------FASE PROPULSADA-----------"  
    
    if(efp == 0) write(1,'( A/)') " (considerando aceleração da gravidade e massa variável)"
    if(efp == 1) write(1,'( A/)') " (considerando aceleração da gravidade, força de arrasto e massa média da queima)"
    write(1,'(F25.13,A)') hq, "  ..... hq  = altura no final da queima (m)"
    write(1,'(F25.13,A)') vq, "  ..... vq  = velocidade no final da queima (m/s)"
    write(1,'(F25.13,A)') vq*3.6d0, "  ..... vq  = velocidade no final da queima (km/h)"
    
    write(1,'(/A )') "-----FASE BALÍSTICA ASCENDENTE------"  
    
    write(1,'( A/)') " (admitindo força peso e força de arrasto)"
    write(1,'(F25.13,A)') tH, "  ..... tH  = tempo até o apogeu (s)"
    write(1,'(F25.13,A)') hH, "  ..... hH  = altura do apogeu (m)"
    
    write(1,'(/A )') "-----FASE BALÍSTICA DESCENDENTE-----"  
    
    write(1,'( A/)') " (admitindo força peso e força de arrasto)"
    write(1,'(F25.13,A)') tI, "  ..... tI  = tempo até o impacto (s)"
    write(1,'(F25.13,A)') vI, "  ..... vI  = velocidade no impacto (m/s)"    
    write(1,'(F25.13,A)') vI*3.6d0, "  ..... vI  = velocidade no impacto (km/h)"    
    
    
close(1)

i = system("Notepad " // './resultados/'//trim(adjustl(caso)) // '_saida.out')

end subroutine solucao
!--------------------------------------------------------------
real*8 function hp(Emed, tq, mf, mp, t)

real*8 Emed, tq, mf, mp, t

hp= Emed*tq*tq/mp*(mf/mp+1-t/tq)*log(1-mp/tq/(mf+mp)*t)-0.5d0*g*t**2+Emed*tq/mp*t

end function hp
!--------------------------------------------------------------
real*8 function vp(Emed, tq, mf, mp, t)

real*8 Emed, tq, mf, mp, t

vp = -Emed*tq/mp*log(1-mp*t/tq/(mf+mp))-g*t

end function vp
!--------------------------------------------------------------
real*8 function ap(Emed, tq, mf, mp, t)

real*8 Emed, tq, mf, mp, t

ap = (tq*(Emed-g*(mf+mp))+g*t*mp)/(tq*(mf+mp)-mp*t)
end function ap
!--------------------------------------------------------------
real*8 function hp2(Emed, m, Theta, t)

real*8 Emed, m, Theta, t


hp2 = 0.5d0*log(cosh(sqrt(theta*(Emed-g*m)/m)*t)**2)/theta

end function hp2
!--------------------------------------------------------------
real*8 function vp2(Emed, m, Theta, t)

real*8 Emed, m, Theta, t

vp2 = sqrt((Emed-g*m)/m/theta)*tanh(sqrt(theta*(Emed-g*m)/m)*t)

end function vp2
!--------------------------------------------------------------
real*8 function ap2(Emed, m, Theta, t)

real*8 Emed, m, Theta, t


ap2 = (Emed-g*m)/m/(cosh(sqrt(theta*(Emed-g*m)/m)*t)**2)

end function ap2
!--------------------------------------------------------------
real*8 function hBA(hq, theta, g, tq, vq, t)

real*8 hq, theta, g, tq, vq, t

hBA = 0.5d0*log(-exp(2*hq*Theta)*(sin(sqrt(g)*sqrt(Theta)*t)*sqrt(g)*&
    cos(sqrt(g)*sqrt(Theta)*tq)**2-sin(sqrt(g)*sqrt(Theta)*t)*cos(sqrt(g)*&
    sqrt(Theta)*tq)*sin(sqrt(g)*sqrt(Theta)*tq)*vq*sqrt(Theta)-sin(sqrt(g)*&
    sqrt(Theta)*t)*sqrt(g)-cos(sqrt(g)*sqrt(Theta)*t)*sin(sqrt(g)*sqrt(Theta)*tq)*&
    sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)+cos(sqrt(g)*sqrt(Theta)*t)*sin(sqrt(g)*&
    sqrt(Theta)*tq)**2*vq*sqrt(Theta))**2/g/(cos(sqrt(g)*sqrt(Theta)*tq)-1)/(cos(sqrt(g)*sqrt(Theta)*tq)+1))/Theta

end function hBA
!--------------------------------------------------------------
real*8 function vBA(hq, theta, g, tq, vq, t)

real*8 hq, theta, g, tq, vq, t

vBA = (cos(sqrt(g)*sqrt(Theta)*t)*g*cos(sqrt(g)*sqrt(Theta)*tq)**2-&
    sqrt(g)*cos(sqrt(g)*sqrt(Theta)*t)*cos(sqrt(g)*sqrt(Theta)*tq)*vq*&
    sqrt(Theta)*sin(sqrt(g)*sqrt(Theta)*tq)-cos(sqrt(g)*sqrt(Theta)*t)*g+&
    sin(sqrt(g)*sqrt(Theta)*t)*sin(sqrt(g)*sqrt(Theta)*tq)*g*cos(sqrt(g)*sqrt(Theta)*tq)-&
    sqrt(g)*sin(sqrt(g)*sqrt(Theta)*t)*vq*sqrt(Theta)+sqrt(g)*sin(sqrt(g)*sqrt(Theta)*t)*&
    vq*sqrt(Theta)*cos(sqrt(g)*sqrt(Theta)*tq)**2)/(sqrt(Theta)*sin(sqrt(g)*sqrt(Theta)*t)*&
    sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)**2-sin(sqrt(g)*sqrt(Theta)*t)*&
    cos(sqrt(g)*sqrt(Theta)*tq)*vq*Theta*sin(sqrt(g)*sqrt(Theta)*tq)-sqrt(Theta)*&
    sin(sqrt(g)*sqrt(Theta)*t)*sqrt(g)-sqrt(Theta)*cos(sqrt(g)*sqrt(Theta)*t)*sin(sqrt(g)*&
    sqrt(Theta)*tq)*sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)+cos(sqrt(g)*sqrt(Theta)*t)*vq*&
    Theta-cos(sqrt(g)*sqrt(Theta)*t)*vq*Theta*cos(sqrt(g)*sqrt(Theta)*tq)**2)

end function vBA
!--------------------------------------------------------------
real*8 function aBA(hq, theta, g, tq, vq, t)

real*8 hq, theta, g, tq, vq, t

aBA = (g+vq**2*Theta)*g/(-4*sqrt(Theta)*sin(sqrt(g)*sqrt(Theta)*t)*sqrt(g)*&
    cos(sqrt(g)*sqrt(Theta)*tq)**2*cos(sqrt(g)*sqrt(Theta)*t)*vq+g*cos(sqrt(g)*&
    sqrt(Theta)*tq)**2-cos(sqrt(g)*sqrt(Theta)*tq)**2*vq**2*Theta+2*cos(sqrt(g)*&
    sqrt(Theta)*t)**2*Theta*vq**2*cos(sqrt(g)*sqrt(Theta)*tq)**2-2*cos(sqrt(g)*&
    sqrt(Theta)*t)**2*g*cos(sqrt(g)*sqrt(Theta)*tq)**2-2*sin(sqrt(g)*sqrt(Theta)*t)*&
    g*cos(sqrt(g)*sqrt(Theta)*t)*sin(sqrt(g)*sqrt(Theta)*tq)*cos(sqrt(g)*sqrt(Theta)*tq)+2*&
    sin(sqrt(g)*sqrt(Theta)*t)*cos(sqrt(g)*sqrt(Theta)*tq)*sin(sqrt(g)*sqrt(Theta)*&
    tq)*vq**2*Theta*cos(sqrt(g)*sqrt(Theta)*t)-2*cos(sqrt(g)*sqrt(Theta)*tq)*sin(sqrt(g)*&
    sqrt(Theta)*tq)*vq*sqrt(Theta)*sqrt(g)+4*sqrt(Theta)*cos(sqrt(g)*sqrt(Theta)*t)**2*&
    sin(sqrt(g)*sqrt(Theta)*tq)*sqrt(g)*cos(sqrt(g)*sqrt(Theta)*tq)*vq+cos(sqrt(g)*&
    sqrt(Theta)*t)**2*g-cos(sqrt(g)*sqrt(Theta)*t)**2*Theta*vq**2+2*sin(sqrt(g)*&
    sqrt(Theta)*t)*sqrt(g)*sqrt(Theta)*cos(sqrt(g)*sqrt(Theta)*t)*vq-g)

end function aBA
!--------------------------------------------------------------
real*8 function hBD( theta, g, tH, xH, t)

real*8  theta, g, tH, vq, t, xH

hBD = 0.5d0*(2*sqrt(g)*sqrt(Theta)*t+2*log(2d0)-log(exp(-2*sqrt(g)*sqrt(Theta)*tH-2*xH*Theta)*&
    (exp(2*sqrt(g)*sqrt(Theta)*t)+exp(2*sqrt(g)*sqrt(Theta)*tH))**2))/Theta

end function hBD
!--------------------------------------------------------------
real*8 function vBD(theta, g, tH, xH, t)

real*8 theta, g, tH, t, xH

vBD = -sqrt(g)/sqrt(Theta)*(exp(2*sqrt(g)*sqrt(Theta)*t)-exp(2*sqrt(g)*sqrt(Theta)*tH))/&
    (exp(2*sqrt(g)*sqrt(Theta)*t)+exp(2*sqrt(g)*sqrt(Theta)*tH))

end function vBD
!--------------------------------------------------------------
real*8 function aBD( theta, g, tH, xH, t)

real*8 theta, g, tH, xH, t

aBD = -4*g*exp(2*sqrt(g)*sqrt(Theta)*(t+tH))/(exp(2*sqrt(g)*sqrt(Theta)*t)+exp(2*sqrt(g)*sqrt(Theta)*tH))**2

end function aBD
!--------------------------------------------------------------
end program trajeto_analitico    