MODULE konstanten
!In diesem Modul befinden sich alle verwendeten Konstanten. a,b,s1 und s2 sind die Konstanten
!des Potentials/der Kraft. gamma ist ein Maß für die Stärke der Reibung. kb ist die Boltzmann-
!konstante in eV, pi die Kreiszahl und N die Anzahl der Teilchen
IMPLICIT NONE
SAVE
DOUBLE PRECISION, PARAMETER :: a=5.d0,b=0.27d0,s1=1.3d0,s2=3.0d0,masse=28.d-4,gamma=10.d0&
			       &,kb=8.6173324d-5,pi=3.14159265359d0
INTEGER, PARAMETER :: N = 6*6*6
END MODULE
!*********************************************************************************************
PROGRAM newtonrk4
USE konstanten
!Programm zum Lösen der Newtonschen Bewegungsgleichung für ein thermodynamisches Problem.
!In dem Vektor y(1,.....,6*N) stehen jeweils für beispielsweise dem ersten Teilchen:
!y(1)=x	y(2)=y y(3)=z y(4)=p_x y(5)=p_y y(6)=p_z
!sodass im Abstand von 6 Einträgen von jedem Teilchen die Koordinaten und Impulse stehen.
!In diesem Programm werden die Impulse, nachdem der Grundzustand erreicht wurde per Monte 
!Carlo über die Maxwell-Boltzmann Verteilung verteilt. Die Längeneinheit beträgt Angström (A),
!die Energieeinheit eV, Kraft eV/A, Zeit ps. die Einheit der Masse wird in eV ps^2/A^2 angege-
!ben, wobei 1u =1.036426866d-4 eV ps^2/A^2. Die Teilchen sollen eine Masse von 28u besitzen.
!Benutzt wird eine Masse von masse = 28.d-4 eV ps^2/A^2
IMPLICIT NONE
INTEGER :: i,j,k,anzahl,iostat
DOUBLE PRECISION :: zeit,temperatur,dt,zeitende,ort,energie,e0,r1,kinetische, potentielle&
		    &,abweichung,ortsbetrag,nenner2,rw,phiw,thetaw
!dt ist die Schrittweite,ort gibt die Ausdehnung des Systems, in energie steckt die kin. Ener-
!gie nach Ausführung von montecarlo. e0 gibt nach E_kin=3/2 kb*T die bei einem idealen Gas zu
!erwartete kinetische Energie. Mit r1 wird der Abstand der Teilchen zueinander ausgerechnet. 
!abweichung gibt die gemittelte,betragliche Abweichung vom Grundzustand. ortsbetrag ist
!ein dummy für die Berechnung der Ausdehnung des Systems. rw,phiw und thetaw sind die Kugelko-
!ordinaten der Winkelgeschwindigkeit des Systems.
DOUBLE PRECISION, EXTERNAL :: kraft 
!Subroutine zur berechnung der Ableitungen der Koordinaten
DOUBLE PRECISION, DIMENSION(6*N) :: y = 0.d0, r0=0.d0,r01=0.d0
!y beinhaltet die Koordinaten und Impulse.In r0 ist der Grundzustand zu finden.r01 ist dummy.
DOUBLE PRECISION, DIMENSION(3) :: schwerp = 0.d0,drehimp=0.d0,mitte=0.d0,winkelg,winkelg1=0.d0
!Mit schwerp wird der Schwerpunktsimpuls berechnet, mit mitte der Massenschwerpunkt. winkelg1
!ist für Zwischenrechnungen da.
REAL :: zeit1,zeit2
DOUBLE PRECISION, DIMENSION(3,3) :: drehz,drehy
!Drehmatrizen zur Rotation des Grundzustands
DOUBLE PRECISION, DIMENSION(3*N) :: l=0.d0

!Legen der Teilchen auf ein kubisches Gitter. Nur zu Beginn nötig gewesen.
!do i=0,5		!legt die Teilchen auf ein kubisches Gitter zwischen
!  do j=0,5		!den Koordinaten (0,0,0) und (25.8,25.8,25.8)
!    do k=0,5	
!	   y(6*k+1+6*6*j+6*6*6*i)=dble(k)*4.3d0
!      y(6*k+2+6*6*j+6*6*6*i)=dble(j)*4.3d0
!      y(6*k+3+6*6*j+6*6*6*i)=dble(i)*4.3d0
!    end do
!  end do
!end do

open(unit=7,file="grundzustand.txt",status="old",iostat=iostat,action="read")
if (iostat == 0) then
  do i=1,6*N
    read(7,*) y(i)
  end do
else
  write(*,*) 'es ist fehler',iostat,'aufgetreten,Grundzustand'
  stop
end if
close(unit=7,iostat=iostat)
write(*,*) iostat
r0 = y
kinetische = 0.d0
potentielle = 0.d0
energie = 0.d0
temperatur = 4700.d0			
zeit = 0.d0
dt = 3.d-3 				!1.d-3 bei 8000K noch gut
zeitende = dt*20000.d0
anzahl = int((zeitende-zeit)/dt)
!Verteilen der Impulse
call montecarlo(temperatur,y)
do i=0,N-1
  energie = energie + (y(6*i+4)*y(6*i+4)+y(6*i+5)*y(6*i+5)+y(6*i+6)*y(6*i+6))/(2.d0*masse)
end do
e0 = 3.0d0/2.0d0 *kb*temperatur
write(*,*) energie/(e0*N)
!Bestimmung Schwerpunktsimpuls
do i=0,N-1
  schwerp(1) = schwerp(1) + y(6*i+4)
  schwerp(2) = schwerp(2) + y(6*i+5)
  schwerp(3) = schwerp(3) + y(6*i+6)
end do
schwerp = schwerp/N
write(*,*) "Schwerpunktsimpuls",schwerp
!Transformation ins Schwerpunktsystem
do i=0,N-1
  y(6*i+4) = y(6*i+4) - schwerp(1)
  y(6*i+5) = y(6*i+5) - schwerp(2)
  y(6*i+6) = y(6*i+6) - schwerp(3)
end do
!Massenmittelpunkt
do i=0,N-1
  mitte(1) = mitte(1) + y(6*i+1)
  mitte(2) = mitte(2) + y(6*i+2)
  mitte(3) = mitte(3) + y(6*i+3)
end do
mitte = mitte/N
!Verschieben des Schwerpunkts auf 0 ->Relativkoordinaten
do i=0,N-1
  r0(6*i+1) = r0(6*i+1)-mitte(1)
  r0(6*i+2) = r0(6*i+2)-mitte(2)
  r0(6*i+3) = r0(6*i+3)-mitte(3)
end do
call drehimpuls(y,winkelg)
!Winkelgeschwindigkeit in Kugelkoordinaten
rw = Sqrt(winkelg(1)**2 +winkelg(2)**2 +winkelg(3)**2)
phiw = atan2(winkelg(2),winkelg(1))
thetaw = acos(winkelg(3)/rw)
!Ausgabe der Kugelkoordinaten der Winkelgeschwindigkeit und der Schrittweite
open(unit=33,file="winkel.txt",status="old",iostat=iostat,action="write")
if (iostat == 0) then
  write(33,*) rw
  write(33,*) phiw
  write(33,*) thetaw
  write(33,*) dt
else
  write(*,*) 'es ist fehler',iostat,'aufgetreten, winkel.txt'
  stop
end if
close(unit=33,iostat=iostat)

!Drehmatrizen
drehz=reshape((/cos(-phiw),sin(-phiw),0.d0,-sin(-phiw),cos(-phiw),0.d0,0.d0,0.d0,1.d0 /)&
      &,(/3,3/))
drehy=reshape((/cos(-thetaw),0.d0,-sin(-thetaw),0.d0,1.d0,0.d0,sin(-thetaw),0.d0&
      &,cos(-thetaw)/),(/3,3/))
!Test, ob Rotationen korrekt, wenn korrekt: so steht w in z-Richtung; x- u. y-Komponente = 0
do j=1,3
  do k=1,3
    winkelg1(j) = winkelg1(j)+drehz(j,k)*winkelg(k)
  end do
end do
write(*,*) "winkelg1",winkelg1
winkelg = 0.d0
do j=1,3
  do k=1,3
    winkelg(j) = winkelg(j)+drehy(j,k)*winkelg1(k)
  end do
end do
write(*,*) "Winkelgeschwindigkeit transformiert", winkelg
!Rotation des Grundzustandsvektors
do i=0,N-1
  do j=1,3
    do k=1,3
      r01(6*i+j) = r01(6*i+j)+drehz(j,k)*r0(6*i+k)
    end do
  end do
end do
r0 = 0.d0
do i=0,N-1
  do j=1,3
    do k=1,3
      r0(6*i+j) = r0(6*i+j)+drehy(j,k)*r01(6*i+k)
    end do
  end do
end do
!Rotation der Orts- und Impulskoordinaten von y
r01 = 0.d0
do i=0,N-1
  do j=1,3
    do k=1,3
      r01(6*i+j) = r01(6*i+j)+drehz(j,k)*y(6*i+k)
      r01(6*i+j+3) = r01(6*i+j+3)+drehz(j,k)*y(6*i+k+3)
    end do
  end do
end do
y = 0.d0
do i=0,N-1
  do j=1,3
    do k=1,3
      y(6*i+j) = y(6*i+j)+drehy(j,k)*r01(6*i+k)
      y(6*i+j+3) = y(6*i+j+3)+drehy(j,k)*r01(6*i+k+3)
    end do
  end do
end do

!l enthält die Kugelkoordinaten des gedrehten Grundzustands (r,phi,theta),
!sodass man jetzt, da w in r-Richtung steht, auf Phi nach jedem Iterationsschritt
!rw*dt aufaddieren kann, sich der Grundzustand also mitdreht.
do i=0,N-1
  l(3*i+1) = sqrt(r0(6*i+1)**2 +r0(6*i+2)**2 +r0(6*i+3)**2)	!Abstand zu Mittelpunkt
  l(3*i+2) = atan2(r0(6*i+2),r0(6*i+1))				!Winkel Phi
  l(3*i+3) = acos(r0(6*i+3)/l(3*i+1))				!Winkel Theta
end do

call cpu_time(zeit1)
open(unit=44,file="drehimpuls.txt",status="old",iostat=iostat,action="write")
open(unit=8,file="test.txt",status="old",iostat=iostat,action="write")
open(unit=23,file="trajektorien.txt",status="old",iostat=iostat,action="write")
if (iostat == 0) then
do i=1, anzahl
    zeit = dble(i-1)*dt
    kinetische = 0.d0
    potentielle = 0.d0
    ort = 0.d0
    nenner2 = 0.d0
    abweichung = 0.d0
    !Berechnung der potentiellen Energie
    do j=0,N-2
      do k=j+1,N-1
        r1 = sqrt((y(6*j+1)-y(6*k+1))**2 +(y(6*j+2)-y(6*k+2))**2 +(y(6*j+3)-y(6*k+3))**2)
        potentielle = potentielle + a*exp(-0.5d0*(r1/s1)**2)-b*exp(-0.5d0*(r1/s2)**2) 
      end do
    end do
    !Berechnung der Ausdehnung des Arrangements
    do j=0,N-1
      ort = ort + sqrt(y(6*j+1)**2 +y(6*j+2)**2 +y(6*j+3)**2)
    end do
    !Berechnung der kinetischen Energie
    do j=0,N-1
      kinetische = kinetische + (y(6*j+4)*y(6*j+4)+y(6*j+5)*y(6*j+5)+&
      		   &y(6*j+6)*y(6*j+6))/(2.d0*masse)
    end do
    !Berechnung der Abweichung von dem Grundzustand
    do j=0,N-1
      ortsbetrag = Sqrt((y(6*j+1)-r0(6*j+1))**2 &
      		   &+ (y(6*j+2)-r0(6*j+2))**2 +(y(6*j+3)-r0(6*j+3))**2)
      if (ortsbetrag < 15.d0) then
        abweichung = abweichung + ortsbetrag
        nenner2 = nenner2 + 1.d0
      end if
    end do
    write(8,'(1x,g11.6,6(1x,G22.15))') zeit,kinetische,potentielle&
    		&,potentielle+kinetische,ort/n,abweichung/nenner2
    call rungekutta4(dt,zeit,y,kraft)
!**********************Korrektur des Drehimpulses*********************************************
    !Addition des durch die Winkelgeschwindigkeit veränderten Winkels in Phi
    do j=0,N-1
      l(3*j+2) = l(3*j+2) + rw*dt
    end do
    !Umschreiben in kartesische Koordinaten
    do j=0,N-1
      r0(6*j+1) = l(3*j+1)*sin(l(3*j+3))*cos(l(3*j+2))
      r0(6*j+2) = l(3*j+1)*sin(l(3*j+3))*sin(l(3*j+2))
      r0(6*j+3) = l(3*j+1)*cos(l(3*j+3))
    end do    
!Nach 50 Iterationsschritten sollen die Trajektorien von 5 Teilchen ausgegeben werden
    if (modulo(i,50) == 0) then
      write(23,'(3(1x,G22.15))') y(99*6+1),y(99*6+2),y(99*6+3)
      write(23,'(3(1x,G22.15))') y(98*6+1),y(98*6+2),y(98*6+3)
      write(23,'(3(1x,G22.15))') y(97*6+1),y(97*6+2),y(97*6+3)
      write(23,'(3(1x,G22.15))') y(93*6+1),y(93*6+2),y(93*6+3)
      write(23,'(3(1x,G22.15))') y(92*6+1),y(92*6+2),y(92*6+3)
    end if   
!Nach 200 Iterationsschritten soll der Zustandsvektor ausgegeben werden    
    if (modulo(i,200)== 0) then
      open(unit=9,file="endzustand.txt",status="old",iostat=iostat,action="write")
      if (iostat == 0) then
  	do j=1,6*N
          write(9,*) y(j)
        end do
      else
        write(*,*) 'es ist fehler',iostat,'aufgetreten,Zustandsvektor'
        stop
      end if
      close(unit=9,iostat=iostat)
      write(*,*) iostat
      schwerp = 0.d0
      do j=0,N-1
        schwerp(1) = schwerp(1) + y(6*j+4)
        schwerp(2) = schwerp(2) + y(6*j+5)
        schwerp(3) = schwerp(3) + y(6*j+6)
      end do
      write(*,*) 'Dies ist ein Test, ob der Schwerpunktsimpuls verschwunden ist'
      write(*,*) schwerp/N
      drehimp = 0.d0
      do j=0,N-1
        drehimp(1) = drehimp(1)+(y(6*j+2))*y(6*j+6)-(y(6*j+3))*y(6*j+5)
        drehimp(2) = drehimp(2)+(y(6*j+3))*y(6*J+4)-(y(6*j+1))*y(6*j+6)
        drehimp(3) = drehimp(3)+(y(6*j+1))*y(6*j+5)-(y(6*j+2))*y(6*j+4)
      end do		
      write(*,*) "Drehimpuls",drehimp
      write(44,*) drehimp
    end if
!Ende der Ausgabe
  end do
else
  write(*,*) 'es ist fehler',iostat,'aufgetreten,test und trajektorie'
  stop
end if
close(unit=8,iostat=iostat)
write(*,*) iostat

open(unit=9,file="endzustand.txt",status="old",iostat=iostat,action="write")
if (iostat == 0) then
  do i=1,6*N
    write(9,*) y(i)
  end do
else
  write(*,*) 'es ist fehler',iostat,'aufgetreten,endzustand'
  stop
end if
close(unit=9,iostat=iostat)
write(*,*) iostat

call cpu_time(zeit2)
write(*,*) zeit2-zeit1

END PROGRAM newtonrk4
!*********************************************************************************************
!Subroutine für die Integration der DGL
!In F müssen an der i-ten Stelle die Ableitungen von y(i) nach T stehen.Die Integration 
!erfolgt über die Variable T. H ist die Schrittweite in T.
SUBROUTINE RUNGEKUTTA4(H,T,Y,F)
USE konstanten
IMPLICIT NONE
DOUBLE PRECISION, EXTERNAL :: F
DOUBLE PRECISION, DIMENSION(6*N), intent(inout)	:: y
DOUBLE PRECISION, INTENT(INOUT)	 :: T
DOUBLE PRECISION, INTENT(IN)	 :: H
DOUBLE PRECISION, DIMENSION(6*N) :: k1,k2,k3,HF
DOUBLE PRECISION		 :: T1,hh,h6
INTEGER				 :: i
hh = h*0.5d0
h6 = h/6.0d0
call F(T,Y,k1)
T1 = T + hh
HF = y + hh*k1
call F(T1,HF,k2)
HF = y + hh*k2
call F(T1,HF,k3)
HF = y + k3*H
k2 = k2 + k3
T = T + H
call F(T,HF,k3)
y = y + h6 *(k1+2.0d0*k2+k3)
!Falls ein Teilchen nach einem Iterationsschritt in einer Dimension außerhalb des Intervalls 
![-25,25] liegt, so wird es an der anderen Seite wieder eingeführt. Die Ausgabe "sprung" zeigt
!während der Berechnungen an, ob ein Teilchen außerhalb des Kastens lag und kann mit etwaigen 
!Sprüngen der Gesamtenergie identifiziert werden.
do i=0,N-1
  if (y(6*i+1)>25.d0) then
    y(6*i+1) = -25.d0
    write(*,*) "sprung"
  end if
  if (y(6*i+2)>25.d0) then
    y(6*i+2) = -25.d0
    write(*,*) "sprung"
  end if
  if (y(6*i+3)>25.d0) then
    y(6*i+3) = -25.d0
    write(*,*) "sprung"
  end if
  if (y(6*i+1)<-25.d0) then
    y(6*i+1) = 25.d0
    write(*,*) "sprung"
  end if
  if (y(6*i+2)<-25.d0) then
    y(6*i+2) = 25.d0
    write(*,*) "sprung"
  end if
  if (y(6*i+3)<-25.d0) then
    y(6*i+3) = 25.d0
    write(*,*) "sprung"
  end if
end do
RETURN
END SUBROUTINE RUNGEKUTTA4
!*********************************************************************************************
!Subroutine enthält die Bewegungsgleichungen/ die Ableitungen für den Vektor y(i)
SUBROUTINE KRAFT(T,Y,HF)
USE konstanten
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(6*N), INTENT(IN) :: Y
DOUBLE PRECISION, INTENT(IN) :: T
DOUBLE PRECISION, DIMENSION(6*N), INTENT(INOUT) :: HF
INTEGER	:: i,j
DOUBLE PRECISION, PARAMETER :: a1 = a/(s1*s1),a2 = b/(s2*s2),a3=-0.5d0/(s1*s1)&
			       &,a4=-0.5d0/(s2*s2)
DOUBLE PRECISION :: r, temp !r beinhaltet |r-r'|, temp wie im Abschnitt 
HF = 0.d0		    !"Potential und effektive Wechselwirkung"
!Gleichung für dx/dt = p/m, Ableitungen der Ortskoordinaten
do i=0,N-1
  HF(6*i+1) = y(6*i+4)/masse
  HF(6*i+2) = y(6*i+5)/masse
  HF(6*i+3) = y(6*i+6)/masse
end do
!Gleichung für dp/dt= F, Ableitungen der Impulskoordinaten
do i=0,(N-2)
  do j=i+1,(N-1)
    r = sqrt((y(6*i+1)-y(6*j+1))**2 + (y(6*i+2)-y(6*j+2))**2 + (y(6*i+3)-y(6*j+3))**2)
    temp= a1*exp(a3*r*r)-a2 *exp(a4*(r)*(r))
    HF(6*i+4) = HF(6*i+4) + temp*(y(6*i+1)-y(6*j+1))
    HF(6*j+4) = HF(6*j+4) - temp*(y(6*i+1)-y(6*j+1))
    HF(6*i+5) = HF(6*i+5) + temp*(y(6*i+2)-y(6*j+2))
    HF(6*j+5) = HF(6*j+5) - temp*(y(6*i+2)-y(6*j+2))
    HF(6*i+6) = HF(6*i+6) + temp*(y(6*i+3)-y(6*j+3))
    HF(6*j+6) = HF(6*j+6) - temp*(y(6*i+3)-y(6*j+3))
  end do
end do
!Reibungsterm: dp/dt = F "-gamma p". Nur am Anfang nötig gewesen
!do i=0,N-1
!  HF(6*i+4) = HF(6*i+4) - gamma*y(6*i+4)
!  HF(6*i+5) = HF(6*i+5) - gamma*y(6*i+5)
!  HF(6*i+6) = HF(6*i+6) - gamma*y(6*i+6)
!end do
RETURN
END SUBROUTINE KRAFT

!*********************************************************************************************
!Diese Subroutine dient der Impulsvergabe nach der Maxwell Boltzmann Verteilung mit Hilfe 
!einer Monte Carlo Methode.
SUBROUTINE MONTECARLO(T,y)
USE konstanten
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: T
DOUBLE PRECISION, DIMENSION(6*N), INTENT(INOUT)	:: y
DOUBLE PRECISION	:: vert,test,theta,phi,betrag,random@,vertmax,betragmax
INTEGER	:: i,j
!T ist die Temperatur, vert ist der zufällig Maxwell Boltzmann Verteilungswert, theta ist der
!zufällig erhaltene Winkel theta des Impulses, gleiches gilt für phi und betrag.
!vertmax gibt eine Funktion für das Maximum der Verteilung. Diese Funktion ist ein Fit der
!Verteilungsmaxima bei verschiedenen Temperaturen und wurde mit Mathematica durchgeführt.
!betragmax ist der Impulsbetrag bei dem die Verteilung einen Wert von 10^-5 annimmt.

vertmax = 181.69723638904276d0- 49036.1596250d0/(T*T) + &   !Gleichungen gelten nur bis 10000K
	&12007.293973803038d0/T - 0.2326018089098843d0 *T +& 
 	&0.00020143492718019792d0 *T*T - 9.889627385975435d-8 *T*T*T + &
 	&2.8465004602038517d-11 *T*T*T*T - 4.898019128141887d-15 *T*T*T*T*T + &
 	&4.950284303906514d-19 *T*T*T*T*T*T - 2.7071514724946854d-23 *T*T*T*T*T*T*T + &
 	&6.176463667836648d-28 *T*T*T*T*T*T*T*T					
betragmax = 0.0017642832638790238d0+ 0.0012577681693344063d0*Sqrt(T)+0.002d0
write(*,*) vertmax,betragmax
do i=0,N-1
  j = 0
100 j = j+1
  vert = random@()*vertmax	!random@() generiert eine Zufallszahl aus dem Intervall (0,1]
  theta = pi*random@()
  phi = 2.d0*pi*random@()
  betrag = betragmax*random@()
  test = betrag*betrag/(2.d0*masse*kb*T)
  if (test > 397.d0) then	!Überprüfung, ob Argument in E-Funktion nicht zu negativ
    goto 100			!Exp(-397)=10^-173=0, also vert größer -> neue Zufallszahlen!
  else if (vert > (2.d0*kb*T*masse*pi)**(-3/2) *4.d0*pi*betrag*betrag*exp(-test)) then
    if (j > 2.d6) then
      write(*,*) "Mehr als 2.d6 Überschreitungen!"
      stop	!Falls über 2*10^6 Impulsbeträge eines einzigen Teilchens nicht der Verteilung
    end if	!entsprechen, so wird das Programm abgebrochen. (Kann höher gesetzt werden)
    goto 100
  end if
  y(6*i+4) = betrag*sin(theta)*cos(phi)	!Überführung in kartesische Koordinaten
  y(6*i+5) = betrag*sin(theta)*sin(phi)
  y(6*i+6) = betrag*cos(theta)
end do
RETURN
END SUBROUTINE MONTECARLO
!*********************************************************************************************
!Diese Subroutine berechnet den Drehimpuls und anschließend über das Gaußsche Eliminations-
!verfahren die Winkelgeschwindigkeit die der Körper besitzt. Diese ist nötig um den Grund-
!zustand des Körpers entsprechend mitzudrehen, um die Abweichung von diesem auszurechnen,
!außerdem kann mit der bekannten Winkelgeschwindigkeit das Koordinatensystem zur Darstellung 
!der Trajektorien der 5 Teilchen entsprechend mitgedreht werden.
SUBROUTINE Drehimpuls(y,winkelg)
USE konstanten
IMPLICIT NONE
DOUBLE PRECISION, DIMENSION(6*N),INTENT(INOUT) :: y
DOUBLE PRECISION, DIMENSION(3), INTENT(INOUT)  :: winkelg
DOUBLE PRECISION, DIMENSION(3)	 :: mitte=0.d0, drehimp=0.d0,vektor=0.d0
DOUBLE PRECISION, DIMENSION(3,3) :: traeg=0.d0
DOUBLE PRECISION, DIMENSION(3,4) :: temp1,temp2
INTEGER :: i,j

winkelg = 0.d0
!Drehimpulskorrektur, zunächst Bestimmung des Schwerpunkts
do i=0,N-1
  mitte(1) = mitte(1) + y(6*i+1)
  mitte(2) = mitte(2) + y(6*i+2)
  mitte(3) = mitte(3) + y(6*i+3)
end do
mitte = mitte/N
write(*,*) mitte
!Verschiebung des Schwerpunkts auf (0,0,0), die Ortsanteile von y sind nun die Relativkoord.
do i=0,N-1
  y(6*i+1) = y(6*i+1)-mitte(1)
  y(6*i+2) = y(6*i+2)-mitte(2)
  y(6*i+3) = Y(6*i+3)-mitte(3)
end do
!Berechnung des Drehimpulses des zusammengesetzten Körpers in Bezug zum Schwerpunkt
do i=0,N-1
  drehimp(1) = drehimp(1)+y(6*i+2)*y(6*i+6)-y(6*i+3)*y(6*i+5)
  drehimp(2) = drehimp(2)+y(6*i+3)*y(6*i+4)-y(6*i+1)*y(6*i+6)
  drehimp(3) = drehimp(3)+y(6*i+1)*y(6*i+5)-y(6*i+2)*y(6*i+4)
end do			
write(*,*) drehimp
!Trägheitstensor
do i=0,N-1
  traeg(1,1) = traeg(1,1) + y(6*i+2)**2 +y(6*i+3)**2
  traeg(2,2) = traeg(2,2) + y(6*i+1)**2 +y(6*i+3)**2
  traeg(3,3) = traeg(3,3) + y(6*i+1)**2 +y(6*i+2)**2
  traeg(1,2) = traeg(2,1) - y(6*i+1)*y(6*i+2)
  traeg(1,3) = traeg(1,3) - y(6*i+1)*y(6*i+3)
  traeg(2,3) = traeg(2,3) - y(6*i+2)*y(6*i+3) 
end do
traeg(2,1) = traeg(1,2)
traeg(3,1) = traeg(1,3)
traeg(3,2) = traeg(2,3)
traeg = masse*traeg
!Gaußsches Eliminationsverfahren
do i=1,3
  do j=1,3
    temp1(i,j) = traeg(i,j)	!Es wird zunächst das lineare Gleichungssystem ine ine neue
    temp1(i,4) = drehimp(i)	!umgeschrieben I*w=L, I ist Matrixwert, w und L Vektor vertig
  end do			!Die erste Spalte des Trägheitstensors bildet die Koeffizien-
end do				!ten von w_1 usw. die vierte Spalte der neuen Matrix bildet  
temp2 = temp1			!der Drehimpuls
do i=1,4
  temp1(3,i) = temp2(3,i) - temp2(3,1)/temp2(2,1) *temp2(2,i)!Eliminierung der ersten
  temp1(2,i) = temp2(2,i) - temp2(2,1)/temp2(1,1) *temp2(1,i)!Einträge in Zeile 2 und 3
end do
temp2 = temp1
do i=1,4
  temp1(3,i) = temp2(3,i) - temp2(3,2)/temp2(2,2) *temp2(2,i)!Eliminierung des zweiten 
end do							     !Eintrags in Zeile 2
!Lösung der Winkelgeschwindigkeit
winkelg(3) = temp1(3,4)/temp1(3,3)		
winkelg(2) = (temp1(2,4)- winkelg(3)*temp1(2,3))/temp1(2,2)
winkelg(1) = (temp1(1,4)- winkelg(3)*temp1(1,3)- winkelg(2)*temp1(1,2))&
	     &/temp1(1,1)
write(*,*) "Drehimpuls",drehimp
!Test,ob mit errechneter Winkelgeschwindigkeit durch Multiplikation mit Trägheitstensor wieder
!der Bahndrehimpuls erhalten wird.
do i=1,3
  do j=1,3
      vektor(i) = vektor(i) + traeg(i,j)*winkelg(j)
  end do
end do
write(*,*) "berechnet", vektor
write(*,*) winkelg
RETURN
END SUBROUTINE drehimpuls