! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2013
!
! ejemplo 6a
!
! ecuacion de onda (hiperbolica)
! d2f/dt2 = nu2 * d2f/du2
!
! donde nu2 es el cuadrado de una constante con dimensiones
! de velocidad siendo su raiz cuadrada el cociente entre
! la tension de la cuerda (tau) y la densidad lineal (rho).
!
! la funcion f(x,t) esta definida para 0 <= x <= 1, 
! donde f(x,0)=1.25*x para x <= 0.8 y f(x,0)=5*(1-x) para x > 0.8,
! f(0,t)=f(1,t)=0 y df/dt (t=0) = 0
!
!
program eje6a
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer :: ii,jj,kk
  real(p) :: tau,rho,nu,dt,dx,nup
  real(p),dimension(101) :: fun,fun1,fun2 
  
! parametros iniciales
!
  write(*,*)'Tension cuerda (N): '
  read(*,*)tau
  write(*,*)'Densidad lineal (kg/m): '
  read(*,*)rho

  nu=sqrt(tau/rho)

  write(*,*)'Nu_p / Nu (>= 1) : '
  read(*,*)nup
  nup=nup*nu
  dx=0.01_p
  dt=dx/nup

! valores iniciales de la funcion
!
  fun=0._p

! pulso unico en x=0.8 
!
  do ii=1,81
    fun(ii)=1.25*float(ii-1)*dx
  enddo
  do ii=82,101
    fun(ii)=1._p-(5._p*float(ii-81)*dx)
  enddo

! archivo de salida
!
  open(7,file='eje6a.out')

  ii=0

! guarda valores iniciales
!
  do jj=1,101,2
    write(7,*)float(ii),float(jj-1),fun(jj)
  enddo

! primer paso
!
  fun1=0._p
  do jj=2,100
    fun1(jj)=0.5_p*(nu/nup)**2*(fun(jj-1)-2._p*fun(jj)+fun(jj+1))
  enddo
  fun1=fun+fun1
  fun2=fun1

! bucle de calculo desde el segundo paso
!
  ii=1
  kk=0
  do while(ii <= 100)

! escribe en archivo
!
    if(kk == 4)then
      kk=0
      do jj=1,101,2
        write(7,*)float(ii),float(jj-1),fun2(jj)
      enddo
    else
      kk=kk+1
    endif

! avanza un paso temporal
!
    fun=fun1
    fun1=fun2
    fun2=0._p
    do jj=2,100
      fun2(jj)=(nu/nup)**2*(fun1(jj-1)-2._p*fun1(jj)+fun1(jj+1))
    enddo
    fun2=2._p*fun1-fun+fun2
    ii=ii+1
  enddo
  close(7)

end program eje6a
