2378 lines
67 KiB
Fortran
2378 lines
67 KiB
Fortran
c*******************************************************************************
|
|
c - bei MUTRACK kann Graphik auch im Batchjob erstellt werden!
|
|
c
|
|
c - NTP_DEBUG einrichten.
|
|
c
|
|
c - bei StartFlaeche.EQ.1 immer E0==0 setzen (gegebenenfalls mit Kommentarausgabe)
|
|
c
|
|
c - S1G1, S1G2, G1G2 einfuehren. (bei der 1. Ueberquerung der Gitter mit t<>0.
|
|
c zu fuellen.
|
|
c*******************************************************************************
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
c===============================================================================
|
|
program ACCEL
|
|
c===============================================================================
|
|
|
|
c Deklarationen:
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
INCLUDE 'accel$sourcedirectory:COM_DIRS.INC'
|
|
INCLUDE 'accel$sourcedirectory:COM_GEO.INC'
|
|
INCLUDE 'accel$sourcedirectory:COM_HVs.INC'
|
|
INCLUDE 'accel$sourcedirectory:COM_Bs.INC'
|
|
|
|
|
|
c die SCHLEIFENVARIABLEN und damit zusammenhaengendes (in der Reihenfolge, wie
|
|
c sie in den 'do 200 ...'-Schleifen auftreten).
|
|
c (Common-Bloecke werden fuer die NTupel-Ausgabe benoetigt)
|
|
|
|
c - Masse und Ladung:
|
|
|
|
real m, m_ ! Masse, Laufvariable fuer Massen-Schleife
|
|
real q, q_ ! Ladung, Laufvariable fuer Ladungs-Schleife
|
|
|
|
|
|
c - Startparameter:
|
|
|
|
c integer randomloop_ ! Laufvariable fuer zufallsverteilte Starts
|
|
real E0_ ! Laufvariable fuer Startenergie_Schleife
|
|
real theta0_ ! Laufvarialbe fuer Startwinkel-Schleife
|
|
real Sin_theta0, Cos_theta0 ! Startwinkel gegen x-Achse
|
|
real phi0_ ! Laufvariable fuer Startwinkel-Schleife
|
|
real Sin_phi0, Cos_phi0 ! azimuthaler Startwinkel (phi0=0: y-Achse)
|
|
real y0_ ! Laufvariable fuer Startpositions_Schleife
|
|
real z0_ ! Laufvariable fuer Startpositions_Schleife
|
|
real r0 ! Radius beim Wuerfeln der Startposition
|
|
real phi_r0 ! Winkel beim Wuerfeln der Startposition
|
|
|
|
! x0(3),v0(3),E0,theta0,phi0 werden in 'COM_ACCEL.INC' declariert
|
|
|
|
|
|
c allgemeine Trajektoriengroessen
|
|
|
|
real dt ! zeitl. Aenderung
|
|
real v_xy ! Geschwindigkeit in x/y-Ebene
|
|
real v0_Betrag
|
|
real Ekin ! kinetische Energie
|
|
COMMON /Ekin/ Ekin
|
|
real radiusQuad ! RadiusQuadrat
|
|
real rQuadHeShield ! RadiusQuadrat (innen)
|
|
real x40(2:3) ! auf x = 40 mm extrapolierter Ort (FUER NTP_40mm)
|
|
real t40 ! auf x = 40 mm extrapolierte Zeit (FUER NTP_40mm)
|
|
common/xt40/ x40,t40 ! fuer NTupel
|
|
|
|
c Variablen fuer den allgemeinen Programmablauf:
|
|
|
|
character uhrzeit*8
|
|
|
|
integer percent_done
|
|
integer zaehler ! Zaehler fuer Monitoring der Trajektorie mittels
|
|
! GRAPHICS und DEBUG
|
|
logical flag_ok ! logische Hilfsvariable
|
|
integer okStepsCounter ! fuer die Berechnung der durchscnittlichen Anzahl
|
|
! Rechenschritte bis zum Ziel
|
|
|
|
logical alreadyTested, wireHit ! Fuer Ueberpruefung von Gitter-Treffern
|
|
real distToWire(2)
|
|
|
|
integer i ! integer-Hilfsvariable
|
|
real help1,help2,help3,help4 ! real-Hilfsvariablen
|
|
|
|
character helpChar*7, ant*1
|
|
character HistogramTitle*19 /'Startort-Verteilung'/
|
|
|
|
real yAbs,zAbs ! = abs(x(2)),abs(x(3))
|
|
real yWire ! y-Position von Gitterstaeben
|
|
|
|
real xMarke,yMarke,zMarke ! Grenzen des Bereiches mit entschaerfter
|
|
! Fehlerbetrachtung
|
|
logical reducedAccur
|
|
COMMON /reducedAccur/ reducedAccur
|
|
|
|
integer lunZwi(0:1) /lunZwi1,lunZwi2/, lunIndx /1/
|
|
integer mappe0,mappe
|
|
|
|
logical calcMap(6) ! gibt an, welche Mappen noch zu rechnen sind
|
|
DATA calcMap / .false.,.false.,.false.,.false.,.false.,.false. /
|
|
|
|
|
|
c Variablen fuer Test, ob Targethalter getroffen wurde:
|
|
|
|
real TgtFactorY, TgtFactorZ
|
|
real TgtConstY , TgtConstZ
|
|
|
|
|
|
c Variablen fuer die Graphikausgabe:
|
|
|
|
real xKoord(1000) ! Koordinatenfelder fuer die
|
|
real yKoord(1000) ! Graphikausgabe
|
|
real zKoord(1000) !
|
|
integer nKoord ! Anzahl der Koordinaten
|
|
|
|
COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord ! fuer Graphikaufruf
|
|
|
|
|
|
c Variablen fuer HBOOK und PAW:
|
|
|
|
integer istat ! fuer HBOOK-Fehlermeldungen
|
|
|
|
integer HB_memsize
|
|
parameter(HB_memsize=1000000)
|
|
real memory(HB_memsize)
|
|
|
|
common /pawc/ memory ! Der Arbeitsbereich fuer HBOOK
|
|
|
|
|
|
c Konstanten:
|
|
|
|
real c ! Lichtgeschwindigkeit in mm/ns
|
|
real meanLifeTime ! mittlere Myon-Lebensdauer in ns
|
|
|
|
parameter (c = 299.7925, meanLifeTime = 2197)
|
|
|
|
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
c HIER GEHT DER PROGRAMMTEXT RICHTIG LOS
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
|
|
c Initialisierungen:
|
|
|
|
INCLUDE 'accel$sourcedirectory:INITIALIZE.INC'
|
|
|
|
|
|
c Einlesen der Parameter aus 'ACCEL.INPUT' und Setzen der entsprechenden
|
|
c Voreinstellungen (liest auch die INFO-files der Potentialmappen):
|
|
|
|
call read_inputFile
|
|
rQuadHeShield = rHeShield*rHeShield
|
|
|
|
|
|
c Grenzen des Bereiches mit entschaerfter Fehlerbetrachtung definieren:
|
|
|
|
xMarke = xEnd_TgtHolder + 1
|
|
yMarke = Dy_Foil - 1
|
|
zMarke = Dz_Foil - 1
|
|
|
|
|
|
c falls Bildschirmausgabe verlangt ist, etwas Abstand schaffen:
|
|
|
|
write(*,*)
|
|
if (n_outWhere.GE.2) then
|
|
do i = 1, 5
|
|
write(*,*)
|
|
enddo
|
|
endif
|
|
|
|
|
|
c Anfangs-Informationen ausgeben und abfragen, ob Einstellungen ok sind. Alle
|
|
c Ausgabefiles initialisieren:
|
|
|
|
call initialize_output
|
|
|
|
|
|
c CERN-Pakete initialisieren (Groesse des COMMONblocks /PAWC/ uebermitteln):
|
|
|
|
call HLIMIT(HB_memsize)
|
|
|
|
|
|
c Graphikausgabe initialisieren:
|
|
|
|
if (GRAPHICS) then
|
|
CALL MZEBRA(-3)
|
|
CALL MZPAW (HB_memsize,' ')
|
|
call masstab_setzen
|
|
CALL HPLSET ('VSIZ',.6) ! AXIS VALUES SIZE
|
|
CALL HPLSET ('TSIZ',.7) ! HISTOGRAM TITLE SIZE
|
|
if (random_pos) then
|
|
if (random_y0z0_equal.OR.random_y0z0_Gauss) then
|
|
CALL HBOOK2 (50,HistogramTitle,
|
|
+ 100, -StartBreite/2., StartBreite/2.,
|
|
+ 100, -StartHoehe /2., StartHoehe /2., 20.)
|
|
else
|
|
CALL HBOOK2 (50,HistogramTitle,
|
|
+ 100, -StartRadius, StartRadius,
|
|
+ 100, -StartRadius, StartRadius, 20.)
|
|
endif
|
|
else
|
|
CALL HBOOK2 (50,HistogramTitle,
|
|
+ 100, max(par(1,yPos),5.), max(par(2,yPos),5.),
|
|
+ 100, max(par(1,zPos),5.), max(par(2,zPos),5.), 20.)
|
|
endif
|
|
|
|
c Textfond waehlen! (-> keine Proportionalschrift!)
|
|
endif
|
|
|
|
|
|
c NTP-relevante Befehle:
|
|
|
|
!c Datei fuer NTupelausgabe oeffnen:
|
|
call HROPEN(lunNTP,'ACCEL',outDir//':'//filename//'.NTP',
|
|
+ 'N',1024,istat)
|
|
if (istat.NE.0) then
|
|
write(*,*)
|
|
write(*,*)'error ',istat,' opening HBOOK-file'
|
|
write(*,*)
|
|
STOP
|
|
endif
|
|
call HBNT(idNTP,filename,'D') ! Disk resident CWN buchen
|
|
|
|
!c die Bloecke des CWN definieren:
|
|
|
|
if (NTP_MISC) then
|
|
call HBNAME(idNTP,'LOOP',schleifenNr,'loop[1,1000]:u,'//
|
|
+ 'startNr[1,100000]:u,Mappe[1,6]:u,steps[0,100000]:u')
|
|
endif
|
|
|
|
if (NTP_Start.OR.random_pos) call HBNAME(idNTP,'x0',x0,
|
|
+ 'x0,y0,z0')
|
|
c + 'x0:r*4:20:[-10,30],y0:r*4:20:[-40,40],z0:r*4:20:[-40,40]')
|
|
if (NTP_Start) call HBNAME(idNTP,'v0',v0,'vx0,vy0,vz0')
|
|
if (NTP_Start.OR.random_E0) call HBNAME(idNTP,'E0',E0,'E0')
|
|
if (NTP_Start.OR.random_angle) call HBNAME(idNTP,'angle0',theta0,
|
|
+ 'theta0,phi0')
|
|
c + 'theta0:r*4:20:[0,180],phi0:r*4:20:[0,360]') ! theta0,phi0
|
|
|
|
if (useDecay) call HBNAME(idNTP,'lifetime',lifetime,'lifetime:r')
|
|
|
|
call HBNAME(idNTP,'dest',gebiet,'Gebiet[0,4]:u,dest[-10,10]:i')
|
|
|
|
call HBNAME(idNTP,'Traj',t,'t,x,y,z,vx,vy,vz')
|
|
if (NTP_Stop) call HBNAME(idNTP,'Ekin',Ekin,'Ekin')
|
|
if (NTP_40mm) call HBNAME(idNTP,'x=40mm',x40,
|
|
+ 'y40,z40,t40')
|
|
|
|
|
|
c - - - - - - - - - -
|
|
c folgendes noch einrichten: (durch verlegen aus 'LOOP')
|
|
c
|
|
c if (NTP_debug) call HBNAME(idNTP,'DEBUG',startNr[1,100000]:u,'//
|
|
c + 'Gebiet[0,Gebiete_Anzahl]:u,Mappe[1,6]:u,'//
|
|
c + 'steps[0,MaxSteps+100]:u')
|
|
c
|
|
c if (NTP_Koord) then ! alle Koordinaten in NTP aufnehmen
|
|
c Call HBNAME(idNTP,'koord',Koord_NTP,
|
|
c + 'tTgt ,xTgt ,yTgt ,zTgt ,vxTgt ,vyTgt ,vzTgt ,ETgt ,'//
|
|
c + 'tGr1 ,xGr1 ,yGr1 ,zGr1 ,vxGr1 ,vyGr1 ,vzGr1 ,EGr1 ,'//
|
|
c + 'tGr1 ,xGr1 ,yGr1 ,zGr1 ,vxGr1 ,vyGr1 ,vzGr1 ,EGr1 '//
|
|
c endif
|
|
c - - - - - - - - - -
|
|
|
|
|
|
c die Einsprungposition fuer den Beginn der Trajektorienberechnungen setzen:
|
|
|
|
if (x0(1).LT.xEndMap1) then
|
|
mappe0 = 1
|
|
elseif (x0(1).LT.xEndMap2) then
|
|
mappe0 = 2
|
|
elseif (x0(1).LT.xEndMap3) then
|
|
mappe0 = 3
|
|
elseif (x0(1).LT.xEndMap4) then
|
|
mappe0 = 4
|
|
elseif (x0(1).LT.xEndMap5) then
|
|
mappe0 = 5
|
|
else
|
|
mappe0 = 6
|
|
endif
|
|
|
|
|
|
c Erstellen des .INFO-files:
|
|
|
|
call make_INFOfile
|
|
|
|
|
|
c Erstellen der Hilfsgroesen fuer den Test, ob Targethalter getoffen wurde:
|
|
|
|
TgtFactorY = (innerDy2_TgtHolder-innerDy1_TgtHolder)/(xEnd_TgtHolder-xFoil)
|
|
TgtConstY = innerDy1_TgtHolder - xFoil*TgtFactorY
|
|
|
|
TgtFactorZ = (innerDz2_TgtHolder-innerDz1_TgtHolder)/(xEnd_TgtHolder-xFoil)
|
|
TgtConstZ = innerDz1_TgtHolder - xFoil*TgtFactorZ
|
|
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c ab hier beginnen die Schleifen. (N.B.: die Laufvariable darf kein Feldelement
|
|
c sein)
|
|
c
|
|
c Besonderheit der Massen- und der Ladungsschleife:
|
|
c Wurde im INPUT-File in der Variablen 'artList' eine Teilchenart spezifi-
|
|
c ziert (-> 'artList_defined'), so werden die Parameter Masse und Ladung nicht
|
|
c entsprechend den Inhalten von par(n,mass) bzw. par(n,charge) eingestellt,
|
|
c sondern entsprechend den zu den Teilchenarten gehoerenden Werten fuer diese
|
|
c Groessen. In diesem Fall besteht die Massenschleife aus genau einem (Leer-)
|
|
c Durchlauf, waehrend die Ladungsschleife fuer jede Teilchenart einen Durchlauf
|
|
c macht, in welcher dann die Einstellung von Ladung UND Masse stattfindet.
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
|
|
c Masse und Ladung:
|
|
c -----------------
|
|
|
|
do 200 m_ = par(1,mass),par(2,mass),par(3,mass) ! MASSE
|
|
if (.NOT.artList_defined) then
|
|
m = m_
|
|
parWert(mass) = m
|
|
endif
|
|
|
|
do 200 q_ = par(1,charge),par(2,charge),par(3,charge) ! LADUNG
|
|
if (.NOT.artList_defined) then
|
|
q = q_
|
|
parWert(charge) = q
|
|
else
|
|
ArtNr = Art_Nr(q_)
|
|
m = Art_Masse(ArtNr)
|
|
q = Art_Ladung(ArtNr)
|
|
parWert(mass) = m
|
|
parWert(charge) = q
|
|
endif
|
|
Energie_Faktor = m / (2.*c*c)
|
|
Beschl_Faktor = q / m * c*c
|
|
|
|
if (useDecay) then ! 'useDecay' setzt 'artList_defined' voraus!
|
|
if (ArtNr.LE.4) then ! es ist ein Myon involviert
|
|
useDecay_ = .true.
|
|
else ! kein Myon involviert
|
|
useDecay_ = .false.
|
|
endif
|
|
endif
|
|
|
|
|
|
c Spannungen:
|
|
c -----------
|
|
|
|
do 200 UTgt = par(1,UTarget),par(2,UTarget),par(3,UTarget) ! U(TARGET)
|
|
parWert(UTarget) = UTgt
|
|
|
|
do 200 UGua = par(1,UGuard),par(2,UGuard),par(3,UGuard) ! U(GUARD)
|
|
parWert(UGuard) = UGua
|
|
|
|
do 200 UG1 = par(1,UGi1),par(2,UGi1),par(3,UGi1) ! U(GITTER)
|
|
parWert(UGi1) = UG1
|
|
|
|
c Magnetfelder:
|
|
c -------------
|
|
|
|
do 200 B_Helm = par(1,BHelm),par(2,BHelm),par(3,BHelm) ! Helmholtzs
|
|
parWert(BHelm) = B_Helm
|
|
|
|
do 200 B_TD = par(1,BTD),par(2,BTD),par(3,BTD) ! TD-Spule
|
|
parWert(BTD) = B_TD
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c Die in den ab hier beginnenden Startparameter-Schleifen eingestellten Start-
|
|
c werte werden bei Verwendung entsprechender Zufallsverteilungen als Offset fuer
|
|
c die letztendlichen Startwerte verwendet:
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Startparameter:
|
|
c ---------------
|
|
|
|
do 200 E0_ = par(1,ener),par(2,ener),par(3,ener) ! E0
|
|
c if (.NOT.random_E0) E0 = E0_
|
|
E0 = E0_
|
|
|
|
c das hier folgende wurde eingefuegt um fuer die Auswertung der Startenergie-
|
|
c verteilung variabler in der Aufteilung der Startenergieintervalle zu sein.
|
|
c (Zuvor musste man um variable Binbreiten zu haben mehrere Accel- und
|
|
c Mutrack-Simulationen durchlaufen lassen).
|
|
|
|
if (E0InterFromFile) then
|
|
lowerE0 = E0Low(nInt(E0_))
|
|
upperE0 = E0Low(nint(E0_+1))
|
|
endif
|
|
|
|
c das hier folgende ist eine spezielle Erweiterung fuer die Auswertung der
|
|
c Energieverluste der hoeherenergetischen Myonen. Hier wird aus sigma(tof)
|
|
c (= 1.6 ns) in Abhaengigkeit der Energie ein sigma(E0) berechnet. Dies ist
|
|
c gedacht fuer Flugzeitemessungen M3S1 bei Fluglaengen von 1 m unter
|
|
c Vernachlaessigung der Auswirkungen der Felder in den Linsen und im Spiegel:
|
|
|
|
if (adjustSigmaE0) then
|
|
help1 = (E0_ + UTgt)
|
|
sigmaE0 = 0.0118 * help1**1.5
|
|
write(lunLOG,*) 'UTgt,E0_,sigmaE0 = ',UTgt,E0_,sigmaE0
|
|
endif
|
|
|
|
do 200 theta0_ = par(1,thetAng),par(2,thetAng),par(3,thetAng) ! theta0
|
|
if (.NOT.random_angle) then
|
|
theta0 = theta0_
|
|
Cos_theta0 = cosd(theta0)
|
|
Sin_theta0 = sind(theta0)
|
|
endif
|
|
|
|
do 200 phi0_ = par(1,phiAng),par(2,phiAng),par(3,phiAng) ! phi0
|
|
if (.NOT.random_angle) then
|
|
phi0 = phi0_
|
|
Cos_phi0 = cosd(phi0)
|
|
Sin_phi0 = sind(phi0)
|
|
endif
|
|
|
|
do 200 y0_ = par(1,yPos),par(2,yPos),par(3,yPos) ! y0
|
|
if (.NOT.random_pos) x0(2) = y0_
|
|
|
|
do 200 z0_ = par(1,zPos),par(2,zPos),par(3,zPos) ! z0
|
|
if (.NOT.random_pos) x0(3) = z0_
|
|
|
|
c die folgenden parWert(n) werden u.U. in der 'Zufallsschleife' weiter unten
|
|
c abgeaendert. Hier werden sie in jedem Fall fuer Tabellenausgaben, Debug-
|
|
c angelegenheiten, Erstellen des .INFO-files u.s.w. erst einmal mit den
|
|
c aktuellen Werten der entsprechenden Schleifen gefuellt:
|
|
|
|
parWert(ener) = E0_
|
|
parWert(thetAng) = theta0_
|
|
parWert(phiAng) = phi0_
|
|
parWert(yPos) = y0_
|
|
parWert(zPos) = z0_
|
|
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c Hier folgen die Befehle, die zu Beginn jeder neuen Schleife faellig sind:
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
SchleifenNr = SchleifenNr + 1 ! Schleifen zaehlen
|
|
okStepsCounter = 0 ! Zaehler fuer Ermittlung der durchschnittlichen
|
|
! Anzahl Schritte bis zum Ziel
|
|
|
|
|
|
c Die Statistikspeicher resetten:
|
|
c Falls nur ein Teilchenstart pro Schleife erfolgt, nimm die Statistik ueber
|
|
c alle Schleifen. Dann erfolgt der Reset nur bei der ersten Schleife:
|
|
|
|
flag_ok = (.NOT.(OneStartPerLoop .AND. SchleifenNr.GT.1))
|
|
|
|
if (flag_ok) call reset_statistics
|
|
|
|
|
|
c Die Kammer zeichnen:
|
|
c Wird pro Schleife nur ein Teilchen gestartet ('OneStartPerLoop'; d.h. kein
|
|
c oder genau ein 'Zufallsstart'), so trage alle Trajektorien in die gleiche
|
|
c Graphik ein. Zeichne die Kammer dann also nur bei der ersten Schleife.
|
|
|
|
if (GRAPHICS .AND. flag_ok) then
|
|
CALL IZPICT ('CHAMBER','M') ! ERZEUGEN VON BILDERN IM PAWC-COMMON-BLOCK
|
|
CALL IZPICT ('HISTO','M')
|
|
CALL IZPICT ('TEXT','M')
|
|
call plot_chamber
|
|
call Graphics_Text ! Text fuer Textwindow erstellen
|
|
call text_plot ! Ausgabe des Textes
|
|
endif
|
|
|
|
|
|
c Ausgabe der aktuellen Settings:
|
|
c Auch dies im Falle von 'OneStartPerLoop' nur bei der ersten Schleife:
|
|
|
|
if ((n_outWhere.NE.0 .OR. smallLogFile) .AND. flag_ok) then
|
|
call output_new_loop
|
|
endif
|
|
|
|
|
|
c Oeffnen der temporaeren files:
|
|
|
|
c - das Startparameterfile:
|
|
c Speicherbedarf pro record in longwords:
|
|
c x0,v0,E0,theta0,phi0 =
|
|
c 3 + 3+ 1+ 1 + 1 = 9
|
|
|
|
open(lunStart,file=filename//'_start.tmp',defaultfile=TMPDir,
|
|
+ form='unformatted',recl=9,status='scratch')
|
|
|
|
|
|
c - das File fuer die Lebensdauern:
|
|
|
|
if (UseDecay_) open (lunDecay,file=filename//'_decay.tmp',defaultfile=TMPDir,
|
|
+ form='unformatted',status='scratch',recl=1)
|
|
|
|
|
|
c - die Zwischenspeicherungsfiles:
|
|
c Speicherbedarf pro record in longwords:
|
|
c steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt =
|
|
c .5 +1+3+3+ .5 + .5 + .5 + .5 + .5 + 1 = 11
|
|
|
|
open (lunZwi1,file=filename//'_zwi1.tmp',defaultfile=TMPDir,
|
|
+ form='unformatted',status='scratch',recl=11)
|
|
open (lunZwi2,file=filename//'_zwi2.tmp',defaultfile=TMPDir,
|
|
+ form='unformatted',status='scratch',recl=11)
|
|
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c Wuerfeln der Startparameter fuer die Teilchen dieser Schleife:
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
write (*,*) 'Wuerfeln der Startparameter ...'
|
|
|
|
c DEBUG beruecksichtigen:
|
|
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> Erstellen der Startparameter <<<<<<<'
|
|
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c gegebenenfalls werden hier die in den Startparameterschleifen eingestellten
|
|
c Startwerte mit einem zufallsverteilten Anteil ueberlagert:
|
|
c Sind keine Zufallsverteilungen verlangt macht die 'Zufallsschleife' genau
|
|
c einen Durchlauf. (-> n_par(0) == 1)
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
do 100 start_nr = 1, n_par(0)
|
|
|
|
if (random_E0) then ! random_ENERGIE
|
|
if (random_E0_equal) then ! -> gleichverteilt
|
|
276 if (E0InterFromFile) then
|
|
E0 = lowerE0 + (upperE0 - lowerE0)*ran(seed)
|
|
else
|
|
E0 = E0_ + lowerE0 + (upperE0 - lowerE0)*ran(seed)
|
|
endif
|
|
if (E0.LT.0) goto 276
|
|
elseif (random_E0_gauss) then ! -> gaussverteilt
|
|
277 call Gauss_Verteilung(sigmaE0,help1)
|
|
E0 = E0_ + help1
|
|
if (E0.LT.0) goto 277
|
|
endif
|
|
parWert(ener) = E0
|
|
endif
|
|
|
|
if (random_pos) then ! random_POSITION
|
|
if (random_y0z0_equal) then ! -> rechteckig, gleichverteilt
|
|
x0(2) = StartBreite * (ran(seed)-.5)
|
|
x0(3) = StartHoehe * (ran(seed)-.5)
|
|
elseif (random_y0z0_Gauss) then ! -> rechteckig, Gaussverteilt
|
|
278 r0 = abs(sigmaPosition*sqrt(-2.*log(1.-ran(seed))))
|
|
phi_r0= 360.*ran(seed)
|
|
x0(2) = r0 * cosd(phi_r0)
|
|
if (abs(x0(2)).GT.StartBreite/2.) goto 278
|
|
x0(3) = r0 * sind(phi_r0)
|
|
if (abs(x0(3)).GT.StartHoehe/2.) goto 278
|
|
elseif (random_r0_equal) then ! -> rund, gleichverteilt
|
|
r0 = StartRadius * sqrt(ran(seed))
|
|
phi_r0= 360. * ran(seed)
|
|
x0(2) = r0 * cosd(phi_r0)
|
|
x0(3) = r0 * sind(phi_r0)
|
|
elseif (random_r0_Gauss) then ! -> rund, Gaussverteilt
|
|
279 r0 = abs(sigmaPosition*sqrt(-2.*log(1.-ran(seed))))
|
|
if (r0.GT.StartRadius) goto 279
|
|
phi_r0= 360.*ran(seed)
|
|
x0(2) = r0 * cosd(phi_r0)
|
|
x0(3) = r0 * sind(phi_r0)
|
|
endif
|
|
x0(2) = y0_ + x0(2)
|
|
x0(3) = z0_ + x0(3)
|
|
parWert(yPos) = x0(2)
|
|
parWert(zPos) = x0(3)
|
|
endif
|
|
|
|
if (StartFlaeche.EQ.1) then
|
|
! Ziehe alle Teilchen auf jeweils 'naechstem' Gitterstab zusammen.
|
|
! Wuerfle Startposition auf der Oberflaeche des Gitterstabes in
|
|
! der x-y-Ebene gleichverteilt. Lasse Teilchen 1e-5 Radien ueber
|
|
! der Drahtoberflaeche starten.
|
|
if (.NOT.random_pos) x0(2) = y0_
|
|
if (abs(x0(2)).LE.innerDy_Gridframe1 .AND. abs(x0(3)).LE.innerDz_Gridframe1) then
|
|
ywire = nint((x0(2)-y_Pos_firstWire1)/distance_wires1)*distance_wires1 + y_Pos_firstWire1
|
|
help1 = 360.*ran(seed)
|
|
x0(1) = xPosition_Grid1 + 1.00001*rWires1*sind(help1)
|
|
x0(2) = ywire + 1.00001*rWires1*cosd(help1)
|
|
! (z-Komponente bleibt unveraendert)
|
|
endif
|
|
parWert(yPos) = x0(2)
|
|
parWert(zPos) = x0(3)
|
|
endif
|
|
|
|
if (StartFlaeche.EQ.2) then
|
|
! wie bei StartFlaeche.EQ.1
|
|
if (.NOT.random_pos) x0(2) = y0_
|
|
if (abs(x0(2)).LE.innerDy_Gridframe2 .AND. abs(x0(3)).LE.innerDz_Gridframe2) then
|
|
ywire = nint((x0(2)-y_Pos_firstWire2)/distance_wires2)*distance_wires2 + y_Pos_firstWire2
|
|
help1 = 360.*ran(seed)
|
|
x0(1) = xPosition_Grid2 + 1.00001*rWires2*sind(help1)
|
|
x0(2) = ywire + 1.00001*rWires2*cosd(help1)
|
|
! (z-Komponente bleibt unveraendert)
|
|
endif
|
|
parWert(yPos) = x0(2)
|
|
parWert(zPos) = x0(3)
|
|
endif
|
|
|
|
if (graphics) then
|
|
if (abs(x0(2)).LE.20 .AND. abs(x0(3)).LE.20) then
|
|
CALL HFILL (50,x0(2),x0(3),1.)
|
|
endif
|
|
endif
|
|
|
|
if (random_angle) then ! random_WINKEL
|
|
4462 if (random_lambert) then ! -> Lambert-verteilt
|
|
call lambert_verteilung(StartLambertOrd,
|
|
+ Cos_theta0,Sin_theta0)
|
|
theta0 = acosd(Cos_theta0)
|
|
elseif (random_gauss) then
|
|
call Gauss_Verteilung_theta(sigmaWinkel,theta0)
|
|
Cos_theta0 = cosd(theta0)
|
|
Sin_theta0 = sind(theta0)
|
|
endif
|
|
|
|
phi0 = 360.*ran(seed)
|
|
Cos_phi0 = cosd(phi0)
|
|
Sin_phi0 = sind(phi0)
|
|
|
|
if (angle_offset) then
|
|
|
|
c -> Es soll aus gewuerfelter Startrichtung (theta0,phi0) und durch die Winkel-
|
|
c schleifen vorgegebenen Startrichtung (theta0_,phi0_) die tatsaechliche
|
|
c Startrichtung berechnet werden. Dafuer werden die gewuerfelten Winkel als
|
|
c 'Streuwinkel' betrachtet.
|
|
c Vorgehensweise:
|
|
c Es werden die Komponenten eines Geschwindigkeitsvektors mit Betrag=1 und durch
|
|
c theta0_,phi0_ bestimmter Richtung berechnet. Danach werden die Komponenten des
|
|
c mit theta0,phi0 gestreuten Geschwindigkeitsvektors und die zugehoerigen Winkel
|
|
c gewonnen, die dann als neudefinierte theta0 und phi0 fuer die tatsaechlichen
|
|
c Startwinkel verwendet werden. Das alles geschieht vollkommen analog zur
|
|
c Winkelaufstreuung in der Triggerfolie im Programm 'MUTRACK'.
|
|
c ('v' wird als Hilfsvariable missbraucht).
|
|
|
|
! Berechnen der 'Geschwindigkeitskomponenten':
|
|
v(1) = cosd(theta0_)
|
|
help1 = sind(theta0_)
|
|
v(2) = help1 * cosd(phi0_)
|
|
v(3) = help1 * sind(phi0_)
|
|
! v_xy ist stets groesser 0 ausser wenn die Zentralrichtung
|
|
! senkrecht nach oben oder unten gerichtet ist. Diese Wahl ist
|
|
! aber sowieso wenig sinnvoll:
|
|
v_xy = SQRT(v(1)*v(1) + v(2)*v(2))
|
|
if (v_xy.EQ.0.) then
|
|
write(*,*)
|
|
write(*,*)' Bei Zufallsverteilung fuer Startwinkel darf die durch die Winkelschleifen'
|
|
write(*,*)' vorgegebene Zentralrichtung nicht senkrecht nach oben oder nach unten weisen!'
|
|
write(*,*)' -> STOP'
|
|
STOP
|
|
endif
|
|
! berechne neue 'Geschwindigkeitskomponenten':
|
|
help1 = v(1)
|
|
help2 = v(2)
|
|
help3 = Sin_theta0*Cos_phi0/v_xy
|
|
help4 = Sin_theta0*Sin_phi0
|
|
v(1) = Cos_theta0*help1 - help3*help2 - help4*help1*v(3)/v_xy
|
|
if (v(1).LT.0.) goto 4462
|
|
v(2) = Cos_theta0*help2 + help3*help1 - help4*help2*v(3)/v_xy
|
|
v(3) = Cos_theta0*v(3) + help4*v_xy
|
|
! Berechne tatsaechlichen Startwinkel:
|
|
if (v(2).EQ.0. .AND. v(3).EQ.0.) then
|
|
if (v(1).GE.0) then
|
|
theta0 = 0.
|
|
else
|
|
theta0 = 180.
|
|
endif
|
|
phi0 = 0.
|
|
else
|
|
theta0 = acosd(v(1))
|
|
phi0 = atan2d(v(3),v(2))
|
|
if (phi0.LT.0) phi0 = phi0+360.
|
|
endif
|
|
Cos_theta0 = cosd(theta0)
|
|
Sin_theta0 = sind(theta0)
|
|
Cos_phi0 = cosd(phi0)
|
|
Sin_phi0 = sind(phi0)
|
|
endif
|
|
|
|
if (theta0.GT.90.) goto 4462
|
|
|
|
parWert(thetAng) = theta0
|
|
parWert(phiAng) = phi0
|
|
endif
|
|
|
|
c Bei Starts von Gitter 1 Startwinkel umrechnen (vielleicht auf
|
|
c Oberflaechennormale am Startort beziehen?)
|
|
|
|
if (StartFlaeche.EQ.1) then
|
|
!
|
|
!
|
|
!
|
|
endif
|
|
|
|
|
|
c Berechnen der Start-Geschwindigkeitskomponenten:
|
|
|
|
v0_Betrag = sqrt(E0/Energie_Faktor)
|
|
|
|
v0(1) = v0_Betrag * Cos_theta0
|
|
v0(2) = v0_Betrag * Sin_theta0 * Cos_phi0
|
|
v0(3) = v0_Betrag * Sin_theta0 * Sin_phi0
|
|
|
|
if (v0(1).LT.0) then
|
|
write(*,*)
|
|
write(*,*) ' >>>> v0(x) negativ!'
|
|
write(*,*)
|
|
STOP
|
|
endif
|
|
|
|
|
|
c Schreiben der Startparameter in das (temporaere) Startparameterfile:
|
|
|
|
write(lunStart) x0,v0,E0,theta0,phi0
|
|
|
|
|
|
c gegebenenfalls die Lebensdauer wuerfeln und in das entsprechende File schreiben:
|
|
|
|
if (UseDecay_) then
|
|
2453 lifeTime = -meanlifeTime * Log(Ran(seed) + 1.0E-37)
|
|
if (lifeTime.LE.0.) goto 2453
|
|
write(lunDecay) lifetime
|
|
endif
|
|
|
|
|
|
c die DEBUG-Daten ausgeben:
|
|
|
|
if (Debug .AND. start_nr.LE.DEBUG_Anzahl) call output_new_particle
|
|
|
|
|
|
c Festlegen des 1. Zeitschritts (beim ersten Versuch 0.1 mm in x-Richtung):
|
|
|
|
if (v0(1).GT.0) then
|
|
dt = 0.1/v0(1)
|
|
if (dt.LT.dtsmall) dt = .01
|
|
else
|
|
dt = .01
|
|
endif
|
|
|
|
|
|
c die auf x=40 mm extrapolierten Koordinaten resetten:
|
|
|
|
if (NTP_40mm) then
|
|
x40(2) = 0.
|
|
x40(3) = 0.
|
|
t40 = 0.
|
|
endif
|
|
|
|
|
|
c schreibe die relevanten Daten in das entsprechende Zwischenspeicherungsfile:
|
|
|
|
! Startortverteilung in y-z-Ebene muss auf Mappenverhaeltnisse skaliert werden:
|
|
help2 = x0(2) / scaleFactor
|
|
help3 = x0(3) / scaleFactor
|
|
write(lunZwi(lunIndx)) 0,0.,x0(1),help2,help3,v0,100+Mappe0, 0, 0,Gebiet0,Mappe0,dt
|
|
! = steps,t,x ,v ,destiny,zaehler,n_dtsmall,Gebiet ,lastMap ,dt
|
|
|
|
100 continue ! naechstes Startparameterset festlegen
|
|
|
|
if (graphics) call schnitt_plot ! Ausgabe der Graphik der Startverteilung
|
|
|
|
if (Debug) then
|
|
write(lunLOG,'(x,79(''-''))')
|
|
write(lunLOG,1001) 'STEP','T','X','Y','Z','Vx','Vy','Vz','E'
|
|
write(lunLOG,'(x,79(''-''))')
|
|
endif
|
|
1001 format (T2,A,T17,A,T25,A,T32,A,T39,A,T47,A,T54,A,
|
|
+ T61,A,T71,A)
|
|
|
|
|
|
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
c jetzt starten die Projektile:
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
|
|
mappe = mappe0
|
|
goto (1,2,3,4,5,6) mappe
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 1. Mappe (dx = 0.125 mm, dy = 0.125 mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
1 calcMap(1) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_1
|
|
else
|
|
call add_map_1
|
|
endif
|
|
start_Nr = 0
|
|
|
|
|
|
c Zurueckspulen der Zwischenspeicherungsfiles:
|
|
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
|
|
|
|
c Ausgabe der Prozentzahl schon gerechneter Trajektorien vorbereiten:
|
|
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0.
|
|
write(*,2001)'map1: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
2001 format ($,x,A,A,A)
|
|
|
|
|
|
c DEBUG beruecksichtigen:
|
|
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 1 <<<<<<<'
|
|
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c hier starten die Projektile in die 1. Potentialmappe:
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Einlesen der Trajektorienparameter:
|
|
|
|
11 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
|
|
|
|
c Ausgabe der Prozentzahl schon gerechneter Trajektorien:
|
|
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
2002 format ($,'+',I3)
|
|
|
|
|
|
c beruecksichtige nur die Teilchen, die sich gerade in der aktuellen Mappe auf-
|
|
c halten:
|
|
|
|
if (destiny.NE.100+mappe) goto 561
|
|
destiny = code_ok
|
|
reachedEndOfMap = .false.
|
|
|
|
|
|
c die DEBUG-Daten ausgeben:
|
|
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
|
|
|
|
c Graphikausgabe beruecksichtigen:
|
|
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0 ! Koordinatenzaehler resetten
|
|
call Save_Graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
|
|
|
|
c Falls Koordinaten im NTupel verlangt sind: Speicher resetten und Startkoordinaten
|
|
c sichern:
|
|
c
|
|
c if (NTP_Koord) then
|
|
c do k = 0, Gebiete_Anzahl ! Trajektorienspeicher resetten
|
|
c do i = 1, 8 ! (1 bis 8: x,y,z,t,vx,vy,vz,E)
|
|
c Koord_NTP(i,k) = 0
|
|
c enddo
|
|
c enddo
|
|
c call Save_NTP_Koord
|
|
c endif
|
|
|
|
|
|
c...............................................................................
|
|
c Integration:
|
|
c Im Bereich des Folienrandes ist es wegen der teilweise extremen
|
|
c Potentialverhaeltnisse kaum moeglich, eine etwas restriktivere Genauigkeits-
|
|
c forderung zu erfuellen. Daher wird im Fall von relativer Fehlerbetrachtung
|
|
c in diesem Bereich grundsaetzlich eine absolute Fehlerbetrachtung mit eps_x,
|
|
c eps_v = 1e-6 durchgefuehrt.
|
|
c (-> logical 'reducedAccur')
|
|
|
|
do while (.NOT.reachedEndOfMap)
|
|
if (log_relativ.AND.x(1).LE.xMarke .AND.
|
|
+ (Abs(x(2)).GE.yMarke .OR. abs(x(3)).GE.zMarke) ) then
|
|
reducedAccur = .true.
|
|
else
|
|
reducedAccur = .false.
|
|
endif
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_1(dt) ! setzt u.U. 'destiny'
|
|
if (destiny.NE.code_ok) goto 551
|
|
if (UseDecay_) call Decay_Test(*551)
|
|
|
|
yAbs = abs(x(2))
|
|
zAbs = abs(x(3))
|
|
|
|
c - Targethalter getroffen?:
|
|
if (x(1).LE.xEnd_TgtHolder .AND. (yAbs.GT.TgtConstY+x(1)*TgtFactorY
|
|
+ .OR. zAbs.GT.TgtConstZ+x(1)*TgtFactorZ)) then
|
|
destiny = code_hit_TgtHolder
|
|
goto 551
|
|
endif
|
|
|
|
if (Steps.GE.MaxStep) then
|
|
destiny = code_lost
|
|
goto 551
|
|
endif
|
|
|
|
|
|
c verarbeite alle 'imonitor' Schritte die Koordinaten fuer GRAPHICS und DEBUG:
|
|
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (GRAPHICS_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
|
|
enddo
|
|
|
|
|
|
c...............................................................................
|
|
551 continue
|
|
Gebiet = upToGrid1
|
|
if (destiny.NE.code_ok) then
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LT.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
elseif (reachedEndOfMap) then
|
|
if (x(1).NE.xEndMap1) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xEndMap1: x(1) = ',x(1)
|
|
x(1) = xEndMap1
|
|
endif
|
|
destiny = 100 + mappe+1
|
|
calcMap(mappe+1) = .true.
|
|
else
|
|
write(*,*)' main: fatal error marker 1'
|
|
STOP
|
|
endif
|
|
4456 format(x,' >> ',A,:,T40,'(',A,')')
|
|
|
|
|
|
c schreibe die Trajektoriendaten fuer GRAPHICS weg:
|
|
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
|
|
|
|
c schreibe aktuelle Daten ins Zwischenspeicherungsfile:
|
|
|
|
561 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
|
|
|
|
c -> das naechste Projektil kann kommen:
|
|
|
|
goto 11
|
|
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 2. Mappe (dx = 0.025 mm, dy = 0.025 mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
2 calcMap(2) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_2
|
|
else
|
|
call add_map_2
|
|
endif
|
|
start_Nr = 0
|
|
|
|
|
|
c Zurueckspulen der Zwischenspeicherungsfiles:
|
|
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
|
|
|
|
c Ausgabe der Prozentzahl schon gerechneter Trajektorien vorbereiten:
|
|
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0
|
|
write(*,2001)'map2: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
|
|
|
|
c DEBUG beruecksichtigen:
|
|
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 2 <<<<<<<'
|
|
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c hier starten die Projektile in die 2. Potentialmappe:
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Einlesen der aktuellen Trajektorienparameter:
|
|
|
|
22 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
|
|
|
|
c Ausgabe der Prozentzahl schon gerechneter Trajektorien:
|
|
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
|
|
|
|
c beruecksichtige nur die Teilchen, die sich gerade in der aktuellen Mappe auf-
|
|
c halten:
|
|
|
|
if (destiny.NE.100+mappe) goto 562
|
|
destiny = code_ok
|
|
if (startFlaeche.EQ.1 .AND. steps.EQ.0 .AND. x(1).GE.xPosition_Grid1) then
|
|
alreadyTested = .true.
|
|
else
|
|
alreadyTested = .false. ! Gittertreffer noch zu pruefen
|
|
endif
|
|
reachedEndOfMap = .false.
|
|
backOneMap = .false.
|
|
|
|
|
|
c die DEBUG-Daten ausgeben:
|
|
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
|
|
|
|
c Graphikausgabe beruecksichtigen:
|
|
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0
|
|
call Save_Graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
|
|
|
|
c...............................................................................
|
|
c Integration:
|
|
|
|
do while (.NOT.(backOneMap.OR.reachedEndOfMap))
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_2(dt)
|
|
if (destiny.NE.code_ok) goto 552
|
|
if (UseDecay_) call Decay_Test(*552)
|
|
|
|
yAbs = abs(x(2))
|
|
zAbs = abs(x(3))
|
|
|
|
c - wires of first grid: auf 2 Methoden testen:
|
|
c 1.: liegt aktueller Raumpunkt innerhalb eines Gitterstabes?
|
|
c 2.: Falls gerade die Position des Drahtgitters passiert wurde: schneidet die
|
|
c aktuelle Bahntangente die Oberflaeche eines Gitterstabes?
|
|
c (Die 2. Methode wurde eingebaut, um speziell bei Rechnungen mir kleinerer
|
|
c Genauigkeit sicherzustellen, dass Drahttreffer moeglichst auch dann als
|
|
c solche erkannt werden, wenn die Schrittweite im Gitterbereich groesser
|
|
c als der Drahtdurchmesser sein sollte).
|
|
help1 = abs(x(1)-xPosition_Grid1)
|
|
if (help1.LE.rWires1 .AND. yAbs.LE.y_Pos_lastWire1+rWires1 .AND.
|
|
+ zAbs.LE.innerDz_Gridframe1) then
|
|
ywire = nint((yAbs-y_Pos_firstWire1)/distance_wires1)*distance_wires1 + y_Pos_firstWire1
|
|
if ( (help1*help1 + (yAbs-ywire)*(yAbs-ywire)) .LE.
|
|
+ rQuadWires1) then
|
|
destiny = code_hit_grid1
|
|
goto 552
|
|
endif
|
|
endif
|
|
c Das jetzt folgende besser einmal aendern:
|
|
c Wenn entweder jetztige oder vorige Position oder beide im Bereich des Gitters
|
|
c sind: pruefe, ob Verbindungslinie der Punkte die Drahtoberflaeche schneidet.
|
|
c Auch Treffer der Rahmen etc ueber Schnitte mit den Verbindungslinien machen.
|
|
c (Dann wird ausgeschaltet, dass Teilchen von vor Hindernis bis hinter Hinderniss
|
|
c springt ohne dass Treffer bemerkt wird)
|
|
if (.NOT.alreadyTested .AND. x(1).GE.xPosition_Grid1 .AND.
|
|
+ yAbs.LE.y_Pos_lastWire1+rWires1 .AND. zAbs.LE.innerDz_Gridframe1) then
|
|
ywire = nint((x(2)-y_Pos_firstWire1)/distance_wires1)*distance_wires1 + y_Pos_firstWire1
|
|
distToWire(2) = x(2) - yWire
|
|
distToWire(1) = x(1) - xPosition_Grid1
|
|
call Test_WireHit(distToWire,rQuadWires1,v(1),v(2),WireHit)
|
|
if (WireHit) then
|
|
destiny = code_hit_grid1
|
|
goto 552
|
|
endif
|
|
alreadyTested = .true.
|
|
endif
|
|
c - Querbalken bei Gitter1:
|
|
if (x(1).GE.xStart_Balken .AND. x(1).LE.xEnd_Balken .AND.
|
|
+ yAbs.LE.Dy_Balken .AND.
|
|
+ zAbs.GE.innerDz_Balken .AND. zAbs.LE.outerDz_Balken) then
|
|
destiny = code_wand
|
|
goto 552
|
|
endif
|
|
c - frame of first grid:
|
|
if (x(1).GE.xStart_Gridframe1 .AND. x(1).LE.xEnd_Gridframe1 .AND.
|
|
+ (yAbs.GE.innerDy_Gridframe1 .OR. zAbs.GE.innerDz_Gridframe1) ) then
|
|
destiny = code_wand
|
|
goto 552
|
|
endif
|
|
|
|
if (Steps.GE.MaxStep) then
|
|
destiny = code_lost
|
|
goto 552
|
|
endif
|
|
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (Graphics_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
c...............................................................................
|
|
552 continue
|
|
if (x(1).LT.xPosition_Grid1) then
|
|
Gebiet = upToGrid1
|
|
else
|
|
Gebiet = upToGrid2
|
|
endif
|
|
if (destiny.NE.code_ok) then
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LT.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
elseif (BackOneMap) then
|
|
if (x(1).NE.xStartMap2) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xStartMap2: x(1) = ',x(1)
|
|
x(1) = xStartMap2
|
|
endif
|
|
destiny = 100 + mappe-1
|
|
calcMap(mappe-1) = .true.
|
|
elseif (reachedEndOfMap) then
|
|
if (x(1).NE.xEndMap2) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xEndMap2: x(1) = ',x(1)
|
|
x(1) = xEndMap2
|
|
endif
|
|
destiny = 100 + mappe+1
|
|
calcMap(mappe+1) = .true.
|
|
else
|
|
write(*,*)' main: fatal error marker 2'
|
|
STOP
|
|
endif
|
|
|
|
|
|
c schreib die Trajektoriendaten fuer GRAPHICS weg:
|
|
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
|
|
|
|
c schreibe aktuelle Daten ins Zwischenspeicherungsfile:
|
|
|
|
562 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
|
|
|
|
c -> das naechste Projektil kann kommen:
|
|
|
|
goto 22
|
|
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 3. Mappe (dx = 0.125mm, dy = 0.125mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
3 calcMap(3) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_3
|
|
else
|
|
call add_map_3
|
|
endif
|
|
start_Nr = 0
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0
|
|
write(*,2001)'map3: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 3 <<<<<<<'
|
|
c...............................................................................
|
|
33 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
if (destiny.NE.100+mappe) goto 563
|
|
destiny = code_ok
|
|
reachedEndOfMap = .false.
|
|
backOneMap = .false.
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0
|
|
call save_graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
c...............................................................................
|
|
do while (.NOT.(backOneMap.OR.reachedEndOfMap))
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_3(dt)
|
|
cd write(lunLog,*)'t,x,v = ',t,x,v
|
|
if (destiny.NE.code_ok) goto 553
|
|
if (UseDecay_) call Decay_Test(*553)
|
|
|
|
yAbs = abs(x(2))
|
|
zAbs = abs(x(3))
|
|
c - frame of first grid:
|
|
if (x(1).GE.xStart_Gridframe1 .AND. x(1).LE.xEnd_Gridframe1 .AND.
|
|
+ (yAbs.GE.innerDy_Gridframe1 .OR. zAbs.GE.innerDz_Gridframe1) ) then
|
|
destiny = code_wand
|
|
goto 553
|
|
endif
|
|
|
|
if (Steps.GE.MaxStep) then
|
|
destiny = code_lost
|
|
goto 553
|
|
endif
|
|
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (Graphics_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
enddo
|
|
c...............................................................................
|
|
553 continue
|
|
Gebiet = upToGrid2
|
|
if (destiny.NE.code_ok) then
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LT.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
elseif (BackOneMap) then
|
|
if (x(1).NE.xStartMap3) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xStartMap3: x(1) = ',x(1)
|
|
x(1) = xStartMap3
|
|
endif
|
|
destiny = 100 + mappe-1
|
|
calcMap(mappe-1) = .true.
|
|
elseif (reachedEndOfMap) then
|
|
if (x(1).NE.xEndMap3) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xEndMap3: x(1) = ',x(1)
|
|
x(1) = xEndMap3
|
|
endif
|
|
destiny = 100 + mappe+1
|
|
calcMap(mappe+1) = .true.
|
|
else
|
|
write(*,*)' main: fatal error marker 3'
|
|
STOP
|
|
endif
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
563 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
goto 33
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 4. Mappe (dx = 0.050 mm, dy = 0.050 mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
4 calcMap(4) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_4
|
|
else
|
|
call add_map_4
|
|
endif
|
|
start_Nr = 0
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0
|
|
write(*,2001)'map4: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 4 <<<<<<<'
|
|
c...............................................................................
|
|
44 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
if (destiny.NE.100+mappe) goto 564
|
|
destiny = code_ok
|
|
if (startFlaeche.EQ.2 .AND. steps.EQ.0 .AND. x(1).GE.xPosition_Grid2) then
|
|
alreadyTested = .true.
|
|
else
|
|
alreadyTested = .false. ! Gittertreffer noch zu pruefen
|
|
endif
|
|
reachedEndOfMap = .false.
|
|
backOneMap = .false.
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0
|
|
call save_graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
c...............................................................................
|
|
do while (.NOT.(backOneMap.OR.reachedEndOfMap))
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_4(dt)
|
|
if (destiny.NE.code_ok) goto 554
|
|
if (UseDecay_) call Decay_Test(*554)
|
|
|
|
yAbs = abs(x(2))
|
|
zAbs = abs(x(3))
|
|
c - wires of second grid:
|
|
help1 = abs(x(1)-xPosition_Grid2)
|
|
if (help1.LE.rWires2 .AND. yAbs.LE.y_Pos_lastWire2+rWires2 .AND.
|
|
+ zAbs.LE.innerDz_Gridframe2) then
|
|
ywire = nint((yAbs-y_Pos_firstWire2)/distance_wires2)*distance_wires2 + y_Pos_firstWire2
|
|
if ( (help1*help1 + (yAbs-ywire)*(yAbs-ywire)) .LE.
|
|
+ rQuadWires2) then
|
|
destiny = code_hit_grid2
|
|
goto 554
|
|
endif
|
|
endif
|
|
if (.NOT.alreadyTested .AND. x(1).GE.xPosition_Grid2 .AND.
|
|
+ yAbs.LE.y_Pos_lastWire2+rWires2 .AND. zAbs.LE.innerDz_Gridframe2) then
|
|
ywire = nint((x(2)-y_Pos_firstWire2)/distance_wires2)*distance_wires2 + y_Pos_firstWire2
|
|
distToWire(2) = x(2) - yWire
|
|
distToWire(1) = x(1) - xPosition_Grid2
|
|
call Test_WireHit(distToWire,rQuadWires2,v(1),v(2),WireHit)
|
|
if (WireHit) then
|
|
destiny = code_hit_grid2
|
|
goto 554
|
|
endif
|
|
alreadyTested = .true.
|
|
endif
|
|
c - frame of second grid:
|
|
if (x(1).GE.xStart_Gridframe2 .AND. x(1).LE.xEnd_Gridframe2 .AND.
|
|
+ (yAbs.GE.innerDy_Gridframe2 .OR. zAbs.GE.innerDz_Gridframe2) ) then
|
|
destiny = code_wand
|
|
goto 554
|
|
endif
|
|
|
|
if (Steps.GE.MaxStep) then
|
|
destiny = code_lost
|
|
goto 554
|
|
endif
|
|
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (Graphics_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
enddo
|
|
c...............................................................................
|
|
554 continue
|
|
if (x(1).LT.xPosition_Grid2) then
|
|
Gebiet = upToGrid2
|
|
else
|
|
Gebiet = upToHeShield
|
|
endif
|
|
if (destiny.NE.code_ok) then
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LT.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
elseif (BackOneMap) then
|
|
if (x(1).NE.xStartMap4) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xStartMap4: x(1) = ',x(1)
|
|
x(1) = xStartMap4
|
|
endif
|
|
destiny = 100 + mappe-1
|
|
calcMap(mappe-1) = .true.
|
|
elseif (reachedEndOfMap) then
|
|
if (x(1).NE.xEndMap4) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xEndMap4: x(1) = ',x(1)
|
|
x(1) = xEndMap4
|
|
endif
|
|
destiny = 100 + mappe+1
|
|
calcMap(mappe+1) = .true.
|
|
else
|
|
write(*,*)' main: fatal error marker 4'
|
|
STOP
|
|
endif
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
564 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
goto 44
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 5. Mappe (dx = 0.125 mm, dy = 0.125 mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
5 calcMap(5) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_5
|
|
else
|
|
call add_map_5
|
|
endif
|
|
start_Nr = 0
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0
|
|
write(*,2001)'map5: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 5 <<<<<<<'
|
|
c...............................................................................
|
|
55 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
if (destiny.NE.100+mappe) goto 565
|
|
destiny = code_ok
|
|
reachedEndOfMap = .false.
|
|
backOneMap = .false.
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0
|
|
call save_graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
c...............................................................................
|
|
do while (.NOT.(backOneMap.OR.reachedEndOfMap))
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_5(dt)
|
|
if (Steps.GE.MaxStep) destiny = code_lost
|
|
if (destiny.NE.code_ok) goto 555
|
|
if (UseDecay_) call Decay_Test(*555)
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (Graphics_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
enddo
|
|
c...............................................................................
|
|
555 continue
|
|
Gebiet = upToHeShield
|
|
if (destiny.NE.code_ok) then
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LT.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
elseif (BackOneMap) then
|
|
if (x(1).NE.xStartMap5) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xStartMap5: x(1) = ',x(1)
|
|
x(1) = xStartMap5
|
|
endif
|
|
destiny = 100 + mappe-1
|
|
calcMap(mappe-1) = .true.
|
|
elseif (reachedEndOfMap) then
|
|
if (x(1).NE.xEndMap5) then
|
|
write(*,*)' %%%%% ACCEL: x(1).NE.xEndMap5: x(1) = ',x(1)
|
|
x(1) = xEndMap5
|
|
endif
|
|
destiny = 100 + mappe+1
|
|
calcMap(mappe+1) = .true.
|
|
else
|
|
write(*,*)' main: fatal error marker 5'
|
|
STOP
|
|
endif
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
565 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
goto 55
|
|
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
c innerhalb der 6. Mappe (dx = 0.250 mm, dy = 0.250 mm, dz = 0.500 mm)
|
|
c:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
|
|
|
|
6 calcMap(6) = .false.
|
|
if (HVs_from_map) then
|
|
call read_map_6
|
|
else
|
|
call add_map_6
|
|
endif
|
|
start_Nr = 0
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
if (log_percent) then
|
|
call time(uhrzeit)
|
|
percent_done = 0
|
|
write(*,2001)'map6: ',Uhrzeit,' %: 0'
|
|
else
|
|
write(*,*) 'integrating ...'
|
|
endif
|
|
if (DEBUG) write (lunLOG,*) '>>>>>>> MAPPE 6 <<<<<<<'
|
|
c...............................................................................
|
|
66 read(lunZwi(lunIndx),END=99999) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
start_nr = start_nr + 1
|
|
if (log_percent) then
|
|
if (100.*real(start_nr)/real(n_par(0)).GE.percent_done+5) then
|
|
percent_done = percent_done + 5
|
|
write(*,2002) percent_done
|
|
endif
|
|
endif
|
|
if (destiny.NE.100+mappe) goto 566
|
|
destiny = code_ok
|
|
if (Debug .AND. start_Nr.LE.DEBUG_Anzahl) then
|
|
Debug_ = .true.
|
|
write(lunLOG,'(x,A,I4)') 'Teilchen Nr.',start_nr
|
|
call Output_Debug
|
|
else
|
|
Debug_ = .false.
|
|
endif
|
|
if (graphics .AND. start_Nr.LE.graphics_Anzahl) then
|
|
graphics_ = .true.
|
|
nKoord = 0
|
|
call save_graphics_Koord
|
|
else
|
|
graphics_ = .false.
|
|
endif
|
|
c...............................................................................
|
|
67 continue
|
|
call INTEGRATIONSSTEP_RUNGE_KUTTA_6(dt)
|
|
if (destiny.NE.code_ok) goto 556
|
|
if (UseDecay_) call Decay_Test(*556)
|
|
|
|
c - He-Shield getroffen bzw. He-Fenster erreicht?
|
|
radiusQuad = x(1)*x(1) + x(2)*x(2)
|
|
if (radiusQuad.GE.rQuadHeShield) then ! aufgeschlagen?
|
|
help1 = v(1)*v(1)+v(2)*v(2)
|
|
help2 = 2*(x(1)*v(1)+x(2)*v(2))
|
|
help3 = radiusQuad - rQuadHeShield
|
|
dt = (-help2+SQRT(help2*help2-4*help1*help3))/(2*help1)
|
|
t = t + dt
|
|
x(1) = x(1) + dt*v(1) ! den Ort berechnen, an dem
|
|
x(2) = x(2) + dt*v(2) ! das Teilchen auf das Rohr
|
|
x(3) = x(3) + dt*v(3) ! aufschlaegt
|
|
if (abs(x(2)).GE.dy_HeWindow .OR.
|
|
+ abs(x(3)).GE.dz_HeWindow) then
|
|
destiny = code_wand
|
|
endif
|
|
goto 556
|
|
endif
|
|
|
|
if (Steps.GE.MaxStep) then
|
|
destiny = code_lost
|
|
goto 556
|
|
endif
|
|
|
|
if (GRAPHICS_.or.Debug_) then
|
|
zaehler = zaehler + 1
|
|
if (zaehler.EQ.iMonitor) then
|
|
if (Graphics_) call Save_Graphics_Koord
|
|
if (Debug_) call Output_Debug
|
|
zaehler = 0
|
|
endif
|
|
endif
|
|
goto 67
|
|
c...............................................................................
|
|
556 continue
|
|
Gebiet = upToHeShield
|
|
if (v(1).LE.0. .and. StartFlaeche.EQ.0) destiny = code_reflektiert
|
|
c if (v(1).LE.0.) destiny = code_reflektiert
|
|
lastMap = mappe
|
|
if (debug_) then
|
|
if (destiny.LE.0) then
|
|
write(lunLog,4456) code_text(destiny)
|
|
else
|
|
write(lunLog,4456) code_text(destiny),Gebiet_text(Gebiet)
|
|
endif
|
|
endif
|
|
if (Graphics_) then
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
566 write(lunZwi(1-lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
goto 66
|
|
|
|
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
c Dieser Teil des Programmes bestimmt, in welcher Mappe jeweils weitergerechnet
|
|
c werden soll:
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
|
|
99999 continue
|
|
|
|
c Vertauschen der Indizierung fuer zu lesendes und zu schreibendes Zwischen-
|
|
c speicherungsfile (0->1, 1->0):
|
|
|
|
lunIndx = 1-lunIndx
|
|
|
|
|
|
c Falls immer noch Teilchen unterwegs sind, bestimme naechste abzuarbeitende
|
|
c Mappe und fahre fort zu integrieren:
|
|
|
|
do mappe = 1,6
|
|
if (calcMap(mappe)) goto (1,2,3,4,5,6) mappe
|
|
enddo
|
|
|
|
|
|
c falls nein:
|
|
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
c HIER IST DER PROGRAMMKODE FUER DIE INTEGRATION DER TRAJEKTORIEN
|
|
c BEENDET! -> SCHREIBE NTUPEL UND ACCEL-FILE:
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
|
|
rewind (lunStart) ! fuer Ausgabe in .ACCEL-file
|
|
rewind (lunZwi1)
|
|
rewind (lunZwi2)
|
|
if (UseDecay_) rewind (lunDecay)
|
|
|
|
do start_nr = 1, n_par(0)
|
|
|
|
c lies Daten ein:
|
|
|
|
if (UseDecay_) read(lunDecay) lifetime
|
|
read(lunZwi(lunIndx)) steps,t,x,v,destiny,zaehler,n_dtsmall,Gebiet,lastMap,dt
|
|
|
|
|
|
c fuer NTupel- und Graphikausgabe: falls das Teilchen noch lebt, interpoliere
|
|
c die Trajektorie bis x = 40 mm (eventuelle Treffer des LN-Schildes bleiben
|
|
c dabei unberuecksichtigt!):
|
|
|
|
if (destiny.EQ.code_ok) then
|
|
gebiet = upToLnShield
|
|
dt = ( 40-x(1) ) / v(1)
|
|
t40 = t + dt
|
|
x40(2) = x(2) + v(2)*dt
|
|
x40(3) = x(3) + v(3)*dt
|
|
endif
|
|
|
|
|
|
c Schreibe das Integrationsergebnis ins .ACCEL-file:
|
|
|
|
if (NTP_stop) Ekin=(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))*Energie_Faktor
|
|
read (lunStart) x0,v0,E0,theta0,phi0
|
|
if (scaleFactor.EQ.1) then
|
|
call HFNT(idNTP)
|
|
else
|
|
! Ort und Zeit muessen entsprechend skaliert werden:
|
|
t = t * scaleFactor
|
|
x(1) = x(1) * scaleFactor
|
|
x(2) = x(2) * scaleFactor
|
|
x(3) = x(3) * scaleFactor
|
|
call HFNT(idNTP)
|
|
! skaliere fuer Plot bzw. Debug-Zwecke wieder zurueck:
|
|
t = t / scaleFactor
|
|
x(1) = x(1) / scaleFactor
|
|
x(2) = x(2) / scaleFactor
|
|
x(3) = x(3) / scaleFactor
|
|
endif
|
|
|
|
|
|
c Plotte letztes Teilstueck der Trajektorie:
|
|
|
|
if (destiny.EQ.code_ok .AND. Graphics .AND.
|
|
+ start_Nr.LE.graphics_Anzahl) then
|
|
x(1) = 40
|
|
x(2) = x40(2)
|
|
x(3) = x40(3)
|
|
t = t40
|
|
nKoord = 0
|
|
call save_graphics_Koord
|
|
call Save_Graphics_Koord
|
|
call plot_trajectory
|
|
endif
|
|
|
|
|
|
c gib die Debug-Information aus:
|
|
|
|
if (debug_) call Output_Debug
|
|
|
|
|
|
c Zaehle mit, bei wie vielen Teilchen trotz dt<dtsmall der tolerierte Fehler
|
|
c ueberschritten wurde. Notiere auch, wieoft dies beim gleichen Teilchen maximal
|
|
c vorkam:
|
|
|
|
if (n_dtsmall.GT.0) then
|
|
dtsmall_counter = dtsmall_counter + 1
|
|
if (n_dtsmall.gt.n_dtsmall_Max) then
|
|
n_dtsmall_Max = n_dtsmall
|
|
endif
|
|
endif
|
|
|
|
|
|
c Zaehle mit, wie viele Trajektorien wegen steps>MaxStep abgebrochen werden:
|
|
|
|
if (destiny.EQ.code_lost) lost_counter = lost_counter + 1
|
|
|
|
|
|
c Zaehle zustaendigen 'destiny'-scaler hoch:
|
|
|
|
if (destiny.GT.0) destiny = destiny + (Gebiet-1)*highest_code_Nr
|
|
statDestiny(destiny) = statDestiny(destiny) + 1
|
|
|
|
if (destiny.EQ.code_ok) okStepsCounter = okStepsCounter + steps
|
|
|
|
enddo
|
|
|
|
|
|
c Schliesse die noch offenen und nicht mehr benoetigten files:
|
|
|
|
close (lunZwi1)
|
|
close (lunZwi2)
|
|
close (lunStart)
|
|
if (UseDecay_) close (lunDecay)
|
|
|
|
|
|
c Mittlere Anzahl an Integrationsschritten fuer Trajektorien mit destiny=0
|
|
c ausgeben:
|
|
|
|
if (statDestiny(code_ok).NE.0) then
|
|
write(*,'(xA,F7.2)')'Mittlere Anzahl an Integrationsschritten bis zum Ziel: ',
|
|
+ real(okStepsCounter)/real(statDestiny(code_ok))
|
|
endif
|
|
|
|
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
c ES FOLGEN SCHREIBEN DES LOGFILES, ERSTELLEN DES POSTSKRIPTS
|
|
C UND SPRUNG IN NEUE SCHLEIFE BZW. PROGRAMMENDE
|
|
czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz
|
|
|
|
|
|
c das Summary ausgeben und die Werte in die Tabellen schreiben:
|
|
c Falls nur ein Teilchenstart pro Schleife erfolgt, werte die Statistiken
|
|
c erst nach der letzten Schleife aus.
|
|
|
|
NotLastLoop = .NOT.(SchleifenNr.EQ.SchleifenZahl)
|
|
flag_ok = .NOT.(OneStartPerLoop.AND.NotLastLoop)
|
|
|
|
if (flag_ok) then
|
|
call eval_statistics
|
|
if (n_outWhere.GT.0 .OR. smallLogFile) call Summary
|
|
if (createTabellen .or. createPhysTab) call output_tabellen
|
|
endif
|
|
|
|
|
|
c das PostScript-file erstellen:
|
|
|
|
c Wird pro Schleife nur ein Teilchen gestartet ('OneStartPerLoop'; d.h. kein
|
|
c oder genau ein 'Zufallsstart'), so trage alle Trajektorien in die gleiche
|
|
c Graphik ein. Das Postskript braucht dann also erst bei der letzten Schleife
|
|
c erstellt zu werden:
|
|
|
|
if (GRAPHICS .AND. flag_ok) then
|
|
c call schnitt_plot ! Ausgabe der Graphik der Schnittebene
|
|
|
|
if (n_postSkript.LE.0) then
|
|
goto 396
|
|
elseif (n_postSkript.EQ.1) then
|
|
if (n_outWhere.LT.2) then
|
|
write(*,*)'.....................................'//
|
|
+ '.........................................'
|
|
write(*,'(2X,A18,I3,A,I3)')'Schleife ',
|
|
+ SchleifenNr,' von ',SchleifenZahl
|
|
endif
|
|
write(*,394)'(P) Ps-file erstellen',
|
|
+ '(R) Restliche ps-files erstellen'
|
|
write(*,394)'(N) ps-file Nicht erstellen',
|
|
+ '(K) Keine ps-files mehr erstellen'
|
|
write(*,394)'(G) Graphikausgabe beenden',
|
|
+ '(A) programm Abbrechen'
|
|
394 format(T6,A,T40,A)
|
|
|
|
27 write(*,395)' [RETURN] = (N) -> '
|
|
395 format($,x,A)
|
|
read(*,'(A)') helpChar
|
|
|
|
28 do i = 1,7 ! bis zu sechs blanks werden akzeptiert
|
|
ant = helpChar(i:i)
|
|
if (ant.NE.' ') goto 29
|
|
enddo
|
|
ant = 'n'
|
|
|
|
29 write(*,*)'==========================='//
|
|
+ '====================================================='
|
|
|
|
if (ant.EQ.'n' .OR. ant.EQ.'N') then
|
|
goto 396
|
|
elseif (ant.EQ.'r' .OR. ant.EQ.'R') then
|
|
n_postSkript = 2
|
|
elseif (ant.EQ.'k' .OR. ant.EQ.'K') then
|
|
n_postSkript = 0
|
|
goto 396
|
|
elseif (ant.EQ.'g' .OR. ant.EQ.'G') then
|
|
call HPLEND
|
|
GRAPHICS = .false.
|
|
goto 200
|
|
elseif (ant.EQ.'a' .OR. ant.EQ.'A') then
|
|
call HPLEND
|
|
call TERMINATE_OUTPUT
|
|
STOP
|
|
elseif (ant.NE.'p' .AND. ant.NE.'P') then
|
|
goto 27
|
|
endif
|
|
endif
|
|
|
|
write (helpChar(1:7),'(''_'',I6)') SchleifenNr
|
|
if (filename.NE.' ') then
|
|
call MAKE_PS(filename//helpChar)
|
|
else
|
|
call MAKE_PS('ACCEL'//helpChar)
|
|
endif
|
|
|
|
|
|
396 continue
|
|
|
|
CALL IZPICT ('CHAMBER','S') ! LOESCHEN DER BILDER IM PAWC-COMMON-BLOCK
|
|
CALL IZPICT ('HISTO','S')
|
|
CALL IZPICT ('TEXT','S')
|
|
|
|
call iclrwk (1,0) ! CLEAREN DER PAW-'WORKSTATIONS'
|
|
call iclrwk (4,0)
|
|
call iclrwk (5,0)
|
|
|
|
CALL HRESET (50,' ') ! RESETTEN DES HISTOGRAMMS
|
|
|
|
endif
|
|
|
|
c -> das gleiche von vorne mit neuen Settings (d.h. neue Schleife)
|
|
|
|
200 continue
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
|
|
c Jetzt sind alle Schleifen abgearbeitet -> fertigmachen zum Programmende:
|
|
|
|
if (Graphics) call HPLEND ! Beende HIGZ Graphikbibliothek
|
|
|
|
call TERMINATE_OUTPUT
|
|
|
|
|
|
END
|
|
|
|
|
|
C===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE Lambert_Verteilung(n_Lambert,cos_theta,sin_theta)
|
|
c ============================================================
|
|
|
|
IMPLICIT NONE
|
|
|
|
real cos_theta,sin_theta
|
|
|
|
real n_Lambert ! Ordnung der Lambert-Verteilung
|
|
real dummy
|
|
integer seed
|
|
common /seed/ seed
|
|
|
|
dummy = ran(seed)
|
|
|
|
if (n_Lambert.EQ.0.) then
|
|
cos_theta = (1.-dummy)
|
|
sin_theta = sqrt(1.-cos_theta*cos_theta)
|
|
elseif (n_Lambert.EQ.1.) then
|
|
cos_theta = sqrt(1.-dummy)
|
|
sin_theta = sqrt(dummy)
|
|
else
|
|
cos_theta = (1.-dummy)**(1./(n_Lambert + 1))
|
|
sin_theta = sqrt(1.-cos_theta*cos_theta)
|
|
endif
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE Gauss_Verteilung(sigma,wert)
|
|
c =======================================
|
|
|
|
IMPLICIT NONE
|
|
|
|
real sigma ! Breite der Gaussverteilung
|
|
real wert ! gewuerfelte Returnvariable
|
|
real radius,phi
|
|
|
|
integer seed
|
|
common /seed/ seed
|
|
|
|
c Da die eindimensionale Gaussfunktion nicht integrierbar ist, wird erst
|
|
c ein Punkt in der Ebene mit der entsprechenden zweidimensionalen Gaussfunktion
|
|
c gewuerfelt. Von diesem Punkt wird dann die x-Komponente zurueckgegeben, die
|
|
c eindimensional Gaussverteilt ist:
|
|
|
|
radius = sigma*Sqrt(-2.*log(1.-ran(seed)))
|
|
phi = 360.*ran(seed)
|
|
wert = radius * cosd(phi)
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE Gauss_Verteilung_theta(sigma,theta)
|
|
c ==============================================
|
|
|
|
IMPLICIT NONE
|
|
|
|
real sigma,theta
|
|
real radius,phi,ratio
|
|
|
|
integer i, seed
|
|
common /seed/ seed
|
|
|
|
c Man beachte, dass hier Winkel gewuerfelt werden! D.h., dass die Variable
|
|
c 'radius' einen Radius in einer 2dimensionalen 'Winkel'-Ebene darstellt.
|
|
c Es wird angenommen, dass sigma in degree angegeben wird (daher die sind()-
|
|
c Funktion in der Zuweisung fuer 'ratio' anstelle der sin()-Fkt.).
|
|
|
|
i = 1
|
|
|
|
10 radius = sigma*Sqrt(-2.*log(1.-ran(seed)))
|
|
phi = 360.*ran(seed)
|
|
theta = abs(radius * cosd(phi))
|
|
! nur theta zwischen 0 und 90 deg sollen eine Chance haben:
|
|
if (theta.GT.90) then
|
|
i = i + 1
|
|
if (i.LE.10000) then
|
|
goto 10
|
|
else
|
|
write(*,*)
|
|
write(*,*) 'SUBROUTINE Gauss_Verteilung_theta:'
|
|
write(*,*) ' Nach 10000 Versuchen noch keinen Winkel < 90 degree gewuerfelt.'
|
|
write(*,*) ' Vorgegebenes Sigma der Winkelverteilung: ',sigma
|
|
write(*,*)
|
|
STOP
|
|
endif
|
|
endif
|
|
|
|
c Zitat aus TP's 'TESTSEED.FOR', aus welchem diese Routine abgeschrieben
|
|
c ist:
|
|
c
|
|
c Now we habe a GAUSSIAN, but we need for multiple scattering
|
|
c GAUSSIAN*SIN(x) =: g(x). This is not integrateable analytically, but
|
|
c we can choose the VON NEUMANN REJECTION to get what we want.
|
|
c As auxiliary function we choose the GAUSSIAN =: f(x), because it
|
|
c satisfies g(x) <= f(x) for all x.
|
|
c We must build the ratio g(x)/f(x) = sin(x) and compare it to
|
|
c another random number:
|
|
|
|
ratio = sind(theta)
|
|
if (ran(seed).GT.ratio) goto 10 ! Verteilung zurechtbiegen
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE reset_statistics
|
|
c ===========================
|
|
|
|
IMPLICIT NONE
|
|
|
|
integer Nr,k
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c der allgemeine Statistikspeicher: (*) : braucht nicht resettet zu werden
|
|
c ---------------------------------
|
|
c
|
|
c statMem(1,Nr): 1. Wert: x(1) (*)
|
|
c statMem(2,Nr): Summe_ueber_i( x(i)-x(1) )
|
|
c statMem(3,Nr): Summe_ueber_i( (x(i)-x(1))**2. )
|
|
c statMem(4,Nr): kleinster Wert
|
|
c statMem(5,Nr): groesster Wert
|
|
c statMem(6,Nr): Mittelwert (*)
|
|
c statMem(7,Nr): Varianz (*)
|
|
c statMem(8,Nr): Anzahl der Werte
|
|
c statMem(9,Nr): Anzahl der Werte in Prozent von 'StartsProSchleife' (*)
|
|
c ('StartsProSchleife' == n_par(0))
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c Ergebnis-Statistik-Speicher resetten:
|
|
|
|
do Nr = 1, stat_Anzahl
|
|
statMem(2,Nr) = 0. ! Summe der Werte
|
|
statMem(3,Nr) = 0. ! Summe der Quadrate
|
|
statMem(4,Nr) = 1.e10 ! Minimalwert
|
|
statMem(5,Nr) = -1.e10 ! Maximalwert
|
|
statMem(8,Nr) = 0. ! Anzahl
|
|
enddo
|
|
|
|
|
|
c StartZaehler resetten:
|
|
|
|
start_nr = 0
|
|
|
|
|
|
c der Statistikspeicher fuer das Teilchen-Schicksal:
|
|
|
|
do k = smallest_code_Nr, Gebiete_Anzahl*highest_code_Nr
|
|
statDestiny(k) = 0
|
|
enddo
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE fill_statMem(wert,Nr)
|
|
c ================================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
real wert
|
|
integer Nr
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
c Wird die Varianz der Verteilung einer Groesse x gemaess der Formel
|
|
c
|
|
c Var(x) = SQRT( <x**2> - <x>**2 ) , < > -> Erwartungswert
|
|
c
|
|
c mit
|
|
c <x> = 1/n * Summe_ueber_i( x(i) )
|
|
c <x**2> = 1/n * Summe_ueber_i( x(i)**2 )
|
|
c
|
|
c berechnet, so tritt manchmal aufgrund der beschraenkten Genauigkeit der
|
|
c numerischen Speicher das Problem auf, dass bei grossen Werten x(i) und
|
|
c kleiner Streuung der Ausdruck unter der Wurzel negativ wird, was erstens
|
|
c unphysikalisch ist und zweitens zum Programmabbruch fuehrt.
|
|
c
|
|
c Dieses Problem liesse sich vermeiden, wenn man die Groessen x(i) relativ
|
|
c zu ihrem Erwartungswert angeben wuerde, der aber erst im nachhinein bekannt
|
|
c ist.
|
|
c
|
|
c Als Naeherungsloesung verwende ich daher fuer die Berechnung der Varianz die
|
|
c x(i) relativ zu x(1), also zum ersten Wert gemessen, der gerade bei kleiner
|
|
c Streuung, bei der das numerische Problem auftritt, nahe am Erwartungswert
|
|
c liegen sollte.
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
c 1. Zaehle mit.
|
|
c 2. Speichere den ersten Wert.
|
|
c 3. Summiere die Abweichungen vom ersten Wert:
|
|
c 4. Summiere die Quadratischen Abweichungen vom ersten Wert:
|
|
c 5. Speichere den kleinsten Wert (wurde noch kein Wert aufgenommen, so ist
|
|
c statMem(4,Nr) = 1.e10):
|
|
c 6. Speichere den groessten Wert (wurde noch kein Wert aufgenommen, so ist
|
|
c statMem(5,Nr) = -1.e10):
|
|
|
|
statMem(8,Nr) = statMem(8,Nr) + 1. ! 1.
|
|
if (statMem(8,Nr).EQ.1) statMem(1,Nr) = wert ! 2.
|
|
statMem(2,Nr) = statMem(2,Nr) + (wert-statMem(1,Nr)) ! 3.
|
|
statMem(3,Nr) = statMem(3,Nr) + (wert-statMem(1,Nr))**2. ! 4.
|
|
if (statMem(4,Nr).GT.wert) statMem(4,Nr) = wert ! 5.
|
|
if (statMem(5,Nr).LT.wert) statMem(5,Nr) = wert ! 6.
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE eval_statistics
|
|
c ==========================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
real n ! Anzahl der Werte, == statMem(8,Nr)
|
|
real radiant
|
|
|
|
integer Nr
|
|
|
|
do Nr = 1, Stat_Anzahl
|
|
if (statNeeded(Nr)) then
|
|
n = statMem(8,Nr)
|
|
if (n.ne.0.) then
|
|
! Berechne Mittelwert:
|
|
statMem(6,Nr) = statMem(2,Nr)/n + statMem(1,Nr)
|
|
! Berechne Varianz:
|
|
radiant = ( statMem(3,Nr) - (statMem(2,Nr)**2. )/n)/n
|
|
statMem(7,Nr) = sqrt(radiant)
|
|
! Berechne Anteil an allen Gestarteten in Prozent:
|
|
statMem(9,Nr) = 100.*n/real(n_par(0))
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE SAVE_NTP_KOORD
|
|
c =========================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
real Ekin
|
|
|
|
Ekin = (v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) * Energie_Faktor
|
|
|
|
Koord_NTP(1,Gebiet) = x(1)
|
|
Koord_NTP(2,Gebiet) = x(2)
|
|
Koord_NTP(3,Gebiet) = x(3)
|
|
|
|
Koord_NTP(4,Gebiet) = t
|
|
|
|
Koord_NTP(5,Gebiet) = v(1)
|
|
Koord_NTP(6,Gebiet) = v(2)
|
|
Koord_NTP(7,Gebiet) = v(3)
|
|
|
|
Koord_NTP(8,Gebiet) = Ekin
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE SAVE_GRAPHICS_KOORD
|
|
c ==============================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
c Variablen fuer die Graphikausgabe:
|
|
|
|
real xKoord(1000) ! Koordinatenfelder fuer die
|
|
real yKoord(1000) ! Graphikausgabe
|
|
real zKoord(1000) !
|
|
integer nKoord ! Anzahl der Koordinaten
|
|
|
|
c dieser Commonblock ist auch im Hauptprogramm sowie in den Routinen der Datei
|
|
c 'SUB_MUPIC.FOR' definiert
|
|
|
|
COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord ! fuer Graphikaufruf
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
nKoord = nKoord + 1
|
|
|
|
xKoord(nKoord) = x(1)*scaleFactor
|
|
yKoord(nKoord) = x(2)*scaleFactor
|
|
zKoord(nKoord) = x(3)*scaleFactor
|
|
|
|
if (nKoord.EQ.1000) then
|
|
call plot_trajectory
|
|
xKoord(1) = xKoord(1000)
|
|
yKoord(1) = yKoord(1000)
|
|
zKoord(1) = zKoord(1000)
|
|
nKoord = 1
|
|
endif
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE Output_Debug
|
|
c =======================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
INCLUDE 'accel$sourcedirectory:COM_GEO.INC'
|
|
|
|
|
|
real Ekin
|
|
integer iostat
|
|
|
|
Ekin = (v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) * Energie_Faktor
|
|
write(lun(1),1,iostat=iostat) steps,t,x,v,Ekin
|
|
|
|
1 format(X,I4,7X,F6.1,2X,F7.2,X,F6.2,X,F6.2,2X,F6.2,X,
|
|
+ F6.2,X,F6.2,2X,G13.6)
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE Decay_Test(*)
|
|
c ========================
|
|
|
|
IMPLICIT NONE
|
|
|
|
INCLUDE 'accel$sourcedirectory:COM_ACCEL.INC'
|
|
|
|
real dt
|
|
|
|
if (scaleFactor*t.GT.lifeTime) then ! Teilchen zerfallen
|
|
dt = t - lifeTime/scaleFactor
|
|
t = lifeTime/scaleFactor
|
|
x(1) = x(1) - dt*v(1)
|
|
x(2) = x(2) - dt*v(2)
|
|
x(3) = x(3) - dt*v(3)
|
|
destiny = code_decay
|
|
RETURN 1
|
|
endif
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|
|
|
|
|
|
OPTIONS /EXTEND_SOURCE
|
|
|
|
SUBROUTINE test_wireHit(distToWire,WireRadiusQuad,v_x,v_y,WireHit)
|
|
c ==================================================================
|
|
|
|
c Diese Routine ueberprueft, ob bei gegebenem Abstandsvektor 'distToWire'
|
|
c zwischen Drahtposition und Teilchenort und gegebener Geschwindigkeit v eines
|
|
c Teilchens bei geradliniger Bewegung der Draht mit Radius WireRadius getroffen
|
|
c wird.
|
|
|
|
IMPLICIT NONE
|
|
|
|
real DistToWire(2),WireRadiusQuad,v_x,v_y
|
|
logical WireHit
|
|
|
|
real steigung, help, radiant
|
|
|
|
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
|
|
if (v_x.GT.v_y) then
|
|
steigung = v_y/v_x
|
|
help = distToWire(2) - distToWire(1) * steigung
|
|
radiant = (1+steigung*steigung)*WireRadiusQuad - help*help
|
|
else
|
|
steigung = v_x/v_y
|
|
help = distToWire(1) - distToWire(2) * steigung
|
|
radiant = (1+steigung*steigung)*WireRadiusQuad - help*help
|
|
endif
|
|
|
|
if (radiant.ge.0) then
|
|
wireHit = .true.
|
|
else
|
|
wireHit = .false.
|
|
endif
|
|
|
|
|
|
END
|
|
|
|
|
|
c===============================================================================
|