PROGRAM Lyapunov
  IMPLICIT NONE !boa prática (se não colocar o fortran aceita qualquer variável)
  CHARACTER(LEN=16), PARAMETER :: arquivoSaida = 'EspectroLyapunov'
  INTEGER, PARAMETER :: transiente = 1.d5, pontos = 1.d5! Valores a serem calculados
  DOUBLE PRECISION, PARAMETER :: sigma = 1.d1, r = 28.d0, b = 8.d0/3.d0
  DOUBLE PRECISION, PARAMETER :: passoRK4 = 1.d-2! passoRK4 de integração
  DOUBLE PRECISION, DIMENSION(3) :: v, expoenteLyapunov ! x,y,z do Lorenz
  DOUBLE PRECISION, DIMENSION(3,3) :: vv! Variáveis Linearizadas
  INTEGER :: ii
  
  expoenteLyapunov = 0.d0
  v(1) = 0.d0! condicao inicial 
  v(2) = 1.d0! condicao inicial
  v(3) = 0.d0! condicao inicial
!
  vv = 0.d0
  FORALL(ii = 1:3) vv(ii,ii) = 1.d0
!
  DO ii = 1, pontos+transiente! loop de integração
    CALL RK4(v, vv)
    CALL OrtoGRS(ii > transiente, vv, expoenteLyapunov)
  END DO  

  expoenteLyapunov = expoenteLyapunov/(DFLOAT(pontos)*passoRK4)
  WRITE (*,*) expoenteLyapunov, SUM(expoenteLyapunov)
!  CALL SalvaSolucao  
!
CONTAINS! Subrotinas e funções independentes
  SUBROUTINE SistemaLorenz(x, xx, x_ponto, xx_ponto)!Sistema de Lorenz
    IMPLICIT NONE 
    DOUBLE PRECISION, INTENT(IN), DIMENSION(3) :: x! não altera o valor
    DOUBLE PRECISION, INTENT(IN), DIMENSION(3,3) :: xx! Matriz das eqs linearizadas
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3) :: x_ponto! calcula o valor da derivada
    DOUBLE PRECISION, INTENT(OUT), DIMENSION(3,3) :: xx_ponto! Derivada d para as eqs linearizadas
    
    x_ponto(1) = sigma*(x(2) - x(1)) !Espaço de fase
    x_ponto(2) = r*x(1) - x(2) - x(1)*x(3) ! dy/dt=rx-y-xz
    x_ponto(3) = x(1)*x(2) - b*x(3)  
!
    xx_ponto(1,:) = sigma*(xx(2,:) - xx(1,:))! Equações linearizadas (Espaço tangente)
    xx_ponto(2,:) = r*xx(1,:) - xx(2,:) - xx(1,:)*x(3) - x(1)*xx(3,:) ! Equações linearizadas
    xx_ponto(3,:) = xx(1,:)*x(2) + x(1)*xx(2,:) - b*xx(3,:)  ! Equações linearizadas
    
  END SUBROUTINE SistemaLorenz
  SUBROUTINE RK4(x, xx)! Calcula novos valores de x, y, z a partir dos velhos
    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3) :: x
    DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3,3) :: xx
    DOUBLE PRECISION, DIMENSION(3) :: k1, k2, k3, k4
    DOUBLE PRECISION, DIMENSION(3,3) :: kk1, kk2, kk3, kk4

    CALL SistemaLorenz(x, xx, k1, kk1)
    CALL SistemaLorenz(x + k1*passoRK4/2.d0, xx + kk1*passoRK4/2.d0, k2, kk2)
    CALL SistemaLorenz(x + k2*passoRK4/2.d0, xx + kk2*passoRK4/2.d0, k3, kk3)
    CALL SistemaLorenz(x + k3*passoRK4, xx + kk3*passoRK4, k4, kk4)
    
    x = x + passoRK4*(k1 + 2.d0*(k2 + k3) + k4)/6.d0
    xx = xx + passoRK4*(kk1 + 2.d0*(kk2 + kk3) + kk4)/6.d0
  END SUBROUTINE RK4
 
  SUBROUTINE OrtoGRS(acumulaLyapunov, xx, expoenteAcumulado)
    IMPLICIT NONE
    LOGICAL, INTENT(IN) :: acumulaLyapunov
    DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3) :: expoenteAcumulado
    DOUBLE PRECISION, INTENT(INOUT), DIMENSION(3,3) :: xx
    DOUBLE PRECISION, DIMENSION(3) :: norma
    DOUBLE PRECISION, DIMENSION(2) :: projecao
    INTEGER :: i, j
!    
    norma(1) = DSQRT(SUM(xx(:,1)*xx(:,1)))
    xx(:,1) = xx(:,1)/norma(1)
!
    DO i = 2, 3
      DO j = 1, i - 1
        projecao(j) = SUM(xx(:,i)*xx(:,j))
      END DO
!
      DO j = 1, i - 1
        xx(:,i) = xx(:,i) - projecao(j)*xx(:,j)
      END DO 
!
      norma(i) = DSQRT(SUM(xx(:,i)*xx(:,i)))
      xx(:,i) = xx(:,i)/norma(i)      
    END DO
!
    IF (acumulaLyapunov) WHERE (norma /= 0.d0) expoenteAcumulado  = expoenteAcumulado + DLOG(norma)     
  END SUBROUTINE OrtoGRS
 
  SUBROUTINE SalvaSolucao! Salva solução do Lorenz
    IMPLICIT NONE
    INTEGER :: i
  
    OPEN(101, FILE = arquivoSaida, STATUS = 'unknown')
 !   DO i = 0, tempo
!      WRITE(101, '(4F10.4)') DFLOAT(i)*passoRK4, xSolucao(i), ySolucao(i), zSolucao(i)
!    END DO
    CLOSE(101)    
  END SUBROUTINE SalvaSolucao
END PROGRAM Lyapunov