630 lines
18 KiB
Fortran
630 lines
18 KiB
Fortran
c===============================================================================
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE READ_INFOS
|
|
c =====================
|
|
|
|
INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC'
|
|
|
|
logical map_error /.false./
|
|
COMMON /map_error/ map_error
|
|
|
|
c if (.NOT.Use_MUTRACK .OR. write_geo) then
|
|
CALL READ_INFO_L1
|
|
if (.NOT.idealMirror) then
|
|
call read_INFO_SP_1
|
|
call read_INFO_SP_2
|
|
call read_INFO_SP_3
|
|
endif
|
|
if (lense2) CALL READ_INFO_L2andFo
|
|
if (TriggerInBeam) CALL READ_INFO_FO
|
|
c endif
|
|
if (.NOT.upTOTDFoilOnly) then
|
|
CALL READ_INFO_L3
|
|
CALL READ_INFO_M2
|
|
endif
|
|
|
|
if (map_error) then
|
|
write(*,*)'----------------------------------------'//
|
|
+ '----------------------------------------'
|
|
call exit
|
|
endif
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE READ_INFO_L1
|
|
c =======================
|
|
|
|
IMPLICIT NONE
|
|
|
|
character*(*) Nr
|
|
parameter (Nr='L1')
|
|
|
|
INCLUDE 'mutrack$sourcedirectory:COM_LUNS.INC'
|
|
INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC'
|
|
INCLUDE 'mutrack$sourcedirectory:MAP_DEF_L1.INC'
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
namelist /geometry/
|
|
+ DistanceCyl,iRadiusCyl,
|
|
+ LengthOuterCyl,oRadiusOuterCyl,
|
|
+ LengthInnerCyl,oRadiusInnerCyl,
|
|
+ RadiusVacTube
|
|
|
|
namelist /grid_info/
|
|
+ Dx,Dr, imax,jmax, xmax,rmax
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
logical map_error
|
|
COMMON /map_error/ map_error
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Einlesen der Mappen-Informationen:
|
|
|
|
open (lunREAD,file=mappenName//'.INFO',defaultfile=mappenDir,
|
|
+ readonly,status='old')
|
|
read (lunREAD,nml=grid_info)
|
|
rewind (lunREAD)
|
|
read (lunREAD,nml=geometry)
|
|
close (lunREAD)
|
|
|
|
|
|
c Pruefen, ob die eingelesene Mappe den Anforderungen genuegt:
|
|
|
|
if (iRadiusCyl.NE.rmax) then
|
|
write(*,*)' L1-Mappe: rmax der Mappe stimmt nicht mit Innenradius der Zylinder ueberein'
|
|
write(*,*)' -> STOP'
|
|
STOP
|
|
endif
|
|
|
|
|
|
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:
|
|
|
|
xEnterMap = xCenterOfLense - xmax
|
|
xLeaveMap = xCenterOfLense + xmax
|
|
|
|
|
|
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_MAP_L1
|
|
c ======================
|
|
|
|
IMPLICIT NONE
|
|
|
|
character*(*) Nr
|
|
parameter (Nr='L1')
|
|
|
|
INCLUDE 'mutrack$sourcedirectory:MAP_DEF_L1.INC'
|
|
INCLUDE 'mutrack$sourcedirectory:READ_MAP.INC'
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE INTEGRATIONSSTEP_RUNGE_KUTTA_L1(dt)
|
|
c ==============================================
|
|
|
|
IMPLICIT NONE
|
|
SAVE
|
|
|
|
character*(*) Nr
|
|
parameter (Nr='L1')
|
|
|
|
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'
|
|
INCLUDE 'mutrack$sourcedirectory:MAP_DEF_L1.INC'
|
|
|
|
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_L1(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_L1(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_L1(x1,EFeld1,*999)
|
|
|
|
|
|
c mache zweiten dt/2 - Schritt:
|
|
|
|
call SINGLESTEP_RUNGE_KUTTA_L1(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_L1(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, neuen so gross wie sinnvoller weise moeglich
|
|
c machen:
|
|
|
|
t = t + dt
|
|
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.2) then
|
|
destiny = code_wand
|
|
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_L1(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_L1/ 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_L1(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_L1(xTest,E2,*999)
|
|
|
|
do i = 1, 3
|
|
xTest(i) = x0(i) + v2(i) * dt
|
|
v3(i) = v0(i) + E2(i) * help
|
|
enddo
|
|
call EFeld_L1(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_L1(x,E,*)
|
|
c ==========================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'mutrack$sourcedirectory:MAP_DEF_L1.INC'
|
|
|
|
real real_i,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 radius ! Betrag des Radiusvektors in y-z-Ebene
|
|
real x(3),E(3) ! Ort und Feldstaerke
|
|
real E_(2) ! Hilfsspeicher fuer Feldberechnung
|
|
real Erad ! radiale 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 Rechne in Gittereinheiten um:
|
|
|
|
radius = sqrt(x(2)*x(2)+x(3)*x(3))
|
|
|
|
real_i = (xmax - abs(x(1)-xCenterOfLense)) / Dx
|
|
real_j = radius / Dr
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Mache die Tests und berechne die Feldstaerke:
|
|
|
|
c Teste, ob Raumpunkt innerhalb der Potentialmappe liegt:
|
|
|
|
if (real_j.GT.jmax) then
|
|
returnCode_EFeld = 2
|
|
RETURN 1
|
|
elseif (real_i.LT.0.) then
|
|
E(1) = 0.
|
|
E(2) = 0.
|
|
E(3) = 0.
|
|
RETURN
|
|
c else
|
|
c returnCode_EFeld = 0
|
|
endif
|
|
|
|
|
|
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
|
|
|
|
|
|
c...............................................................................
|
|
c Berechnen des elektrischen Feldes:
|
|
c ----------------------------------
|
|
c
|
|
c Potentialverlauf ist symmetrisch zu real_i = imax (Mittelebene des inneren
|
|
c Zylinders) und zu j=0 (Linsenachse):
|
|
c
|
|
c map(imax+i,-j) == map(imax-1,j),
|
|
c map(i,-j) == map(i,j).
|
|
c
|
|
c Entlang j=0 ist also Erad=0 und damit E(2)=0 und E(3)=0, bei i=imax ist E(1)=0
|
|
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) = 2. * Abstand_i *(map(ihelp)-map(ihelp-1)) ! i=imax -> map(imax+1,j) = map(imax-1,j)
|
|
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
|
|
if (x(1).GT.xCenterOfLense) E(1) = - E(1)
|
|
|
|
|
|
c Berechne die radiale Komponente der Feldstaerke:
|
|
|
|
if (real_j.LT.1e-10) then
|
|
E(2) = 0.
|
|
E(3) = 0.
|
|
RETURN
|
|
endif
|
|
|
|
j = stuetzstelle_j(1)
|
|
|
|
do n = 1, 2
|
|
i = stuetzstelle_i(n)
|
|
ihelp = j*(imax+1) + i
|
|
if (j.EQ.jmax) then
|
|
E_(n) = map(ihelp-(imax+1)) - map(ihelp)
|
|
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
|
|
Erad = E_(1) + Abstand_i_Betrag*(E_(2)-E_(1))
|
|
Erad = Erad / Dr ! Reskalierung entsprechend r-Gitterkonstanten
|
|
|
|
|
|
c Berechne E(2) und E(3) aus Erad:
|
|
|
|
E(2) = Erad * x(2) / radius
|
|
E(3) = Erad * x(3) / radius
|
|
|
|
|
|
cd write(18,*)'x,E = ',x,E
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|