C		Programme de simulation d'un paquet d'ondes
C		Dominique Lefebvre 
C		Scientific Programming
C		Juin 2009
		
		PROGRAM PaquetOndes
		IMPLICIT NONE

C		Declaration des constantes
		DOUBLE PRECISION pi, Omega0, a
		DOUBLE COMPLEX ci, u, Enveloppe
		INTEGER n, nt , i, it
C		Declaration des parametres		
		PARAMETER (pi = 3.141592654)
		PARAMETER (n = 1000)		 	! pas d'integration sur Omega
		PARAMETER (nt = 2000)			! pas d'integration sur le temps
C    	Declaration des variables
		DOUBLE PRECISION Omega, dOmega,k,z,t,t0,dt
		REAL TTAB(nt), UTAB(nt)

C		Definition des caracteristiques du paquet 'onde
C		Omega0 : frequence du signal (en Hz)
C		a :      intervalle entre les noeuds (en m)
		Omega0 = 50.
        a=2.
C		Initialisation des constantes et parametres
		ci = (0.d0, 1.d0)
		dt = 0.1/Omega0
		Omega0 = 2.*pi*Omega0
        dOmega = Omega0/(n-1)

C		calcul et representation de U(z,t) en fonction de z
10		CONTINUE
		WRITE(*,'("z en cm (-1 pour sortir) : ",$)')
		READ (*,*) z
		IF (z . EQ. -1) GOTO 20

C		Boucle sur t
		DO it = 1, nt
			u = (0.d0, 0.d0)
			t = (it - nt/2)*dt
C			Boucle sur Omega
			DO i = 1, n
				Omega = (i - 1)*dOmega
C				Calcul de k avec l'equation de dispersion
				k = 2./a*ASIN(Omega/Omega0)
C				Integration de F(Omega)
				u = u + Enveloppe(Omega)*EXP(ci*(Omega*t - k*z))
			ENDDO
 			TTAB(it) = REAL(t)
			UTAB(it) = REAL(DBLE(u))	! extraction de la partie reelle
		ENDDO

C 		Dessin du paque d'ondes
		CALL TraceCourbe(TTAB,UTAB,nt)


C		Condition de sortie (simule le REPEAT WHILE)
		GOTO 10

C		fin du programme
20		WRITE(*,'("Fin du programme PaquetOndes",$)')
		END

C       ---------------------------------------------------------------
C		Routine de calcul de l'enveloppe gaussienne
C       ---------------------------------------------------------------
		DOUBLE COMPLEX FUNCTION Enveloppe (Omega)
		IMPLICIT NONE
		DOUBLE PRECISION FrequCentrale, Omega, Largeur, pi
		PARAMETER (pi = 3.141592654)
C		Definition des parametres de la courbe enveloppe
			FrequCentrale = 50.0		! frequence centrale du paquet
			Largeur = 5.0				! largeur du paquet
			FrequCentrale = 2.*pi*FrequCentrale
        	Largeur = 2.*pi*Largeur

C		Definition de la gaussienne enveloppe		
			Enveloppe = EXP(-((Omega - FrequCentrale)/Largeur)**2)
		END

C       ---------------------------------------------------------------
C		Routine de calcul de l'element MAX d'un vecteur de REAL
C		T : vecteur de REAL
C		n : nombre  (INTEGER)
C		---------------------------------------------------------------
		REAL FUNCTION MaxArray (T,n)
		IMPLICIT NONE
		REAL T(*), MaxTab
		INTEGER i,n
			MaxTab = T(1)
			DO i = 2,n 
				Maxtab = MAX(T(i),MaxTab)
			ENDDO
			MaxArray = MaxTab
		END

C       ---------------------------------------------------------------
C		Routine de calcul de l'element MIN d'un vecteur de REAL
C		T : vecteur de REAL
C		n : nombre  (INTEGER)
C		---------------------------------------------------------------
		REAL FUNCTION MinArray (T,n)
		IMPLICIT NONE
		REAL T(*), MinTab
		INTEGER i,n
			MinTab = T(1)
			DO i = 2,n 
				Mintab = MIN(T(i),MinTab)
			ENDDO
			MinArray = MinTab
		END

C		---------------------------------------------------------------
C		Routine de trace de courbe
C		X, Y : vecteurs REAL
C       n : nombre de points
C		---------------------------------------------------------------
		SUBROUTINE TraceCourbe(X,Y,n)
		IMPLICIT NONE
		REAL X(*), Y(*)
		INTEGER n
		REAL Ymin, Ymax, Xmin, XMax
		REAL MinArray, MaxArray

C			Initilaisation de la fenetre graphique
        	CALL METAFL('XWIN')
			CALL DISINI()
			CALL PAGFLL(255)
			CALL COLOR('BLACK')
			CALL WINFNT ('ARIAL')
       		CALL PAGHDR('Dominique LEFEBVRE - ',' - Scientific Programming',4,0)
C			Definition des legendes du graphe
			CALL NAME ('Distance', 'X')
       		CALL NAME ('Amplitude', 'Y')
        	CALL TITLIN ('Paquet d''ondes', 2)
 			CALL CENTER()
C			Definition de la fenetre de trace du graphe
			Xmin = MinArray(X,n)
			Xmax = MaxArray(X,n)
			Ymin = -100.
			Ymax = 100.
			CALL GRAF(Xmin,Xmax,Xmin,1.,Ymin,Ymax,Ymin,50.)
			CALL BOX2D()
			CALL AXGIT()
			CALL TITLE()		
			CALL SOLID()
C			Trace de la courbe                                          
        	CALL CURVE (X, Y, n)
C			Fin		
 			CALL DISFIN()
		END