PROGRAM Lyapunov
  IMPLICIT NONE !boa prática (se não colocar o fortran aceita qualquer variável)
  CHARACTER(LEN=17), PARAMETER :: arquivoSaida = 'EspectroLyapunov2'
  INTEGER, PARAMETER :: transiente = 1.d6, pontos = 1.d6, passos = 1.d3+1! Valores a serem calculados
  DOUBLE PRECISION, PARAMETER :: pMinimo = 0.d0, pMaximo = 1.d2
  DOUBLE PRECISION, PARAMETER :: c = 5.d0, b = 1.d2
  DOUBLE PRECISION, PARAMETER :: passoRK4 = 1.d-3! passoRK4 de integração
  DOUBLE PRECISION, DIMENSION(3,passos) :: expoenteLyapunov ! 
  DOUBLE PRECISION, DIMENSION(3) :: v! x,y,z do Vallis
  DOUBLE PRECISION, DIMENSION(3,3) :: vv! Variáveis Linearizadas
  DOUBLE PRECISION :: p, pPasso
  INTEGER :: ii, jj
  
  pPasso = (pMaximo - pMinimo)/DFLOAT(passos - 1)
  
  expoenteLyapunov = 0.d0
  p = pMinimo - pPasso
  DO ii = 1, passos
    WRITE(*,*) ii
    v(1) = 0.d0! condicao inicial 
    v(2) = 1.d0! condicao inicial
    v(3) = 0.d0! condicao inicial
!
    vv = 0.d0
    FORALL(jj = 1:3) vv(jj,jj) = 1.d0
!
    p = p + pPasso
    DO jj = 1, pontos+transiente! loop de integração
      CALL RK4(v, vv)
      CALL OrtoGRS(jj > transiente, vv, expoenteLyapunov(:,ii))
    END DO  
  END DO  

  expoenteLyapunov = expoenteLyapunov/(DFLOAT(pontos)*passoRK4)
  CALL SalvaEspectro  
!
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) = b*x(2) - c*(x(1) + p) !Espaço de fase
    x_ponto(2) = x(1)*x(3) - x(2) !
    x_ponto(3) = -x(1)*x(2) - x(3)  + 1.d0
!
    xx_ponto(1,:) = b*(xx(2,:)) -c*xx(1,:) ! Equações linearizadas (Espaço tangente)
    xx_ponto(2,:) = xx(1,:)*x(3) + x(1)*xx(3,:) - xx(2,:) ! Equações linearizadas
    xx_ponto(3,:) = -xx(1,:)*x(2) - x(1)*xx(2,:) - 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 SalvaEspectro! Salva solução do Lorenz
    IMPLICIT NONE
    INTEGER :: i
  
    OPEN(101, FILE = arquivoSaida, STATUS = 'unknown')
    p = pMinimo - pPasso
    DO i = 1, passos
      p = p + pPasso
      WRITE(101, '(4F12.6)') p, expoenteLyapunov(:,i)
    END DO
    CLOSE(101)    
  END SUBROUTINE SalvaEspectro
END PROGRAM Lyapunov