SHARE
TWEET

Untitled

a guest Oct 8th, 2015 230 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. module constantes
  3. implicit none
  4. real(8), parameter :: g = 9.81, l = 0.1, m = 6.e-3, R = 8.32, &
  5. theta0 = 1., S = 1.e-4, P0 = 200., gam=5./3., T0=300.
  6. real(8), parameter :: V0 = S*l*theta0
  7. real(8), parameter :: pi=acos(-1.)
  8. real(8) :: nmol=P0*V0/(R*T0)
  9. real(8) :: Temp=100, Temp_final=300, tau= 100., E0=0.1, Emax=15., alpha=0.5, tmax=10
  10. end module constantes
  11.  
  12. subroutine nomfichier(nom) !pour entrer manuellement les nom des fichiers de resultats
  13.  
  14.     integer :: l
  15.     character (len=40) :: nom
  16.     character (len=20) :: nom2
  17.     write(*,'("Entrer un nom de fichier : ",$)')
  18.     read(*,*) nom2
  19.     l = len_trim(nom2)
  20.     nom=nom2(1:l)//'.res'
  21.     write(*,*) "Fichier resultat : ", nom
  22.  
  23. end subroutine nomfichier
  24.  
  25.  
  26. program main
  27.  
  28.     use constantes
  29.     implicit none
  30.  
  31.     real(8), dimension(2)     ::      theta
  32.     real(8), dimension(2)     ::      e
  33.     integer                   ::      nsteps, i, n,j, nsteps_j, jmax
  34.     real(8)                   ::      t, dt, dE
  35.     external                  ::      force, force_champ
  36.     character (len=40)        ::      nom !pour rentrer a la main le nom d'un fichier
  37.    character (len=4)         ::      ET
  38.    theta(1)=0.5             !condition initiale sur la position
  39.    theta(2)=0.              !C.I sur la vitesse
  40.    dt=0.001
  41.    dE=1.                   ! variation de E0 dans la boucle
  42.    nsteps = nint(tmax/dt)
  43.    n=2                     !dimension du tableau theta
  44.  
  45.    nom='res.dat'
  46.    !
  47.    !open (2,file=nom)
  48.  
  49.    call initialisation(dt, ET,jmax)
  50.    call nomfichier(nom)
  51.  
  52.    open(1,file=nom)
  53.  
  54.    Var_loop: do j=0, jmax
  55.  
  56.        time_loop: do i=0, nsteps
  57.            t=i*dt
  58.            call energie(t,theta,e,n)
  59.  
  60.            Select case (ET) !Choix du type d'etude
  61.                 case ("T" , "D")
  62.                     call rk4(t, theta , dt , n ,force) !etude de l'effet de la temperature (sans champ)
  63.  
  64.                case ("E")
  65.                    call rk4(t, theta , dt , n ,force_champ) ! etude de l'effet du champ a temperature fixe
  66.  
  67.                 case ("ET")
  68.                     call rk4(t, theta , dt , n ,force_champ) ! temperature et champ qui changent
  69.  
  70.             end select
  71.  
  72.             write(1,*) t, theta(1) ,theta(2), e(1), e(2), e(1)+e(2)
  73.         enddo time_loop
  74.  
  75.         write(1,*) " "
  76.         write(1,*) " "
  77.  
  78.         Select case (ET)
  79.             case ("T")
  80.                 Temp=Temp+50
  81.             case ("E")
  82.                 E0=E0+0.05
  83.             case default
  84.  
  85.         end select
  86.  
  87.     enddo Var_loop
  88.     close(1)
  89.     write(*,*) "jmax vaut: ", jmax, "Temp vaut ", Temp
  90.     write(*,*) 'Les résultats sont dans le fichier: ', nom
  91.  
  92. end program main
  93.  
  94.  
  95.  
  96. subroutine initialisation(dt, ET,jmax) !choix des parametres d'etude
  97.    use constantes
  98.    implicit none
  99.    integer                    :: jmax
  100.    real(8)                    :: dt, pas
  101.    character (len=4)          :: ET
  102.  
  103.    write(*,*) "Choisir le parametre a étudier: Temperature sans champ (T), champ variable, temperature fixe (E)&
  104.                ou temperature variable et champ fixe (ET) ? ou par defaut (D) ?"
  105.    read(*,*) ET
  106.  
  107.    select case (ET)
  108.  
  109.    case ('E')
  110.        write(*,*) "Choisir la valeur de depart de E0, sa valeur maximal, le pas et le coefficient de frottement alpha "
  111.        read(*,*) E0, Emax, pas, alpha
  112.        jmax= floor(abs((E0-Emax)/pas))
  113.    case('T')
  114.        write(*,*) "Choisir la valeur de depart de T initiale, sa valeur maximal, et le pas"
  115.        read(*,*) Temp, Tmax, pas
  116.        jmax= floor(abs((Temp-Tmax)/pas))
  117.    case('ET')
  118.        write(*,*) "Choisir la valeur de depart de T initiale, sa valeur maximal,&
  119.                    le pas et E0 et le coefficient de frottement alpha"
  120.        read(*,*) Temp, Tmax, pas, E0, alpha
  121.        jmax= floor(abs((Temp-Tmax)/pas))
  122.    case('D')
  123.        jmax=1
  124.        E0=1
  125.        alpha=1
  126.        Temp=200
  127.    end select
  128.  
  129. end subroutine initialisation
  130.  
  131. subroutine mean_and_variance(n,x,mean,variance)
  132.    !Computes the mean and variance of desired array x
  133.    implicit none
  134.    integer                             :: i
  135.    integer, intent(in)                 :: n
  136.    real(8), dimension(n), intent(in)   :: x
  137.    real(8), intent(inout)              :: mean, variance
  138.  
  139.    do i=1, n
  140.        mean = mean + x(i)
  141.    enddo
  142.    mean = mean/n
  143.  
  144.    do i=1, n   !somme des carrés
  145.        variance = variance + x(i)**2
  146.    enddo
  147.    variance = variance/n - mean**2
  148.    !formule de König-Huygens, plus pratique
  149.  
  150. end subroutine
  151.  
  152.  
  153.  
  154.  
  155.  
  156. subroutine energie(t,theta,e,n) ! Calcul des energies pot et cinetique
  157.    use constantes
  158.    implicit none
  159.    integer ,intent(in)                 ::   n
  160.    real(8) ,intent(in)                 ::   t
  161.    real(8) ,dimension(n) , intent(in)  ::   theta
  162.    real(8) ,dimension(n) , intent(inout) ::   e
  163.  
  164.    e(1)=0.5*m*l**2*theta(2)**2 !energie cinetique
  165.    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))
  166.  
  167. end subroutine energie
  168.  
  169.  
  170. subroutine force(t,theta,dtheta,n) !calcul des derivees
  171.    use constantes
  172.    implicit none
  173.    integer ,intent(in)                 ::   n
  174.    real(8) ,intent(in)                 ::   t
  175.    real(8) ,dimension(n) , intent(in)  ::   theta
  176.    real(8) ,dimension(n) , intent(out) ::   dtheta
  177.  
  178.    dtheta(2)=g/l*sin(theta(1)) + ((nmol*R*Temp)/(m*l**2*theta0))*&
  179.    ((1.-theta(1)/theta0)**gam-(1+theta(1)/theta0)**gam)/(1.-theta(1)**2/theta0**2)**gam
  180.  
  181.    dtheta(1)=theta(2)
  182.  
  183. end subroutine force
  184.  
  185.  
  186. subroutine force_champ(t,theta,dtheta,n) ! diviser par m*l le champ et la dissipation ?
  187.    use constantes
  188.    implicit none
  189.    integer ,intent(in)                 ::   n
  190.    real(8) ,intent(in)                 ::   t
  191.    real(8) ,dimension(n) , intent(in)  ::   theta
  192.    real(8) ,dimension(n) , intent(out) ::   dtheta
  193.  
  194.         dtheta(1)=theta(2)
  195.  
  196.    dtheta(2)=g/l*sin(theta(1)) + ((nmol*R*Temp)/(m*l**2*theta0))*&
  197.    ((1.-theta(1)/theta0)**gam-(1+theta(1)/theta0)**gam)/(1.-theta(1)**2/theta0**2)**gam &
  198.    +(E0*sin(2*pi*t/tau) -alpha*theta(2))
  199.  
  200. end subroutine force_champ
  201.  
  202.  
  203.  
  204. subroutine rk4(x,y,dx,n, deriv)             ! ici x est le temps, y = theta, dx=dt
  205.  
  206.    implicit none
  207.    integer                 ,intent(in)     :: n
  208.    real(8)                 ,intent(in)     :: x, dx
  209.    real(8)                                 :: ddx
  210.    real(8), dimension(n)   ,intent(inout)  :: y
  211.    real(8), dimension(n)                   :: yp, k1, k2, k3, k4
  212.  
  213.    ddx = 0.5*dx
  214.    call deriv(x,y,k1,n)           ;   yp = y + ddx*k1
  215.    call deriv(x + ddx,yp,k2,n)    ;   yp = y + ddx*k2
  216.    call deriv(x + ddx,yp,k3,n)    ;   yp = y +  dx*k3
  217.    call deriv(x + dx, yp,k4,n)    ;   y  = y +  dx*(k1 + 2.*k2 + 2.*k3 + k4 )/6.
  218.  
  219. end subroutine rk4
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top