c c------------------------------------------------------------------------------ c c Changes starting on 17-Oct-2000, TP, PSI c c - add position of muons at TD foil and position of c foil electrons when hitting MCP3 to NTuple c - start changes for migration to NT and Unix; avoid using c logicals or environment variables; cancel OPTIONS/EXTEND_SOURCE; c use lower case filenames always c c****************************************************************************** c* ... MUTRACK.FOR (Stand: Februar '96) * c* * c* Dieses Programm integriert Teilchenbahnen in der UHV-Kammer der NEMU- * c* Apparatur. Startpunkte koennen zwischen der Moderatorfolie und dem MCP2 * c* frei gewaehlt werden, Endpunkt der Berechnungen ist (sofern die Teilchen * c* nicht vorher schon ausscheiden) die Ebene des MCP2. Bis jetzt koennen * c* also nur Bewegungen in Strahlrichtung, nicht entgegen derselben berechnet * c* werden. * c* Das Programm selbst rechnet den zweistufigen Beschleuniger als ideal, * c* bietet aber die Moeglichkeit Simulationen von TP oder AH (Programm 'Accel')* c* mit realem Beschleuniger einzulesen. Die Integration der Teilchenbahnen * c* erstreckt sich bei diesen Simulationen bis etwa zum He-Schild, MUTRACK * c* rechnet dann von dort aus weiter. * c* Verschiedene Einstellungen koennen in ineinandergreifenden Schleifen in * c* aequidistanten Schritten variiert werden (z.B. Spannungen des Transport- * c* Systems, Startgroessen der Teilchen, Masse und Ladung ...). Ein Teil dieser* c* Groessen kann aber auch alternativ nach verschiedenen frei waehlbaren * c* Zufallsverteilungen gewurfelt werden. * c* Die Integrationsergebnisse koennen in der Form von NTupeln abgespeichert * c* werden, was sie der Darstellung und Auswertung mit dem CERN-Programm PAW * c* zugaenglich macht. * c* Neben der reinen Integrationsarbeit fuehrt Mutrack Statistiken ueber * c* verschiedene Groessen (z.Z. verschiedene Flugzeiten und Ortsverteilungen) * c* die Mittelwerte und Standandartabweichungen sowie Minimal- und Maximalwerte* c* umfassen. * c* Diese Groessen koennen einfach ausgegeben oder in einem Tabellenfile abge- * c* speichert werden, welches von PHYSICA mittels der Fortran-Routine * c* 'READ_DATA' eingelesen werden kann. Verschiedene PHYSICA-Makros * c* (.PCM-files) ermoeglichen dann die Darstellung dieser statistischen * c* Groessen in Form von 2D- und 3D-Graphiken. (z.B. Abhaengigkeit der Trans- * c* mission von HV-Settings des Transportsystems). * c* Die momentan vorhandenen Routinen heissen * c* * c* MUINIT.PCM * c* HELP.PCM * c* MUPLOT_1DIM.PCM * c* MUPLOT_2DIM.PCM * c* TYPE_LOGHEADER.PCM * c* TYPE_PARAMS_GRAPHIC.PCM * c* TYPE_PARAMS_TEXT.PCM * c* * c* Nach dem Start (von dem Directory aus, in dem obige Routinen abgelegt sind)* c* muss PHYSICA mit dem Befehl '@MUINIT' initialisiert werden. Danach koennen * c* obige Routinen ueber Aliasse angesprochen werden. Weitere Informationen * c* hierzu erhaelt man, indem man in PHYSICA nach der Initialisierung 'MUHELP' * c* eingibt. * c* Der Sourcecode fuer Mutrack ist ueber verschiedene .FOR-Dateien verteilt, * c* die jeweils zu einem Problembereich gehoerige Subroutinen enthalten. Die * c* zur Zeit vorhandenen Files und die darin enthaltenen Routinen sind: c* c* MUTRACK.FOR c* SUB_ARTLIST.FOR c* SUB_OUTPUT.FOR c* SUB_INPUT.FOR c* SUB_INTEGR_FO.FOR c* SUB_INTEGR_L1.FOR c* SUB_INTEGR_L3.FOR c* SUB_INTEGR_M2.FOR c* SUB_PICTURE.FOR c* SUB_TRIGGER.FOR c* c* c* Includefiles mit COMMON-Blöcken: c* c* COM_DIRS.INC c* COM_KAMMER.INC c* COM_LUNS.INC c* COM_MUTRACK.INC c* COM_OUTPUT.INC c* COM_TD_EXT.INC c* COM_TD_INT.INC c* COM_WINKEL.INC c* GEO_KAMMER.INPUT c* GEO_TRIGGER.INC c* c* c* Icludefile mit Defaultwerten fuer eine Reihe benutzerdefinierbarer und Programm- c* interner Groessen: c* c* INITIALIZE.INC c* c* c* Includefiles fuer die Potentialmappen: c* c* MAP_DEF_FO.INC c* MAP_DEF_L1.INC c* MAP_DEF_L3.INC c* MAP_DEF_M2.INC c* c* READ_MAP.INC c* c* c* Benoetigte Eingabfiles: c* c* MUTRACK.INPUT (fuer die Integrationen zu verwendende Einstellungen) c* kammer_geo.INPUT (Spezifizierung der Kammergeometrie) c* mappenName.INFO (Dateien mit Angaben ueber zugehoerige Potentialmappen) c* mappenName.MAPPE (die Potentialmappen) c* MUTRACK_NR.DAT (zuletzt vergebene Nummern der Ausgabedateien, wird c* von Mutrack verwaltet). c* c* c* Fuer die Erstellung der Potentialmappen mit dem Triumf-Programm stehen folgende c* Hilfsmittel zur Verfuegung: c* c* BESCHL-INIT.FOR c* LENSE-INIT.FOR c* c* Diese Boundary-Routinen stellen folgende Moeglichkeiten zur Verfuegung: c* c* Initialisierung von Scratch, von 2D und von 3D-Mappen. Kontrollmoeglichkeiten c* ueber die Ausgabe der Potentialbereiche. c* c* Die Mappen koennen von PHYSICA aus mittels der FORTRAN-Routine ' ' c* und den .PCM-Makros ' ' ... angeschaut und ausgegeben werden. c* c* c* c* Liste der moeglichen Ausgabefiles: c* c* MU_nnnn.LOG c* MU_nnnn.GEO c* MU_nnnn.PHYSICA c* MU_nnnn.NTP c* MU_nnnn._tab c* c* Diese Version von MUTRACK enthaelt nur noch rudimentaere Anteile des ursprueng- c* lichen Programmes von Thomas Wutzke. Hauptunterschiede und Erweiterungen sind: c* c* # Ersetzen der Euler-Integration durch ein schrittweitenkontrolliertes c* Runge-Kutta Verfahren. Der dieser Implementation zugrundeliegende Algo- c* rythmus entstammt dabei dem Buch 'NUMERICAL RECIPES, The Art of Scientific c* Computing' (Fortran Version) von Press, Flannery, Teukolsky und Vetterling, c* Cambridge University Press (1989). c* c* # Verbesserter Algorythmus zur Berechnung der Feldstaerken aus den Potential- c* Mappen. c* c* # Implementierung des gesamten Statistikaparates. (Zuvor waren PAW-Ntupel die c* einzige Ausgabeform abgesehen von den Debuginformationen). c* c* # Uebersichtlichere Gestalltung der Ein- und Ausgabe, sowie der Debug-Infos. c* c* # Implementierung der Moeglichkeit, verschiedenen Parameter in Schleifen zu c* durchlaufen. c* c* # Implementierung der fuer die graphische Darstellung mit PHYSICA notwendigen c* Routinen. c* c* # Implementierung des Triggerdetektors. c* c* # Implementierung der Graphikausgabe der Teilchenbahnen (diese Routinen wurden c* in ihrer ersten Fassung von Michael Birke geschrieben). c* c* # Umstellen der Potentialmappen auf 'unformattiert' und Einschraenken der c* Mappen auf den wirklich benoetigten Bereich (d.h. z.B. Ausnutzen der c* Symmetrie der Linsen, wodurch die Mappengroesse bei den Linsen mehr als c* halbiert werden konnte. c* c* # Implementierung der Moeglichkeit, die Kammergeometrie (d.h. die Positionen c* der verwendeten Elemente) sowie die Potentialmappen (z.B. fuer unter- c* schiedliche Linsenkonfigurationen) ueber ein .INPUT-Eingabefile ohne c* Umschreiben des Sourcecodes aendern zu koennen. c* c* Das Programm verwendet fuer Graphikdarstellung und NTupel-Erzeugung Routinen der c* zum PAW-Komplex gehoerenden CERN-Bibliotheken 'HPLOT' und 'HBOOK'. c* c* Am Anfang der Deatei 'COM_MUTRACK.INC' findet sich eine Liste der wichtigsten c* Aenderungen ueber die verschiedenen Versionen ab 1.4.1. c* c* Gruss, Anselm Hofer c****************************************************************************** c C =============== program MUTRACK C =============== c Deklarationen: Implicit None INCLUDE 'com_mutrack.inc' INCLUDE 'com_dirs.inc' INCLUDE 'com_td_ext.inc' INCLUDE 'com_winkel.inc' INCLUDE 'com_kammer.inc' INCLUDE 'geo_trigger.inc' c die SCHLEIFENVARIABLEN fuer die 'do 200 ...'-Schleifenund und damit c zusammenhaengendes (Common-Bloecke werden fuer die NTupel-Ausgabe benoetigt): c - 'virtuelle' Flugstreckenverlaengerungen: real delta_L1,delta_L2 c - Energieverlust in der Triggerfolie und Dicke derselben: real E_loss c - Drehwinkel: c (alfaTgt, alfaSp, alfaTD und ihre Winkelfunktionen werden in 'COM_WINKEL.INC' c erledigt: COMMON /ANGELS/) real y_intersectSP ! Benoetigt fuer Schnittpkt. der Trajektorie real yUppLeft, yLowLeft ! mit forderer Spiegelebene real x_intersectTD ! Benoetigt fuer Schnittpkt. der Trajektorie ! mit TD-Folie real x_intersectTDMap ! ... mit TD-Mappe common /x_intersectTD/ x_intersectTD,x_intersectTDMap c - Masse und Ladung: real m, m_ ! Masse, Laufvariable fuer Massen-Schleife real q, q_ ! Ladung, Laufvariable fuer Ladungs-Schleife integer qInt COMMON /charge/ qInt ! fuer 'NTP_charge' integer nNeutral,nCharged ! fuer Ausgabe des gewuerfelten neutralen anteils COMMON /nNeutral/ nNeutral,nCharged c - MCP2: real U_MCP2 ! Spannung am MCP2 c - Triggerdetektor: U_F, U_V, U_H und U_MCP3 werden in 'COM_TD_EXT.INC' c erledigt. (COMMON /TRIGGERSETTINGS/) c - Transportsystem: real U_Tgt ! Target-Spannung real U_Gua ! Spannung am Guardring real U_G1 ! Spannung am ersten Gitter real U_L1 ! Spannung an Linse 1 real U_Sp ! Spiegelspannung real U_L2 ! Spannung an Linse 2 real U_L3 ! Spannung an Linse 3 COMMON /U_L2/ U_L2 ! fuer die Addition der 'L2andFo'-Mappe real last_U_L2 / -1.E10 / ! fuer die Addition der 'L2andFo'-Mappe real last_U_F / -1.E10 / c - Magnetfeldstaerken: real B_Helm ! Magnetfeld der Helmholtzspulen real B_TD ! Magnetfeld der Kompensationsspule am TD c - Startparameter: 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_MUTRACK.INC' declariert c allgemeine Trajektoriengroessen real dt ! zeitl. Aenderung real v_xy ! Geschwindigkeit in x/y-Ebene real v_square, v0_Betrag, v_Betrag real Ekin ! kinetische Energie real a1,a2 ! Beschleunigung in 1. bzw. 2. Beschl.Stufe real aFoil ! Beschleunigung zwischen Massegitter und Folie real radiusQuad ! RadiusQuadrat real radiusQuad_ ! RadiusQuadrat real radius real S1xM2 ! Zeit vom Start bis zur MCP2-Ebene real S1M2 ! Zeit vom Start bis zum MCP2 (Treffer voarausgesetzt) real S1Fo ! Zeit vom Start bis zur Folie real S1FoOnly ! Zeit vom Start bis zur Folie real FoM2 ! Zeit zwischen Folie und MCP2 real FoM2Only ! wie FoM2, falls keine anderen TOFs verlangt real S1M3 ! Zeit vom Start bis Eintreffen der FE auf MCP3 real M3M2 ! Zeit vom Eintreffen der FE auf MCP3 bis MCP2 real alfa ! Bahnwinkel gegen die Triggerfolienebene real E_Verlust /0./ ! Energieverlust in der Folie real delta_E_Verlust ! Streuung des Energieverlustes in der Folie real thetaAufstreu ! Ablenkung aus vorheriger Richtung in der Folie real phiAufstreu ! azimuthaler Winkel der Ablenkung gegenueber Horiz. COMMON /FOLIE/ E_Verlust,thetaAufstreu,phiAufstreu real Beschl_Faktor ! Faktor bei Berechn. der Beschleunigung im EFeld COMMON /BESCHL_FAKTOR/ Beschl_Faktor real length1 ! = d_Folie_Achse + MappenLaenge_FO real length2 ! = xTD - d_Folie_Achse - MappenLaenge_FO ! = xTD-length1 c Groessen der Folienelektronen ('FE'): integer nFE ! jeweilige Anzahl an FE (2 <= nFE <= 5) real E0FE ! Startenergie der Folienelektronen real ct0,st0,cf0,sf0 ! die Winkelfunktionen der Startwinkel der FE real f0 ! 'phi0' fuer die FE real x0FE(3) ! Startort der Folienelektronen auf der TD-Folie real xFE(3),vFE(3) ! Ort und Geschw. der FE real tFE ! Zeit real tFE_min ! kuerzeste gueltige FE-Flugzeit je Projektil integer tFE_(5) /-1,-1,-1,-1,-1/ ! Flugzeit jedes FE in ps (fuer NTP) c c---------------- c c-TP-10/2000 add variables to have position information of muons at c TD and FE at MCP3 in NTuple; up to 5 electrons possible c real xFE_MCP(5), yFE_MCP(5), zFE_MCP(5) common /TrigDet/ x0FE, xFE_MCP, yFE_MCP, zFE_MCP c c---------------- c COMMON /S1xM2/ S1xM2 ! fuer NTupel COMMON /TIMES/ S1M2,S1Fo,FoM2,S1M3,M3M2,tFE_ ! fuer NTupel common /FoM2Only/ FoM2Only COMMON /S1FoOnly/ S1FoOnly c Variablen fuer den allgemeinen Programmablauf: integer qIndxMu common /qIndxMu/ qIndxMu integer ntpid(1) ! fuer das Einlesen des NTupels von ACCEL oder von integer ntpnr ! FoilFile integer firstEventNr external firstEventNr logical NTPalreadyWritten real Spiegel_Faktor ! Faktor bei Berechn. der Reflektionszeit im Spiegel integer bis_Spiegel ! verschiedene Label integer bis_L3_Mappe, bis_MCP2_Mappe, MCP2_Mappe character uhrzeit*8 integer percent_done logical fill_NTP real radiusQuad_HeShield real radiusQuad_LNShield real radiusQuad_L1 real radiusQuad_L2 real radiusQuad_L3 real radiusQuad_Blende real radiusQuad_Rohr real radiusQuad_MCP2 ! Radiusquadrat des MCP2 real radiusQuad_MCP2active ! Radiusquadrat der aktiven Flaeche des MCP2 real radiusQuad_Sp ! Radiusquadrat der Spiegeldraehte real rWires_Sp ! Radius der Spiegeldraehte logical check_Blende /.false./ real xChangeKoord ! legt den Ort nach dem Spiegel fest, bei parameter (xChangeKoord = 75.) ! dem das Koordinatensystem gewechselt wird integer n_return ! die Returnvariable fuer Aufruf von 'TD_CALC' integer zaehler ! Zaehler fuer Monitoring der Trajektorie in den ! Gebieten, in denen stepwise integriert werden ! muss logical flag, flag_ok integer okStepsCounter integer i, k ! integer-Hilfsvariablen real help1, help2 ! real-Hilfsvariablen real help3, help4 ! real-Hilfsvariablen real YieldPlus,YieldNeutral ! Ladungsanteile nach TD-Foliendurchgang integer startLabel ! das Einsprunglabel beim Teilchenstart character helpChar*7, ant*1 character HistogramTitle*32 /'Schnitt bei x = (i. Teil)'/ d real dtmin_L1, dtmin_Sp, dtmin_L2andFo, dtmin_FO, dtmin_L3, dtmin_M2 d real dtmax_L1, dtmax_Sp, dtmax_L2andFo, dtmax_FO, dtmax_L3, dtmax_M2 d real x_dtmin_L1(3), x_dtmax_L1(3), x_dtmin_FO(3), x_dtmax_FO(3) d real x_dtmin_L2andFo(3), x_dtmax_L2andFo(3) d real x_dtmin_L3(3), x_dtmax_L3(3), x_dtmin_M2(3), x_dtmax_M2(3) d real x_dtmin_Sp(3), x_dtmax_Sp(3) d d ! /ntp_steps/ enthaelt auch 'steps' (ueber COM-MUTRACK.INC) d COMMON /ntp_steps/ dtmin_L1, x_dtmin_L1, dtmax_L1, x_dtmax_L1, d + dtmin_Sp, x_dtmin_Sp, dtmax_Sp, x_dtmax_Sp, d + dtmin_L2andFo, x_dtmin_L2andFo, dtmax_L2andFo, x_dtmax_L2andFo, d + dtmin_FO, x_dtmin_FO, dtmax_FO, x_dtmax_FO, d + dtmin_L3, x_dtmin_L3, dtmax_L3, x_dtmax_L3, d + dtmin_M2, x_dtmin_M2, dtmax_M2, x_dtmax_M2 real x40(2:3),v40(3),t40,E40 ! Speicher fuer Trajektoriengroessen bei x=40mm COMMON /NTP_40mm/ x40,v40,t40,E40 cMBc logical writeTraj2File cMBc common /writeTraj2File/ writeTraj2File c Variablen fuer Test ob Draht getroffen wurde: real distToWire(2) integer DrahtNr logical WireHit real WireRadiusQuad_G1,WireRadiusQuad_G2 real WireRadiusQuad_Sp c Variablen fuer die Graphikausgabe: real xKoord(1000),xKoord_(1000) ! Koordinatenfelder fuer die real yKoord(1000),yKoord_(1000) ! Graphikausgabe real zKoord(1000),zKoord_(1000) ! cMBc real tKoord(1000),tKoord_(1000) ! integer nKoord,nKoordSave ! Anzahl der Koordinaten cMBc COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord,tKoord COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord 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) c------------------------------------------------------------------------------- c Konstanten und Variable fuer Berechnung der Winkelaufstreuung in Triggerfolie c mittels Meyer-Formel (L.Meyer, phys.stat.sol. (b) 44, 253 (1971)): real g1, g2 ! Tabellierte Funktionen der Referenz real effRedThick ! effektive reduzierte Dicke ('tau' der Referenz) c - Parameter: real Z1, Z2 ! die atomaren Nummern von Projektil und Target real a0 ! Bohrscher Radius in cm real screeningPar ! Screeningparameter 'a' in cm fuer Teilchen der ! Kernladungszahl Z1=1 in Kohlenstoff (Z2 = 6) ! bei Streichung von Z1 (vgl. Referenz, S. 268) real r0Meyer ! r0(C) berechnet aus dem screeningParameter 'a' ! und dem ebenfalls bei Meyer angegebenem ! Verhaeltnis a/r0=0.26 (vgl. Referenz, S. 263 oben) real eSquare ! elektrische Ladung zum Quadrat in keV*cm real HWHM2sigma ! Umrechnungsfaktor von (halber!) Halbwertsbreite ! nach Sigma der Gaussfunktion real Na ! die Avogadrokonstante real mMolC ! molare Masse von C in ug real Pi ! die Kreiszahl parameter (Z1 = 1, Z2 = 6, a0 = 5.29E-9, ScreeningPar = 2.5764E-9) parameter (r0Meyer = 9.909E-9, eSquare = 1.44E-10, HWHM2sigma = 1./1.17741) parameter (Na = 6.022e23, mMolC = 12.011e6, Pi = 3.141592654) c - Bei der Berechnung von Sigma auftretende Vorfaktoren. c (Meyer_faktor 1 wird benoetigt fuer Berechnung der reduzierten Dicke aus der c 'ug/cm2'-Angabe der Foliendicke. Meyer_faktor2 und Meyer_faktor3 werden c direkt fuer die Berechnung von sigma aus den beiden tabellierten Funktionen c g1 und g2 verwendet): real Meyer_Faktor1, Meyer_Faktor2, Meyer_Faktor3 parameter (Meyer_faktor1 = Pi*screeningPar*screeningPar * Na/mMolC) ! Na/mMolC = 1/m(C-Atom) parameter (Meyer_faktor2 = (2*Z1*Z2 * eSquare)/ScreeningPar * 180./Pi + * HWHM2sigma) parameter (Meyer_faktor3 = (screeningPar/r0Meyer) * (screeningPar/r0Meyer)) c------------------------------------------------------------------------------- c Kommentar zur Berechnung der Winkelaufstreuung nach Meyer: c c Als Bedingung fuer die Gueltigkeit der Rechnung wird verlangt, dass c c (1) die Anzahl n der Stoesse >> 20*(a/r0)^(4/3) sein muss. Fuer Protonen auf c Graphit ist laut Referenz a/r0 gleich 0.26 (mit Dichte von 3.5 g/ccm habe c ich einen Wert von 0.29 abgeschaetzt). Fuer Myonen hat man den selben c Wert zu nehmen. Damit ergibt sich die Forderung, dass n >> 3.5 sein muss. c c (2) unabhaengig von (1) n >> 5 sein muss, was (1) also mit einschliesst. c c Mit n = Pi*r0*r0*Teilchen/Flaeche ergibt sich fuer eine Foliendicke von c 3 ug/cm^2 als Abschaetzung fuer n ein Wert von 37. (r0 ueber r0 = 0.5 N^(1/3) c und 3.5 g/ccm zu 8.9e-9 cm abgeschaetzt). D.h., dass die Bedingungen in c unserem Fall gut erfuellt sind. c In dem Paper wird eine Formel fuer Halbwertsbreiten angegeben. Ich habe nicht c kontrolliert, in wie weit die Form der Verteilung tatsaechlich einer Gauss- c verteilung entspricht. Zumindest im Bereich der Vorwaertsstreuung sollte c die in diesem Programm verwendete Gaussverteilung aber eine sehr gute c Naeherung abgeben. Abweichungen bei groesseren Winkeln koennten jedoch u. U. c die absolute Streuintensitaet in Vorwaertsrichtung verfaelschen. czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c HIER GEHT DER PROGRAMMTEXT RICHTIG LOS czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c Initialisierungen: INCLUDE 'initialize.inc' c Einlesen der Parameter aus 'MUTRACK.INPUT' und Setzen der entsprechenden c Voreinstellungen. Einlesen der Kammergeometrie sowie der INFO-files der c Feldmappen: call read_inputFile c Berechnen der RadiusQuadrate: radiusQuad_HeShield = rHeShield*rHeShield radiusQuad_LNShield = rLNShield*rLNShield radiusQuad_Rohr = radius_Rohr*radius_Rohr radiusQuad_L1 = iRadiusCyl_L1*iRadiusCyl_L1 radiusQuad_L2 = iRadiusCyl_L2*iRadiusCyl_L2 radiusQuad_L3 = iRadiusCyl_L3*iRadiusCyl_L3 radiusQuad_Blende = radius_Blende*radius_Blende radiusQuad_MCP2 = radius_MCP2*radius_MCP2 radiusQuad_MCP2active = radius_MCP2active*radius_MCP2active WireRadiusQuad_G1 = dWires_G1/2. * dWires_G1/2. WireRadiusQuad_G2 = dWires_G2/2. * dWires_G2/2. WireRadiusQuad_Sp = dWires_Sp/2. * dWires_Sp/2. rWires_Sp = dWires_Sp/2. radiusQuad_Sp = rWires_Sp * rWires_Sp c Einlesen der Feldmappen: write(*,*)'----------------------------------------'// + '----------------------------------------' if (.NOT.(par(1,UL1).EQ.0. .AND. n_par(UL1).LE.1)) call READ_MAP_L1 if (.NOT.(idealMirror .OR. (par(1,USp).EQ.0. .AND. n_par(USp).LE.1))) then call read_Map_SP_1 call read_Map_SP_2 call read_Map_SP_3 endif if (TriggerInBeam .AND. .NOT.lense2 .AND. ! 'lense2' muss noch in sub_input richtig gesetzt werden! (-> foilfile) + .NOT.(par(1,UFolie).EQ.0. .AND. n_par(UFolie).LE.1) ) then call READ_MAP_FO endif if (.NOT.(par(1,UL3).EQ.0. .AND. n_par(UL3).LE.1)) then if (.NOT.(par(1,UMCP2).EQ.0. .AND. n_par(UMCP2).LE.1)) then if (xLeaveMap_L3.GT.xEnterMap_M2) then write(*,*) write(*,*)' Potentialmappen von Linse 3 und MCP2 ueberlappen!' write(*,*)' Dies ist in der aktuellen Implementierung des Programmes' write(*,*)' nicht vorgesehen!' write(*,*) write(*,*)' -> STOP' write(*,*) STOP endif endif call READ_MAP_L3 endif if (.NOT.(par(1,UMCP2).EQ.0. .AND. n_par(UMCP2).LE.1)) call READ_MAP_M2 c Eingelesene Simulationsparameter auf Schirm geben und bestaetigen lassen. c Die Ausgabefiles initialisieren: call initialize_output c falls ein 'FoilFile' erstellt werden soll, schreibe das .INFO-files: if (createFoilFile) call make_INFOFile if (Use_MUTRACK) Use_ACCEL = .false. c Defaultwert fuer 'fill_NTP' setzen (wird weiter unten ueberschrieben, falls c fuer das Fuellen des NTupels spezielle Triggerbedingung verlangt ist): if (createNTP) then fill_NTP = .true. else fill_NTP = .false. endif c CERN-Pakete initialisieren (Groesse des COMMONblocks PAWC uebermitteln): if (.NOT.fromScratch.OR.Graphics.OR.createNTP.OR.createFoilFile) call HLIMIT(HB_memsize) c Graphikausgabe initialisieren: if (GRAPHICS) then call masstab_setzen CALL HPLSET ('VSIZ',.6) ! AXIS VALUES SIZE write(HistogramTitle(17:22),'(F6.1)') schnitt_x write(HistogramTitle(25:25),'(I1)') schnitt_p CALL HPLSET ('TSIZ',.7) ! HISTOGRAM TITLE SIZE CALL HBOOK2 (50,HistogramTitle,100,-30.,30.,100,-30.,30.,20.) endif c falls fruehere Simulation fortgefuehrt werden soll, oeffne entsprechende Datei: if (.NOT.fromScratch) then if (use_ACCEL) then call HROPEN(lunREAD,'ACCEL',ACCEL_Dir//':'//fileName_ACCEL//'.NTP', + ' ',1024,istat) else call HROPEN(lunREAD,'MUread',outDir//':'//fileName_MUTRACK//'.NTP', + ' ',1024,istat) endif call HRIN(0,99999,0) call HIDALL(ntpid,ntpNr) call HDELET(ntpid(1)) i = NTP_read - ntpid(1) call HRIN(NTP_read-i,9999,i) ! NTP_read = NTP_write+1 call HBNAME(NTP_read,' ',0,'$CLEAR') ! alles resetten ! fuer die benoetigten Bloecke des CWN die entsprechenden Speicher- ! lokalisationen uebermitteln: if (random_E0) call HBNAME(NTP_read,'E0',E0,'$SET') if (random_pos) call HBNAME(NTP_read,'x0',x0,'$SET') if (random_angle) call HBNAME(NTP_read,'angle0',theta0,'$SET') ! theta0,phi0 if (UseDecay_prevSim) call HBNAME(NTP_read,'lifetime',lifetime,'$SET') if (smearS1Fo .AND. use_MUTRACK) then call HBNAME(NTP_read,'S1FoS',S1FoOnly,'$SET') endif call HBNAME(NTP_read,'dest',gebiet,'$SET') ! gebiet,destiny call HBNAME(NTP_read,'Traj',t,'$SET') ! t,x,v endif c NTP-relevante Befehle: c BAD LUCK!!! Das Packen der Real-Variablen im folgenden hat KEINERLEI VER- c KLEINERUNG DER FILEGROESSE bewirkt!!!! (fuer die Integers habe ich noch c keinen Test gemacht). -> wohl besser wieder herausnehmen. Ich verliere c u.U. nur Genauigkeit und habe nur einen eingeschraenkten Wertebereich zur c Verfuegung! if (createNtp.OR.createFoilFile) then !c Datei fuer NTupelausgabe oeffnen: call HROPEN(lunNTP,'MUwrite',outDir//':'//filename//'.NTP', + 'N',1024,istat) if (istat.NE.0) then write(*,*) write(*,*)'error ',istat,' opening HBOOK-file' write(*,*) STOP endif call HBNT(NTP_write,filename,'D') ! D: Disk resident CWN buchen !c die Bloecke des CWN definieren: if (.NOT.OneLoop) call HBNAME(NTP_write,'LOOP',schleifenNr,'loop[1,1000]:u') if (M2_triggered .OR. Fo_triggered.AND.upToTDFoilOnly) then ! -> Gebiet und Destiny stehen hier sowieso fest, nimm ! diese Groessen daher erst gar nicht mehr in das NTupel auf! else call HBNAME(NTP_write,'DEST',gebiet,'Gebiet[0,20]:u,dest[-10,10]:i') endif if (NTP_Start .OR. createFoilFile.AND.random_pos) then call HBNAME(NTP_write,'X0',x0,'x0,y0,z0') endif if (NTP_Start) call HBNAME(NTP_write,'V0',v0,'vx0,vy0,vz0') if (NTP_Start .OR. createFoilFile.AND.random_E0) then call HBNAME(NTP_write,'E0',E0,'E0') endif if (NTP_Start .OR. createFoilFile.AND.random_angle) then call HBNAME(NTP_write,'ANGLE0',theta0,'theta0,phi0') endif if (NTP_lifetime .OR. createFoilFile.AND.UseDecay) then call HBNAME(NTP_write,'LIFETIME',lifetime,'lifetime:r') endif if (NTP_40mm) call HBNAME(NTP_write,'X=40MM',x40, + 'y40,z40,vx40,vy40,vz40,t40,E40') if (NTP_S1xM2) call HBNAME(NTP_write,'S1xM2',S1xM2,'S1xM2') if (NTP_Times) then if (TriggerInBeam) then if (generate_FE) then call HBNAME(NTP_write,'TIMES',S1M2, + 'S1M2,S1Fo,FoM2,S1M3,M3M2:r,TFE(5):i') else call HBNAME(NTP_write,'TIMES',S1M2, + 'S1M2,S1Fo,FoM2') endif else call HBNAME(NTP_write,'TIMES',S1M2, + 'S1M2') endif endif if (NTP_FoM2Only) then call HBNAME(NTP_write,'FoM2',FoM2Only,'FoM2') endif if (NTP_Folie) then call HBNAME(NTP_write,'FOLIE',E_Verlust, + 'ELoss,thetStreu,phiStreu') c c-------------------------- c c-TP-10/2000 add position at foil and MCP3 (FE) c call HBNAME(NTP_write, 'TrigDet', x0FE, + 'x0FE,y0FE,z0FE,xFE(5),yFE(5),zFE(5)') c c-------------------------- c endif if (NTP_charge) call HBNAME(NTP_write,'CHARGE',qInt,'q[-5,5]:i') if (NTP_stop.OR.createFoilFile) then call HBNAME(NTP_write,'TRAJ',t,'t,x,y,z,vx,vy,vz') endif c if (createFoilFile .AND. smearS1Fo .AND. .NOT.NTP_times) then if (smearS1Fo) then call HBNAME(NTP_write,'S1FoS',S1FoOnly,'S1FoS') endif if (NTP_stop) then call HBNAME(NTP_write,'EKIN',Ekin,'Ekin') endif d if (NTP_steps) then d call HBNAME(NTP_write,'STEP',steps,'steps[1,100000]:u,'// d + 'dtminL1, xdtminL1, ydtminL1, zdtminL1,'// d + 'dtmaxL1, xdtmaxL1, ydtmaxL1, zdtmaxL1,'// d + 'dtminL2, xdtminL2, ydtminL2, zdtminL2,'// d + 'dtmaxL2, xdtmaxL2, ydtmaxL2, zdtmaxL2,'// d + 'dtminFO, xdtminFO, ydtminFO, zdtminFO,'// d + 'dtmaxFO, xdtmaxFO, ydtmaxFO, zdtmaxFO,'// d + 'dtminL3, xdtminL3, ydtminL3, zdtminL3,'// d + 'dtmaxL3, xdtmaxL3, ydtmaxL3, zdtmaxL3,'// d + 'dtminM2, xdtminM2, ydtminM2, zdtminM2,'// d + 'dtmaxM2, xdtmaxM2, ydtmaxM2, zdtmaxM2') d endif endif c die Label definieren: assign 7 to bis_Spiegel assign 14 to bis_L3_Mappe assign 16 to bis_MCP2_Mappe assign 17 to MCP2_Mappe c die Einsprungposition fuer den Beginn der Trajektorienberechnungen setzen: if (Use_MUTRACK) then assign 113 to startLabel elseif (Use_ACCEL) then assign 3 to startLabel elseif (Gebiet0.EQ.target .OR. Gebiet0.EQ.upToGrid1) then assign 1 to startLabel elseif (Gebiet0.EQ.upToGrid2) then assign 2 to startLabel elseif (Gebiet0.EQ.upToHeShield) then assign 3 to startLabel elseif (Gebiet0.EQ.upToLNShield) then assign 4 to startLabel elseif (Gebiet0.EQ.upToL1Map) then assign 5 to startLabel elseif (Gebiet0.EQ.upToExL1) then assign 6 to startLabel elseif (Gebiet0.EQ.upToEnSp) then assign 7 to startLabel elseif (Gebiet0.EQ.upToExSp) then assign 8 to startLabel elseif (Gebiet0.EQ.upToChKoord) then assign 9 to startLabel elseif (Gebiet0.EQ.upToEnTD) then assign 10 to startLabel elseif (Gebiet0.EQ.upToExTD) then if (log_alpha0_KS) then assign 111 to startLabel else assign 112 to startLabel endif elseif (Gebiet0.EQ.upToL2andFoMap) then c assign 12 to startLabel elseif (Gebiet0.EQ.upToExL2) then c assign 13 to startLabel elseif (Gebiet0.EQ.upToL3Map) then assign 12 to startLabel elseif (Gebiet0.EQ.upToExL3) then assign 13 to startLabel elseif (Gebiet0.EQ.upToM2Map) then assign 14 to startLabel elseif (Gebiet0.EQ.upToMCP2) then assign 15 to startLabel endif c Abkuerzungen 'Length1' und 'length2' setzen: length1 = d_Folie_Achse + MappenLaenge_FO length2 = xTD - d_Folie_Achse - MappenLaenge_FO c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c ab hier beginnen die Schleifen: c (Bemerkung: eine Laufvariable darf kein Feldelement 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 Bei Aenderungen in der Abfolge der Schleifen muss die Anweisungszeile c 'DATA reihenfolge /.../' in 'INITIALIZE.INC' entsprechend editiert werden! c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c zusaetliche Flugstrecken vor TD und MCP2 (gehen NUR in t, NICHT in x ein!!): c ---------------------------------------------------------------------------- do 200 Delta_L1 = par(1,DeltaL1),par(2,DeltaL1),par(3,DeltaL1) parWert(DeltaL1) = Delta_L1 do 200 Delta_L2 = par(1,DeltaL2),par(2,DeltaL2),par(3,DeltaL2) parWert(DeltaL2) = Delta_L2 c Foliendicke und Energieverlust: c ------------------------------- do 200 E_loss = par(1,Eloss),par(2,Eloss),par(3,Eloss) ! Eloss parWert(Eloss) = E_loss mean_E_Verlust = E_loss do 200 Thickness = par(1,Thickn),par(2,Thickn),par(3,Thickn)! Thickness parWert(Thickn) = Thickness c MCP2: c ----- do 200 U_MCP2 = par(1,UMCP2),par(2,UMCP2),par(3,UMCP2) ! U(MCP2) parWert(UMCP2) = U_MCP2 c Winkel: c ------- do 200 alfaTgt = par(1,alfTgt),par(2,alfTgt),par(3,alfTgt) ! ALPHA(TARGET) parWert(alfTgt) = alfaTgt Sin_alfaTgt= sind(alfaTgt) Cos_alfaTgt= cosd(alfaTgt) do 200 alfaSp = par(1,alfSp),par(2,alfSp),par(3,alfSp) ! ALPHA(SPIEGEL) parWert(alfSp) = alfaSp Sin_alfaSp = sind(alfaSp) Cos_alfaSp = cosd(alfaSp) Tan_alfaSp = tand(alfaSp) help1 = dSpiegel/2.+DreharmLaenge ! Berechne die y-Werte der 'oberen linken' (yUppLeft) und der ! 'unteren linken' (yLowLeft) Spiegelecke: if (idealMirror) then yUppLeft = + bSpiegel/2. * Sin_alfaSp + + help1 * Cos_alfaSp yLowLeft = - bSpiegel/2. * Sin_alfaSp + + help1 * Cos_alfaSp endif ! Berechne Schnittpunkt y_intersectSp der vorderen Spiegelebene bzw. ! der vorderen Mappenkante mit der Geraden x = xSpiegel: if (.NOT.idealMirror) help1 = help1 + xSpGrid1 y_intersectSp = help1/Cos_alfaSp do 200 alfaTD = par(1,alfTD),par(2,alfTD),par(3,alfTD) ! ALPHA(TRIGGERDETEKTOR) parWert(alfTD) = alfaTD Sin_alfaTD = sind(alfaTD) Cos_alfaTD = cosd(alfaTD) Tan_alfaTD = tand(alfaTD) ! Berechne Schnittpunkt 'x_intersectTD' der x-Achse mit der Folien- ! ebene bzw im Fall von 'GridInFrontOfFoil' mit dem Gitter vor der ! Triggerfolie: help1 = d_Folie_Achse if (gridInFrontOfFoil) help1 = help1 + d_Grid_Folie x_intersectTD = xTD - help1/Cos_alfaTD help1 = d_Folie_Achse + mappenLaenge_Fo x_intersectTDMap = xTD - help1/Cos_alfaTD c TriggerDetektor: c ---------------- do 200 U_V = par(1,UVorne),par(2,UVorne),par(3,UVorne) ! U(VORNE) parWert(UVorne) = U_V do 200 U_H = par(1,UHinten),par(2,UHinten),par(3,UHinten) ! U(HINTEN) parWert(UHinten) = U_H do 200 U_MCP3 = par(1,UMCP3),par(2,UMCP3),par(3,UMCP3) ! U(MCP3) parWert(UMCP3) = U_MCP3 do 200 U_F = par(1,UFolie),par(2,UFolie),par(3,UFolie) ! U(FOLIE) parWert(UFolie) = U_F c Transportsystem: c ---------------- do 200 U_L2 = par(1,UL2),par(2,UL2),par(3,UL2) ! U(Linse 2) parWert(UL2) = U_L2 c gegebenenfalls die Mappe 'L2andFo' zusammenbauen: if (lense2) then if ( .NOT.(par(1,UL2).EQ.0. .AND. n_par(UL2).LE.1) .OR. + .NOT.(par(1,UFolie).EQ.0. .AND. n_par(UFolie).LE.1) ) then ! Addiere die Mappen nur erneut, falls die jetztige Konfiguration ! nicht mit der letzten uebereinstimmt: if (U_L2.NE.last_U_L2 .OR. U_F.NE.last_U_F) then call ADD_MAP_L2andFo last_U_L2 = U_L2 last_U_F = U_F endif endif endif do 200 U_Sp = par(1,USp),par(2,USp),par(3,USp) ! U(SPIEGEL) parWert(USp) = U_Sp do 200 U_L1 = par(1,UL1),par(2,UL1),par(3,UL1) ! U(Linse 1) parWert(UL1) = U_L1 do 200 U_L3 = par(1,UL3),par(2,UL3),par(3,UL3) ! U(Linse 3) parWert(UL3) = U_L3 c die Magnetfelder: c ----------------- do 200 B_Helm = par(1,BHelm),par(2,BHelm),par(3,BHelm) ! Helmholtzsp. parWert(BHelm) = B_Helm do 200 B_TD = par(1,BTD),par(2,BTD),par(3,BTD) ! TD-Spule parWert(BTD) = B_TD 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 qIndxMu = q_ ! fuer Verwendung in function firstEventNr! ArtNr = Art_Nr(q_) m = Art_Masse(ArtNr) q = Art_Ladung(ArtNr) parWert(mass) = m parWert(charge) = q endif ! gegebenenfalls ein Flag fuer die Beruecksichtigung des Myonen- ! zerfalles setzen: 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 Beschleuniger: c -------------- do 200 U_Tgt = par(1,UTgt),par(2,UTgt),par(3,UTgt) ! U(TARGET) parWert(UTgt) = U_Tgt do 200 U_Gua = par(1,UGua),par(2,UGua),par(3,UGua) ! U(GUARD) parWert(UGua) = U_Gua do 200 U_G1 = par(1,UG1),par(2,UG1),par(3,UG1) ! U(GITTER) parWert(UG1) = U_G1 parIndx(5) = parIndx(5) + 1 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c haeufig benoetigte Faktoren, die von der aktuellen Masse, Ladung und Hoch- c spannungen abhaengen: c (bei Linse 2 wird die Spannung direkt auf die Potentialmappe aufmultipliziert. c Daher wird dort 'Beschl_Faktor' verwendet und kein 'Beschl_Faktor_L2' benoetigt) Energie_Faktor = m / (2.*c*c) Beschl_Faktor = q / m * c*c Beschl_Faktor_L1 = Beschl_Faktor * U_L1 Beschl_Faktor_Sp = Beschl_Faktor * U_Sp Beschl_Faktor_FO = Beschl_Faktor * U_F Beschl_Faktor_L3 = Beschl_Faktor * U_L3 Beschl_Faktor_M2 = Beschl_Faktor * U_MCP2 aFoil = - Beschl_Faktor * U_F / d_Grid_Folie if (U_Sp.EQ.0. .OR. q.EQ.0.) then Spiegel_Faktor = 0 else Spiegel_Faktor = 2.*dspiegel / (Beschl_Faktor * U_Sp) !<-- pruefen! endif ! Die Beschleunigungen in den beiden (idealen) Beschleunigerstufen: a1 = Beschl_Faktor * (U_Tgt - U_G1) / (XGrid1 - XTarget) a2 = Beschl_Faktor * U_G1 / (XGrid2 - xGrid1) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Falls 'fromScratch': c Die in den ab hier beginnenden Startparameter-Schleifen gesetzten Werte c werden gegebenenfalls weiter unten durch zufallsverteilte Offsets modi- c fiziert. (-> 'Zufallschleife': 'do 100 randomloop_ = 1, n_par(0)) c Andernfalls: c Wurden waehrend ACCEL oder 'foilfile' fuer die Startparameter Zufalls- c verteilungen verwendet, so werden die entsprechenden Groessen aus dem c betreffenden NTupel eingelesen. c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Startparameter: c --------------- do 200 E0_ = par(1,ener),par(2,ener),par(3,ener) ! E0 if (.NOT.random_E0) then E0 = E0_ v0_Betrag = sqrt(E0/Energie_Faktor) endif if (E0InterFromFile) then lowerE0 = E0Low(nInt(E0_)) upperE0 = E0Low(nint(E0_+1)) endif c falls Energieverlustberechnung aus ICRU-Tabelle verlangt ist und mittlerer c Energieverlust nicht fuer jedes Teilchen extra berechnet werden soll (sinnvoll c wenn alle Teilchen gleiche Startenergie haben oder Streuung der Startenergien c klein ist, so dass die Streuung des mittleren Energieverlustes vernachlaessigt c werden kann): if (log_E_Verlust_ICRU .AND. .NOT.calculate_each) then if (random_E0_equal) then Ekin = E0_ + (upperE0+lowerE0)/2. else Ekin = E0_ endif if (Gebiet0.EQ.target .OR. Gebiet0.EQ.upToGrid1) then Ekin = Ekin + q*(U_Tgt - U_F) elseif (Gebiet0.EQ.upToGrid2) then Ekin = Ekin + q*(U_G1 - U_F) endif call CALC_ELOSS_ICRU(Ekin,q,m,Thickness,mean_E_Verlust) endif if (log_Meyer_F_Function) then if (random_E0_equal) then Ekin = E0_ + (upperE0+lowerE0)/2. else Ekin = E0_ endif if (Gebiet0.EQ.target .OR. Gebiet0.EQ.upToGrid1) then Ekin = Ekin + q*(U_Tgt - U_F) elseif (Gebiet0.EQ.upToGrid2) then Ekin = Ekin + q*(U_G1 - U_F) endif effRedThick = Meyer_Faktor1 * Thickness call Get_F_Function_Meyer(effRedThick,Ekin) 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) then x0(2) = y0_ endif do 200 z0_ = par(1,zPos),par(2,zPos),par(3,zPos) ! z0 if (.NOT.random_pos) then x0(3) = z0_ endif 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 u.s.w. erst einmal mit den aktuellen Werten der c entsprechenden Schleifen gefuellt: parWert(ener) = E0_ parWert(thetAng) = theta0_ parWert(phiAng) = phi0_ parWert(yPos) = y0_ parWert(zPos) = z0_ c falls fruehere Simulation fortgefuehrt wird: c Berechne diejenige Eventnummer in NTP_read, ab welcher die relevanten c Simulationsparameter von ACCEL bzw. des 'FoilFiles' mit den gegenwaertigen c MUTRACK-(Schleifen)-Parametern uebereinstimmen: if (.NOT.fromScratch) eventNr = firstEventNr() c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Hier folgen die Befehle, die zu Beginn jeder neuen Schleife faellig sind: c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SchleifenNr = SchleifenNr + 1 ! Schleifen zaehlen okStepsCounter = 0 ! 'okStepsCounter' dient der Bestimmung ! der mittleren Anzahl von Integrations- ! schritten bis zum Ziel nNeutral = 0 ! noch wurden keine Teilchen in der TD-Folie nCharged = 0 ! neutralisiert 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 ('CHAM_1','M') ! ERZEUGEN VON BILDERN IM PAWC-COMM-BLOCK CALL IZPICT ('CHAM_2','M') CALL IZPICT ('HISTO','M') CALL IZPICT ('TEXT','M') call plot_chamber(schnitt_p) 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(q_) ! (q_) wegen der neutral fractions endif c Ausgabe der Prozentzahl schon gerechneten Trajektorien vorbereiten: if (log_percent) then call time(uhrzeit) percent_done = 0 write(*,1001)Uhrzeit,' %: 0' endif 1001 format ($,6x,A,A) c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c bei 'fromScratch': c Hier wird gegebenenfalls bei Zufallsverteilung von Startparametern ein ent- c sprechend gewuerfelter Offset auf den aktuellen Wert aufgeschlagen. c Ansonsten: c Lies bei Zufallsverteilungen die entsprechenden Startwerte aus dem NTupel c ein. c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - do 100 randomloop_ = 1, n_par(0) if (.NOT.fromScratch) then eventNr = eventNr + 1 ! Eventnummer im NTP if (smearS1Fo.AND.use_MUTRACK) call HGNTB(NTP_read,'S1FoS',eventNr,istat) if (istat.NE.0) then write(*,*) write(*,*)' error executing ''call HGNTB(',NTP_read,',''S1FoS'',eventNr,istat)''' write(*,*)' eventNr = ',eventNr write(*,*)' -> STOP' write(*,*) call exit endif ! Einlesen von 'Gebiet' und 'destiny': call HGNTB(NTP_read,'dest',eventNr,istat) ! gegebenenfallsvon freuher verwendete Gebietskodierung ! (ohne Linse2) aktualisieren if (gebiet.GE.10 .AND. MutrackVersionIndx.LE.1) gebiet = gebiet+2 c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(',NTP_read,',''dest'',eventNr,istat)''' c write(*,*)' eventNr = ',eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif ! Einlesen der Trajektoriendaten 't,x(3),v(3)': call HGNTB(NTP_read,'Traj',eventNr,istat) c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(',NTP_read,',''Traj'',eventNr,istat)''' c write(*,*)' eventNr = ',eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif if (Use_Accel) then ! Uebersetzen der von ACCEL verwendeten code_Nummern fuer die ! moeglichen Teilchenschicksale in von MUTRACK verwendete ! code_Nummern: if (destiny.EQ.-5) then destiny = code_frontOfMapAc elseif (destiny.EQ.-4) then destiny = code_leftMapAc elseif (destiny.EQ.-3) then gebiet = upToGrid2 destiny = code_grid elseif (destiny.EQ.-2) then gebiet = upToGrid1 destiny = code_grid elseif (destiny.EQ.-1) then destiny = code_hit_TgtHolder elseif (destiny.EQ.code_ok) then Gebiet = upToHeShield elseif (destiny.EQ.+1) then destiny = code_decay elseif (destiny.EQ.+2) then destiny = code_reflektiert elseif (destiny.EQ.+3) then destiny = code_wand elseif (destiny.EQ.+4) then destiny = code_lost elseif (destiny.EQ.+5) then destiny = code_dtsmall else write(*,*)'UNKNOWN ACCEL-CODE-NR: destiny = ',destiny call exit endif ! Auf xGrid2 zurueckrechnen, damit unabhaengiger Test auf ! Treffer des He-Fensters gemacht werden kann (nur, falls ! Teilchen nicht schon anderweitig gestorben ist). Auch ! notwendig fuer Graphikausgabe. if (destiny.EQ.0) then dt = (xGrid2-x(1))/v(1) ! < 0. t = t + dt x(1) = xGrid2 x(2) = x(2)+v(2)*dt x(3) = x(3)+v(3)*dt endif ! falls Kryo verdreht ist, rechne in Kammerkoordinaten um: if (alfaTgt.NE.0.) then if (alfaTgtVertically) then help1 = x(1) x(1) = help1 * Cos_alfaTgt - x(3) * Sin_alfaTgt x(3) = help1 * Sin_alfaTgt + x(3) * Cos_alfaTgt help1 = v(1) v(1) = help1 * Cos_alfaTgt - v(3) * Sin_alfaTgt v(3) = help1 * Sin_alfaTgt + v(3) * Cos_alfaTgt else help1 = x(1) x(1) = help1 * Cos_alfaTgt - x(2) * Sin_alfaTgt x(2) = help1 * Sin_alfaTgt + x(2) * Cos_alfaTgt help1 = v(1) v(1) = help1 * Cos_alfaTgt - v(2) * Sin_alfaTgt v(2) = help1 * Sin_alfaTgt + v(2) * Cos_alfaTgt endif endif endif endif if (random_E0) then ! random_ENERGIE if (fromScratch) then if (random_E0_equal) then ! -> gleichverteilt 300 E0 = E0_ + lowerE0 + (upperE0 - lowerE0)*ran(seed) if (E0.LT.0) goto 300 elseif (random_E0_gauss) then ! -> gaussverteilt 310 call Gauss_Verteilung(sigmaE0,help1) E0 = E0_ + help1 if (E0.LT.0) goto 310 endif else ! Einlesen von 'E0' aus NTP: call HGNTB(NTP_read,'E0',eventNr,istat) c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(NTP_read,''E0'',eventNr,istat)''' c write(*,*)' eventNr = ',eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif endif parWert(ener) = E0 v0_Betrag = sqrt(E0/Energie_Faktor) endif if (random_pos) then ! random_POSITION if (fromScratch) then 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 320 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 320 x0(3) = r0 * sind(phi_r0) if (abs(x0(3)).GT.StartHoehe/2.) goto 320 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 330 r0 = abs(sigmaPosition*sqrt(-2.*log(1.-ran(seed)))) if (r0.GT.StartRadius) goto 330 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) else ! Einlesen von 'x0(3)' aus NTP: call HGNTB(NTP_read,'x0',eventNr,istat) c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(',NTP_read,',''x0'',eventNr,istat)''' c write(*,*)' eventNr = ',eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif endif parWert(yPos) = x0(2) parWert(zPos) = x0(3) endif if (random_angle) then ! random_WINKEL if (fromScratch) then 340 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 neuetheta0,phi0 als die tatsaechlichen Startwinkel c verwendet werden. Das alles geschieht vollkommen analog zur Winkelaufstreuung c in der Triggerfolie. c v wird als Hilfsvariable missbraucht. ! Berechnung 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 340 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 340 else ! Einlesen von 'theta0' und 'phi0' aus NTP: call HGNTB(NTP_read,'angle0',eventNr,istat) c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(',NTP_read,',''angle0'',eventNr,istat)''' c write(*,*)' eventNr = ',eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif endif parWert(thetAng) = theta0 parWert(phiAng) = phi0 endif ! Berechnung der Start-Geschwindigkeitskomponenten: v0(1) = v0_Betrag * Cos_theta0 v0(2) = v0_Betrag * Sin_theta0 * Cos_phi0 v0(3) = v0_Betrag * Sin_theta0 * Sin_phi0 if (fromScratch) then ! den Zeit-Speicher resetten: t = 0. ! Startparameter in Koordinatenspeicher uebergeben: do i = 1, 3 x(i) = x0(i) v(i) = v0(i) enddo endif c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Hier folgen die restl. Vorbereitungen zum Start des individuellen Projektils: c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c n_dtsmall resetten: n_dtsmall = 0 c Aufstreuwinkel resetten: thetaAufstreu = 0. phiAufstreu = 0. c x-Komponente der Startgeschwindigkeit ueberpruefen: if (v0(1).LT.0) then write(*,*) write(*,*) ' >>>> v(x) beim Start negativ!' write(*,*) call exit endif c die Lebensdauer wuerfeln: c (wird eine fruehere Simulation fortgefuehrt und wurde dort bereits der Myonen- c zerfall beruecksichtigt, so verwende die dort gewuerfelten Lebenszeiten) if (UseDecay_) then if (.NOT.UseDecay_prevSim) then 350 lifeTime = -meanlifeTime * Log(Ran(seed) + 1.0E-37) if (lifeTime.LE.0.) goto 350 elseif (.NOT.fromScratch) then call HGNTB(NTP_read,'lifetime',eventNr,istat) c if (istat.NE.0) then c write(*,*) c write(*,*)' error executing ''call HGNTB(',NTP_read,',''lifetime'',eventNr,istat)''' c write(*,*)' eventNr = ', eventNr c write(*,*)' -> STOP' c write(*,*) c call exit c endif endif endif c die Ladung resetten (falls in der Folie Neutralisierung stattgefunden hat): c ('qInt' wird fuer 'NTP_charge' benoetigt) q = parWert(charge) qInt = int(q) c Ausgabe der Prozentzahl schon gerechneter Trajektorien: if (log_percent) then if (100.*real(start_nr(1))/real(n_par(0)) + .GE.percent_done+5) then percent_done = percent_done + 5 write(*,1002) percent_done endif endif 1002 format ($,'+',I3) c andere Variablen auf den richtigen Stand bringen: if (fromScratch) then destiny = code_ok ! bis jetzt ist dem Teilchen noch nichts zugestossen Gebiet = Gebiet0 endif start_nr(1) = start_nr(1) + 1 ! Projektil-Startnummer erhoehen steps = 0 ! es wurden noch keine Integrationsschritte durchgefuehrt NTPalreadyWritten = .false. ! fuer 'createFoilFile' c die DEBUG-Daten ausgeben: if (Debug .AND. start_Nr(1).LE.DEBUG_Anzahl) then Debug_ = .true. call output_new_particle call Output_Debug else Debug_ = .false. endif c StartKoordinaten fuer Graphikausgabe sichern: if (graphics .AND. (start_Nr(1).LE.graphics_Anzahl .OR. OneStartPerLoop)) then graphics_ = .true. if (Use_ACCEL) then nKoord = 1 xKoord(1) = x0(1) yKoord(1) = x0(2) zKoord(1) = x0(3) else nKoord = 0 endif if (.NOT.(Use_MUTRACK.OR.Gebiet0.EQ.upToExTD)) call Save_Graphics_Koord else graphics_ = .false. endif c gegebenenfalls 'fill_NTP' resetten: if (Fo_triggered.OR.M2_triggered.OR.xM2_triggered) fill_NTP = .false. c Falls Schrittweiteninformationen im NTupel verlangt sind: Speicher resetten c und Startkoordinaten sichern: d if (NTP_steps) then d dtmin_L1 = +1.e10 d x_dtmin_L1(1) = 0 d x_dtmin_L1(2) = 0 d x_dtmin_L1(3) = 0 d dtmax_L1 = -1.e10 d x_dtmax_L1(1) = 0 d x_dtmax_L1(2) = 0 d x_dtmax_L1(3) = 0 d d dtmin_L2andFo = +1.e10 d x_dtmin_L2andFo(1) = 0 d x_dtmin_L2andFo(2) = 0 d x_dtmin_L2andFo(3) = 0 d dtmax_L2andFo = -1.e10 d x_dtmax_L2andFo(1) = 0 d x_dtmax_L2andFo(2) = 0 d x_dtmax_L2andFo(3) = 0 d d dtmin_FO = +1.e10 d x_dtmin_FO(1) = 0 d x_dtmin_FO(2) = 0 d x_dtmin_FO(3) = 0 d dtmax_FO = -1.e10 d x_dtmax_FO(1) = 0 d x_dtmax_FO(2) = 0 d x_dtmax_FO(3) = 0 d d dtmin_L3 = +1.e10 d x_dtmin_L3(1) = 0 d x_dtmin_L3(2) = 0 d x_dtmin_L3(3) = 0 d dtmax_L3 = -1.e10 d x_dtmax_L3(1) = 0 d x_dtmax_L3(2) = 0 d x_dtmax_L3(3) = 0 d d dtmin_M2 = +1.e10 d x_dtmin_M2(1) = 0 d x_dtmin_M2(2) = 0 d x_dtmin_M2(3) = 0 d dtmax_M2 = -1.e10 d x_dtmax_M2(1) = 0 d x_dtmax_M2(2) = 0 d x_dtmax_M2(3) = 0 d endif if (NTP_40mm) then x40(2) = 0. x40(3) = 0. v40(1) = 0. v40(2) = 0. v40(3) = 0. t40 = 0. E40 = 0. endif c Die Flugzeiten resetten: S1xM2 = 0. S1M2 = 0. S1Fo = 0. FoM2 = 0. S1M3 = 0. M3M2 = 0. c Falls das Teilchen schon nicht mehr existiert, gehe gleich zur Ausgabe: if (destiny.NE.code_ok) goto 555 ! (nur bei '.NOT.fromScratch' moeglich) czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c hier starten die Projektile: czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz goto startLabel ! StartLabel = Gebiet0 als label c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c erste Beschleunigerstufe: (homogenes Feld) 1 Gebiet = upToGrid1 steps = Steps + 1 if (a1.NE.0.) then help1 = v(1)*v(1) + 2.*a1*(xGrid1-x(1)) if (help1.LT.0) then ! Reflektion noch vor 1. Gitter dt = -2*v(1)/a1 t = t + dt !x(1) bleibt unveraendert x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt v(1) = -v(1) !v(2) bleibt unveraendert !v(3) bleibt unveraendert destiny = code_reflektiert goto 555 endif dt = (sqrt(help1) - v(1))/a1 ! (ergibt sich aus x=v*t+1/2*a*t**2 mit richtiger V.Z.-Wahl ('+')) v(1) = v(1) + a1*dt else if (v(1).EQ.0) then write(*,*) write(*,*)'ERROR: v(x) beim Start = 0. und '// + 'Beschleunigung = 0' write(*,*) STOP endif dt = (xGrid1-xTarget) / v(1) endif t = t + dt !v(2) bleibt unveraendert !v(3) bleibt unveraendert x(1) = xGrid1 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt c - Aufgeschlagen? if (Abs(x(2)).gt.dygrid1/2. .OR. + Abs(x(3)).gt.dzgrid1/2.) then flag = .true. destiny = code_wand else flag = .false. endif c - Gitterstab getroffen? if (testOnWireHit) then DrahtNr = nInt(x(2)/dist_Wires_G1) distToWire(1) = 0. distToWire(2) = x(2) - DrahtNr * dist_Wires_G1 call Test_WireHit(distToWire,WireRadiusQuad_G1,v(1),v(2),WireHit) if (WireHit) then flag = .true. destiny = code_grid endif endif c - Koordinatentransformation in Kammersystem: if (alfaTgt.NE.0.) then if (alfaTgtVertically) then help1 = x(3) help2 = v(1) help3 = v(3) x(1) = xgrid1 * Cos_alfaTgt - help1 * Sin_alfaTgt x(3) = xgrid1 * Sin_alfaTgt + help1 * Cos_alfaTgt v(1) = help2 * Cos_alfaTgt - help3 * Sin_alfaTgt v(3) = help2 * Sin_alfaTgt + help3 * Cos_alfaTgt else help1 = x(2) help2 = v(1) help3 = v(2) x(1) = xgrid1 * Cos_alfaTgt - help1 * Sin_alfaTgt x(2) = xgrid1 * Sin_alfaTgt + help1 * Cos_alfaTgt v(1) = help2 * Cos_alfaTgt - help3 * Sin_alfaTgt v(2) = help2 * Sin_alfaTgt + help3 * Cos_alfaTgt endif endif c - zerfallen? if (useDecay_) call Decay_Test(*555) c - falls aufgeschlagen: if (flag) goto 555 c - Koordinatentransformation zurueck in Beschleunigersystem: if (alfaTgt.NE.0.) then if (alfaTgtVertically) then x(1) = xGrid1 x(3) = help1 v(1) = help2 v(3) = help3 else x(1) = xGrid1 x(2) = help1 v(1) = help2 v(2) = help3 endif endif c - Datenausgabe: if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zweite Beschleunigerstufe: (homogenes Feld) 2 Gebiet = upToGrid2 steps = Steps + 1 if (a2.NE.0.) then help1 = v(1)*v(1) + 2.*a2*(XGrid2-XGrid1) if (help1.LT.0) then ! Reflektion noch vor 2. Gitter dt = -2*v(1)/a2 t = t + dt !x(1) bleibt unveraendert x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt v(1) = -v(1) !v(2) bleibt unveraendert !v(3) bleibt unveraendert destiny = code_reflektiert goto 555 endif dt = (sqrt(help1) - v(1))/a2 v(1) = v(1) + a2*dt else if (v(1).EQ.0) then ! (kann nur bei Start in 2. Stufe passieren) write(*,*) write(*,*)'ERROR: v(x) beim Start = 0. und '// + 'Beschleunigung = 0' write(*,*) STOP endif dt = (XGrid2-XGrid1) / v(1) endif t = t + dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt c - Aufgeschlagen? if (Abs(x(2)).gt.dygrid2/2. .OR. + Abs(x(3)).gt.dzgrid2/2.) then flag = .true. destiny = code_wand else flag = .false. ! <- noetig, falls Start auf 1. Gitter endif c - Gitterstab getroffen? if (testOnWireHit) then DrahtNr = nInt(x(2)/dist_Wires_G2) distToWire(1) = 0 distToWire(2) = x(2) - DrahtNr * dist_Wires_G2 call Test_WireHit(distToWire,WireRadiusQuad_G2,v(1),v(2),WireHit) if (WireHit) then flag = .true. destiny = code_grid endif endif c - Koordinatentransformation in Kammersystem: if (alfaTgt.NE.0.) then if (alfaTgtVertically) then x(1) = xgrid2 * Cos_alfaTgt - x(3) * Sin_alfaTgt x(3) = xgrid2 * Sin_alfaTgt + x(3) * Cos_alfaTgt help1 = v(1) v(1) = help1 * Cos_alfaTgt - v(3) * Sin_alfaTgt v(3) = help1 * Sin_alfaTgt + v(3) * Cos_alfaTgt else x(1) = xgrid2 * Cos_alfaTgt - x(2) * Sin_alfaTgt x(2) = xgrid2 * Sin_alfaTgt + x(2) * Cos_alfaTgt help1 = v(1) v(1) = help1 * Cos_alfaTgt - v(2) * Sin_alfaTgt v(2) = help1 * Sin_alfaTgt + v(2) * Cos_alfaTgt endif else x(1) = xgrid2 endif c - zerfallen? if (useDecay_) call Decay_Test(*555) c - falls aufgeschlagen: if (flag) goto 555 c - Datenausgabe: if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen zweitem Gitter und He-Shield: (feldfrei) 3 Gebiet = upToHeShield Steps = Steps + 1 radiusQuad = x(1)*x(1) + x(2)*x(2) help1 = v(1)*v(1)+v(2)*v(2) help2 = x(1)*v(1)+x(2)*v(2) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_HeShield))-help2) + /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 das Schild x(3) = x(3) + dt*v(3) ! durchquert if (useDecay_) call Decay_Test(*555) if (Abs(x(2)).gt.DYHESHIELD/2. .OR. + Abs(x(3)).gt.DZHESHIELD/2.) then destiny = code_wand goto 555 endif if (Debug_) call Output_Debug c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c Groessen bei x=40 mm berechnen: if (NTP_40mm) then dt = (40-x(1))/v(1) x40(2) = x(2)+v(2)*dt x40(3) = x(3)+v(3)*dt v40(1) = v(1) v40(2) = v(2) v40(3) = v(3) t40 = t + dt ! help1 = v(1)*v(1)+v(2)*v(2) noch bekannt von 'upToHeShield' v_square = help1 + v(3)*v(3) E40 = v_square * Energie_Faktor endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen He-Shield und LN-Shield: (feldfrei) 4 Gebiet = upToLNShield Steps = Steps + 1 radiusQuad = x(1)*x(1) + x(2)*x(2) ! help1 = v(1)*v(1)+v(2)*v(2) ! noch bekannt von 'upToHeShield' help2 = x(1)*v(1)+x(2)*v(2) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_LNShield))-help2) + /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 das Schild x(3) = x(3) + dt*v(3) ! durchquert if (useDecay_) call Decay_Test(*555) if (abs(x(2)).gt.dyLNShield/2. .OR. + Abs(x(3)).gt.dzLNShield/2.) then destiny = code_wand goto 555 endif if (Debug_) call Output_Debug c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen LN-Shield und Beginn der L1-Mappe: (feldfrei) 5 Gebiet = upToL1Map Steps = Steps + 1 dt = (xEnterMap_L1 - x(1)) / v(1) t = t + dt x(1) = xEnterMap_L1 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug if (radiusQuad.GT.radiusQuad_L1) then ! Teilchen fliegt an L1 vorbei destiny = code_vorbei goto 555 endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c innerhalb L1: (inhom. Felder -> Integrationen) 6 Gebiet = upToExL1 ! GebietsNummer fuer L1 setzen c Teste, ob das Teilchen ueberhaupt eine Beschleunigung erfaehrt (Spannung=0?, c Ladung=0?). Falls nicht, steppe gleich bis zum Mappenende: if (Beschl_Faktor_L1.EQ.0) then d dtmax_L1 = 0. d dtmin_L1 = 0. dt = (xLeaveMap_L1 - x(1)) / v(1) ! Zeit bis zum Mappenende x(1) = xLeaveMap_L1 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt goto 5106 endif c............................................................................... c Das Teilchen spuert eine Beschleunigung, es muss also integriert werden. c Gehe als ersten Versuch 0.5 mm in das Gebiet hinein: dt = .5/v(1) zaehler = 0 c............................................................................... c hierher wird zurueckgesprungen, solange die Integration in der L1 bleibt 5006 call INTEGRATIONSSTEP_RUNGE_KUTTA_L1(dt) d if (NTP_steps) then d if (dt.LT.dtmin_L1) then d dtmin_L1 = dt d x_dtmin_L1(1) = x(1) d x_dtmin_L1(2) = x(2) d x_dtmin_L1(3) = x(3) d endif d if (dt.GT.dtmax_L1) then d dtmax_L1 = dt d x_dtmax_L1(1) = x(1) d x_dtmax_L1(2) = x(2) d x_dtmax_L1(3) = x(3) d endif d endif c............................................................................... 5106 Steps = Steps + 1 ! neuer Ort, Zeit und Geschwindigkeit sind festgelegt c do some tests: if (destiny.EQ.code_dtSmall) then ! n_dtsmall>maxBelowDtSmall goto 555 elseif (destiny.EQ.code_wand) then ! aufgeschlagen radiusQuad = x(2)*x(2) + x(3)*x(3) ! -> den Ort berechnen, an help1 = v(2)*v(2)+v(3)*v(3) ! dem das Teilchen auf- help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_L1))-help2) + /help1 t = t + dt x(1) = x(1) + dt*v(1) x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? goto 555 c elseif (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'L1-1: ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP elseif (useDecay_) then ! zerfallen? call Decay_Test(*555) endif if (x(1).LT.xEnterMap_L1) then if (v(1).LT.0) then ! reflektiert? destiny = code_reflektiert goto 555 else ! darf nicht sein! write(*,*) write(*,*)' L1: x(1).LT.xEnterMap .AND. v(1).GE.0. -> STOP' write(*,*) STOP endif elseif (Steps.GE.MaxStep) then ! Teilchen verloren destiny = code_lost goto 555 elseif (x(1).GE.xLeaveMap_L1) then ! Verlasse L1 dt = (xLeaveMap_L1 - x(1))/v(1) ! rechne zurueck auf exaktes t = t + dt ! Mappenende x(1) = xLeaveMap_L1 x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug goto bis_Spiegel ! -> Mache bei upToEnSp weiter c elseif (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'L1-2: ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP 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 goto 5006 ! naechster Integrationsschritt in L1-Mappe c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen Linse 1 und Spiegel: (feldfrei) 7 Gebiet = upToEnSp Steps = Steps + 1 c - berechne Schnittpunkt mit forderer Spiegelebene: help2 = v(2)/v(1) ! Steigung der Bahn in der x-y-Ebene if (help2.GE.Tan_alfaSp) then ! Teilchen fliegt am Spiegel vorbei dt = (600-x(1))/v(1) t = t + dt x(1) = x(1) + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt if (useDecay_) call Decay_Test(*555) destiny = code_vorbei goto 555 else ! help1 == neues x(1) help1 = (x(2) - y_intersectSP + xSpiegel*Tan_alfaSp + - xLeaveMap_L1*help2) / (Tan_alfaSp - help2) dt = (help1-x(1)) / v(1) t = t + dt x(1) = help1 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt endif if (useDecay_) call Decay_Test(*555) if (Debug_) call Output_Debug c Berechnung der Trajektorie bei idealem Spiegel: if (idealMirror) then ! ~~~ 40: if ~~~~~~~~~~~ c - pruefe, ob das Teilchen die ForderSEITE des Spiegels trifft: if ( x(2).GT.yUppLeft .OR. x(2).LT.yLowLeft .OR. + abs(x(3)).GT.HSpiegel/2.) then ! -> Teilchen fliegt am Spiegel vorbei dt = (600-x(1))/v(1) t = t + dt x(1) = x(1) + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt destiny = code_vorbei goto 555 endif c - pruefe, ob das Teilchen einen Gitterstab des Spiegels trifft: if (testOnWireHit) then help1 = x(2)-yLowLeft ! Abstand zum Bezugspunkt DrahtNr = nInt(help1/(Sin_alfaSp*dist_Wires_Sp)) distToWire(2) = help1 - DrahtNr * Sin_alfaSp*dist_Wires_Sp distToWire(1) = distToWire(2)/Tan_alfaSp call Test_WireHit(distToWire,WireRadiusQuad_Sp,v(1),v(2),WireHit) if (WireHit) then destiny = code_grid goto 555 endif endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c im Spiegel: (homogenes Feld) 8 Gebiet = upToExSp Steps = Steps + 1 c - pruefe, ob Teilchen nicht zuviel Energie senkrecht zum Spiegel hat: if (Spiegel_Faktor.EQ.0.) then ! Spannung=0. oder q=0 dt = (600-x(1))/v(1) t = t + dt x(1) = x(1) + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt destiny = code_durchSpiegel goto 555 endif ! help1 == Winkel in xy-Ebene zwischen Bewegungsrichtung und Spiegelfront help1 = alfaSp - atand(v(2)/v(1)) ! help2 = Geschw.Komponente senkrecht auf den Spiegel gerichtet ! help3 = Geschw.Komponente parallel zum Spiegel, zu positiven y hin v_xy = sqrt( v(1)*v(1) + v(2)*v(2) ) help2 = sind(help1) * v_xy help3 = cosd(help1) * v_xy if (help2*help2*Energie_Faktor.GE.q*U_Sp) then ! Teilchen tritt durch Spiegel durch dt = (600-x(1))/v(1) t = t + dt x(1) = x(1) + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt destiny = code_durchSpiegel goto 555 endif if (Graphics_) call Save_Graphics_Koord c - berechne Zeit, bis Teilchen wieder auf Spiegelforderseite ist: dt = help2 * Spiegel_Faktor ! Spiegel_Faktor == 2 / a t = t + dt c - berechne Versetzung in xy-Ebene, parallel zur Spiegelebene, c in 'positiver y-Richtung' (speichere in 'help1'): help1 = help3*dt c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c falls Graphikausgabe verlangt ist: c Um die Teilchenbahn im Innern des Spiegels anzudeuten, berechne die Orte bei c t+dt/4, t+td/2 und t+3dt/4. Bestimme dafuer erst die jeweilige Versetzung c senkrecht zur Spiegelebene aus dx = vx * t + 1/2 * a * t**2. c (speichere in help4): if (Graphics_) then help4 = help2*dt*.25 - (dt*dt*.0625)/Spiegel_faktor nKoord = nKoord + 1 xKoord(nKoord) = x(1)+help4*Sin_alfaSp+help1*.25*Cos_alfaSp yKoord(nKoord) = x(2)-help4*Cos_alfaSp+help1*.25*Sin_alfaSp zKoord(nKoord) = x(3) + v(3)*dt*.25 help4 = help2*dt*.50 - (dt*dt*.2500)/Spiegel_faktor nKoord = nKoord + 1 xKoord(nKoord) = x(1)+help4*Sin_alfaSp+help1*.50*Cos_alfaSp yKoord(nKoord) = x(2)-help4*Cos_alfaSp+help1*.50*Sin_alfaSp zKoord(nKoord) = x(3)+v(3)*dt*.50 help4 = help2*dt*.75 - (dt*dt*.5625)/Spiegel_faktor nKoord = nKoord + 1 xKoord(nKoord) = x(1)+help4*Sin_alfaSp+help1*.75*Cos_alfaSp yKoord(nKoord) = x(2)-help4*Cos_alfaSp+help1*.75*Sin_alfaSp zKoord(nKoord) = x(3)+v(3)*dt*.75 endif c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c - berechne Austrittsort: x(1) = x(1) + help1 * Cos_alfaSp x(2) = x(2) + help1 * Sin_alfaSp x(3) = x(3) + v(3)*dt c - berechne Austrittsgeschwindigkeit (help2 geht bei Spiegelung in -help2 ueber): v(1) = help3 * Cos_alfaSp - help2 * Sin_alfaSp v(2) = help2 * Cos_alfaSp + help3 * Sin_alfaSp if (v(2).LE.0) then write(*,*) write(*,*)'ERROR: nach Spiegel ist v(y) <= 0.' write(*,*) STOP endif if (useDecay_) call Decay_Test(*555) c - pruefe, ob Austrittspunkt auf forderer Spiegelflaeche liegt: if (x(2).GT.yUppLeft .OR. x(2).LT.yLowLeft .OR. + abs(x(3)).GT.hSpiegel/2.) then ! Teilchen trifft auf Spiegelwand destiny = code_wand goto 555 endif c - pruefe, ob das Teilchen einen Gitterstab des Spiegels trifft: if (testOnWireHit) then help1 = x(2)-yLowLeft ! Abstand zum Bezugspunkt DrahtNr = nInt(help1/(Sin_alfaSp*dist_Wires_Sp)) distToWire(2) = help1 - DrahtNr * Sin_alfaSp*dist_Wires_Sp distToWire(1) = distToWire(2)/Tan_alfaSp call Test_WireHit(distToWire,WireRadiusQuad_Sp,v(1),v(2),WireHit) if (WireHit) then destiny = code_grid goto 555 endif endif if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug goto 9 endif ! ~~~ 40: endif ~~~~ c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c innerhalb der Spiegelmappe (dx = 0.050 mm, dy = 0.050 mm) Gebiet = upToExSp nKoordSave = nKoord c Die Anweisungen dieses Abschnitts verlaufen analog zu denen c von Linse 1. -> Fuer Kommentare siehe dort! if (Beschl_Faktor_Sp.EQ.0. .OR. q.EQ.0) then d dtmax_Sp = 0. d dtmin_Sp = 0. dt = (600-x(1))/v(1) t = t + dt x(1) = x(1) + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt destiny = code_durchSpiegel goto 555 endif dt = 0.5/v(1) reachedEndOfMap = .false. zaehler = 0 c Rechne in Spiegelmappen-Koordinaten um: c Im Spiegelmappensystem: x-Achse verlaueft entlang der forderen Mappenkante, c y-Achse aus dem Spiegel heraus. (entgegen der Richtung zunehmender Mappen- c j-indizierung!) 5008 help1= x(1) - xSpiegel x(1) = - x(2)*Cos_alfaSp + help1*Sin_alfaSp + + (dSpiegel/2.+DreharmLaenge+xSpGrid1) x(2) = x(2)*Sin_alfaSP + help1*Cos_alfaSP help1= v(1) v(1) = - v(2)*Cos_alfaSp + help1*Sin_alfaSp v(2) = v(2)*Sin_alfaSP + help1*Cos_alfaSP c mache Integrationsschritt: call INTEGRATIONSSTEP_RUNGE_KUTTA_Sp(dt) ! setzt u.U. auch 'destiny' Steps = Steps + 1 c do some tests: if (Steps.GE.MaxStep) destiny = code_lost ! Teilchen verloren c - Potentialmappe nach Reflektion wieder verlasssen? if (x(1).LT.0) then reachedEndOfMap = .true. c - Spiegelrahmen getroffen? elseif (x(1).GE.xSpGrid1 .AND. + (abs(x(2)).GT.bSpiegel/2. .OR. abs(x(3)).GT.hSpiegel/2.)) then destiny = code_wand c - Gitterstab getroffen? else help1 = min(abs(x(1)-xSpGrid1),abs(x(1)-xSpGrid1)) if (help1.LE.rWires_Sp) then DrahtNr = nInt(x(2)/dist_Wires_Sp) distToWire(2) = x(2) - DrahtNr * dist_Wires_Sp if ( (help1*help1 + distToWire(2)*distToWire(2)).LE. + radiusQuad_Sp) destiny = code_grid endif endif c if (destiny.NE.code_ok) then c if (x(1).LT.xSpGrid1) then c if (v(1).GT.0) then c gebiet = UpToGrid c else c gebiet = upToExMap c endif c else c gebiet = RetToGrid c endif c endif c Rechne in Kammerkoordinaten zurueck: help1= x(1)-(dSpiegel/2.+DreharmLaenge+xSpGrid1) x(1) = help1*Sin_alfaSP + x(2)*Cos_alfaSP + xSpiegel x(2) = - help1*Cos_alfaSP + x(2)*Sin_alfaSP help1= v(1) v(1) = help1*Sin_alfaSP + v(2)*Cos_alfaSP v(2) = - help1*Cos_alfaSP + v(2)*Sin_alfaSP d if (NTP_steps) then d if (dt.LT.dtmin_Sp) then d dtmin_Sp = dt d x_dtmin_Sp(1) = x(1) d x_dtmin_Sp(2) = x(2) d x_dtmin_Sp(3) = x(3) d endif d if (dt.GT.dtmax_Sp) then d dtmax_Sp = dt d x_dtmax_Sp(1) = x(1) d x_dtmax_Sp(2) = x(2) d x_dtmax_Sp(3) = x(3) d endif d endif c zerfallen? if (useDecay_) call Decay_Test(*555) c Bahnberechnung abgebrochen? if (destiny.NE.code_ok) goto 555 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 if (.NOT.reachedEndOfMap) goto 5008 ! naechster Integrationsschritt c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen Spiegel und Koordinatenwechsel-Ebene y=xChangeKoord: (feldfrei) 9 Gebiet = upToChKoord Steps = Steps + 1 if (x(2).LT.xChangeKoord) then ! gegebenenfalls flag fuer Graphikausgabe des Punktes setzen flag = .true. else flag = .false. endif dt = (xChangeKoord - x(2)) / v(2) t = t + dt x(1) = x(1) + v(1)*dt x(2) = xChangeKoord x(3) = x(3) + v(3)*dt help4 = x(1)-xSpiegel radiusQuad = help4*help4 + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(1)*v(1)+v(3)*v(3) help2 = help4*v(1)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_.AND.flag) call Save_Graphics_Koord if (Debug_) call Output_Debug c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c falls Graphikausgabe verlangt ist: Gib jetzt die Trajektorie im 'horizontalen' c Teil der Kammer aus und resette nKoord: if (Graphics_) then call plot_horizontal if (schnitt_p.eq.1) call schnitt ! Schnittebene ! die letzten Koordinaten fuer Plot der Trajektorie im 2. Kammerteil ! uebernehmen (in neues KoordinatenSystem transformiert): k = nKoord if (idealMirror) then nKoord = 7 else if (nKoord.LT.nKoordSave) then ! => ein 'turn over' fand statt waehrend das Teilchen in der ! Spiegelmappe war => x(999) -> x(1), x(1000) -> x(2) nKoord = nKoord + (999-nKoordSave) else nKoord = nKoord - nKoordSave + 1 endif nKoord = nKoord-2 endif do i = nKoord, 1, -1 xKoord_(i) = yKoord(k) yKoord_(i) = xSpiegel - xKoord(k) zKoord_(i) = zKoord(k) k = k - 1 if (k.EQ.0) then k = 998 endif enddo do i = 1, nKoord xKoord(i) = xKoord_(i) yKoord(i) = yKoord_(i) zKoord(i) = zKoord_(i) enddo endif c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c - Vollziehe Koordinatenwechsel: neuer Ursprung in der Spiegelaufhaengung, c x-Richtung in bisherige y-Richtung (also wiederum entlang Strahlachse), c y-Richtung in bisheriger negativer x-Richtung. z-Richtung bleibt gleich. help1 = x(2) x(2) = xSpiegel - x(1) x(1) = help1 help1 = v(1) v(1) = v(2) v(2) = - help1 if (Debug_) then write (lun(1),*) 'KOORDINATENWECHSEL:' call Output_Debug endif c Beruecksichtige gegebenenfalls die Flugzeit in 'delta_L1', welches 'vor dem c Triggerdetektor' eingeschoben werden kann: dt = Delta_L1 / v(1) x(1) = x(1)+v(1)*dt x(2) = x(2)+v(2)*dt x(3) = x(3)+v(3)*dt t = t + dt c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: if (lense2) then ! ~~~~~~~~~~~~~~ ******* ~~~~~~~~~~~~~~~ c Bei 'lense2' wird fuer das Feld der Linse 2 und das Feld der TD-Folie eine c gemeinsame Mappe verwendet. Hierbei ist allerdings der Triggerwinkel auf 0 c Grad festgelegt. Da es in Zukunft in der Praxis wohl kaum noch vorkommen wird, c dass der Triggerdetektor verdreht wird, sollte diese Einschraenkung jedoch c keine grossen Auswirkungen haben. c Ist der Triggerdetektor nicht im Strahl, so wird der Anteil der Triggerfolie c gleich Null gesetzt. c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c zwischen KOORDINATENWECHSEL und Beginn der L2andFo-Mappe: (feldfrei) 10 Gebiet = upToL2andFoMap Steps = Steps + 1 dt = (xEnterMap_L2andFo - x(1)) / v(1) t = t + dt x(1) = xEnterMap_L2andFo x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug if (radiusQuad.GT.radiusQuad_L2) then ! Teilchen fliegt an L2 vorbei destiny = code_vorbei goto 555 endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c innerhalb der gemeinsamen Mappe von Linse 2 und dem Feld vor der Trigger- c Detektor-Folie: 11 Gebiet = upToExL2 ! Gebietsnummer fuer L2 setzen c Die Anweisungen dieses Abschnitts verlaufen analog zu denen c von Linse 1. -> Fuer Kommentare siehe dort! if (Beschl_Faktor.EQ.0. .OR. (U_L2.EQ.0. AND. U_F.EQ.0.)) then c WRITE(*,*) 'HALLOHALLO!' d dtmax_L2andFo = 0. d dtmin_L2andFo = 0. dt = (xEndLense_L2 - x(1)) / v(1) ! Zeit bis zum Linsenende x(1) = xEndLense_L2 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_L2) then destiny = code_wand radiusQuad_ = radiusQuad_L2 goto 5111 endif if (TriggerInBeam) then Gebiet = upToEnTD ! Gebietsnummer fuer upToTD setzen ! Zeit bis zum Mappenende (falls TD im Strahl: bis Triggerfolie) dt = (xLeaveMap_L2andFo - xEndLense_L2) / v(1) x(1) = xLeaveMap_L2andFo x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then destiny = code_wand radiusQuad_ = radiusQuad_Rohr endif endif goto 5111 endif dt = .5/v(1) zaehler = 0 reachedEndOfMap = .false. c die Integrationsroutine will x bereits relativ zum Mappenanfang geliefert c bekommen: 5011 x(1) = x(1) - xEnterMap_L2andFo call INTEGRATIONSSTEP_RUNGE_KUTTA_L2(dt) x(1) = x(1) + xEnterMap_L2andFo d if (NTP_steps) then d if (dt.LT.dtmin_L2andFo) then d dtmin_L2andFo = dt d x_dtmin_L2andFo(1) = x(1) d x_dtmin_L2andFo(2) = x(2) d x_dtmin_L2andFo(3) = x(3) d endif d if (dt.GT.dtmax_L2andFo) then d dtmax_L2andFo = dt d x_dtmax_L2andFo(1) = x(1) d x_dtmax_L2andFo(2) = x(2) d x_dtmax_L2andFo(3) = x(3) d endif d endif 5111 Steps = Steps + 1 if (destiny.EQ.code_dtSmall) then ! n_dtsmall>maxBelowDtSmall goto 555 elseif (destiny.EQ.code_wand) then ! aufgeschlagen radiusQuad = x(2)*x(2) + x(3)*x(3) ! -> den Ort berechnen, an help1 = v(2)*v(2)+v(3)*v(3) ! dem das Teilchen auf- help2 = x(2)*v(2)+x(3)*v(3) ! schlaegt dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_))-help2)/help1 t = t + dt x(1) = x(1) + dt*v(1) x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? goto 555 elseif (useDecay_) then ! zerfallen? call Decay_Test(*555) endif if (Gebiet.EQ.upToExL2) then ! ----> noch innerhalb von Linse 2 radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_L2) then destiny = code_wand radiusQuad_ = radiusQuad_L2 goto 5111 endif if (x(1).LT.xEnterMap_L2andFo) then if (v(1).LT.0) then ! reflektiert destiny = code_reflektiert goto 555 else ! darf nicht sein! write(*,*) write(*,*)' L2: x(1).LT.xEnterMap .AND. v(1).GE.0. -> STOP' write(*,*) call exit endif elseif (x(1).GT.xEndLense_L2) then ! Verlasse L2 Gebiet = upToEnTD endif else ! ----> zw. Linse 2 und TD-Folie: c if (x(1).EQ.xLeaveMap_L2andFo) then ! Verlasse Mappe if (reachedEndOfMap) then ! Verlasse Mappe c WRITE(*,*) 'HALLO: x(1).EQ.xLeaveMap_L2andFo !!' ! ==================================================== ! muss in Integrationsroutine richtig abgestimmt sein! ! ==================================================== if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug if (TriggerInBeam) then ! rechne in Triggerkoordinaten um (Folie == x=0) x(1) = 0 goto 112 else goto bis_L3_Mappe endif endif if (radiusQuad.GT.radiusQuad_Rohr) then destiny = code_wand radiusQuad_ = radiusQuad_Rohr goto 5111 endif endif if (Steps.GE.MaxStep) then ! Teilchen verloren destiny = code_lost goto 555 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 5011 ! naechster Integrationsschritt in gleicher Feldmappe endif ! if (lense2) then.... ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: if (.NOT.TriggerInBeam) goto bis_L3_Mappe c zwischen Koordinatenwechselebene und Triggerfolie: (feldfrei) 12 Gebiet = upToEnTD c Die Anweisungen dieses Abschnitts verlaufen in weiten Teilen parallel zu denen c von Linse 1. -> Fuer Kommentare zu diesen Bereichen siehe dort! if (Beschl_Faktor_FO.EQ.0. .OR. gridInFrontOfFoil) then ! => keine Integration in der Folienmappe Steps = Steps + 1 help1 = v(2)/v(1) ! Steigung der Bahn in der x-y-Ebene d dtmin_FO = 0. d dtmax_FO = 0. c - berechne Schnittpunkt der Trajektorie mit Ebene der Triggerfolie bzw. bei c 'GridInFrontOfFoil' mit Ebene des Gitters vor der Triggerfolie: c Folienebene : y'= (x_intersectTD - x') / Tan_alfaTD c Trajektorie : y'= y + v(2)/v(1)*(x'-x) = y + help1*(x'-x) c => Schnittpunkt: x'= (x_intersectTD/Tan_alfaTD - y + help1*x)/(help1 + 1/Tan_alfaTD) c = (x_intersectTD + Tan_alfaTD*(help1*x-y))/(1+help1*Tan_alfaTD) c (erste Gleichung hat Probleme bei Tan_alfaTD = 0!) if (atand(help1).EQ.alfaTD-90) then ! ueberpruefen<<<<<<<<<<<<<<<<<< ! Teilchen fliegt parallel zur Folie => fliegt an TD vorbei destiny = code_vorbei goto 555 else ! help2 == neues x(1) if (Tan_alfaTD.EQ.0) then dt = (x_intersectTD-x(1)) / v(1) x(1) = x_intersectTD else help2 = (x_intersectTD+Tan_alfaTD* + (help1*xChangeKoord-x(2)))/(1+help1*Tan_alfaTD) if (help2.LT.xChangeKoord) then ! Teilchen fliegt 'steiler' als Folienebene ! -> kein Schnittpunkt mit dt.gt.0 => fliegt an TD vorbei destiny = code_vorbei goto 555 else ! Bahntangente kreuzt Folienebene dt = (help2-x(1)) / v(1) x(1) = help2 endif endif endif c -> Teilchenort in Folienebene bzw. bei 'GridInFrontOfFoil' in Ebene des c geerdeten Gitters vor der Triggerfolie: x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug c Koordinatentransformation vom Kammersystem in das System des Triggerdetektors: c (Ursprung in Folienmitte, x-Achse senkrecht zur Folie, y-Achse parallel zur c Folie. Wenn der TD nicht verdreht ist, verlaufen die Achsen parallel zu c denen des Kammersystems): if (alfaTD.NE.0) then x(2) = (xTD-x(1))*Sin_alfaTD + x(2)*Cos_alfaTD help1= v(1) ! zur Zwischenspeicherung v(1) = help1*Cos_alfaTD + v(2)*Sin_alfaTD v(2) = -help1*Sin_alfaTD + v(2)*Cos_alfaTD endif if (.NOT.GridInFrontOfFoil) then x(1) = 0 else ! -> berechne Schnittpunkt der Trajektorie mit Folienebene unter ! der Annahme einer idealen Potentialrampe: if (aFoil.NE.0.) then help1 = v(1)*v(1) + 2.*aFoil*(d_Grid_Folie) if (help1.LT.0) then ! Reflektion noch vor Folie dt = -2*v(1)/aFoil t = t + dt x(1) = - d_Grid_Folie x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt v(1) = -v(1) destiny = code_reflektiert goto 555 endif dt = (sqrt(help1) - v(1))/aFoil ! (ergibt sich aus x=v*t+1/2*a*t**2 mit richtiger V.Z.-Wahl ('+')) v(1) = v(1) + aFoil*dt else dt = d_Grid_Folie / v(1) endif t = t + dt x(1) = 0 ! im Triggersystem x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt endif ! if (GridInFrontOfFoil) ... goto 112 c............................................................................... else ! (if Beschl_Faktor_FO.EQ.0 .OR. GridInFrontOfFoil) then ... ! => Integration in der Folienmappe: ! alte Version: ab xChangeKoord wurde integriert, wobei das EFeld im ! Bereich vor der Mappe als 0 zurueckgegeben wurde. ! Ab Version 1.2.9: Bis Schnittpunkt der Trajektorie mit Beginn der ! Potentialmappe wird extrapoliert, dann erst integriert: c Einschub ab Version 1.2.9: *************************************************** Steps = Steps + 1 help1 = v(2)/v(1) ! Steigung der Bahn in der x-y-Ebene c - berechne Schnittpunkt der Trajektorie mit Beginn der Potentialmappe: c Mappenebene : y'= (x_intersectTDMap - x') / Tan_alfaTD c Trajektorie : y'= y + v(2)/v(1)*(x'-x) = y + help1*(x'-x) c => Schnittpunkt: x'= (x_intersectTDMap/Tan_alfaTD - y + help1*x)/(help1 + 1/Tan_alfaTD) c = (x_intersectTDMap + Tan_alfaTD*(help1*x-y))/(1+help1*Tan_alfaTD) c (erste Gleichung hat Probleme bei Tan_alfaTD = 0!) if (atand(help1).EQ.alfaTD-90) then ! ueberpruefen<<<<<<<<<<<<<<<<<< ! Teilchen fliegt parallel zur Mappe => fliegt an TD vorbei destiny = code_vorbei goto 555 ! stimmt so u.U. noch nicht ganz. Kommt aber eigentlich eh nie vor! ! (stimmt bis jetzt wohl nur fuer positive alpha(TD) else if (Tan_alfaTD.EQ.0) then dt = (x_intersectTDMap-x(1)) / v(1) x(1) = x_intersectTDMap else ! help2 == neues x(1): help2 = (x_intersectTDMap+Tan_alfaTD* + (help1*xChangeKoord-x(2)))/(1+help1*Tan_alfaTD) ! folgendes herauskommentiert, da es teilweise passierte, dass ! der Mappenanfang ueber xChangekoord hinausreichte und die ! Trajektorien dann faelschlicherweise abgebrochen worden sind. c if (help2.LT.xChangeKoord) then c ! Teilchen fliegt 'steiler' als Mappenebene c ! -> kein Schnittpunkt mit dt.gt.0 => fliegt an TD vorbei c destiny = code_vorbei c goto 555 c else ! Bahntangente kreuzt Mappenebene dt = (help2-x(1)) / v(1) x(1) = help2 c endif endif endif c -> Teilchenort in Mappenebene: x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug c Ende des Einschubes ab Version 1.2.9: **************************************** c => Jetzt erfolgt Start in die Folienmappe: reachedEndOfMap = .false. ! Folienebene wurde noch nicht erreicht dt = .5/v(1) ! 1. Testschritt 0.5 mm in x-Richtung zaehler = 0 c Rechne in Folienmappen-Koordinaten um: 5012 if (alfaTD.NE.0) then help1= x(1)-xTD x(1) = help1*Cos_alfaTD + x(2)*Sin_alfaTD + length1 x(2) = -help1*Sin_alfaTD + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD + v(2)*Sin_alfaTd v(2) = -help1*Sin_alfaTD + v(2)*Cos_alfaTd else x(1) = x(1) - length2 endif c mache Integrationssschritt: call INTEGRATIONSSTEP_RUNGE_KUTTA_FO(dt) d if (NTP_steps) then d if (dt.LT.dtmin_FO) then d dtmin_FO = dt d x_dtmin_FO(1) = x(1) d x_dtmin_FO(2) = x(2) d x_dtmin_FO(3) = x(3) d endif d if (dt.GT.dtmax_FO) then d dtmax_FO = dt d x_dtmax_FO(1) = x(1) d x_dtmax_FO(2) = x(2) d x_dtmax_FO(3) = x(3) d endif d endif c Rechne in Kammerkoordinaten zurueck: if (alfaTD.NE.0) then help1= x(1)-length1 x(1) = help1*Cos_alfaTD - x(2)*Sin_alfaTD + xTD x(2) = help1*Sin_alfaTD + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD - v(2)*Sin_alfaTD v(2) = help1*Sin_alfaTD + v(2)*Cos_alfaTD else x(1) = x(1) + length2 endif Steps = Steps + 1 ! neuer Ort, Zeit und Geschwindigkeit sind festgelegt c do some tests: if (destiny.EQ.code_dtSmall) then ! n_dtSmall>maxBelowDtSmall goto 555 endif radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then ! aufgeschlagen help1 = v(2)*v(2)+v(3)*v(3) ! -> den Ort berechnen, an help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /help1 t = t + dt x(1) = x(1) + dt*v(1) x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 elseif (useDecay_) then ! zerfallen? call Decay_Test(*555) endif if (destiny.EQ.code_reflektiert) then ! reflektiert goto 555 c elseif (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'FO-1: ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP elseif (Steps.GE.MaxStep) then ! Teilchen verloren destiny = code_lost goto 555 endif if (reachedEndOfMap) then ! Folienebene erreicht ! rechne in Triggerkoordinaten um (Folie == x=0) if (alfaTD.NE.0) then x(2) = (xTD-x(1))*Sin_alfaTD + x(2)*Cos_alfaTD help1= v(1) ! zur Zwischenspeicherung v(1) = help1*Cos_alfaTD + v(2)*Sin_alfaTD v(2) = -help1*Sin_alfaTD + v(2)*Cos_alfaTD endif x(1) = 0 goto 112 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 goto 5012 ! naechster Integrationsschritt in FELD VOR FOLIE c............................................................................... endif ! (if Beschl_Faktor_FO.EQ.0) then ... c Einsprunglabel fuer Starts auf der Triggerfolie mit Startwinkelangaben c im Kammersystem => transformiere Geschwindigkeitsvektor in das Triggersystem: 111 if (alfaTD.NE.0) then help1= v(1) ! zur Zwischenspeicherung v(1) = help1*Cos_alfaTD + v(2)*Sin_alfaTD v(2) = -help1*Sin_alfaTD + v(2)*Cos_alfaTD endif c - pruefe, ob das Projektil die Folie trifft: 112 radiusQuad = x(2)*x(2) + x(3)*x(3) If (radiusQuad.GT.radiusQuad_Folie) then ! zurueckrechnen in das Kammersystem: if (alfaTD.NE.0) then help1= x(1) x(1) = (help1-d_Folie_Achse)*Cos_alfaTD - + x(2)*Sin_alfaTD + xTD x(2) = (help1-d_Folie_Achse)*Sin_alfaTD + + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD - v(2)*Sin_alfaTD v(2) = help1*Sin_alfaTD + v(2)*Cos_alfaTD else x(1) = x(1) + xTD - d_Folie_Achse endif destiny = code_vorbei goto 555 endif c So verlangt, schreibe die aktuellen Trajektoriengroessen in das 'FoilFile': c (hier ist sichergestellt, dass die Folie getroffen worden ist, Wechsel- c wirkungen mit der Folie wurden aber noch nicht beruecksichtigt). c HIER WERDEN 'X' UND 'V' IM TRIGGERSYSTEM ABGESPEICHERT! if (createFoilFile) then ! falls die Flugzeit bis zur Triggerfolie verschmiert in das ! NTupel aufgenommen werden soll: if (smearS1Fo) then call Gauss_Verteilung(sigmaS1Fo,help4) S1FoOnly = t + help4 endif if (NTP_stop) then Ekin=(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))*Energie_Faktor endif call HFNT(NTP_write) NTPalreadyWritten = .true. endif c - Zeitpunkt bei Erreichen der Folie sichern: 113 S1Fo = t if (createNTP.AND.Fo_triggered) fill_NTP = .true. if (statNeeded(Nr_S1Fo)) call fill_statMem(S1Fo,Nr_S1Fo) c - Speichern der Koordinaten fuer die Statistiken: if (statNeeded(Nr_y_Fo)) then call fill_statMem( x(2),Nr_y_Fo) endif if (statNeeded(Nr_z_Fo)) then call fill_statMem( x(3),Nr_z_Fo) endif if (statNeeded(Nr_r_Fo)) then radius = SQRT(x(2)*x(2) + x(3)*x(3)) call fill_statMem(radius,Nr_r_Fo) endif c - speichere Auftreffort des Projektils fuer die Berechnung der Folienelektronen: if (generate_FE) then x0FE(1) = x(1) x0FE(2) = x(2) x0FE(3) = x(3) endif c - falls nur bis zur Folie gerechnet werden soll, beende hier die Integration: if (upToTDFoilOnly) then ! zurueckrechnen in das Kammersystem: if (alfaTD.NE.0) then help1= x(1) x(1) = (help1-d_Folie_Achse)*Cos_alfaTD - + x(2)*Sin_alfaTD + xTD x(2) = (help1-d_Folie_Achse)*Sin_alfaTD + + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD - v(2)*Sin_alfaTD v(2) = help1*Sin_alfaTD + v(2)*Cos_alfaTD else x(1) = x(1) + xTD - d_Folie_Achse endif if (generate_FE) Gebiet = UpToExTD goto 555 endif c - pruefe, ob das Projektil auf das Stuetzgitter aufschlaegt: if (testOnWireHit .AND. ran(seed).GT.TransTDFoil) then destiny = code_Stuetzgitter ! zurueckrechnen in das Kammersystem: if (alfaTD.NE.0) then help1= x(1) x(1) = (help1-d_Folie_Achse)*Cos_alfaTD - + x(2)*Sin_alfaTD + xTD x(2) = (help1-d_Folie_Achse)*Sin_alfaTD + + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD - v(2)*Sin_alfaTD v(2) = help1*Sin_alfaTD + v(2)*Cos_alfaTD else x(1) = x(1) + xTD - d_Folie_Achse endif goto 555 endif c - Energieverlust und Winkelaufstreuung: if (log_E_Verlust .OR. log_Aufstreu) then if (Debug_) then Steps = Steps + 1 call Output_Debug endif v_square = v(1)*v(1) + v(2)*v(2) + v(3)*v(3) v_Betrag = SQRT(v_square) Ekin = v_square * Energie_Faktor endif c -- Energieverlust (vorerst nur Gaussverteilt): if (log_E_Verlust_defined.OR.log_Meyer_Gauss) then ! Berechne Bahnwinkel relativ zur Folienebene fuer effektive Folien- ! dicke: alfa = atand(SQRT(v(2)*v(2)+v(3)*v(3))/v(1)) endif if (log_E_Verlust) then if (calculate_each) then call CALC_ELOSS_ICRU(Ekin,q,m,Thickness,E_Verlust) else E_Verlust = mean_E_Verlust endif if (log_E_Verlust_defined) E_Verlust = E_Verlust / cosd(alfa) if (debug_) write (lunLOG,*) ' mittlerer Energieverlust: ',E_Verlust ! Now we have the mean energy loss. We still have to modify it ! according to the distribution of energy losses, i.e. ! E_Verlust -> E_Verlust + delta_E_Verlust: delta_E_Verlust = 0. if (log_E_Straggling_sigma) then 400 call Gauss_Verteilung(sigmaE,delta_E_Verlust) if (debug_) write (lunLOG,*) ' sigmaE,delta_E_Verlust: ',sigmaE,delta_E_Verlust if (E_Verlust+delta_E_Verlust.LT.0.) goto 400 elseif (log_E_Straggling_equal) then 410 delta_E_Verlust = lowerE + (upperE - lowerE)*ran(seed) if (E_Verlust+delta_E_Verlust.LT.0) goto 410 elseif (log_E_Straggling_Lindhard) then ! Streuung in Abhaengigkeit von mittlerer Energie in Folie: call E_Straggling_Lindhard(Ekin-0.5*E_Verlust,m,sigmaE) 420 call Gauss_Verteilung(sigmaE,delta_E_Verlust) if (debug_) write (lunLOG,*) ' sigmaE,delta_E_Verlust: ',sigmaE,delta_E_Verlust if (E_Verlust+delta_E_Verlust.LT.0.) goto 420 elseif (log_E_Straggling_Yang) then ! Streuung in Abhaengigkeit von mittlerer Energie in Folie! call E_Straggling_Yang(Ekin-0.5*E_Verlust,m,sigmaE) 430 call Gauss_Verteilung(sigmaE,delta_E_Verlust) if (debug_) write (lunLOG,*) ' sigmaE,delta_E_Verlust: ',sigmaE,delta_E_Verlust if (E_Verlust+delta_E_Verlust.LT.0.) goto 430 endif if (E_Verlust+delta_E_Verlust.GE.Ekin) then destiny = code_stopped_in_foil goto 555 endif E_Verlust = E_Verlust + delta_E_Verlust ! help1 == Reduzierungsfaktor fuer Geschw.Betrag help1 = sqrt( (Ekin - E_Verlust)/Ekin ) v(1) = help1 * v(1) v(2) = help1 * v(2) v(3) = help1 * v(3) v_Betrag = help1 * v_Betrag if (debug_) write (lunLOG,*) ' Energieverlust: ',E_Verlust endif c -- Winkelaufstreuung (vorerst nur Gaussverteilt): if (log_aufstreu) then if (log_Meyer_F_Function) then call throwMeyerAngle(thetaAufstreu) else if (log_Meyer_Gauss) then if (log_E_Verlust) Ekin = Ekin - .5 * E_Verlust ! mittlere Energie effRedThick = Meyer_Faktor1 * Thickness / cosd(alfa) call g_Functions(g1,g2,effRedThick) sigmaAufstreu = Meyer_Faktor2 / Ekin * (g1 + Meyer_Faktor3*g2) if (debug_) then write (lunLOG,*) ' effekt. red. Dicke: ',effRedThick write (lunLOG,*) ' Sigma(Streuwinkel): ',sigmaAufstreu endif endif call Gauss_Verteilung_theta(sigmaAufstreu,thetaAufstreu) endif st0 = sind(thetaAufstreu) ct0 = cosd(thetaAufstreu) phiAufstreu = 360.*ran(seed) v_xy = SQRT(v(1)*v(1) + v(2)*v(2)) ! v_xy stets groesser 0 ! wegen v(1)>0 help1 = v(1) help2 = v(2) help3 = v_Betrag*st0*cosd(phiAufstreu)/v_xy help4 = st0*sind(phiAufstreu) v(1) = ct0*help1 - help3*help2 - help4*help1*v(3)/v_xy v(2) = ct0*help2 + help3*help1 - help4*help2*v(3)/v_xy v(3) = ct0*v(3) + help4*v_xy if (debug_) write (lunLOG,*) ' Aufstreuung: theta, phi =', + thetaAufstreu,phiAufstreu endif if (Debug_ .AND. (log_E_Verlust .OR. log_Aufstreu)) then call Output_Debug endif c - Neutralisierung in der Folie? if (log_neutralize) then if (neutral_fract(q_).EQ.-1.0) then v_square = v(1)*v(1) + v(2)*v(2) + v(3)*v(3) Ekin = v_square * Energie_Faktor call chargeStateYields(Ekin,m,YieldPlus,YieldNeutral) YieldNeutral = 100. * YieldNeutral else YieldNeutral = neutral_fract(q_) endif if (100.*ran(seed).LE.YieldNeutral) then q = 0. qInt = 0 if (debug_) then write (lunLOG,*) ' Teilchen wurde neutralisiert' endif nNeutral = nNeutral + 1 else nCharged = nCharged + 1 endif endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c im TriggerDetektor: (homogene Felder) 13 Gebiet = upToExTD Steps = Steps + 1 c der Weg des Projektils innerhalb der Triggerkammer: call TRIGGER(m,q,t,x,v,Debug_,graphics_,n_return) c Koordinatentransformation vom System des Triggerdetektors in das Kammersystem: c ('d_Achse_Ground' == Abstand zwischen TD-Achse und 'Ground'-Gitter) if (alfaTD.NE.0) then help1= x(1) x(1) = (help1-d_Folie_Achse)*Cos_alfaTD - + x(2)*Sin_alfaTD + xTD x(2) = (help1-d_Folie_Achse)*Sin_alfaTD + + x(2)*Cos_alfaTD help1= v(1) v(1) = help1*Cos_alfaTD - v(2)*Sin_alfaTD v(2) = help1*Sin_alfaTD + v(2)*Cos_alfaTD else x(1) = x(1) + xTD - d_Folie_Achse endif c Was ist im TD mit dem Teilchen passiert? if (n_return.NE.0) then ! -->> das Projektil kam nicht bei GROUND an if (n_return.GT.100 .AND. n_return.LE.120) then ! -> abgebrochen statTD(1,n_return-100) = statTD(1,n_return-100)+1 ! Grund notieren destiny = code_lostInTD ! im TD verloren elseif (n_return.GT.0..AND.n_return.LE.75) then ! -> pfosten getroffen! pfostenHit(n_return,1) = pfostenHit(n_return,1)+1 destiny = code_wand elseif (n_return.EQ.-5) then ! -> im TD auf Gitterstab statTD(1,17) = statTD(1,17)+1 ! destiny = code_grid elseif (n_return.EQ.-9) then ! -> NICHT im MCP3 registriert statTD(1,18) = statTD(1,18)+1 ! destiny = code_notRegInM3 elseif (n_return.EQ.-10) then ! -> im MCP3 registriert statTD(1,16) = statTD(1,16)+1 ! '16' zaehlt MCP3-Treffer destiny = code_wand endif goto 555 ! naechstes Projektil else ! -->> das Projektil kam bei GROUND an statTD(1,15) = statTD(1,15)+1 ! '15' zaehlt GROUND-Treffer endif if (useDecay_) call Decay_Test(*555) c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen KOORDINATENWECHSEL BZW. GROUND-GITTER und Beginn der L3-Mappe: c (feldfrei) 14 Gebiet = upToL3Map Steps = Steps + 1 dt = (xEnterMap_L3 - x(1)) / V(1) t = t + dt x(1) = xEnterMap_L3 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug if (radiusQuad.GT.radiusQuad_L3) then ! Teilchen fliegt an L3 vorbei destiny = code_vorbei goto 555 endif c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c innerhalb L3: (inhom. Felder -> Integrationen) 15 Gebiet = upToExL3 ! Gebietsnummer fuer L3 setzen c Die Anweisungen dieses Abschnitts verlaufen analog zu denen c von Linse 1. -> Fuer Kommentare siehe dort! if (Beschl_Faktor_L3.EQ.0. .OR. q.EQ.0) then ! q=0 -> in Folie neutralisiert d dtmax_L3 = 0. d dtmin_L3 = 0. dt = (xLeaveMap_L3 - x(1)) / v(1) ! Zeit bis zum Mappenende x(1) = xLeaveMap_L3 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_L3) destiny = code_wand goto 5115 endif dt = .5/v(1) zaehler = 0 5015 call INTEGRATIONSSTEP_RUNGE_KUTTA_L3(dt) d if (NTP_steps) then d if (dt.LT.dtmin_L3) then d dtmin_L3 = dt d x_dtmin_L3(1) = x(1) d x_dtmin_L3(2) = x(2) d x_dtmin_L3(3) = x(3) d endif d if (dt.GT.dtmax_L3) then d dtmax_L3 = dt d x_dtmax_L3(1) = x(1) d x_dtmax_L3(2) = x(2) d x_dtmax_L3(3) = x(3) d endif d endif 5115 Steps = Steps + 1 if (destiny.EQ.code_dtSmall) then ! n_dtsmall>maxBelowDtSmall goto 555 elseif (destiny.EQ.code_wand) then ! aufgeschlagen radiusQuad = x(2)*x(2) + x(3)*x(3) ! -> den Ort berechnen, an help1 = v(2)*v(2)+v(3)*v(3) ! dem das Teilchen auf- help2 = x(2)*v(2)+x(3)*v(3) ! schlaegt dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_L3))-help2) + /help1 t = t + dt x(1) = x(1) + dt*v(1) x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? goto 555 c elseif (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'L3-1: ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP elseif (useDecay_) then ! zerfallen? call Decay_Test(*555) endif if (x(1).LT.xEnterMap_L3) then if (v(1).LT.0) then ! reflektiert? destiny = code_reflektiert goto 555 else ! darf nicht sein! write(*,*) write(*,*)' L3: x(1).LT.xEnterMap .AND. v(1).GE.0. -> STOP' write(*,*) STOP endif elseif (Steps.GE.MaxStep) then ! Teilchen verloren destiny = code_lost goto 555 elseif (x(1).GE.xLeaveMap_L3) then ! Verlasse L3 dt = (xLeaveMap_L3 - x(1))/v(1) ! rechne zurueck auf exaktes t = t + dt ! Mappenende x(1) = xLeaveMap_L3 x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug goto bis_MCP2_Mappe c elseif (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'L3-2: ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP 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 5015 ! naechster Integrationsschritt in gleicher Feldmappe c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c zwischen L3-Mappe und MCP2-Mappe (feldfrei) 16 Gebiet = upToM2Map if (x(1).EQ.xEnterMap_M2) goto MCP2_Mappe Steps = Steps + 1 dt = (xEnterMap_M2 - x(1)) / v(1) t = t + dt x(1) = xEnterMap_M2 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then help1 = v(2)*v(2)+v(3)*v(3) help2 = x(2)*v(2)+x(3)*v(3) dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /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 (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_wand goto 555 endif if (useDecay_) call Decay_Test(*555) if (Graphics_) call Save_Graphics_Koord if (Debug_) call Output_Debug c::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: c vor MCP2: (inhom. Felder -> Integrationen) 17 Gebiet = upToMCP2 c Beruecksichtige gegebenenfalls die Flugzeit in 'delta_L2', welches 'vor dem c MCP2' eingeschoben werden kann. Addiert wird vorerst nur die Flugzeit in c dieser zusaetzlichen Flugstrecke. Korrekterweise muessten alle nachfolgenden c Positionen um 'delta_L2' verschoben werden. Dies zu implementieren ist c allerdings momentan aus Zeitgruenden nicht moeglich. dt = Delta_L2 / v(1) t = t + dt c Die Anweisungen dieses Abschnitts verlaufen analog zu denen c von Linse 1. -> Fuer Kommentare siehe dort! if (Beschl_Faktor_M2.EQ.0. .OR. q.EQ.0) then ! q=0 -> in Folie neutralisiert d dtmax_M2 = 0. d dtmin_M2 = 0. if (xBlende.GT.x(1)) then dt = (xBlende - x(1)) / v(1) ! Zeit bis zur Blende x(1) = xBlende x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) then destiny = code_wand elseif (radiusQuad.GE.radiusQuad_Blende) then if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? destiny = code_hitBlende goto 555 endif endif dt = (xMCP2 - x(1)) / v(1) ! Zeit bis MCP2 x(1) = xMCP2 x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GT.radiusQuad_Rohr) destiny = code_wand reachedEndOfMap = .true. goto 5117 endif dt = .5/v(1) reachedEndOfMap = .false. zaehler = 0 if (xBlende.GT.0) check_Blende = .true. 5017 call INTEGRATIONSSTEP_RUNGE_KUTTA_M2(dt) d if (NTP_steps) then d if (dt.LT.dtmin_M2) then d dtmin_M2 = dt d x_dtmin_M2(1) = x(1) d x_dtmin_M2(2) = x(2) d x_dtmin_M2(3) = x(3) d endif d if (dt.GT.dtmax_M2) then d dtmax_M2 = dt d x_dtmax_M2(1) = x(1) d x_dtmax_M2(2) = x(2) d x_dtmax_M2(3) = x(3) d endif d endif 5117 Steps = Steps + 1 if (destiny.EQ.code_dtSmall) then ! n_dtsmall>maxBelowDtSmall goto 555 elseif (check_Blende.AND.x(1).GE.xBlende) then dt = (xBlende - x(1)) / v(1) ! zurueck zur Blende ... x(1) = xBlende x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt radiusQuad = x(2)*x(2) + x(3)*x(3) if (radiusQuad.GE.radiusQuad_Blende) then destiny = code_hitBlende goto 555 endif dt = -dt ! ... wieder zum aktuellen Ort x(1) = xBlende + v(1)*dt x(2) = x(2) + v(2)*dt x(3) = x(3) + v(3)*dt t = t + dt check_Blende = .false. elseif (destiny.EQ.code_wand) then radiusQuad = x(2)*x(2) + x(3)*x(3) ! -> den Ort berechnen, an help1 = v(2)*v(2)+v(3)*v(3) ! dem das Teilchen auf- help2 = x(2)*v(2)+x(3)*v(3) ! schlaegt dt = (SQRT(help2*help2-help1*(radiusQuad-radiusQuad_Rohr))-help2) + /help1 t = t + dt x(1) = x(1) + dt*v(1) x(2) = x(2) + dt*v(2) x(3) = x(3) + dt*v(3) if (useDecay_) call Decay_Test(*555) ! vor Aufschlag zerfallen? goto 555 elseif (useDecay_) then ! zerfallen? call Decay_Test(*555) endif if (destiny.EQ.code_reflektiert) then ! reflektiert goto 555 elseif (reachedEndOfMap) then ! MCP2-Ebene erreicht c if (destiny.NE.code_ok) then ! voruebergehend fuer Testzwecke c write(*,*) c write(*,*)'M2 ''destiny.NE.code_ok'' where it should -> STOP' c write(*,*)' destiny = ',destiny,': ',code_text(destiny) c write(*,*) c STOP c endif if (createNTP.AND.xM2_triggered) fill_NTP = .true. S1xM2 = t if (statNeeded(Nr_S1xM2))call fill_statMem(S1xM2,Nr_S1xM2) radiusQuad = x(2)*x(2) + x(3)*x(3) radius = SQRT(radiusQuad) if (statNeeded(Nr_y_xM2)) call fill_statMem( x(2),Nr_y_xM2) if (statNeeded(Nr_z_xM2)) call fill_statMem( x(3),Nr_z_xM2) if (statNeeded(Nr_r_xM2)) call fill_statMem(radius,Nr_r_xM2) if (radiusQuad.LE.radiusQuad_MCP2active) then S1M2 = t ! Zeiten werden sowohl fuer Statistiken FoM2 = t - S1Fo ! als auch fuer NTupel benoetigt if (statNeeded(Nr_S1M2)) call fill_statMem(S1M2,Nr_S1M2) if (statNeeded(Nr_FoM2)) call fill_statMem(FoM2,Nr_FoM2) if (createNTP.AND.M2_triggered) fill_NTP = .true. if (statNeeded(Nr_y_M2)) call fill_statMem( x(2),Nr_y_M2) if (statNeeded(Nr_z_M2)) call fill_statMem( x(3),Nr_z_M2) if (statNeeded(Nr_r_M2)) call fill_statMem(radius,Nr_r_M2) else ! am MCP2 vorbei if (radiusQuad.LE.radiusQuad_MCP2) then destiny = code_hitMCP2inactive else destiny = code_vorbei if (Graphics_) then ! Damit diese Trajektorie 40mm ueber die nKoord = nKoord + 1 ! MCP2-Ebene hinausgezeichnet wird dt = 40./v(1) t = t + dt xKoord(nKoord) = x(1)+40. yKoord(nKoord) = x(2)+v(2)*dt zKoord(nKoord) = x(3)+v(3)*dt goto 556 endif endif endif goto 555 elseif (Steps.GE.MaxStep) then ! Teilchen verloren destiny = code_lost goto 555 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 5017 ! naechster Integrationsschritt im Feld vor MCP2 czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c HIER IST DER PROGRAMMKODE FUER DIE TRAJEKTORIENBERECHNUNG c DER PROJEKTILE BEENDET! czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz 555 if (Graphics_) call Save_Graphics_Koord 556 if (Debug_) call Output_Debug c Gib Trajektorie in Graphikfenster aus: if (Graphics_) then if (Gebiet.LE.upToChKoord) then ! Bahnberechnung wurde vor call plot_horizontal ! Koordinatenwechsel abgebrochen if (schnitt_p.eq.1) call schnitt else call plot_vertikal if (schnitt_p.eq.2) call schnitt endif nKoord = 0 endif c Pruefe, ob Teilchen reflektiert wurde: if ((Gebiet.EQ.upToExL1 .OR. Gebiet.EQ.upToEnTD .OR. + Gebiet.EQ.upToExL3 .OR. Gebiet.EQ.upToMCP2) .AND. + v(1).LE.0.) then destiny = code_reflektiert endif c Zaehle mit, bei wie vielen Teilchen trotz dtMaxStep abgebrochen werden: if (destiny.EQ.code_lostInTD) then lostInTD_counter = lostInTD_counter + 1 elseif (destiny.EQ.code_lost) then lost_counter = lost_counter + 1 endif c bei DEBUG: Ausgabe des Teilchenschicksals und des aktuellen Gebiets: if (debug_) then indx = index(code_text(destiny),':') if (indx.EQ.0) then write(lun(1),*) 'destiny : ',code_text(destiny) else write(lun(1),*) 'destiny : ',code_text(destiny)(1:indx-1) endif indx = index(Gebiet_text(Gebiet),':') if (indx.EQ.0) then write(lun(1),*) 'Gebiet : ',Gebiet_text(Gebiet) else write(lun(1),*) 'Gebiet : ',Gebiet_text(Gebiet)(1:indx-1) endif endif czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c HIER STARTEN DIE FOLIENELEKTRONEN czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz if (generate_FE) then ! ~~~ 1: if ~~~~~~~~~~~~ if (Gebiet.GE.UpToExTD) then ! ~~~ 2: if ~~~~~~~~~~~~ c sekundaerelektronen nFE = int(4.*ran(seed))+2 ! Anzahl wuerfeln: [2,5] tFE_min = 0. ! tFE_min: kuerzeste FE-Flugzeit: ! bekam noch keinen Wert zugewiesen c - die Laufzeiten der Folienelektronen: c c---------------------------------------- c c-TP-10/2000 reset of end positions of electrons c do k = 1, 3 xFE_MCP(k) = 0. yFE_MCP(k) = 0. zFE_MCP(k) = 0. enddo c c---------------------------------------- c do 450, k = 1, nFE xFE(1) = x0FE(1) xFE(2) = x0FE(2) xFE(3) = x0FE(3) E0FE = 0.003*ran(seed) ! Start-Energie wuerfeln: [0,3) eV v_Betrag = sqrt(2.*E0FE/511.015)*c ! Startgeschwindigkeit call Lambert_Verteilung(1.,ct0,st0) ! Startwinkel wuerfeln f0 = 360.*ran(seed) cf0 = cosd(f0) sf0 = sind(f0) vFE(1) = v_Betrag * ct0 ! Geschwindigkeitskomponenten vFE(2) = v_Betrag * st0*cf0 vFE(3) = v_Betrag * st0*sf0 tFE = 0. nKoord = 0 start_nr(2) = start_nr(2) + 1 ! (2): FE call TRIGGER(511.015,-1.,tFE,xFE,vFE,DEBUG_FE.AND.Debug_, + plot_FE,n_return) if (plot_FE) call plot_vertikal if (n_return.NE.-10) then C - das FE kam nicht am MCP3 an -> if (n_return.GT.100 .AND. n_return.LE.120) then ! -> abgebrochen statTD(2,n_return-100) = statTD(2,n_return-100)+1 ! Grund notieren elseif (n_return.GT.0 .AND. n_return.LE.75) then ! -> pfosten getroffen! pfostenHit(n_return,2) = pfostenHit(n_return,2)+1 elseif (n_return.EQ.0) then ! -> GROUND getroffen statTD(2,15) = statTD(2,15)+1 ! '15' zaehlt GROUND-Treffer elseif (n_return.EQ.-5) then ! -> im TD auf Gitterstab statTD(2,17) = statTD(2,17)+1 elseif (n_return.EQ.-9) then ! -> NICHT im MCP3 registriert statTD(2,18) = statTD(2,18)+1 endif tFE_(k) = -1 ! -1: FE kam nicht bei MCP3 an c c--------------------------------------- c c-TP-10/2000 c xFE_MCP(k) = -100. yFE_MCP(k) = -100. zFE_MCP(k) = -100. c c--------------------------------------- c goto 450 ! naechstes FE endif c - das FE kam beim MCP3 an -> statTD(2,16) = statTD(2,16)+1 ! '16' zaehlt MCP3-Treffer tFE_(k)=int(1000.*tFE) ! tFE in ps. (braucht als Integer ! weniger Speicherplatz) c c--------------------------------------- c c-TP-10/2000 c xFE_MCP(k) = xFE(1) yFE_MCP(k) = xFE(2) zFE_MCP(k) = xFE(3) c c--------------------------------------- c c fuer die Statistik: die Flugzeiten saemtlicher das MCP3 erreichender FE abspeichern: if (statNeeded(Nr_t_FE)) call fill_statMem(tFE,Nr_t_FE) c kuerzeste Elektronenflugzeit fuer das aktuelle Projektilteilchen notieren: if (tFE_min.EQ.0. .OR. tFE.LT.tFE_min) tFE_min = tFE 450 continue ! -> naechstes Folienelektron c die Flugzeiten der nicht gestartenen Folienelektronen (nFE+1 bis 5) auf 0. setzen: do while (nFE.LT.5) nFE = nFE + 1 tFE_(nFE) = 0. enddo c Jetzt sind die Folienelektronen durchgelaufen. c Fuelle Statistiken fuer die 'gemessenen' Teilchenflugzeiten (mit Beruecksichti- c gung der Flugzeiten der Folienelektronen). M3M2 aber nur, wenn MCP2 auch c getroffen wurde: if (tFE_min.NE.0.) then S1M3 = S1Fo + tFE_min ! +, da Stop verzoegert if (statNeeded(Nr_S1M3)) then call fill_statMem(S1M3,Nr_S1M3) endif if (destiny.EQ.code_ok) then M3M2 = FoM2 - tFE_min ! -, da Start verzoegert if (statNeeded(Nr_M3M2)) call fill_statMem(M3M2,Nr_M3M2) endif endif else ! ~~~ 2: else ~~~~~~~~~~ do k= 1, 5 tFE_(k) = 0. ! nicht gestartet enddo endif ! ~~~ 2: endif ~~~~~~~~~ endif ! ~~~ 1: endif~~~~~~~~~~ czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c ES FOLGEN DATENAUSGABE, SPRUNG IN NEUE SCHLEIFE UND PROGRAMMENDE czzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzzz c trage das NTupel ein: c So verlangt, schreibe die aktuellen Trajektoriengroessen in das NTupel: c (falls bei 'createFoilFile' 'NTPalreadyWritten' nicht gesetzt ist schied c dieses Teilchen schon vor der Triggerfolie aus. Ist es dagegen gesetzt wurden c die Trajektoriendaten mit dem Erreichen der Triggerfolie abgespeichert um sie c in den im Triggersystem gueltigen Werten zu haben): if (fill_NTP .OR. createFoilFile) then if (NTPalreadyWritten) then NTPalreadyWritten = .false. else if (NTP_stop) then Ekin=(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))*Energie_Faktor endif if (smearS1Fo .AND. .NOT.use_MUTRACK) then if (s1fo.NE.0) then call Gauss_Verteilung(sigmaS1Fo,help4) S1FoOnly = S1Fo + help4 else S1FoOnly = 0. endif endif FoM2Only = FoM2 call HFNT(NTP_write) endif endif c Nimm das Schicksal des Teilchens in den zugehoerigen Statistikspeicher auf: if (destiny.GT.0) destiny = destiny + (Gebiet-1)*highest_code_Nr statDestiny(destiny) = statDestiny(destiny) + 1 if (destiny.EQ.code_ok) okStepsCounter = okStepsCounter + steps c -> das naechste Projektil kann kommen 100 continue c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Jetzt sind alle Projektile dieser Schleife abgearbeitet! c Mittlere Anzahl an Integrationsschritten fuer Trajektorien mit destiny = c 'code_ok' ausgeben: if (statDestiny(code_ok).NE.0) then write(*,'(6x,A,F7.2)')'Mittlere Anzahl an Integrationsschritten bis zum Ziel: ', + real(okStepsCounter)/real(statDestiny(code_ok)) endif 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 call schnitt_plot ! Ausgabe der Graphik der Schnittebene if (n_postSkript.LE.0) then goto 620 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(*,1003)'(P) Ps-file erstellen', + '(R) Restliche ps-files erstellen' write(*,1003)'(N) ps-file Nicht erstellen', + '(K) Keine ps-files mehr erstellen' write(*,1003)'(G) Graphikausgabe beenden', + '(A) programm Abbrechen' 1003 format(T6,A,T40,A) helpChar = 'n' 600 write(*,1004)' [RETURN] = (N) -> ' 1004 format($,x,A) read(*,'(A)') helpChar do i = 1,7 ! bis zu sechs blanks werden akzeptiert ant = helpChar(i:i) if (ant.NE.' ') goto 610 enddo ant = 'N' 610 write(*,*)'==========================='// + '=====================================================' call str$upcase(ant,ant) if (ant.EQ.'N') then goto 620 elseif (ant.EQ.'R') then n_postSkript = 2 elseif (ant.EQ.'K') then n_postSkript = 0 goto 620 elseif (ant.EQ.'G') then call HPLEND GRAPHICS = .false. goto 200 elseif (ant.EQ.'A') then call HPLEND call TERMINATE_OUTPUT STOP elseif (ant.NE.'P') then goto 600 endif endif write (helpChar(1:7),'(''_'',I6)') SchleifenNr if (filename.NE.' ') then call MAKE_PS(filename//helpChar) else call MAKE_PS('MUTRACK'//helpChar) endif 620 continue CALL IZPICT ('CHAM_1','S') ! LOESCHEN DER BILDER IM PAWC-COMMON-BLOCK CALL IZPICT ('CHAM_2','S') CALL IZPICT ('HISTO','S') CALL IZPICT ('TEXT','S') call iclrwk (1,flag_ok) ! CLEAREN DER WORKSTATIONS call iclrwk (3,flag_ok) call iclrwk (4,flag_ok) call iclrwk (5,flag_ok) 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: c - HIGZ-Graphikbibliothek schliessen: if (Graphics) call HPLEND c - HBOOK-Datei schliessen: if (.NOT.fromScratch) then if (use_ACCEL) then call HREND('ACCEL') elseif (Use_MUTRACK) then call HREND('MUread') endif close (lunRead) endif 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 randomVar integer seed common /seed/ seed randomVar = ran(seed) if (n_Lambert.EQ.0.) then cos_theta = (1.-randomVar) sin_theta = sqrt(1.-cos_theta*cos_theta) elseif (n_Lambert.EQ.1.) then cos_theta = sqrt(1.-randomVar) sin_theta = sqrt(randomVar) else cos_theta = (1.-randomVar)**(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 real zwoPi parameter (zwoPi = 2.*3.1415927) 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 = zwoPi * ran(seed) wert = radius * cos(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 1 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 1 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 1 ! Verteilung zurechtbiegen END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE G_Functions(G1,G2,tau) c ================================= c Diese Routine gibt in Abhaengigkeit von der reduzierten Dicke 'tau' c Funktionswerte fuer g1 und g2 zurueck. g1 und g2 sind dabei die von c Meyer angegebenen tabellierten Funktionen fuer die Berechnung von Halbwerts- c breiten von Streuwinkelverteilungen. (L.Meyer, phys.stat.sol. (b) 44, 253 c (1971)) IMPLICIT NONE real tau,g1,g2 real tau_(26),g1_(26),g2_(26) real help integer i DATA tau_ /0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, + 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, + 10.0, 12.0, 14.0, 16.0, 18.0, 20.0 / DATA g1_ /0.050,0.115,0.183,0.245,0.305,0.363,0.419,0.473,0.525,0.575, + 0.689,0.799,0.905,1.010,1.100,1.190,1.370,1.540,1.700,1.850, + 1.990,2.270,2.540,2.800,3.050,3.290 / DATA g2_ / 0.00,1.25,0.91,0.79,0.73,0.69,0.65,0.63,0.61,0.59, + 0.56,0.53,0.50,0.47,0.45,0.43,0.40,0.37,0.34,0.32, + 0.30,0.26,0.22,0.18,0.15,0.13 / if (tau.LT.tau_(1)) then write(*,*) write(*,*)'SUBROUTINE G_Functions:' write(*,*)' Fehler bei Berechnung der g-Funktionen fuer Winkelaufstreuung:' write(*,*)' aktuelles tau ist kleiner als kleinster Tabellenwert:' write(*,*)' tau = ',tau write(*,*)' tau_(1) = ',tau_(1) write(*,*) STOP endif i = 1 10 i = i + 1 if (i.EQ.27) then write(*,*) write(*,*)'SUBROUTINE G_Functions:' write(*,*)' Fehler bei Berechnung der g-Funktionen fuer Winkelaufstreuung:' write(*,*)' aktuelles tau ist groesser als groesster Tabellenwert:' write(*,*)' tau = ',tau write(*,*)' tau_(26) = ',tau_(26) write(*,*) STOP elseif (tau.gt.tau_(i)) then goto 10 endif c lineare Interpolation zwischen Tabellenwerten: help = (tau-tau_(i-1))/(tau_(i)-tau_(i-1)) g1 = g1_(i-1) + help*(g1_(i)-g1_(i-1)) g2 = g2_(i-1) + help*(g2_(i)-g2_(i-1)) END c=============================================================================== options /extend_source subroutine Get_F_Function_Meyer(tau,Ekin) c ========================================= implicit none real tau real Ekin real thetaSchlange,thetaSchlangeMax real theta,thetaMax,thetaStep real f1,f2,F c------------------------------------ c - Parameter: real Z1, Z2 ! die atomaren Nummern von Projektil und Target c real a0 ! Bohrscher Radius in cm real screeningPar ! Screeningparameter 'a' in cm fuer Teilchen der ! Kernladungszahl Z1=1 in Kohlenstoff (Z2 = 6) ! bei Streichung von Z1 (vgl. Referenz, S. 268) real r0Meyer ! r0(C) berechnet aus dem screeningParameter 'a' ! und dem ebenfalls bei Meyer angegebenem ! Verhaeltnis a/r0=0.26 (vgl. Referenz, S. 263 oben) real eSquare ! elektrische Ladung zum Quadrat in keV*cm real Pi ! die Kreiszahl c parameter (a0 = 5.29E-9) parameter (Z1 = 1, Z2 = 6, ScreeningPar = 2.5764E-9) parameter (r0Meyer = 9.909E-9, eSquare = 1.44E-10) parameter (Pi = 3.141592654) real Meyer_Faktor3 real Meyer_Faktor4 real zzz ! 'Hilfsparameter' real Meyer_Faktor5 parameter (Meyer_faktor3 = (screeningPar/r0Meyer) * (screeningPar/r0Meyer)) parameter (Meyer_faktor4 = screeningPar / (2.*Z1*Z2*eSquare) * Pi/180.) parameter (zzz = screeningPar / (2.*Z1*Z2*eSquare)) parameter (Meyer_faktor5 = zzz*zzz / (8*Pi*Pi)) c------------------------------------ integer nBin,nBinMax parameter (nBinMax=201) real value(0:nBinMax) /0.,nBinMax*0./ real area(nBinMax) / nBinMax*0./ real integ(0:nBinMax) /0.,nBinMax*0./ common /MeyerTable/ value,area,integ,thetaStep,nBin integer i real rhelp integer HB_memsize parameter(HB_memsize=500000) real memory(HB_memsize) COMMON /PAWC/ memory c nur noch fuer Testzwecke: real fValues(203) real fValuesFolded(203) integer idh parameter (idh = 50) INCLUDE 'mutrack$sourcedirectory:COM_DIRS.INC' character filename*20 ! Name der Ausgabe-Dateien COMMON /filename/ filename c------------------------------------------------------------------------------- c Festlegen des maximalen Theta-Wertes sowie der Schrittweite: if (tau.LT.0.2) then write(*,*) 'Subroutine ''Get_F_Function_Meyer'':' write(*,*) 'Effektive Dicke ist kleiner als 0.2 => kann ich nicht ... => STOP' call exit elseif (tau.LE.2.) then ! => Tabelle A thetaSchlangeMax = 4.0 elseif (tau.LE.8.) then ! => Tabelle B thetaSchlangeMax = 7.0 elseif (tau.LE.20.) then ! => Tabelle C thetaSchlangeMax = 20.0 else write(*,*) 'Subroutine ''Get_F_Function_Meyer'':' write(*,*) 'Effektive Dicke ist groesser als 20 => kann ich nicht ... => STOP' call exit endif thetaMax = thetaSchlangeMax / Meyer_Faktor4 / Ekin if (thetaMax.GT.50) then thetaStep = .5 elseif (thetaMax.GT.25) then thetaStep = .25 elseif (thetaMax.GT.12.5) then thetaStep = .125 else thetaStep = .0625 endif c Tabelle der F-Werte erstellen: nBin = 0 do theta = thetaStep, thetaMax, thetaStep ! Berechne aus theta das 'reduzierte' thetaSchlange (dabei gleich ! noch von degree bei theta in Radiant bei thetaSchlange umrechnen): thetaSchlange = Meyer_faktor4 * Ekin * theta ! Auslesen der Tabellenwerte fuer die f-Funktionen: call F_Functions_Meyer(tau,thetaSchlange,f1,f2) if (thetaSchlange.EQ.-1) then ! wir sind jenseits von thetaSchlangeMax goto 10 endif ! Berechnen der Streuintensitaet: F = Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2) nBin = nBin + 1 if (nBin.GT.nBinMax) then write(*,*) 'nBin > nBinMax => EXIT' call exit endif value(nBin) = sind(theta)*F fValues(nBin+1) = F ! fuer Testzwecke fValuesFolded(nBin+1) = sind(theta)*F ! fuer Testzwecke enddo c Berechnen der Flaecheninhalte der einzelnen Kanaele sowie der Integrale: 10 do i = 1, nBin area(i) = (value(i)+value(i-1))/2. * thetaStep integ(i) = integ(i-1) + area(i) enddo c Normiere totale Flaeche auf 1: rHelp = integ(nBin) do i = 1, nBin value(i) = value(i) / rHelp area(i) = area(i) / rHelp integ(i) = integ(i) / rHelp enddo c vorerst noch: gib Tabelle in Datei und Histogrammfile aus: ! Berechne die Werte fuer theta=0: call F_Functions_Meyer(tau,0.,f1,f2) F = Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2) fValues(1) = F fValuesFolded(1) = 0. ! Gib die Werte in das Tabellenfile aus: c theta = 0. c open (10,file=outDir//':'//filename//'.TAB',status='NEW') c do i = 1, nBin+1 c write(10,*) theta, fValues(i), fValuesFolded(i) c theta = theta + thetaStep c enddo c close (10) ! Buchen und Fuellen der Histogramme: call HBOOK1(idh,'F',nBin+1,-0.5*thetaStep,(real(nBin)+0.5)*thetaStep,0.) call HPAK(idh,fValues) call HRPUT(idh,outDir//':'//filename//'.RZ','N') call HDELET(idh) call HBOOK1(idh+1,'F*sin([q])',nBin+1,-0.5*thetaStep,(real(nBin)+0.5)*thetaStep,0.) call HPAK(idh+1,fValuesFolded) call HRPUT(idh+1,outDir//':'//filename//'.RZ','U') call HDELET(idh+1) END c=============================================================================== options /extend_source subroutine throwMeyerAngle (theta) c ================================== implicit none real lowerbound,y1,y2,f,root,radiant,fraction integer bin,nBin integer nBinMax parameter (nBinMax=201) real theta,thetaStep real value(0:nBinMax) /0.,nBinMax*0./ real area(nBinMax) / nBinMax*0./ real integ(0:nBinMax) /0.,nBinMax*0./ common /MeyerTable/ value,area,integ,thetaStep,nBin real rhelp real random integer seed common /seed/ seed c bin: Nummer des Bins, innerhalb dessen das Integral den Wert von c random erreicht oder ueberschreitet: random = ran(seed) bin = 1 do while (random.GT.integ(bin)) bin = bin + 1 if (bin.GT.nBin) then write(*,*) 'error 1' call exit endif enddo fraction = (random-integ(bin-1)) / (integ(bin)-integ(bin-1)) y1 = value(bin-1) y2 = value(bin) f = thetaStep / (y2-y1) rHelp = y1*f radiant = rHelp*rHelp + fraction*thetaStep*(y1+y2)*f root = SQRT(radiant) lowerBound = real(bin-1)*thetaStep if (f.GT.0) then theta = lowerBound - rHelp + root else theta = lowerBound - rHelp - root endif END c=============================================================================== options /extend_source subroutine F_Functions_Meyer(tau,thetaSchlange,f1,f2) c ===================================================== implicit none c Diese Routine gibt in Abhaengigkeit von 'thetaSchlange' und 'tau' c Funktionswerte fuer f1 und f2 zurueck. f1 und f2 entsprechen dabei den c bei Meyer angegebenen Funktion gleichen Namens. Die in dieser Routine c verwendeten Tabellen sind eben dieser Referenz entnommen: c L.Meyer, phys.stat.sol. (b) 44, 253 (1971) real tau,thetaSchlange real f1, f2, f1_(2), f2_(2) integer column_,column,row integer iColumn real weightCol, weightRow c------------------------------------------------------------------------------- c die Tabellendaten der Referenz (Tabellen 2 und 3): integer nColumn parameter (nColumn = 25) real tau_(nColumn) / + 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, + 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10., 12., 14., 16., 18., 20. / integer nRowA parameter (nRowA = 25) real thetaSchlangeA(nRowA) / + .00, .05, .10, .15, .20, .25, .30, .35, .40, .45, .50, .60, + .70, .80, .90, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0 / integer nRowB parameter (nRowB = 24) real thetaSchlangeB(nRowB) / + 0.0, 0.2, 0.4, 0.5, 0.6, 0.8, 1.0, 1.2, 1.4, 1.5, 1.6, 1.8, + 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0 / integer nRowC parameter (nRowC = 24) real thetaSchlangeC(nRowC) / + 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, + 7.0, 8.0, 9.0, 10., 11., 12., 13., 14., 15., 16., 18., 20. / real f1_A(9,nRowA) + /1.69E+2,4.55E+1,2.11E+1,1.25E+1,8.48E+0,6.21E+0,4.80E+0,3.86E+0,3.20E+0, + 9.82E+1,3.72E+1,1.97E+1,1.20E+1,8.27E+0,6.11E+0,4.74E+0,3.83E+0,3.17E+0, + 3.96E+1,2.58E+1,1.65E+1,1.09E+1,7.73E+0,5.82E+0,4.58E+0,3.72E+0,3.10E+0, + 1.76E+1,1.58E+1,1.27E+1,9.26E+0,6.93E+0,5.38E+0,4.31E+0,3.55E+0,2.99E+0, + 8.62E+0,1.01E+1,9.45E+0,7.58E+0,6.02E+0,4.85E+0,3.98E+0,3.33E+0,2.84E+0, + 4.65E+0,6.55E+0,6.91E+0,6.06E+0,5.11E+0,4.28E+0,3.62E+0,3.08E+0,2.66E+0, + 2.74E+0,4.45E+0,5.03E+0,4.78E+0,4.27E+0,3.72E+0,3.23E+0,2.82E+0,2.47E+0, + 1.77E+0,3.02E+0,3.71E+0,3.76E+0,3.53E+0,3.20E+0,2.86E+0,2.55E+0,2.27E+0, + 1.22E+0,2.19E+0,2.78E+0,2.96E+0,2.91E+0,2.73E+0,2.51E+0,2.28E+0,2.07E+0, + 8.82E-1,1.59E+0,2.12E+0,2.35E+0,2.39E+0,2.32E+0,2.19E+0,2.03E+0,1.87E+0, + 6.55E-1,1.20E+0,1.64E+0,1.88E+0,1.97E+0,1.96E+0,1.90E+0,1.79E+0,1.68E+0, + 3.80E-1,7.15E-1,1.01E+0,1.22E+0,1.35E+0,1.40E+0,1.41E+0,1.39E+0,1.34E+0, + 2.26E-1,4.45E-1,6.44E-1,8.08E-1,9.28E-1,1.01E+0,1.05E+0,1.06E+0,1.05E+0, + 1.39E-1,2.80E-1,4.21E-1,5.45E-1,6.46E-1,7.22E-1,7.75E-1,8.07E-1,8.21E-1, + 8.22E-2,1.76E-1,2.78E-1,3.71E-1,4.53E-1,5.21E-1,5.74E-1,6.12E-1,6.37E-1, + 5.04E-2,1.11E-1,1.86E-1,2.57E-1,3.22E-1,3.79E-1,4.27E-1,4.65E-1,4.94E-1, + 2.51E-2,5.60E-2,9.24E-2,1.31E-1,1.69E-1,2.02E-1,2.40E-1,2.71E-1,2.97E-1, + 1.52E-2,3.20E-2,5.08E-2,7.23E-2,9.51E-2,1.18E-1,1.41E-1,1.63E-1,1.83E-1, + 1.03E-2,2.05E-2,3.22E-2,4.55E-2,6.01E-2,7.53E-2,9.02E-2,1.05E-1,1.19E-1, + 8.80E-3,1.48E-2,2.25E-2,3.13E-2,4.01E-2,5.03E-2,6.01E-2,7.01E-2,8.01E-2, + 6.10E-3,1.15E-2,1.71E-2,2.28E-2,2.89E-2,3.52E-2,4.18E-2,4.86E-2,5.55E-2, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,1.71E-2,1.98E-2,2.28E-2,2.58E-2, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,8.90E-3,1.02E-2,1.16E-2,1.31E-2, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,4.90E-3,5.70E-3,6.40E-3,7.20E-3, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,2.90E-3,3.40E-3,3.90E-3,4.30E-3/ real f1_B(9,nRowB) + /2.71E+0,1.92E+0,1.46E+0,1.16E+0,9.52E-1,8.03E-1,6.90E-1,5.32E-1,4.28E-1, + 2.45E+0,1.79E+0,1.39E+0,1.12E+0,9.23E-1,7.82E-1,6.75E-1,5.23E-1,4.23E-1, + 1.87E+0,1.48E+0,1.20E+0,9.96E-1,8.42E-1,7.24E-1,6.32E-1,4.98E-1,4.07E-1, + 1.56E+0,1.30E+0,1.09E+0,9.19E-1,7.89E-1,6.86E-1,6.03E-1,4.80E-1,3.95E-1, + 1.28E+0,1.11E+0,9.62E-1,8.33E-1,7.27E-1,6.40E-1,5.69E-1,4.59E-1,3.81E-1, + 8.23E-1,7.90E-1,7.29E-1,6.64E-1,6.01E-1,5.44E-1,4.94E-1,4.12E-1,3.49E-1, + 5.14E-1,5.36E-1,5.29E-1,5.07E-1,4.78E-1,4.47E-1,4.16E-1,3.60E-1,3.13E-1, + 3.19E-1,3.58E-1,3.76E-1,3.78E-1,3.70E-1,3.57E-1,3.45E-1,3.08E-1,2.76E-1, + 2.02E-1,2.40E-1,2.64E-1,2.77E-1,2.82E-1,2.80E-1,2.65E-1,2.59E-1,2.39E-1, + 1.67E-1,1.96E-1,2.20E-1,2.36E-1,2.44E-1,2.47E-1,2.45E-1,2.35E-1,2.21E-1, + 1.33E-1,1.61E-1,1.85E-1,2.02E-1,2.12E-1,2.18E-1,2.18E-1,2.14E-1,2.03E-1, + 8.99E-2,1.12E-1,1.32E-1,1.48E-1,1.59E-1,1.67E-1,1.68E-1,1.75E-1,1.72E-1, + 6.24E-2,7.94E-2,9.50E-2,1.09E-1,1.20E-1,1.29E-1,1.35E-1,1.42E-1,1.43E-1, + 4.55E-2,5.74E-2,6.98E-2,8.11E-2,9.09E-2,9.92E-2,1.06E-1,1.15E-1,1.19E-1, + 3.35E-2,4.22E-2,5.19E-2,6.11E-2,6.95E-2,7.69E-2,8.33E-2,9.28E-2,9.85E-2, + 2.50E-2,3.16E-2,3.92E-2,4.66E-2,5.35E-2,6.00E-2,6.57E-2,7.49E-2,8.13E-2, + 1.90E-2,2.40E-2,2.99E-2,3.58E-2,4.16E-2,4.70E-2,5.20E-2,6.05E-2,6.70E-2, + 1.47E-2,1.86E-2,2.32E-2,2.79E-2,3.25E-2,3.70E-2,4.12E-2,4.89E-2,5.51E-2, + 8.10E-3,1.04E-2,1.30E-2,1.57E-2,1.84E-2,2.12E-2,2.40E-2,2.93E-2,3.42E-2, + 4.80E-3,6.20E-3,7.70E-3,9.30E-3,1.09E-2,1.26E-2,1.44E-2,1.79E-2,2.14E-2, + 2.80E-3,3.80E-3,4.70E-3,5.70E-3,6.70E-3,7.50E-3,8.90E-3,1.13E-2,1.36E-2, + 1.70E-3,2.30E-3,2.90E-3,3.60E-3,4.20E-3,4.90E-3,5.60E-3,7.20E-3,8.80E-3, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,2.00E-3,2.80E-3,3.50E-3, + 0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,0.00 ,8.80E-4,1.20E-3,1.60E-3/ real f1_C(7,nRowC) + /3.65E-1,2.62E-1,2.05E-1,1.67E-1,1.41E-1,1.21E-1,1.05E-1, + 3.33E-1,2.50E-1,1.95E-1,1.61E-1,1.36E-1,1.18E-1,1.03E-1, + 2.75E-1,2.18E-1,1.76E-1,1.48E-1,1.27E-1,1.11E-1,9.80E-2, + 2.04E-1,1.75E-1,1.50E-1,1.29E-1,1.13E-1,1.01E-1,9.00E-2, + 1.41E-1,1.31E-1,1.19E-1,1.08E-1,9.71E-2,8.88E-2,8.01E-2, + 9.32E-2,9.42E-2,9.10E-2,8.75E-2,8.00E-2,7.44E-2,6.91E-2, + 5.98E-2,6.52E-2,6.72E-2,6.62E-2,6.40E-2,6.12E-2,5.82E-2, + 3.83E-2,4.45E-2,4.80E-2,4.96E-2,4.98E-2,4.90E-2,4.77E-2, + 2.46E-2,3.01E-2,3.40E-2,3.65E-2,3.79E-2,3.84E-2,3.83E-2, + 1.59E-2,2.03E-2,2.39E-2,2.66E-2,2.85E-2,2.97E-2,3.04E-2, + 1.04E-2,1.37E-2,1.66E-2,1.92E-2,2.12E-2,2.27E-2,2.37E-2, + 4.39E-3,6.26E-3,8.26E-3,9.96E-3,1.15E-2,1.29E-2,1.41E-2, + 2.06E-3,3.02E-3,4.24E-3,5.28E-3,6.32E-3,7.32E-3,8.26E-3, + 1.21E-3,1.69E-3,2.24E-3,2.85E-3,3.50E-3,4.16E-3,4.82E-3, + 8.50E-4,1.10E-3,1.38E-3,1.65E-3,2.03E-3,2.45E-3,2.88E-3, + 5.90E-4,7.40E-4,8.50E-4,9.90E-4,1.23E-3,1.49E-3,1.71E-3, + 3.90E-4,4.60E-4,5.20E-4,6.30E-4,7.65E-4,9.65E-4,1.12E-3, + 2.40E-4,2.70E-4,3.10E-4,3.98E-4,4.97E-4,6.03E-4,7.18E-4, + 1.50E-4,1.70E-4,2.15E-4,2.70E-4,3.35E-4,4.35E-4,5.00E-4, + 1.00E-4,1.20E-4,1.46E-4,1.90E-4,2.40E-4,2.88E-4,3.43E-4, + 0.00 ,0.00 ,1.04E-4,1.41E-4,1.80E-4,2.10E-4,2.50E-4, + 0.00 ,0.00 ,8.20E-5,1.06E-4,1.38E-4,1.58E-4,1.85E-4, + 0.00 ,0.00 ,5.40E-5,7.00E-5,8.60E-5,1.03E-4,1.20E-4, + 0.00 ,0.00 ,4.20E-5,5.40E-5,6.50E-5,7.70E-5,8.80E-5/ real f2_A(9,nRowA) + / 3.52E+3, 3.27E+2, 9.08E+1, 3.85E+1, 2.00E+1, 1.18E+1, 7.55E+0, 5.16E+0, 3.71E+0, + 2.58E+2, 1.63E+2, 7.30E+1, 3.42E+1, 1.85E+1, 1.11E+1, 7.18E+0, 4.96E+0, 3.59E+0, + -1.12E+2, 4.84E+0, 3.56E+1, 2.34E+1, 1.45E+1, 9.33E+0, 6.37E+0, 4.51E+0, 3.32E+0, + -5.60E+1,-1.12E+1, 9.87E+0, 1.24E+1, 9.59E+0, 7.01E+0, 5.16E+0, 3.83E+0, 2.91E+0, + -2.13E+1,-1.22E+1,-2.23E+0, 3.88E+0, 5.15E+0, 4.65E+0, 3.87E+0, 3.12E+0, 2.45E+0, + -8.25E+0,-9.58E+0,-5.59E+0,-1.40E+0, 1.76E+0, 2.71E+0, 2.71E+0, 2.35E+0, 1.95E+0, + -3.22E+0,-6.12E+0,-5.28E+0,-2.87E+0,-1.92E-1, 1.32E+0, 1.69E+0, 1.74E+0, 1.48E+0, + -1.11E+0,-3.40E+0,-4.12E+0,-3.08E+0,-6.30E-1, 3.60E-1, 9.20E-1, 1.03E+0, 1.04E+0, + -2.27E-1,-2.00E+0,-2.93E+0,-2.69E+0,-1.48E+0,-3.14E-1, 2.69E-1, 5.28E-1, 6.09E-1, + 1.54E-1,-1.09E+0,-2.10E+0,-2.15E+0,-1.47E+0,-6.77E-1,-1.80E-1, 1.08E-1, 2.70E-1, + 3.28E-1,-6.30E-1,-1.50E+0,-1.68E+0,-1.34E+0,-8.43E-1,-4.60E-1,-1.85E-1,-4.67E-3, + 3.32E-1,-2.06E-1,-7.32E-1,-9.90E-1,-9.42E-1,-8.20E-1,-6.06E-1,-4.51E-1,-3.01E-1, + 2.72E-1,-3.34E-2,-3.49E-1,-5.65E-1,-6.03E-1,-5.79E-1,-5.05E-1,-4.31E-1,-3.45E-1, + 2.02E-1, 2.80E-2,-1.54E-1,-3.00E-1,-3.59E-1,-3.76E-1,-4.60E-1,-3.40E-1,-3.08E-1, + 1.38E-1, 4.84E-2,-5.56E-2,-1.44E-1,-2.04E-1,-2.39E-1,-2.54E-1,-2.49E-1,-2.48E-1, + 9.47E-2, 4.86E-2,-1.08E-2,-6.44E-2,-1.02E-1,-1.34E-1,-1.62E-1,-1.79E-1,-1.87E-1, + 5.33E-2, 3.71E-2, 1.85E-2, 1.63E-3,-1.69E-2,-3.69E-2,-5.66E-2,-7.78E-2,-9.33E-2, + 3.38E-2, 2.40E-2, 1.62E-2, 9.90E-3, 3.76E-3,-4.93E-3,-1.66E-2,-3.05E-2,-4.22E-2, + 2.12E-2, 1.56E-2, 1.05E-2, 7.80E-3, 7.92E-3, 6.30E-3, 3.20E-4,-8.50E-3,-1.66E-2, + 1.40E-2, 9.20E-3, 5.30E-3, 4.70E-3, 6.31E-3, 8.40E-3, 5.30E-3, 8.80E-4,-3.30E-3, + 9.20E-3, 4.70E-3, 1.70E-3, 2.60E-3, 4.49E-3, 6.60E-3, 6.00E-3, 4.70E-3, 2.80E-3, + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 / real f2_B(9,nRowB) + / 2.75E+0, 1.94E+0, 9.13E-1, 6.06E-1, 4.26E-1, 3.14E-1, 2.40E-1, 1.51E-1, 1.03E-1, + 1.94E+0, 1.16E+0, 7.56E-1, 5.26E-1, 3.81E-1, 2.87E-1, 2.23E-1, 1.43E-1, 9.78E-2, + 5.85E-1, 5.04E-1, 4.10E-1, 3.30E-1, 2.69E-1, 2.17E-1, 1.78E-1, 1.22E-1, 8.71E-2, + 7.83E-2, 2.00E-1, 2.35E-1, 2.19E-1, 1.97E-1, 1.73E-1, 1.48E-1, 1.08E-1, 7.93E-2, + -1.82E-1, 1.56E-2, 1.04E-1, 1.36E-1, 1.38E-1, 1.31E-1, 1.19E-1, 9.46E-2, 7.19E-2, + -2.71E-1,-1.66E-1,-7.29E-2,-4.74E-3, 3.60E-2, 5.50E-2, 6.28E-2, 5.98E-2, 5.09E-2, + -1.87E-1,-1.58E-1,-1.09E-1,-5.80E-2,-2.03E-2, 2.48E-3, 1.99E-2, 3.36E-2, 3.27E-2, + -1.01E-1,-1.05E-1,-8.95E-2,-6.63E-2,-3.93E-2,-2.38E-2,-9.22E-3, 8.47E-3, 1.52E-2, + -5.19E-2,-6.47E-2,-6.51E-2,-5.62E-2,-4.51E-2,-3.49E-2,-2.45E-2,-8.19E-3, 2.05E-3, + -3.68E-2,-4.89E-2,-5.36E-2,-5.06E-2,-4.27E-2,-3.65E-2,-2.80E-2,-1.33E-2,-3.47E-3, + -2.33E-2,-3.69E-2,-4.41E-2,-4.38E-2,-3.97E-2,-3.50E-2,-2.88E-2,-1.60E-2,-6.68E-3, + -8.76E-3,-2.07E-2,-2.90E-2,-3.17E-2,-3.09E-2,-2.92E-2,-2.63E-2,-1.79E-2,-1.03E-2, + -1.20E-3,-1.11E-2,-1.90E-2,-2.20E-2,-2.32E-2,-2.24E-2,-2.10E-2,-1.66E-2,-1.11E-2, + 1.72E-3,-4.82E-3,-1.02E-2,-1.42E-2,-1.65E-2,-1.66E-2,-1.60E-2,-1.39E-2,-1.09E-2, + 2.68E-3,-1.18E-3,-5.19E-3,-8.30E-5,-1.01E-2,-1.14E-2,-1.16E-2,-1.16E-2,-9.99E-3, + 2.81E-3, 8.21E-4,-1.96E-3,-3.99E-3,-5.89E-3,-7.13E-3,-8.15E-3,-9.05E-3,-8.60E-3, + 2.61E-3, 1.35E-3,-2.99E-4,-1.79E-3,-3.12E-3,-4.44E-3,-5.61E-3,-7.01E-3,-7.27E-3, + 2.06E-3, 1.45E-3, 4.64E-4,-5.97E-4,-1.71E-3,-2.79E-3,-3.84E-3,-5.29E-3,-5.90E-3, + 1.07E-3, 9.39E-4, 8.22E-4, 3.58E-4,-1.15E-4,-6.60E-4,-1.18E-3,-2.15E-3,-2.88E-3, + 4.97E-4, 5.46E-4, 6.15E-4, 5.56E-4, 3.14E-4, 9.80E-5,-1.30E-4,-5.98E-4,-1.07E-4, + 1.85E-4, 3.11E-4, 4.25E-4, 4.08E-4, 3.63E-4, 3.04E-4, 2.24E-4, 2.80E-5,-2.10E-4, + 4.80E-5, 1.48E-4, 2.44E-4, 2.80E-4, 3.01E-4, 3.11E-4, 3.13E-4, 2.40E-4, 1.10E-4, + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 1.39E-4, 1.80E-4, 1.80E-4, + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 4.38E-5, 7.30E-5, 8.40E-5/ real f2_C(7,nRowC) + / 7.36E-2, 4.21E-2, 2.69E-2, 1.83E-2, 1.34E-2, 1.01E-2, 7.88E-3, + 5.79E-2, 3.61E-2, 2.34E-2, 1.64E-2, 1.21E-2, 9.26E-3, 7.28E-3, + 2.94E-2, 2.17E-2, 1.60E-2, 1.23E-2, 9.49E-3, 7.45E-3, 5.95E-3, + 2.30E-3, 7.07E-3, 7.76E-3, 7.02E-3, 6.13E-3, 5.17E-3, 4.34E-3, + -7.50E-3,-2.00E-3, 9.93E-4, 2.36E-3, 2.82E-3, 2.86E-3, 2.72E-3, + -8.27E-3,-5.37E-3,-2.58E-3,-7.96E-4, 3.75E-4, 9.71E-4, 1.28E-3, + -5.79E-3,-5.12E-3,-3.86E-3,-2.46E-3,-1.20E-3,-3.74E-4, 1.74E-4, + -3.26E-3,-3.43E-3,-3.26E-3,-2.68E-3,-1.84E-3,-1.12E-3,-4.54E-4, + -1.46E-3,-1.49E-3,-2.20E-3,-2.18E-3,-1.85E-3,-1.40E-3,-8.15E-4, + -4.29E-4,-9.44E-4,-1.29E-3,-1.50E-3,-1.51E-3,-1.36E-3,-9.57E-4, + -3.30E-5,-3.66E-4,-6.78E-4,-9.38E-4,-1.09E-3,-1.09E-3,-9.56E-4, + 1.50E-4, 3.10E-5,-1.38E-4,-3.06E-4,-4.67E-4,-5.48E-4,-6.08E-4, + 1.00E-4, 8.50E-5, 2.30E-5,-6.60E-5,-1.58E-4,-2.40E-4,-3.05E-4, + 5.40E-5, 6.50E-5, 4.90E-5, 1.20E-5,-3.60E-5,-8.90E-5,-1.31E-4, + 2.90E-5, 4.30E-5, 4.40E-5, 2.90E-5, 5.10E-6,-2.20E-5,-4.80E-5, + 1.40E-5, 2.40E-5, 2.80E-5, 2.60E-5, 1.90E-5, 7.50E-6,-1.10E-5, + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , + 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 , 0.00 / c=============================================================================== c Bestimme, welche Reihen der Tabellen fuer Interpolation benoetigt werden: if (tau.LT.tau_(1)) then write(*,*) 'tau is less than the lowest tabulated value:' write(*,*) 'tau = ',tau write(*,*) 'minimum = ',tau_(1) call exit elseif (tau.GT.tau_(nColumn)) then write(*,*) 'tau is greater than the highest tabulated value:' write(*,*) 'tau = ',tau write(*,*) 'maximum = ',tau_(nColumn) call exit endif column_ = 2 do while (tau.GT.tau_(column_)) column_ = column_ + 1 enddo ! Das Gewicht der Reihe zu groesserem Tau: weightCol = (tau-tau_(column_-1)) / (tau_(column_)-tau_(column_-1)) c Besorge fuer gegebenes 'thetaSchlange' die interpolierten f1- und f2 -Werte c der beiden relevanten Reihen: c iColumn = 1 => Reihe mit hoeherem Index c iColumn = 2 => Reihe mit kleinerem Index iColumn = 1 5 continue if (column_.LE.9) then ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Werte aus 1. Tabelle: 0.2 <= tau <= 1.8 column = column_ if (thetaSchlange.LT.thetaSchlangeA(1)) then write(*,*) 'thetaSchlange is less than the lowest tabulated value in table 1:' write(*,*) 'thetaSchlange = ',thetaSchlange write(*,*) 'minimum = ',thetaSchlangeA(1) call exit elseif (thetaSchlange.GT.thetaSchlangeA(nRowA)) then c write(*,*) 'thetaSchlange is greater than the highest tabulated value in table 1:' c write(*,*) 'thetaSchlange = ',thetaSchlange c write(*,*) 'maximum = ',thetaSchlangeA(nRowA) c call exit thetaSchlange = -1. RETURN endif row = 2 do while (thetaSchlange.GT.thetaSchlangeA(row)) row = row + 1 enddo ! Gewicht des Tabellenwertes zu groesseren ThetaSchlange: weightRow = (thetaSchlange-thetaSchlangeA(row-1)) / + (thetaSchlangeA(row)-thetaSchlangeA(row-1)) f1_(iColumn) = (1.-weightRow) * f1_A(column,row-1) + + weightRow * f1_A(column,row) f2_(iColumn) = (1.-weightRow) * f2_A(column,row-1) + + weightRow * f2_A(column,row) elseif (column_.LE.18) then ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Werte aus 2. Tabelle: 2.0 <= tau <= 7.0 column = column_ - 9 if (thetaSchlange.LT.thetaSchlangeB(1)) then write(*,*) 'thetaSchlange is less than the lowest tabulated value in table 1:' write(*,*) 'thetaSchlange = ',thetaSchlange write(*,*) 'minimum = ',thetaSchlangeB(1) call exit elseif (thetaSchlange.GT.thetaSchlangeB(nRowB)) then c write(*,*) 'thetaSchlange is greater than the highest tabulated value in table 1:' c write(*,*) 'thetaSchlange = ',thetaSchlange c write(*,*) 'maximum = ',thetaSchlangeB(nRowB) c call exit thetaSchlange = -1. RETURN endif row = 2 do while (thetaSchlange.GT.thetaSchlangeB(row)) row = row + 1 enddo ! Gewicht des Tabellenwertes zu groesseren ThetaSchlange: weightRow = (thetaSchlange-thetaSchlangeB(row-1)) / + (thetaSchlangeB(row)-thetaSchlangeB(row-1)) f1_(iColumn) = (1.-weightRow) * f1_B(column,row-1) + + weightRow * f1_B(column,row) f2_(iColumn) = (1.-weightRow) * f2_B(column,row-1) + + weightRow * f2_B(column,row) else ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! Werte aus 3. Tabelle: 8.0 <= tau <= 20. column = column_ - 18 if (thetaSchlange.LT.thetaSchlangeC(1)) then write(*,*) 'thetaSchlange is less than the lowest tabulated value in table 1:' write(*,*) 'thetaSchlange = ',thetaSchlange write(*,*) 'minimum = ',thetaSchlangeC(1) call exit elseif (thetaSchlange.GT.thetaSchlangeC(nRowC)) then c write(*,*) 'thetaSchlange is greater than the highest tabulated value in table 1:' c write(*,*) 'thetaSchlange = ',thetaSchlange c write(*,*) 'maximum = ',thetaSchlangeC(nRowC) c call exit thetaSchlange = -1. RETURN endif row = 2 do while (thetaSchlange.GT.thetaSchlangeC(row)) row = row + 1 enddo ! Gewicht des Tabellenwertes zu groesseren ThetaSchlange: weightRow = (thetaSchlange-thetaSchlangeC(row-1)) / + (thetaSchlangeC(row)-thetaSchlangeC(row-1)) f1_(iColumn) = (1.-weightRow) * f1_C(column,row-1) + + weightRow * f1_C(column,row) f2_(iColumn) = (1.-weightRow) * f2_C(column,row-1) + + weightRow * f2_C(column,row) endif ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if (iColumn.EQ.1) then column_ = column_ - 1 iColumn = 2 goto 5 endif f1 = weightCol*f1_(1) + (1.-weightCol)*f1_(2) f2 = weightCol*f2_(1) + (1.-weightCol)*f2_(2) END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE reset_statistics c =========================== IMPLICIT NONE integer Nr,n,k INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.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 die Scaler fuer den Returncode des TDs und die Pfostenhits sowie die c StartZaehler resetten: do n = 1, 2 ! (1: Projektile, 2: FolienElektronen) start_nr(n) = 0 do k = 1, 18 statTD(n,k) = 0. enddo do k = 1, 75 pfostenHit(k,n) = 0. enddo enddo 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 'mutrack$sourcedirectory:COM_MUTRACK.INC' real wert integer Nr c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Wird die Varianz der Verteilung einer Groesse x gemaess der Formel c c Var(x) = SQRT( - **2 ) , < > -> Erwartungswert c c mit c = 1/n * Summe_ueber_i( x(i) ) c = 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 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 (*): wird im SUB 'eval_statistics' berechnet. c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Zaehle mit: statMem(8,Nr) = statMem(8,Nr) + 1. c Speichere den ersten Wert: if (statMem(8,Nr).EQ.1) statMem(1,Nr) = wert c Summiere die Abweichungen vom ersten Wert: statMem(2,Nr) = statMem(2,Nr) + (wert-statMem(1,Nr)) c Summiere die Quadratischen Abweichungen vom ersten Wert: statMem(3,Nr) = statMem(3,Nr) + (wert-statMem(1,Nr))**2. c Speichere den kleinsten Wert (wurde noch kein Wert aufgenommen, so ist c statMem(4,Nr) = 1.e10): if (statMem(4,Nr).GT.wert) statMem(4,Nr) = wert c Speichere den groessten Wert (wurde noch kein Wert aufgenommen, so ist c statMem(5,Nr) = -1.e10): if (statMem(5,Nr).LT.wert) statMem(5,Nr) = wert END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE eval_statistics c ========================== IMPLICIT NONE 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)) INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC' real n ! Anzahl der Werte, == statMem(8,Nr) real radiant integer Nr,l do Nr = 1, Stat_Anzahl if (statNeeded(Nr)) then n = statMem(8,Nr) if (n.ne.0.) then !c Berechne Mittelwert: statMem(6,Nr) = statMem(2,Nr)/n + statMem(1,Nr) !c Berechne Varianz: radiant = ( statMem(3,Nr) - (statMem(2,Nr)**2. )/n)/n statMem(7,Nr) = sqrt(radiant) !c Berechne Anteil an allen gestarteten in Prozent statMem(9,Nr) = 100.*n/real(n_par(0)) else do l = 1, 9 statMem(l,Nr) = 0. enddo endif endif enddo END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE SAVE_GRAPHICS_KOORD c ============================== IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC' INCLUDE 'mutrack$sourcedirectory:COM_WINKEL.INC' INCLUDE 'mutrack$sourcedirectory:COM_KAMMER.INC' c Variablen fuer die Graphikausgabe: real xKoord(1000) ! Koordinatenfelder fuer die real yKoord(1000) ! Graphikausgabe real zKoord(1000) ! cMBc real tKoord(1000) ! integer nKoord ! Anzahl der Koordinaten cMBc COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord,tKoord ! fuer Graphikaufruf COMMON /GRAPHIX/ xKoord,yKoord,zKoord,nKoord ! fuer Graphikaufruf c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - nKoord = nKoord + 1 xKoord(nKoord) = x(1) yKoord(nKoord) = x(2) zKoord(nKoord) = x(3) cMBc tKoord(nKoord) = t if (nKoord.EQ.1000) then if (Gebiet.LE.upToChKoord) then ! Bahnberechnung wurde vor call plot_horizontal ! Koordinatenwechsel abgebrochen else call plot_vertikal endif xKoord(1) = xKoord( 999) ! die letzten beiden uebernehmen, yKoord(1) = yKoord( 999) ! damit gegebenenfalls der Richtungs- zKoord(1) = zKoord( 999) ! pfeil gezeichnet werden kann. cMBc tKoord(1) = tKoord( 999) xKoord(2) = xKoord(1000) yKoord(2) = yKoord(1000) zKoord(2) = zKoord(1000) cMBc tKoord(2) = tKoord(1000) nKoord = 2 endif END c=============================================================================== OPTIONS /EXTEND_SOURCE SUBROUTINE Output_Debug c ======================= IMPLICIT NONE INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC' INCLUDE 'mutrack$sourcedirectory:COM_WINKEL.INC' INCLUDE 'mutrack$sourcedirectory:COM_KAMMER.INC' real Ekin, temp1, temp2 Ekin = (v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) * Energie_Faktor if (Gebiet.EQ.1 .AND. alfaTgt.NE.0) then if (alfaTgtVertically) then temp1 = xGrid1*Cos_alfaTgt - x(3)*Sin_alfaTgt temp2 = xGrid1*Sin_alfaTgt + x(3)*Cos_alfaTgt write(lun(1),1) steps,Gebiet,t,temp1,x(2),temp2,v,Ekin else temp1 = xGrid1*Cos_alfaTgt - x(2)*Sin_alfaTgt temp2 = xGrid1*Sin_alfaTgt + x(2)*Cos_alfaTgt write(lun(1),1) steps,Gebiet,t,temp1,temp2,x(3),v,Ekin endif else write(lun(1),1) steps,Gebiet,t,x,v,Ekin endif 1 format(X,I4,X,I2,4X,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 'mutrack$sourcedirectory:COM_MUTRACK.INC' real dt if (t.GT.lifeTime) then ! Teilchen zerfallen dt = t - lifeTime t = lifeTime 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 chargeStateYields(E,masse,Yield_plus,Yield_zero) c =========================================================== c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c Die Funktion sowie die Parameter sind uebernommen aus: c c M.Gonin, R.Kallenbach, P.Bochsler: 'Charge exchange of hydrogen atoms c in carbon foils at 0.4 - 120 keV', Rev.Sci.Instrum. 65 (3), March 1994 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IMPLICIT NONE real E ! kinetische Energie in keV real masse ! in keV / c**2 real a_zero,a_minus real k_Fermi,k_zero,k_minus real zwo_k_Fermi real k_Fermi_Quad,k_zero_Quad,k_minus_Quad real vc_minus,vc_plus,v_Bohr,v_rel parameter ( a_zero = 0.953, a_minus = 0.029 ) parameter ( k_Fermi = 1.178 ) ! [v_Bohr] parameter ( k_Fermi_Quad = k_Fermi * k_Fermi ) parameter ( zwo_k_fermi = 2. * k_Fermi ) parameter ( k_zero = 0.991*k_Fermi ) ! [v_Bohr] parameter ( k_zero_Quad = k_zero * k_zero ) parameter ( k_minus = 0.989*k_Fermi ) ! [v_Bohr] parameter ( k_minus_Quad = k_minus * k_minus ) parameter ( vc_minus = 0.284, vc_plus = 0.193 ) ! [v_Bohr] parameter ( v_Bohr = 7.2974E-3 ) ! [c] real Q_zero,Q_minus,D real Yield_minus,Yield_zero,Yield_plus real help1,help2,help3 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if (E.LT.0) then write(*,*) write(*,*) 'error in subroutine ''chargeStateYields'':' write(*,*) 'E = ',E,' < 0!' write(*,*) '-> STOP' write(*,*) STOP endif c Energie in Geschwindigkeit umrechnen (in Einheiten von v_Bohr): c - klassisch: v_rel = SQRT(2.*E/masse) / v_Bohr c - relativistisch: c help1 = 1. + E/masse c v_rel = SQRT(1. - 1./(help1*help1)) / v_Bohr c Die geladenen Anteile berechnen (vgl. obige Referenz): help1 = v_rel*v_rel help2 = zwo_k_Fermi*v_rel Q_zero = 1. + (k_zero_Quad - k_Fermi_Quad - help1) / help2 Q_minus = 1. + (k_minus_Quad - k_Fermi_Quad - help1) / help2 help1 = a_zero * Q_zero help2 = a_minus * Q_minus help3 = (1.-Q_zero)*(1.-Q_minus) D = help1*(help2 + (1.-Q_minus)) + help3 Yield_minus = help1*help2 / D Yield_plus = help3 / D Yield_minus = Yield_minus * exp(-vc_minus/v_rel) Yield_plus = Yield_plus * exp(-vc_plus /v_rel) Yield_zero = 1. - (Yield_minus + Yield_plus) c write(6,*) 'E vrel Neutral Plus Minus' c write(6,*) E, v_rel, yield_zero, yield_plus, yield_minus 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 Teilchen und Draht und gegebener Geschwindigkeit v eines Teilchens c bei geradliniger Bewegung und Drahtradius 'WireRadius' ein Schnittpunkt c von Teilchenbahn und Drahtumfang existiert, ob also der Draht getroffen wird. c Dafuer genuegt es zu pruefen, ob der Radiant der 'Mitternachtsformel' fuer die c entsprechende quadratische Gleichung groesser oder gleich Null ist: IMPLICIT NONE real DistToWire(2),WireRadiusQuad,v_x,v_y logical WireHit real steigung, help, radiant c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if (abs(v_x).GT.abs(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===============================================================================