Warning: Undefined array key "HTTP_REFERER" in /var/www/vhosts/web61970.ssd-space.de/phi.pf-control.de/uni/download.php on line 13
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