!----------------------------------------------------!
! PRÁCTICA 8. Resolución de Ecuaciones Diferenciales  !	   
!					Ordinarias con Valores Iniciales  !                       
! 												      !
!												      !
! Félix Pedrera García							      !  	 
! 								      !
!-----------------------------------------------------!


program EDO
real, parameter:: pi=3.1415926
real:: g0, l, theta0, theta10, h
real::n, t, theta, theta1, F0, omegaf, k !,Amp
real::a, b
integer::mov
character (len=11), fichero

g0=9.8
l=0.5
theta0=pi/4
theta10=0
h=0.02
a=0
b=4
omegaf=2



print*, 'Programa para calcular resolver la ecuación diferencial del péndulo'
pause  
!print*, ' '
!print*, 'Selecciona una de las siguientes opciones:'
!print*, ' ' 
!print*, '(1) Movimiento de poca amplitud, sin(x)-->x'
!print*, '(2) Movimiento general'

!Read*, mov	
!if (mov==2) then 
	!Amp=sin(theta)
	!else 
	!Amp=theta
!end if

print*, ' '
print*, 'Selecciona el tipo de movimiento:'
print*, '  '
print*, '(1) Armonico simple, k=0, F0=0'
print*, '    Condiciones iniciales y parametros:'
print*, '    * Aceleracion de la gravedad: g=9.8 m/s^2.'
print*, '    * Longitud del pendulo: 0.5 m.'
print*, '    * Angulo inicial: pi/4 rad.'
print*, '    * Velocidad angular inicial: 0 rad/s.'
print*, '    * Intervalo de tiempo: 4 s.'
print*, '    * Paso: 0.02 s.'	 
print*, ' '
print*, '(2) Amortiguado y Forzado'
print*, '    Condiciones iniciales y parametros:'
print*, '    * Aceleracion de la gravedad: g=9.8 m/s^2.'
print*, '    * Longitud del pendulo: l=9.8 m.'
print*, '    * Angulo inicial: 0.333 rad.'
print*, '    * Velocidad angular inicial: 0.1 rad/s.'
print*, '    * Intervalo de tiempo: 400 s.'
print*, '    * Paso: 0.1 s.'
print*, '    * Constante de Amortiguamiento: k=0.2 N/m'
print*, '    * Amplitud de la Fuerza externa: 1.25 N.'
print*, '    * Frecuencia de la Fuerza externa: 0.3 rad/s.'
print*, ' '    
print*, '(3) Movimiento personalizado'
print*, ' '

read*, mov
if (mov==1) then 
	k=0
	F0=0
	print*, 'Has seleccionado Movimiento Armónico Simple' 
	pause
end if

if (mov==2) then 
	k=0.2
	F0=1.25
	omegaf=0.3
	l=9.8
	theta0=0.333
	theta10=0.1
	a=0
	b=400
	h=0.1
	print*, 'Has seleccionado Movimiento Amortiguado y Forzado'
	pause
end if 

if (mov==3) then
	print*, 'Has seleccionado Moviento Personalizado'
	print*, ' ' 
	print*, 'Introduce el valor de las siguientes condiciones iniciales y parametros'
	print*, 'Aceleracion de la gravedad, g0'
	read*, g0
	print*, 'Longitud del pendulo, l'
	read*, l
	print*, 'Angulo inicial, theta0'
	read*, theta0
	print*, 'Velocidad angular inicial, theta10'
	read*, theta10
	print*, 'Tiempo inicial, a'
	read*, a
	print*, 'Tiempo final, b'
	read*, b
	print*, 'Paso, h'
	read*, h
	print*, 'Constante de Amortiguamiento, k'
	read*, k
	print*, 'Amplitud de la Fuerza externa, F0'
	read*, F0
	print*, 'Frecuencia de la Fuerza externa, omegaf'
	read*, omegaf
end if

n=(b-a)/h
t=a
theta=theta0
theta1=theta10	

print*, 'Escribe el nombre del fichero de resultados'
read "(A)", fichero

open (1, file=fichero)		

do i=1, n
	d11=G(t, theta, theta1)
	d12=F(t, theta, theta1)
	d21=G(t+0.5*h, theta+0.5*h*d11, theta1+0.5*h*d12)
	d22=F(t+0.5*h, theta+0.5*h*d11, theta1+0.5*h*d12)
	d31=G(t+0.5*h, theta+0.5*h*d21, theta1+0.5*h*d22)
	d32=F(t+0.5*h, theta+0.5*h*d21, theta1+0.5*h*d22)
	d41=G(t+h, theta+h*d31, theta1+h*d32)
	d42=F(t+h, theta+h*d31, theta1+h*d32)
	
	theta=theta+((h/6)*(d11+2*d21+2*d31+d41))
	theta1=theta1+((h/6)*(d12+2*d22+2*d32+d42))
	write (1,*) t, theta, theta1
	print*, t, theta, theta1

	t=a+i*h
end do

close (1)

	print*, 'Aceleracion de la gravedad, g0'
	print*, g0
	print*, 'Longitud del pendulo, l'
	print*, l
	print*, 'Angulo inicial, theta0'
	print*, theta0
	print*, 'Velocidad angular inicial, theta10'
	print*, theta10
	print*, 'Tiempo inicial, a'
	print*, a
	print*, 'Tiempo final, b'
	print*, b
	print*, 'Paso, h'
	print*, h
	print*, 'Constante de Amortiguamiento, k'
	print*, k
	print*, 'Amplitud de la Fuerza externa, F0'
	print*, F0
	print*, 'Frecuencia de la Fuerza externa, omegaf'
	print*, omegaf

			 
contains



function F(t, theta, theta1)
	real::t, theta, theta1, F
	F=((F0*cos(omegaf*t))-(k*theta1)-((g0/l)*theta)) !Amp
end function

function G(t, theta, theta1)
	real::t, theta, theta1, G
	G=theta1
end function 



end program