module constantes
implicit none
real(8), parameter :: g = 9.81, l = 0.1, m = 6.e-3, R = 8.32, &
theta0 = 1., S = 1.e-4, P0 = 200., gam=5./3., T0=300.
real(8), parameter :: V0 = S*l*theta0
real(8), parameter :: pi=acos(-1.)
real(8) :: nmol=P0*V0/(R*T0)
real(8) :: Temp=100, Temp_final=300, tau= 100., E0=0.1, Emax=15., alpha=0.5, tmax=10
end module constantes

subroutine nomfichier(nom) !pour entrer manuellement les nom des fichiers de resultats

    integer :: l
    character (len=40) :: nom
    character (len=20) :: nom2
    write(*,'("Entrer un nom de fichier : ",$)')
    read(*,*) nom2
    l = len_trim(nom2)
    nom=nom2(1:l)//'.res'
    write(*,*) "Fichier resultat : ", nom

end subroutine nomfichier


program main

    use constantes
    implicit none

    real(8), dimension(2)     ::      theta
    real(8), dimension(2)     ::      e
    integer                   ::      nsteps, i, n,j, nsteps_j, jmax
    real(8)                   ::      t, dt, dE
    external                  ::      force, force_champ
    character (len=40)        ::      nom !pour rentrer a la main le nom d'un fichier
    character (len=4)         ::      ET
    theta(1)=0.5             !condition initiale sur la position
    theta(2)=0.              !C.I sur la vitesse
    dt=0.001
    dE=1.                   ! variation de E0 dans la boucle
    nsteps = nint(tmax/dt)
    n=2                     !dimension du tableau theta

    nom='res.dat'
    !
    !open (2,file=nom)

    call initialisation(dt, ET,jmax)
    call nomfichier(nom)

    open(1,file=nom)

    Var_loop: do j=0, jmax

        time_loop: do i=0, nsteps
            t=i*dt
            call energie(t,theta,e,n)

            Select case (ET) !Choix du type d'etude
                case ("T" , "D")
                    call rk4(t, theta , dt , n ,force) !etude de l'effet de la temperature (sans champ)

                case ("E")
                    call rk4(t, theta , dt , n ,force_champ) ! etude de l'effet du champ a temperature fixe

                case ("ET")
                    call rk4(t, theta , dt , n ,force_champ) ! temperature et champ qui changent

            end select

            write(1,*) t, theta(1) ,theta(2), e(1), e(2), e(1)+e(2)
        enddo time_loop

        write(1,*) " "
        write(1,*) " "

        Select case (ET)
            case ("T")
                Temp=Temp+50
            case ("E")
                E0=E0+0.05
            case default

        end select

    enddo Var_loop
    close(1)
    write(*,*) "jmax vaut: ", jmax, "Temp vaut ", Temp
    write(*,*) 'Les résultats sont dans le fichier: ', nom

end program main



subroutine initialisation(dt, ET,jmax) !choix des parametres d'etude
    use constantes
    implicit none
    integer                    :: jmax
    real(8)                    :: dt, pas
    character (len=4)          :: ET

    write(*,*) "Choisir le parametre a étudier: Temperature sans champ (T), champ variable, temperature fixe (E)&
                ou temperature variable et champ fixe (ET) ? ou par defaut (D) ?"
    read(*,*) ET

    select case (ET)

    case ('E')
        write(*,*) "Choisir la valeur de depart de E0, sa valeur maximal, le pas et le coefficient de frottement alpha "
        read(*,*) E0, Emax, pas, alpha
        jmax= floor(abs((E0-Emax)/pas))
    case('T')
        write(*,*) "Choisir la valeur de depart de T initiale, sa valeur maximal, et le pas"
        read(*,*) Temp, Tmax, pas
        jmax= floor(abs((Temp-Tmax)/pas))
    case('ET')
        write(*,*) "Choisir la valeur de depart de T initiale, sa valeur maximal,&
                    le pas et E0 et le coefficient de frottement alpha"
        read(*,*) Temp, Tmax, pas, E0, alpha
        jmax= floor(abs((Temp-Tmax)/pas))
    case('D')
        jmax=1
        E0=1
        alpha=1
        Temp=200
    end select

end subroutine initialisation

subroutine mean_and_variance(n,x,mean,variance)
    !Computes the mean and variance of desired array x
    implicit none
    integer                             :: i
    integer, intent(in)                 :: n
    real(8), dimension(n), intent(in)   :: x
    real(8), intent(inout)              :: mean, variance

    do i=1, n
        mean = mean + x(i)
    enddo
    mean = mean/n

    do i=1, n   !somme des carrés
        variance = variance + x(i)**2
    enddo
    variance = variance/n - mean**2
    !formule de König-Huygens, plus pratique

end subroutine





subroutine energie(t,theta,e,n) ! Calcul des energies pot et cinetique
    use constantes
    implicit none
    integer ,intent(in)                 ::   n
    real(8) ,intent(in)                 ::   t
    real(8) ,dimension(n) , intent(in)  ::   theta
    real(8) ,dimension(n) , intent(inout) ::   e

    e(1)=0.5*m*l**2*theta(2)**2 !energie cinetique
    e(2)=m*g*l*cos(theta(1))+(nmol*R*Temp/(gam-1.))*((1-theta(1)/theta0)**(1-gam)+(1+theta(1)/theta0)**(1-gam))

end subroutine energie


subroutine force(t,theta,dtheta,n) !calcul des derivees
    use constantes
    implicit none
    integer ,intent(in)                 ::   n
    real(8) ,intent(in)                 ::   t
    real(8) ,dimension(n) , intent(in)  ::   theta
    real(8) ,dimension(n) , intent(out) ::   dtheta

    dtheta(2)=g/l*sin(theta(1)) + ((nmol*R*Temp)/(m*l**2*theta0))*&
    ((1.-theta(1)/theta0)**gam-(1+theta(1)/theta0)**gam)/(1.-theta(1)**2/theta0**2)**gam

    dtheta(1)=theta(2)

end subroutine force


subroutine force_champ(t,theta,dtheta,n) ! diviser par m*l le champ et la dissipation ?
    use constantes
    implicit none
    integer ,intent(in)                 ::   n
    real(8) ,intent(in)                 ::   t
    real(8) ,dimension(n) , intent(in)  ::   theta
    real(8) ,dimension(n) , intent(out) ::   dtheta

    dtheta(1)=theta(2)

    dtheta(2)=g/l*sin(theta(1)) + ((nmol*R*Temp)/(m*l**2*theta0))*&
    ((1.-theta(1)/theta0)**gam-(1+theta(1)/theta0)**gam)/(1.-theta(1)**2/theta0**2)**gam &
    +(E0*sin(2*pi*t/tau) -alpha*theta(2))

end subroutine force_champ



subroutine rk4(x,y,dx,n, deriv)             ! ici x est le temps, y = theta, dx=dt

    implicit none
    integer                 ,intent(in)     :: n
    real(8)                 ,intent(in)     :: x, dx
    real(8)                                 :: ddx
    real(8), dimension(n)   ,intent(inout)  :: y
    real(8), dimension(n)                   :: yp, k1, k2, k3, k4

    ddx = 0.5*dx
    call deriv(x,y,k1,n)           ;   yp = y + ddx*k1
    call deriv(x + ddx,yp,k2,n)    ;   yp = y + ddx*k2
    call deriv(x + ddx,yp,k3,n)    ;   yp = y +  dx*k3
    call deriv(x + dx, yp,k4,n)    ;   y  = y +  dx*(k1 + 2.*k2 + 2.*k3 + k4 )/6.

end subroutine rk4