C 		Programme de simulation d'un oscillateur amorti par un frottement fluide
C		Dominique LEFEBVRE - avril 2006

	  PROGRAM OscillateurAmorti

	  IMPLICIT None
C 		Declaration des routines externes
C		ResolutionRK4 : routine d'implementation du schema RK4
C		Derivee       : description du systeme differentiel de l'oscillateur	  
      EXTERNAL ResolutionRK4
	  EXTERNAL Derivee

C		Declaration des parametres RK4
C		Le systeme differentiel comporte deux equations
      INTEGER NbEquation
      PARAMETER (NbEquation = 2)

C      Declaration des constantes
	  REAL PI
	  DATA PI/3.141592654/

C		Declaration des variables
	  INTEGER i, NbPas
	  REAL y(NbEquation), t, PasTemps
      REAL omega0, Q, x0, v0, alpha, beta, Ta, T0

C		Declaration des variables communes
	  COMMON /Oscillateur/alpha, beta

C       Saisie des parametres du mouvement
      WRITE(*,'(a,$)') 'Pulsation propre de l''oscillateur (rd.s-1) : '
      READ(*,*) omega0
      WRITE(*,'(a,$)') 'Facteur de qualite de l''oscillateur : '
      READ(*,*) Q
      WRITE(*,'(a,$)') 'Abscisse initiale x0  (en m) : '
      READ(*,*) x0
      WRITE(*,'(a,$)') 'Vitesse initiale v0 (en m.s-1) : '
      READ(*,*) v0     
 

C		calcul des valeurs intermediaires pour optimiser le traitement
C       dans Derivee
	  alpha = omega0/Q
      beta = omega0*omega0

C		Calcul de la periode propre et la pseudo periode de l'oscillateur
	  T0 = 2*PI/omega0
	  Ta = T0/SQRT(1 - 1/(4*Q*Q))

C		Pour obtenir un trace exploitable, il faut estimer PasTemps
C		en fonction de omega0
	  IF (omega0 .LT. 10.0d0) THEN
		PasTemps = 0.01d0
   	  ELSE
		PasTemps = 0.1d0/omega0
	  ENDIF

C		Calcul du nombre de pas en fonction de Q
	  IF (Q .LE. 0.5d0) THEN
		  NbPas = T0*3/PasTemps
	  ELSE
		  NbPas = (Ta*Q+1)/PasTemps   ! Q*Ta donne approximativement le nb d'oscillations
      ENDIF

C		Ouverture du fichier de stockage des resultats de calcul
C		Ce fichier sera trace par GnuPlot
	  OPEN(10, FILE='OscillateurAmorti.dat')

C		Ecriture des conditions initiales
	  y(1) = x0
      y(2) = v0
      t = 0.0d0
	  WRITE(10,* ) t, y(1), y(2)

C		Boucle de calcul
	  DO i=1, NbPas
		t = i*PasTemps
		CALL ResolutionRK4(t, y, PasTemps, NbEquation, Derivee)
		WRITE(10,*) t, y(1), y(2)
	  ENDDO

C		Fermeture du fichier de donnees
	  CLOSE(10)

	  STOP
	  END