c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_INFO_SP_1 c ========================= IMPLICIT NONE character*4 Nr parameter (Nr='Sp_1') INCLUDE 'mutrack$sourcedirectory:COM_LUNS.INC' INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC' INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_1.INC' c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - namelist /geometry/ xSpGrid1,xSpGrid2,dWires_Sp namelist /grid_info/ + Dx,Dy, imax,jmax, xmin,xmax,ymax c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - logical map_error COMMON /map_error/ map_error c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Einlesen der Mappen-Informationen: open (lunREAD,file=mappenName//'_1.INFO',defaultfile=mappenDir, + readonly,status='old') read (lunREAD,nml=grid_info) rewind (lunREAD) read (lunREAD,nml=geometry) close (lunREAD) c eingelesene imax und jmax um 1 reduzieren, da in 'MUTRACK' die Feldindizes c ab 0 laufen, bei 'RELAX3D' jedoch ab 1: imax = imax-1 jmax = jmax-1 c die geometrischen Daten fuer MUTRACK zusammenstellen: dSpiegel = xSpGrid2 - xSpGrid1 dist_Wires_Sp = 2. * real(jmax) * Dy c checken, ob der reservierte Speicherplatz ausreicht: if ((imax+1)*(jmax+1).NE.maxmem+1) then write(*,*)'----------------------------------------'// + '----------------------------------------' write(*,*) Nr//'-Mappe: ',mappenName write(*,*) ' BENOETIGTER Speicher: (imax+1)*(jmax+1) = ',(imax+1)*(jmax+1) write(*,*) ' RESERVIERTER Speicher: maxmem + 1 = ',maxmem + 1 write(*,*) if ((imax+1)*(jmax+1).GT.maxmem+1) then write(*,*) '=> reservierter Speicherplatz ist ungenuegend.' write(*,*) write(*,*) '=> ''maxmem'' in mutrack$sourcedirectory:MAP_DEF_'//Nr//'.INC angleichen,' write(*,*) ' dann Programm mit ''LINKMUV'' am DCL-Prompt neu kompilieren' write(*,*) ' und linken.' write(*,*) write(*,*) ' Mindestwert fuer ''maxmem'' ist ',(imax+1)*(jmax+1)-1 write(*,*) map_error = .true. endif endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_INFO_SP_2 c ========================= IMPLICIT NONE character*4 Nr parameter (Nr='Sp_2') INCLUDE 'mutrack$sourcedirectory:COM_LUNS.INC' INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC' INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_2.INC' c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - namelist /grid_info/ + Dx,Dy, imax,jmax, xmin,xmax,ymax c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - logical map_error COMMON /map_error/ map_error c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Einlesen der Mappen-Informationen: open (lunREAD,file=mappenName//'_2.INFO',defaultfile=mappenDir, + readonly,status='old') read (lunREAD,nml=grid_info) close (lunREAD) c eingelesene imax und jmax um 1 reduzieren, da in 'MUTRACK' die Feldindizes c ab 0 laufen, bei 'RELAX3D' jedoch ab 1: imax = imax-1 jmax = jmax-1 c checken, ob der reservierte Speicherplatz ausreicht: if ((imax+1)*(jmax+1).NE.maxmem+1) then write(*,*)'----------------------------------------'// + '----------------------------------------' write(*,*) Nr//'-Mappe: ',mappenName write(*,*) ' BENOETIGTER Speicher: (imax+1)*(jmax+1) = ',(imax+1)*(jmax+1) write(*,*) ' RESERVIERTER Speicher: maxmem + 1 = ',maxmem + 1 write(*,*) if ((imax+1)*(jmax+1).GT.maxmem+1) then write(*,*) '=> reservierter Speicherplatz ist ungenuegend.' write(*,*) write(*,*) '=> ''maxmem'' in mutrack$sourcedirectory:MAP_DEF_'//Nr//'.INC angleichen,' write(*,*) ' dann Programm mit ''LINKMUV'' am DCL-Prompt neu kompilieren' write(*,*) ' und linken.' write(*,*) write(*,*) ' Mindestwert fuer ''maxmem'' ist ',(imax+1)*(jmax+1)-1 write(*,*) map_error = .true. endif endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_INFO_SP_3 c ========================= IMPLICIT NONE character*4 Nr parameter (Nr='Sp_3') INCLUDE 'mutrack$sourcedirectory:COM_LUNS.INC' INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC' INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_3.INC' c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - namelist /grid_info/ + Dx,Dy, imax,jmax, xmin,xmax,ymax c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - logical map_error COMMON /map_error/ map_error c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Einlesen der Mappen-Informationen: open (lunREAD,file=mappenName//'_3.INFO',defaultfile=mappenDir, + readonly,status='old') read (lunREAD,nml=grid_info) close (lunREAD) c eingelesene imax und jmax um 1 reduzieren, da in 'MUTRACK' die Feldindizes c ab 0 laufen, bei 'RELAX3D' jedoch ab 1: imax = imax-1 jmax = jmax-1 c checken, ob der reservierte Speicherplatz ausreicht: if ((imax+1)*(jmax+1).NE.maxmem+1) then write(*,*)'----------------------------------------'// + '----------------------------------------' write(*,*) Nr//'-Mappe: ',mappenName write(*,*) ' BENOETIGTER Speicher: (imax+1)*(jmax+1) = ',(imax+1)*(jmax+1) write(*,*) ' RESERVIERTER Speicher: maxmem + 1 = ',maxmem + 1 write(*,*) if ((imax+1)*(jmax+1).GT.maxmem+1) then write(*,*) '=> reservierter Speicherplatz ist ungenuegend.' write(*,*) write(*,*) '=> ''maxmem'' in mutrack$sourcedirectory:MAP_DEF_'//Nr//'.INC angleichen,' write(*,*) ' dann Programm mit ''LINKMUV'' am DCL-Prompt neu kompilieren' write(*,*) ' und linken.' write(*,*) write(*,*) ' Mindestwert fuer ''maxmem'' ist ',(imax+1)*(jmax+1)-1 write(*,*) map_error = .true. endif endif if (map_error) then write(*,*)'----------------------------------------'// + '----------------------------------------' STOP endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_MAP_SP_1 c ======================== IMPLICIT NONE character*(*) Nr parameter (Nr='_1') INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_1.INC' INCLUDE 'mutrack$sourcedirectory:READ_MAP_SP.INC' END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_MAP_SP_2 c ======================== IMPLICIT NONE character*(*) Nr parameter (Nr='_2') INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_2.INC' INCLUDE 'mutrack$sourcedirectory:READ_MAP_SP.INC' END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE READ_MAP_SP_3 c ======================== IMPLICIT NONE character*(*) Nr parameter (Nr='_3') INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_3.INC' INCLUDE 'mutrack$sourcedirectory:READ_MAP_SP.INC' END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE INTEGRATIONSSTEP_RUNGE_KUTTA_SP(dt) c ============================================== IMPLICIT NONE SAVE character*(*) Nr parameter (Nr='Sp') c Diese Subroutine berechnet zu einem vorgegebenen Zeitschritt dt den c Integrationsschritt zweimal: einmal direkt mit dt und einmal ueber zwei c aufeinanderfolgende Schritte mit dt/2. (die beiden dt/2-Schritte werden c zuerst ausgefuehrt). c c Aus der Differenz der beiden Resultate wird eine Abschaetzung fuer den Fehler c des dt-Schrittes gewonnen, die dazu verwendet wird zu entscheiden, ob der c Integrationsschritt mit einem verkuerzten Zeitintervall wiederholt werden c muss, oder ob das Zeitintervall fuer den folgenden ausgedehnt werden kann. c c Die beiden Einzelergebnisse aus dem dt- und den beiden dt/2-Schritten, die c jeweils ueber Runge-Kutta-Rechnung vierter Ordnung erhalten werden, werden c zum Schluss noch zusammengenommen, um ein Resultat mit Genauigkeit fuenfter c Ordnung in dt zu erhalten. c c Der ganze Ablauf entspricht den Ausfuehrungen in Kapitel 15.2 der NUMERICAL c RECIPIES: 'Adaptive Stepsize Control for Runge-Kutta' (vgl. Referenz im c fileheader von 'MUTRACK.FOR') INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC' character*40 MappenName real dl_max COMMON /integration_Sp/ MappenName,dl_max integer i ! Zaehlvariable real dt,dt_half ! zeitl. Aenderung, halbe zeitl. Aenderung real EFeld0(3),EFeld1(3) ! elektr. Felder real x1(3),Dx1(3),Dx2(3) ! fuer Ortsintegration real v1(3),Dv1(3),Dv2(3) ! fuer Geschw.Integration real xDifferenz(3), vDifferenz(3) real maxErr_x,maxErr_v,maxErr ! fuer Fehlerbetrachtung real errCon, safety ! fuer Schrittweitenkontrolle real pShrink, pGrow ! fuer Schrittweitenkontrolle PARAMETER (errCon = 6.e-4, safety = .9) ! vgl. Referenz PARAMETER (pShrink = -.25, pGrow = -.2) ! errCon = (4./safety)**(1/pGrow) logical flag_dtSmall ! wenn bei dt < dtsmall der Fehler ! immer noch zu gross ist. integer returnCode_EFeld COMMON /returnCode_EFeld/ returnCode_EFeld ! 1: TESTORT HINTER DER MAPPE ! 2: Testort neben der Mappe ! 3: TESTORT VOR DER MAPPE c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - flag_dtSmall = .false. ! flag resetten if (dt.LT.dtsmall) dt = dtsmall c berechne EFeld am aktuellen Ort. Speichere in EFeld0, damit sie wiederverwendet c werden kann, falls mit kuerzerem Zeitschritt wiederholt werden muss: call EFeld_SP(x,EFeld0,*999) c............................................................................... 10 continue ! hier gehts wieder von vorne los, falls Zeitschritt dt ! abgeaendert werden muss. dt_half = dt / 2. c mache ersten dt/2 - Schritt: call SINGLESTEP_RUNGE_KUTTA_SP(dt_half,EFeld0,x,v, Dx1,Dv1 ,*999) c berechne EFeld bei x1: x1(1) = x(1) + Dx1(1) x1(2) = x(2) + Dx1(2) x1(3) = x(3) + Dx1(3) v1(1) = v(1) + Dv1(1) v1(2) = v(2) + Dv1(2) v1(3) = v(3) + Dv1(3) call EFeld_SP(x1,EFeld1,*999) c mache zweiten dt/2 - Schritt: call SINGLESTEP_RUNGE_KUTTA_SP(dt_half,EFeld1,x1,v1, Dx2,Dv2 ,*999) c Summiere Ergebnisse der beiden dt/2 -Schritte und speichere in Dx1, Dv1: Dx1(1) = Dx1(1) + Dx2(1) Dx1(2) = Dx1(2) + Dx2(2) Dx1(3) = Dx1(3) + Dx2(3) Dv1(1) = Dv1(1) + Dv2(1) Dv1(2) = Dv1(2) + Dv2(2) Dv1(3) = Dv1(3) + Dv2(3) c mache dt - Schritt: call SINGLESTEP_RUNGE_KUTTA_SP(dt,EFeld0,x,v, Dx2,Dv2 ,*999) c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c Fehlerbetrachtung und Berechnung des endgueltigen Ergebnisses: c Fehlerbetrachtung: c der groesste (absolute bzw. relative) Fehler im Ort soll kleiner als eps_x c sein, der groesste Fehler in der Geschwindigkeit kleiner als eps_v: c -> Bestimme den jeweils groessten Fehler der drei Komponenten des Ortes und c der Geschwindikeit (dh. die groesste Differenz der Aederungen): maxErr_x = 0. maxErr_v = 0. do i = 1, 3 xDifferenz(i) = Dx1(i)-Dx2(i) vDifferenz(i) = Dv1(i)-Dv2(i) if (log_relativ) then if (Dx1(i).NE.0.) + maxErr_x = Max( maxErr_x, Abs( xDifferenz(i)/Dx1(i) ) ) if (Dv1(i).NE.0.) + maxErr_v = Max( maxErr_v, Abs( vDifferenz(i)/Dv1(i) ) ) else maxErr_x = Max( maxErr_x, Abs( xDifferenz(i) ) ) maxErr_v = Max( maxErr_v, Abs( vDifferenz(i) ) ) endif enddo c - Skaliere den jeweils groessten relativen Fehler auf das jeweilige Epsilon: maxErr_x = maxErr_x / eps_x maxErr_v = maxErr_v / eps_v c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c der groessere der beiden reskalierten Fehler bestimmt, ob der Integrations- c schritt mit kleinerem Zeitintervall wiederholt werden muss, bzw. um welchen c Faktor das Zeitintervall fuer den naechsten Schritt vergroessert werden kann: c Liegt der Fehler ausserhalb des Toleranzbereiches und ist dt bereits jetzt c kleiner als dtsmall, so mache keinen neuen Versuch sondern akzeptiere als Not- c loesung den bestehenden Naeherungswert. Setze dt in diesem Fall als Default c fuer den kommenden Integrationsschritt auf dtsmall. Setze aber auch das flag c 'flag_dtsmall', damit gezaehlt werden kann, wie oft dieses Prozedur fuer ein c bestimmtes Teilchen angewendet werden muss. Ist dies zu oft der Fall, so brich c diese Trajektorienberechnung ganz ab (-> destiny = code_dtsmall). c (2. Teil erfolgt weiter unten) maxErr = Max(maxErr_x,maxErr_v) if (maxErr.GT.1.) then if (dt.LT.dtsmall) then ! Fehler immer noch zu gross, obwohl flag_dtsmall = .true. ! dtsmall schon unterschritten ist else !c Bestimme kuerzeren Zeitschritt fuer neuen Versuch (vgl. Referenz): dt = safety * dt * (maxErr**pShrink) goto 10 endif endif c Nimm die Ergebnisse aus dem dt-Schritt und den beiden dt/2-Schritten und c berechne damit den neuen Ort und die neue Geschwindigkeit mit Genauigkeit c fuenfter Ordnung in dt: x(1) = x(1) + Dx1(1) + xDifferenz(1) / 15. x(2) = x(2) + Dx1(2) + xDifferenz(2) / 15. x(3) = x(3) + Dx1(3) + xDifferenz(3) / 15. v(1) = v(1) + Dv1(1) + vDifferenz(1) / 15. v(2) = v(2) + Dv1(2) + vDifferenz(2) / 15. v(3) = v(3) + Dv1(3) + vDifferenz(3) / 15. c alten Zeitschritt addieren: t = t + dt c neuen Zeitschritt so gross wie sinnvoller weise moeglich machen: 3454 if (flag_dtSmall) then if (n_dtsmall.LT.maxBelowDtSmall) then dt = dtSmall ! fuer naechsten RK-Schritt n_dtsmall = n_dtsmall + 1 else destiny = code_dtsmall ! gib Teilchen verloren RETURN endif else if (maxErr.GT.errCon) then dt = safety * dt * (maxErr**pGrow) ! vgl. Referenz else dt = 4. * dt ! <- Vergroesserung des Zeitschritts max. um endif ! Faktor 4! ! pruefen, ob Maximallaenge fuer ersten Testschritt nicht ueberschritten ist: d if (log_confine) dt = min(dt,dl_max/sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))) endif RETURN c hier folgt der Code fuer 'returnCode_EFeld.NE.0': 999 continue if (returnCode_EFeld.EQ.1) then ! Testort hinter der Mappe destiny = code_durchSpiegel RETURN elseif (returnCode_EFeld.NE.0) then write(*,*) write(*,*) 'SINGLESTEP_RUNGE_KUTTA_'//Nr//': ' write(*,*) 'unallowed value of ''returnCode_EFeld'': ',returnCode_Efeld write(*,*) '-> STOP' write(*,*) STOP endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE SINGLESTEP_RUNGE_KUTTA_SP(dt,E0,x0,v0, Dx,Dv, *) c =========================================================== IMPLICIT NONE c Diese Subroutine berechnet bei vorgegebenem Zeitschritt einen einzelnen c Runge-Kutta-Integrationsschritt (4. Ordnung). c Die Vorgehensweise entspricht den Ausfuehrungen in Kapitel 15.1 der c NUMERICAL RECIPIES: 'Runge-Kutta Method'. Die Formeln sind zwar teilweise c etwas umgeschrieben, sind aber mathematisch alle aequivalent zu denen der c Referenz. c Zurueckgegeben werden die errechneten Orts- und Geschwindigkeitsaenderungen c anstatt direkt der neuen Werte, da sonst vor allem bei den Ortskoordinaten c Schwierigkeiten auftreten koennen, wenn in der Subroutine 'INTEGRATIONSSTEP_ c RUNGE_KUTTA' aus der Differenz der neuen Werte aus den beiden dt/2- und dem c dt-Schritt der Fehler abgeschaetzt werden soll (kleine Differenz moeglicher- c weise grosser Werte). real Beschl_Faktor COMMON /Beschl_Faktor_SP/ Beschl_Faktor real E0(3), x0(3), v0(3) ! Eingangsgroessen real E1(3), E2(3), E3(3) ! E-Felder an Testorten real v1(3), v2(3), v3(3) ! Geschwindigkeiten an Testorten real dt,dt_half,dt_sixth ! zeitl. Aenderung, dt/2, dt/6 real help, help_half, help_sixth ! Hilfsvariable, help/2, help/6 real xTest(3) ! Test-Orte real Dx(3), Dv(3) ! Ergebnisspeicher integer i ! Zaehlvariable c = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = dt_half = dt / 2. dt_sixth = dt / 6. help = Beschl_Faktor * dt help_half = help / 2. help_sixth = help / 6. do i = 1, 3 xTest(i) = x0(i) + v0(i) * dt_half v1(i) = v0(i) + E0(i) * help_half enddo call EFeld_SP(xTest,E1,*999) do i = 1, 3 xTest(i) = x0(i) + v1(i) * dt_half v2(i) = v0(i) + E1(i) * help_half enddo call EFeld_SP(xTest,E2,*999) do i = 1, 3 xTest(i) = x0(i) + v2(i) * dt v3(i) = v0(i) + E2(i) * help enddo call EFeld_SP(xTest,E3,*999) do i = 1, 3 Dx(i) = (v0(i) + 2.*(v1(i)+v2(i)) + v3(i)) * dt_sixth Dv(i) = (E0(i) + 2.*(E1(i)+E2(i)) + E3(i)) * help_sixth enddo RETURN 999 RETURN 1 END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE EFeld_SP(x,E,*) c ========================== IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:COM_KAMMER.INC' real x(3),E(3) ! Ort und Feldstaerke c Falls Testort ausserhalb der Mappe liegt oder Anode getroffen hat: integer returnCode_EFeld COMMON /returnCode_EFeld/ returnCode_EFeld ! 1: TESTORT HINTER DER MAPPE ! 2: Testort neben der Mappe ! 3: TESTORT VOR DER MAPPE c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Pruefe, in welchen Mappenbereich der Raumpunkt faellt und rufe entsprechende c Routine auf: c Mappe 1 ist die groebste, Mappe 2 etwas feiner, Mappe 3 ist die feinste. if (x(1).LT.0) then E(1) = 0. E(2) = 0. E(3) = 0. RETURN elseif (x(1).GT.xmax_Sp_1) then returnCode_Efeld = 1 RETURN 1 elseif (x(1).LT.xmin_Sp_2.OR.x(1).GT.xmax_Sp_2) then CALL EFeld_Sp_1(x,E) elseif (x(1).LT.xmin_Sp_3.OR.x(1).GT.xmax_Sp_3) then CALL EFeld_Sp_2(x,E) else CALL EFeld_Sp_3(x,E) endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE EFeld_SP_1(x,E) c ========================== IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_1.INC' real real_i,real_j,real_j_ ! x,r im Mappensystem in Gittereinheiten integer stuetzstelle_i(2) ! naechste Stuetzstellen in x- und integer stuetzstelle_j(2) ! r-Richtung real Abstand_i,Abstand_i_Betrag ! Entfernung zur naechsten Stuetzstelle real Abstand_j,Abstand_j_Betrag ! (in Gittereinheiten!) integer i,j, n, ihelp real x(3),E(3) ! Ort und Feldstaerke real E_(2) ! Hilfsspeicher fuer Feldberechnung integer zaehler_Sp_1 /0/ common /zaehler_Sp_1/ zaehler_Sp_1 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - zaehler_Sp_1 = zaehler_Sp_1 + 1 c Rechne in Gittereinheiten um und teste, ob Raumpunkt innerhalb der Potential- c mappe liegt: real_i = (x(1)-xmin)/Dx ! Umrechnung in Mappensystem findet bereits in real_j_ = x(2)/ Dy ! MUTRACK statt. c Mappe ist zyklisch in y-Richtung. j laeuft von 0 bis jmax vom 0. Gitterstab c zur Mitte zwischen 0. und 1. Gitterstab. c -> 1.) waehle aequivalentes real_j aus Bereich von Mitte zwischen -1. und 0. c Gitterstab bis Mitte zwischen 0. und 1. Gitterstab. c 2.) Nimm Absolutbetrag. n = nInt(real_j_/(2.*real(jmax))) real_j_ = real_j_ - real(n) * (2.*real(jmax)) real_j = abs(real_j_) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Berechne die Feldstaerke: c Bestimme naechstgelegene Stuetzstellen und die Komponenten des Abstands- c Gittervektors zur allernaechsten Stuetzstelle sowie deren Betraege: stuetzstelle_i(1) = nint(real_i) Abstand_i = real_i - stuetzstelle_i(1) ! Abstand zur naeheren Stuetzstelle Abstand_i_Betrag = abs(Abstand_i) if (Abstand_i.gt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) + 1 elseif (Abstand_i.lt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) - 1 else stuetzstelle_i(2) = stuetzstelle_i(1) endif stuetzstelle_j(1) = nint(real_j) Abstand_j = real_j - stuetzstelle_j(1) Abstand_j_Betrag = abs(Abstand_j) if (Abstand_j.gt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) + 1 elseif (Abstand_j.lt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) - 1 else stuetzstelle_j(2) = stuetzstelle_j(1) endif if (real_i.GT.real(imax)) then c if (real_i.GT.real(imax)+.5) then write(*,*)'Sub_integr_SP_1: Fehlermarke 0' STOP c else c real_i = imax c endif endif if (stuetzstelle_j(1).EQ.jmax .AND. Abstand_j.GT.0.) then write(*,*)'Sub_integr_SP_1: Fehlermarke 1' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP elseif (stuetzstelle_j(1).EQ.0 .AND. Abstand_j.LT.0.) then write(*,*)'Sub_integr_SP_1: Fehlermarke 2' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP endif c............................................................................... c Berechnen des elektrischen Feldes: c ---------------------------------- c............................................................................... c Berechne die x-Komponente der Feldstaerke: c Um die Feldstaerke zu bekommen, interpoliere jeweils linear zwischen den c Werten auf den beiden naechstgelegenen j-Ketten: i = stuetzstelle_i(1) do n = 1, 2 j = stuetzstelle_j(n) ihelp = j*(imax+1) + i if (i.EQ.imax) then E_(n) = map(ihelp-1) - map(ihelp) elseif (i.GT.0) then E_(n) = (-0.5+Abstand_i)*(map(ihelp)-map(ihelp-1)) + + ( 0.5+Abstand_i)*(map(ihelp)-map(ihelp+1)) else E_(n) = map(ihelp) - map(ihelp+1) endif enddo E(1) = E_(1) + Abstand_j_Betrag*(E_(2)-E_(1)) E(1) = E(1) / Dx ! Reskalierung entsprechend x-Gitterkonstanten c Berechne die y-Komponente der Feldstaerke: j = stuetzstelle_j(1) do n = 1, 2 i = stuetzstelle_i(n) ihelp = j*(imax+1) + i if (j.EQ.jmax) then E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp-(imax+1))) elseif (j.GT.0) then E_(n) = (-0.5+Abstand_j)*(map(ihelp)-map(ihelp-(imax+1))) + + ( 0.5+Abstand_j)*(map(ihelp)-map(ihelp+(imax+1))) else ! j=0 -> map(i,j-1) = map(i,j+1) == map(i,1) E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp+(imax+1))) endif enddo E(2) = E_(1) + Abstand_i_Betrag*(E_(2)-E_(1)) E(2) = E(2) / Dy ! Reskalierung entsprechend y-Gitterkonstanten if (real_j_.LT.0) E(2) = -E(2) END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE EFeld_SP_2(x,E) c ========================== IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_2.INC' real real_i,real_j,real_j_ ! x,r im Mappensystem in Gittereinheiten integer stuetzstelle_i(2) ! naechste Stuetzstellen in x- und integer stuetzstelle_j(2) ! r-Richtung real Abstand_i,Abstand_i_Betrag ! Entfernung zur naechsten Stuetzstelle real Abstand_j,Abstand_j_Betrag ! (in Gittereinheiten!) integer i,j, n, ihelp real x(3),E(3) ! Ort und Feldstaerke real E_(2) ! Hilfsspeicher fuer Feldberechnung integer zaehler_Sp_2 /0/ common /zaehler_Sp_2/ zaehler_Sp_2 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - zaehler_Sp_2 = zaehler_Sp_2 + 1 c Rechne in Gittereinheiten um und teste, ob Raumpunkt innerhalb der Potential- c mappe liegt: real_i = (x(1)-xmin)/Dx ! Umrechnung in Mappensystem findet bereits in real_j_ = x(2)/ Dy ! MUTRACK statt. c Mappe ist zyklisch in y-Richtung. j laeuft von 0 bis jmax vom 0. Gitterstab c zur Mitte zwischen 0. und 1. Gitterstab. c -> 1.) waehle aequivalentes real_j aus Bereich von Mitte zwischen -1. und 0. c Gitterstab bis Mitte zwischen 0. und 1. Gitterstab. c 2.) Nimm Absolutbetrag. n = nInt(real_j_/(2.*real(jmax))) real_j_ = real_j_ - real(n) * (2.*real(jmax)) real_j = abs(real_j_) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Berechne die Feldstaerke: c Bestimme naechstgelegene Stuetzstellen und die Komponenten des Abstands- c Gittervektors zur allernaechsten Stuetzstelle sowie deren Betraege: stuetzstelle_i(1) = nint(real_i) Abstand_i = real_i - stuetzstelle_i(1) ! Abstand zur naeheren Stuetzstelle Abstand_i_Betrag = abs(Abstand_i) if (Abstand_i.gt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) + 1 elseif (Abstand_i.lt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) - 1 else stuetzstelle_i(2) = stuetzstelle_i(1) endif stuetzstelle_j(1) = nint(real_j) Abstand_j = real_j - stuetzstelle_j(1) Abstand_j_Betrag = abs(Abstand_j) if (Abstand_j.gt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) + 1 elseif (Abstand_j.lt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) - 1 else stuetzstelle_j(2) = stuetzstelle_j(1) endif if (real_i.GT.real(imax)) then c if (real_i.GT.real(imax)+.5) then write(*,*)'Sub_integr_SP_2: Fehlermarke 0' STOP c else c real_i = imax c endif endif if (stuetzstelle_j(1).EQ.jmax .AND. Abstand_j.GT.0.) then write(*,*)'Sub_integr_SP_2: Fehlermarke 1' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP elseif (stuetzstelle_j(1).EQ.0 .AND. Abstand_j.LT.0.) then write(*,*)'Sub_integr_SP_2: Fehlermarke 2' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP endif c............................................................................... c Berechnen des elektrischen Feldes: c ---------------------------------- c............................................................................... c Berechne die x-Komponente der Feldstaerke: c Um die Feldstaerke zu bekommen, interpoliere jeweils linear zwischen den c Werten auf den beiden naechstgelegenen j-Ketten: i = stuetzstelle_i(1) do n = 1, 2 j = stuetzstelle_j(n) ihelp = j*(imax+1) + i if (i.EQ.imax) then E_(n) = map(ihelp-1) - map(ihelp) elseif (i.GT.0) then E_(n) = (-0.5+Abstand_i)*(map(ihelp)-map(ihelp-1)) + + ( 0.5+Abstand_i)*(map(ihelp)-map(ihelp+1)) else E_(n) = map(ihelp) - map(ihelp+1) endif enddo E(1) = E_(1) + Abstand_j_Betrag*(E_(2)-E_(1)) E(1) = E(1) / Dx ! Reskalierung entsprechend x-Gitterkonstanten c Berechne die y-Komponente der Feldstaerke: j = stuetzstelle_j(1) do n = 1, 2 i = stuetzstelle_i(n) ihelp = j*(imax+1) + i if (j.EQ.jmax) then E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp-(imax+1))) elseif (j.GT.0) then E_(n) = (-0.5+Abstand_j)*(map(ihelp)-map(ihelp-(imax+1))) + + ( 0.5+Abstand_j)*(map(ihelp)-map(ihelp+(imax+1))) else ! j=0 -> map(i,j-1) = map(i,j+1) == map(i,1) E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp+(imax+1))) endif enddo E(2) = E_(1) + Abstand_i_Betrag*(E_(2)-E_(1)) E(2) = E(2) / Dy ! Reskalierung entsprechend y-Gitterkonstanten if (real_j_.LT.0) E(2) = -E(2) END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE EFeld_SP_3(x,E) c ========================== IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:MAP_DEF_SP_3.INC' real real_i,real_j,real_j_ ! x,r im Mappensystem in Gittereinheiten integer stuetzstelle_i(2) ! naechste Stuetzstellen in x- und integer stuetzstelle_j(2) ! r-Richtung real Abstand_i,Abstand_i_Betrag ! Entfernung zur naechsten Stuetzstelle real Abstand_j,Abstand_j_Betrag ! (in Gittereinheiten!) integer i,j, n, ihelp real x(3),E(3) ! Ort und Feldstaerke real E_(2) ! Hilfsspeicher fuer Feldberechnung integer zaehler_Sp_3 /0/ common /zaehler_Sp_3/ zaehler_Sp_3 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - zaehler_sp_3 = zaehler_sp_3 + 1 c Rechne in Gittereinheiten um und teste, ob Raumpunkt innerhalb der Potential- c mappe liegt: real_i = (x(1)-xmin)/Dx ! Umrechnung in Mappensystem findet bereits in real_j_ = x(2)/ Dy ! MUTRACK statt. c Mappe ist zyklisch in y-Richtung. j laeuft von 0 bis jmax vom 0. Gitterstab c zur Mitte zwischen 0. und 1. Gitterstab. c -> 1.) waehle aequivalentes real_j aus Bereich von Mitte zwischen -1. und 0. c Gitterstab bis Mitte zwischen 0. und 1. Gitterstab. c 2.) Nimm Absolutbetrag. n = nInt(real_j_/(2.*real(jmax))) real_j_ = real_j_ - real(n) * (2.*real(jmax)) real_j = abs(real_j_) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Berechne die Feldstaerke: c Bestimme naechstgelegene Stuetzstellen und die Komponenten des Abstands- c Gittervektors zur allernaechsten Stuetzstelle sowie deren Betraege: stuetzstelle_i(1) = nint(real_i) Abstand_i = real_i - stuetzstelle_i(1) ! Abstand zur naeheren Stuetzstelle Abstand_i_Betrag = abs(Abstand_i) if (Abstand_i.gt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) + 1 elseif (Abstand_i.lt.0.) then stuetzstelle_i(2) = stuetzstelle_i(1) - 1 else stuetzstelle_i(2) = stuetzstelle_i(1) endif stuetzstelle_j(1) = nint(real_j) Abstand_j = real_j - stuetzstelle_j(1) Abstand_j_Betrag = abs(Abstand_j) if (Abstand_j.gt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) + 1 elseif (Abstand_j.lt.0.) then stuetzstelle_j(2) = stuetzstelle_j(1) - 1 else stuetzstelle_j(2) = stuetzstelle_j(1) endif if (real_i.GT.real(imax)) then c if (real_i.GT.real(imax)+.5) then write(*,*)'Sub_integr_SP_3: Fehlermarke 0' STOP c else c real_i = imax c endif endif if (stuetzstelle_j(1).EQ.jmax .AND. Abstand_j.GT.0.) then write(*,*)'Sub_integr_SP_3: Fehlermarke 1' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP elseif (stuetzstelle_j(1).EQ.0 .AND. Abstand_j.LT.0.) then write(*,*)'Sub_integr_SP_3: Fehlermarke 2' write(*,*)'xmin,xmax,ymax = ',xmin,xmax,ymax write(*,*)'x = ',x write(*,*)'real_j, stuetzstelle_j,abstand_j = ',real_j, stuetzstelle_j,abstand_j STOP endif c............................................................................... c Berechnen des elektrischen Feldes: c ---------------------------------- c............................................................................... c Berechne die x-Komponente der Feldstaerke: c Um die Feldstaerke zu bekommen, interpoliere jeweils linear zwischen den c Werten auf den beiden naechstgelegenen j-Ketten: i = stuetzstelle_i(1) do n = 1, 2 j = stuetzstelle_j(n) ihelp = j*(imax+1) + i if (i.EQ.imax) then E_(n) = map(ihelp-1) - map(ihelp) elseif (i.GT.0) then E_(n) = (-0.5+Abstand_i)*(map(ihelp)-map(ihelp-1)) + + ( 0.5+Abstand_i)*(map(ihelp)-map(ihelp+1)) else E_(n) = map(ihelp) - map(ihelp+1) endif enddo E(1) = E_(1) + Abstand_j_Betrag*(E_(2)-E_(1)) E(1) = E(1) / Dx ! Reskalierung entsprechend x-Gitterkonstanten c Berechne die y-Komponente der Feldstaerke: j = stuetzstelle_j(1) do n = 1, 2 i = stuetzstelle_i(n) ihelp = j*(imax+1) + i if (j.EQ.jmax) then E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp-(imax+1))) elseif (j.GT.0) then E_(n) = (-0.5+Abstand_j)*(map(ihelp)-map(ihelp-(imax+1))) + + ( 0.5+Abstand_j)*(map(ihelp)-map(ihelp+(imax+1))) else ! j=0 -> map(i,j-1) = map(i,j+1) == map(i,1) E_(n) = 2.0*Abstand_j*(map(ihelp)-map(ihelp+(imax+1))) endif enddo E(2) = E_(1) + Abstand_i_Betrag*(E_(2)-E_(1)) E(2) = E(2) / Dy ! Reskalierung entsprechend y-Gitterkonstanten if (real_j_.LT.0) E(2) = -E(2) END c===============================================================================