        PROGRAM MAXWELL
C       Erstellt eine Histogramm fuer die Geschwindigkeit von
C	Punktpartikeln, die durch eindimensionale elastische Stoesse
C	ins thermische Gleichgewicht kommen.
Compiliere mit "g77 maxwell.f"
C
        IMPLICIT NONE
        INTEGER NMAX,NB
        PARAMETER(NMAX=51, NB=45)
C	Geschwindigkeit V, Ort X, 1/Masse Z, Histogramm HISTO,
C	Stosszeit TAU, 
        REAL*8 V(0:NMAX), X(0:NMAX), Z(0:NMAX), DV, TEMPER, SIG
	INTEGER HISTO(-NB:NB), NBIN, M, N, K, KEND
	DATA HISTO/NB*0,0,NB*0/
C
        WRITE(*,*) ' Wieviele Teilchen <', NMAX
        READ (*,*) N
        IF(N.GT.NMAX) N = NMAX
        CALL ANFANG(X,V,Z,N)
C
        WRITE(*,*) ' Wieviele Kollisionen?'
        READ (*,*) KEND
	WRITE(*,*) ' Groesse der "Bins"? [0.1--0.2]'
	READ (*,*) DV
	NBIN = NB
	DO K = 1, KEND
	 DO M = 1, N
          CALL STOSS(X,V,Z,N)
	 END DO
         CALL HISTOGR(V,DV,HISTO,NBIN,N)
	END DO
C
C	Kontrolliert die Konstanz der kinetischen Energie
	WRITE(*,*) ' Temperatur', TEMPER(V(1),Z(1),N)
C	Erzeugt nun einen mit "Gnuplot" lesbare Datei fort.11
	WRITE(11,'(A37)') '# Histogramm fuer die Geschwindigkeit'
	DO M = -NBIN, NBIN
         IF(HISTO(M).GT.0) THEN
         SIG = HISTO(M)
         SIG = DSQRT(SIG)
	 WRITE(11,'(F10.6,I7,F10.2)') DV*M, HISTO(M), SIG
         END IF
	END DO
        WRITE(6,*) ' Geschwindigkeitsverteilung --> fort.11'

        STOP
        END
C       --------------------------
        SUBROUTINE ANFANG(X,V,Z,N)
C
C       Initialisiert Positionen X und Geschwindigkeiten V der N Partikel
C       mit Masse Z und berechnet die kinetische Energie (EKIN ~ Temperatur).
C
        IMPLICIT NONE
        INTEGER N, M
        REAL*8 EKIN, SKAL, DM, X(0:N+1), V(0:N+1), Z(0:N+1)
        REAL RAND
C
         EKIN = 0.0
C	 Mit der Massendiffererenz
         DM   = 1.0D0/3.0D0
C	 sind die Massen Z = 1 +- DM: 4/3, 2/3, 4/3 ...
C	 Linke Wand:
         X(0) = 0.0D0
         V(0) = 0.0D0
         Z(0) = 0.0D0
C
C	Positionen: (M-1)/N < X(M) < M/N
C	Geschwindigkeiten: -1.5 < V(M) < +1.5
        DO M  = 1, N
         X(M) = (M - RAND())/N
         V(M) = 3.0*RAND() - 1.5
         Z(M) = 1.0D0 + DM
         EKIN = EKIN + Z(M)*V(M)**2
         DM   = -DM
        END DO
C
C	 Rechte Wand:
         X(N+1) = 1.0D0
         V(N+1) = 0.0D0
         Z(N+1) = 0.0D0
C       Skalierung der Geschwindigkeiten,
C	damit die Temperatur T = 2*Ekin/N = 1 ist:
         EKIN = EKIN/N
         SKAL = 1.0D0/DSQRT(EKIN)
         EKIN = 0.0D0
        DO M  = 1, N
         V(M) = V(M)*SKAL
         EKIN = EKIN + Z(M)*V(M)**2
	 Z(M) = 1.0D0/Z(M)  ! Reziproke Masse besser fuer die Stossberechnung
        END DO
         EKIN = EKIN/N
C	Kontrolle: Temperatur EKIN sollte 1.0 sein!
        WRITE(6,*) ' 2*Ekin/N =', EKIN
        RETURN
	END
C       -----------------------------
        SUBROUTINE STOSS(X,V,Z,N)
C
C       Berechnet die kleinste Zeit TAU bis zum n"achstem Stoss,
C	dh. TAV = DX/DV und das kleinste positive TAV ist dann TAU.
C	DV und DX sind die Geschwindigkeitsdifferenzen und Abstaende
C	benachbarter Teilchen.
C	Berechnet die Geschwindigskeitaenderung zweier stossender
C	Teilchen und die Bewegung aller Teilchen.
C
        IMPLICIT NONE
        INTEGER N, M, EINS, ZWEI, MSTOSS, M2
        REAL*8 X(0:N+1), V(0:N+1), Z(0:N+1), TAU, TAV, DV, DX, VNEU(2)
C
	 TAU = 1000.0D0
	 MSTOSS = 0
	DO M = 1, N+1
	 DV = V(M-1) - V(M)
	 IF(DV.GT.0.0D0) THEN
	  DX  = X(M) - X(M-1)
	  TAV = DX/DV
	  IF(TAV.LT.TAU) THEN
	   TAU = TAV
	   MSTOSS = M
	  END IF
	 END IF
	END DO
C
C	Bestimmt nun die neuen Positionen der Partikel
	DO M = 1, N
	 X(M) = X(M) + V(M)*TAU
	END DO
	IF(X(MSTOSS-1).GT.X(MSTOSS)) X(MSTOSS) = X(MSTOSS-1)
C
C	und dann die neuen Geschwindigkeiten nach dem Stoss mit
C	v_neu = [2M*V_alt + (m - M)*v_alt]/(M + m)
C	V_neu = [2m*v_alt + (M - m)*V_alt]/(M + m)
C	oder
C	v_neu = [2*V_alt/m + (1/M - 1/m)*v_alt/(1/M + 1/m).
C
C	Die Partikel MSTOSS-1 und MSTOSS kollidieren. Auch die Reflexion
C	an den Waenden ist beruecksichtigt. Denn fuer die Wand ist 1/M =
C	Z(0) = Z(N+1) = 0.0, so dass V(1) --> -V(1) & V(N) --> -V(N), wobei
C	die Geschwindigkeit der Wand V(0) = V(N+1) = 0 sein muss.
C
	EINS =  1
	ZWEI =  1
	DO M = MSTOSS-1, MSTOSS
	 M2 = M + EINS
	 VNEU(ZWEI) = (Z(M)+Z(M))*V(M2) + (Z(M2)-Z(M))*V(M)
	 VNEU(ZWEI) = VNEU(ZWEI)/(Z(M)+Z(M2))
	 EINS = -EINS
	 ZWEI = ZWEI + 1
	END DO
	 V(MSTOSS-1)= VNEU(1)
	 V(MSTOSS)  = VNEU(2)
C
	RETURN
	END
C       ------------------------------------
        SUBROUTINE HISTOGR(V,DV,HISTO,NBIN,N)
C	Erzeugt aus den Geschwindigkeiten V(1), V(2), ... V(N)
C	ein Histogramm von -NBIN bis +NBIN
C
	REAL*8 V(0:N+1), DV
	INTEGER M, N, NBIN, NB, HISTO(-NBIN:NBIN)
C
	DO M = 1, N, 2
	NB = INT(0.5D0 + DABS(V(M))/DV)
	IF(NB.GT.NBIN) NB=NBIN
	IF(V(M).LT.0.0D0) NB = -NB
	HISTO(NB) = HISTO(NB) + 1
	END DO
C
	RETURN
	END
C       ---------------------------------
        FUNCTION TEMPER(V,Z,N)
C	Berechnet die kinetische Energie pro Partikel
	INTEGER M, N
	REAL*8 V(N), Z(N), EKIN
C
	 EKIN = 0.0D0
	DO M = 1, N
	 EKIN = EKIN + V(M)**2/Z(M)
	END DO
	TEMPER = EKIN/N
C
	RETURN
	END
Compiliere mit g77 %
