LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
musrsim/mcv3k/src/mcv3k.f

2265 lines
84 KiB
Fortran
Raw Blame History

C%CREATE MCV3
CC
CC MONTE-CARLO RECHNUNG:
CC DURCHGANG VON SCHWEREN TEILCHEN DURCH MATERIE
CC
CC REV.: ELOSS2 ANSTELLE VON RGEFIN EINBAUT.
CC REV.: UND DOPPELTGENAUE RGE-E-TABELLEN
CC REV.: Einführung eines Parameters für den Befehl 'BUILD_TAB', mit dem die
CC Variable "PREC" und damit die relative Genauigkeit der Interpolation
CC beim Erstellen der Reichweite-Energie-Tabelle bestimmt werden kann.
CC Damit lassen sich zum Programmabbruch führende Fehler in der Routine
CC SPLIN3 vermeiden, die zuvor bei gewissen Materialien wie z. B. reinem
CC Aluminium auftraten. Zusätzlich gibt es hier und in UPLIBK.FOR einige
CC weitere kleine Änderungen. KT 10-JAN-95
CC REV.: Eine und damit vielleicht bereits die alleinige Ursache für Programm-
CC abstürze, die besonders bei hohen Teilchenzahlen auftraten, ist gefun-
CC den und beseitigt. Die Energie-Reichweite-Tabelle wird von einem kubi-
CC schen Spline-Polynom angenähert, das sich aufgrund von Rundungsfehlern
CC unter Umständen numerisch nicht korrekt lösen läßt. In der Funktion
CC ELOSS2 in UPLIBK.FOR wurde bislang in diesem Fall kommentarlos eine
CC Stop-Anweisung ausgeführt. Jetzt wird das betreffende Teilchen aus der
CC weiteren Rechnung genommen und dies in einem Hinweis in der Standard-
CC Ausgabedatei vermerkt. Bei der anzunehmenden Rate von einem bis zwei
CC Teilchen pro 100,000 wird dies den Aussagewert der Rechnungen nicht
CC einschränken. In der Funktion ELOSS2 wurde ferner die Genauigkeit der
CC Konstanten 4*PI/3 und 2*PI/3 von 8 auf 16 Stellen verbessert, wodurch
CC eine Reihe von Fehlberechnungen bereits vermieden wurde. Allen weiter-
CC en Stop-Anweisungen, die zu einem ähnlich unerklärten Programmende
CC führten, wurde die Ausgabe eines Hinweises über die Ursachen des Ab-
CC bruchs und deren mögliche Behebung vorangestellt.
CC In der Standard-Ausgabedatei werden jetzt Datum und Uhrzeit vom Start
CC und Ende des Programms sowie die Dauer der Rechnung festgehalten.
CC KT 30-JAN-95
CC REV.: Der Fehler in der Zeitausgabe bei einer Rechnung über Mitternacht ist
CC korrigiert. Vier neue Eingabedatei-Befehle sind implementiert.
CC Die beiden Befehle IMPULS und WINKEL haben eine analoge Struktur. Sie
CC bewirken die Ausgabe der Impulsbeträge [MeV/c] bzw. des Streuwinkels
CC zur Z-Achse [rad] der Teilchen in die durch den Parameter festgelegte
CC Unit vor bzw. hinter einem Bauteil, das durch die Position des Befehls
CC in der Eingabedatei festgelegt ist. In den Vektor IMPOUT bzw. THEOUT
CC mit jeweils der Länge 201 wird an die Position der Nummer des nachfol-
CC genden Bauteiles die Nummer der Ausgabe-Unit geschrieben. In der Rou-
CC tine COMPUT erfolgt später genau dann eine Ausgabe, wenn der zu dem
CC nachfolgenden Bauteil gehörende Wert des Vektors von Null verschieden
CC ist. Um bei Verwendung des Befehls SKIP keine unerwünschte Ausgabe zu
CC erhalten, sollten 'IMPULS' und 'WINKEL' vor 'GEOMETRY' stehen.
CC Auch der Befehl STOPPED hat die Nummer der Ausgabe-Unit als Parameter
CC mit demselben Prinzip. Er bewirkt die Ausgabe der Stoppverteilung in
CC dem zuletzt definierten Bauteil in Form einer Datei, die soviele Werte
CC enthält wie das Bauteil in Einzellagen aufgeteilt ist.
CC Mit dem Befehl FOCUS läßt sich schließlich der Strahl der einkommenden
CC Teilchen fokussieren. Der erste Parameter gibt die Z-Position [cm] der
CC x=0 Fokusgeraden an, der zweite die Z-Position der y=0 Fokusgeraden,
CC die, keine Ablenkung durch Streuung vorausgesetzt, alle Teilchenbahnen
CC schneiden. Erlaubt sind Eingaben von -1000 bis +1000; ein negativer
CC Wert bewirkt eine Defokussierung, Null beläßt den Strahl parallel.
CC KT 16-FEB-95
CC REV.: Implementierung eines neuen Eingabedatei-Befehls zue Ausgabe der X,Y-
CC Koordinaten ('WHERE_XY'). Die Bearbeitung erfolgt völlig analog zu
CC den Befehlen IMPULS und WINKEL. Ferner die endgültige Lösung des
CC Mitternacht-Problems. KT 21-FEB-95
CC REV.: Änderung der Zuweisungen von Ein- und Ausgabeunits, um den Code auch
CC für AXP-Rechner verwenden zu können. KT 5-MAR-95
CC REV.: Implementierung eines homogenen Magnetfeldes parallel zu Z, Befehl
CC MAGNETIC. Parameter: 1. Z-Koordinate des Anfangs, 2. Stärke in Gauß,
CC wobei das Vorzeichen die Richtung in Bezug auf Z angibt. Ein neuer
CC Befehl für die X,Y-Ausgabe, der nicht nur die Teilchen erfaßt, die
CC das Bauteil treffen: ALL_XY.
CC
CC Alle Änderungen sind mit "c-kt" kommentiert.
cc
cc
cc
CC REV.: Changed the BEAM_XY command a little; now the Beam parameter
cc is also required as input parameter: "1" means isotropic circular
cc "2" gaussian
cc "3" isotropic rectangular
cc TP, 1-jun-1995
CC
CC
cc REV.: New command INIT_OUT to write the start values of location
cc momentum, energy and time to file; usage as the new commands
cc by kt. Also changed the formatted output from the OUTPUT command
cc
cc TP, 1-jun-1995
cc
cc REV.: Option for unformatted output: new command FORMATTED and UFORMATTED;
cc If FORMATTED
cc is set, then formatted output will be made for the 'in front of
cc data', the 'Behind' and the 'init' data, TP, 2-jun-1995
cc
cc
cc TP, 17-Dec-1996: Aenderungen in MCV3K_SUB.FOR um Programm auch auf Alpha
cc Maschinen laufen zu lassen (z.B. neue SPLINE Routine,
cc diverse REAL*8 Deklarationen).
cc AH hat noch einige ueberfluessige Variablendekla.
cc gefunden, die ebenfalls herauskommentiert wurden
cc
cc TP, 12-Sep-2000: some changes to get it run in Windows NT and Unix
cc exchanged SYSOUT and SYSIN, SYSOUT was set
cc to 5 and
cc SYSIN to 6 ? Commented out str$upcase.
cc commented out open and close of
cc open(SYSOUT,...,FILE='SYS$OUTPUT'...)
cc
cc----------------------------------------------------------------------------
cc
cc
BLOCK DATA
COMMON / ATOMIC / RADLEN(92)
C RADIATION LENGTH IN G/CM**2
DATA RADLEN /63.0470,94.3221,82.7559,65.1899,52.6868,
1 42.6983,37.9879,34.2381,32.9303,28.9367,27.7362,
2 25.0387,24.0111,21.8234,21.2053,19.4953,19.2783,
3 19.5489,17.3167,16.1442,16.5455,16.1745,15.8425,
4 14.9444,14.6398,13.8389,13.6174,12.6820,12.8616,
5 12.4269,12.4734,12.2459,11.9401,11.9082,11.4230,
6 11.3722,11.0272,10.7623,10.4101,10.1949,9.9225,
7 9.8029,9.6881,9.4825,9.2654,9.2025,8.9701,8.9945,
8 8.8491,8.8170,8.7244,8.8267,8.4803,8.4819,8.3052,
9 8.3073,8.1381,7.9557,7.7579,7.7051,7.5193,7.5727,
1 7.4377,7.4830,7.3563,7.3199,7.2332,7.1448,7.0318,
2 7.0214,6.9237,6.8907,6.8177,6.7630,6.6897,6.6763,
3 6.5936,6.5433,6.4608,6.4368,6.4176,6.3688,6.2899,
4 6.1907,6.0651,6.2833,6.1868,6.1477,6.0560,6.0726,
5 5.9319,5.9990/
END
C ----------------------------------------------------------------------
C H A U P T P R O G R A M M
C ----------------------------------------------------------------------
CHARACTER*10 CMD, CMDTAB(33), GEOTYP(2), OUTTYP(7), BEAMTY(3)
CHARACTER*20 COMENT
c-kt Uhrzeit, Datum und Systemzeit bei Programmbeginn und -ende sowie
c-kt die Zählvariable für die aus der Rechnung genommenen Teilchen
CHARACTER*8 startzeit, endezeit
CHARACTER*9 startdatum, endedatum
REAL*4 startsec, dauer
INTEGER*4 h, m, s, fehler
c-kt Ende
c-kt Variablen für MAGNETIC: Beginn und Stärke des Magnetfeldes [cm], [Gauss]
REAL*4 MSTART, MGAUSS
REAL*8 RGDATA
INTEGER*4 SEEDS
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY und ALL_XY:
INTEGER*2 IMPOUT, THEOUT, STPOUT, XYOUT, XYALL, LUNINIT
c-kt Deklaration neuer Bezeichner für Ein- und Ausgabeunit: SYSOUT, SYSIN
INTEGER*2 SYSOUT, SYSIN
LOGICAL TEST, RULER, LAYON, LAYGEO, BEAMP,
1 BEAMXY, SEEDOK, MASSOK, TABOK, SCATTR, STRAGG,
2 SKIP, form
COMMON
1 / EXECFL / SCATTR, STRAGG, RELSTR
COMMON
1 / ATOMIC / RADLEN(92)
COMMON
1 / RETAB / NRGMAX, NRGCUR, RGDATA (5,3000)
2 / LAYER / NLAY, ILAYS(200), ILAYE(200),
3 THICK(200), RTHICK(200), RO(200),
4 ZMEAN(200), AMEAN(200), RMEAN(200),
5 ZKOO(200), IGEOT(200), PGEOT(5,200)
6 / LAYERG / NGRP, IGRPS(200), IGRPE(200)
7 / LAYERC / COMENT(200)
8 / PARTIC / BMASS, BMASS2, EKMIN, EKMAX, IBFLAG, IBUNIT
9 / BEAM / PMEAN, PFWHM, PLOW, PHIGH, XFWHM, YFWHM,
1 RMAX,
2 IBTYP, XBEAM, YBEAM, PSIGMA, RMAX2
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
3 , XFOCUS, YFOCUS
c-kt Beginn und Stärke des Magnetfeldes
4 , MSTART, MGAUSS
COMMON
1 / KINEMA / ID, XX, YY, ZZ,
2 PX, PY, PZ, PTOTAL,
c-kt Einführen der Variable TTHETA für den Tangens des Streuwinkels:
3 THETA, PHI, STHETA, CTHETA, TTHETA,
4 SPHI, CPHI,
5 TTOTAL, TLAST,
6 ELOSS
7 / OUTPUT / NLIMIT(200), NCOUNT(200), IOTYPE(200),
8 NLUNIT(200), ELIMIT(200)
9 / COMPON / NZCOMP(20), ACOMP(20), NATOM(20), NELEM, NELEMC
1 / SEEDCO / SEEDS (6)
COMMON
c-kt Deklaration neuer Bezeichner für Ein- und Ausgabeunit: SYSOUT, SYSIN
1 / STATUS / SYSOUT, SYSIN, TEST
3 / INPUT / WPARM(6)
c-kt Zählvariable für die aus der Rechnung genommenen Teilchen:
COMMON
1 / COUNT / fehler
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY und ALL_XY:
COMMON
1 / SGLOUT / IMPOUT(201), THEOUT(201), STPOUT(200), XYOUT(201),
2 XYALL(201), LUNINIT
DATA GEOTYP / 'CIRCLE ','RECTANGLE '/
DATA OUTTYP / 'NOTHING ','IN FRONT ','STOPPED ',
1 'BEHIND ','STP,BEHIND','IN FR.,STP',
2 'ALL '/
DATA BEAMTY / 'ISOTROPIC ','GAUSSIAN ', 'RECTANGULA' /
C ----------------------------------------------------------------------
C TABELLE DER GUELTIGEN KOMMANDOS
C ----------------------------------------------------------------------
DATA CMDTAB(1)/'BEAM_P '/
DATA CMDTAB(2)/'LAYER '/
DATA CMDTAB(3)/'GEOMETRY '/
DATA CMDTAB(4)/'COMPONENT '/
DATA CMDTAB(5)/'BUILD_TAB '/
DATA CMDTAB(6)/'HOLE '/
DATA CMDTAB(7)/'OUTPUT '/
DATA CMDTAB(8)/'COMPUTE '/
DATA CMDTAB(9)/'EXIT '/
DATA CMDTAB(10)/'RULER '/
DATA CMDTAB(11)/'NORULER '/
DATA CMDTAB(12)/'TEST '/
DATA CMDTAB(13)/'NOTEST '/
DATA CMDTAB(14)/'SEEDS '/
DATA CMDTAB(15)/'PARTICLE '/
DATA CMDTAB(16)/'BEAM_XY '/
DATA CMDTAB(17)/'OLD_LAYER '/
DATA CMDTAB(18)/'STRAG_ON '/
DATA CMDTAB(19)/'STRAG_OFF '/
DATA CMDTAB(20)/'SCATT_ON '/
DATA CMDTAB(21)/'SCATT_OFF '/
DATA CMDTAB(22)/'SUMMARY '/
DATA CMDTAB(23)/'SKIP '/
c-kt Die neuen Befehle IMPULS, WINKEL, STOPPED, WHERE_XY, ALL_XY, FOCUS
c-kt und MAGNETIC
DATA CMDTAB(24)/'IMPULS '/
DATA CMDTAB(25)/'WINKEL '/
DATA CMDTAB(26)/'STOPPED '/
DATA CMDTAB(27)/'FOCUS '/
DATA CMDTAB(28)/'WHERE_XY '/
DATA CMDTAB(29)/'MAGNETIC '/
DATA CMDTAB(30)/'ALL_XY '/
c
DATA CMDTAB(31)/'INIT_OUT '/
DATA CMDTAB(32)/'FORMATTED '/
DATA CMDTAB(33)/'UFORMATTED'/
DATA NNAME/33/
c-kt Ende
C ----------------------------------------------------------------------
C INITIALISIERUNG ALLER VARIABLEN
C ----------------------------------------------------------------------
c-kt den Zeitpunkt des Programmbeginns festhalten
CALL DATE(startdatum)
CALL TIME(startzeit)
startsec = SECNDS(0.0)
c-kt Initialisieren des Zählers für die aus der Rechnung genommenen Teilchen
fehler = 0
c-kt Ende
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
XFOCUS = 0.0
YFOCUS = 0.0
c-kt Ende
ZKOMAX = 0.0
NREP = 0
NRGMAX = 3000
NRGCUR = 0
NGRP = 0
IKARD = 0
SYSOUT = 6
SYSIN = 5
NLAY = 0
NELEM = 0
SCATTR = .TRUE.
STRAGG = .TRUE.
RULER = .FALSE.
TEST = .FALSE.
LAYON = .FALSE.
LAYGEO = .FALSE.
TABOK = .FALSE.
BEAMP = .FALSE.
BEAMXY = .FALSE.
SEEDOK = .FALSE.
MASSOK = .FALSE.
SKIP = .FALSE.
FORM = .FALSE.
DO 10 I=1,200
IOTYPE(I) = 1
NLIMIT(I) = 0
NLUNIT(I) = 6
c-kt Integer-Array für die Unit-Nummern zu dem Befehl STOPPED:
STPOUT(I) = 0
10 CONTINUE
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt WHERE_XY und ALL_XY
DO 11 I=1,201
IMPOUT(I) = 0
THEOUT(I) = 0
XYOUT(I) = 0
XYALL(I) = 0
11 CONTINUE
LUNINIT = 0
c-kt Ende
c
c-kt Zuweisen des logischen Namens SYS$OUTPUT an die Unit SYSOUT:
c OPEN(UNIT=SYSOUT,FILE='SYS$OUTPUT',STATUS='NEW')
c-kt Zuweisen des logischen Namens SYS$INPUT an die Unit SYSIN:
c OPEN(UNIT=SYSIN,FILE='SYS$INPUT',READONLY,STATUS='OLD')
c
C ----------------------------------------------------------------------
C LESE DAS NAECHSTE KOMMANDO EIN
C ----------------------------------------------------------------------
50 IKARD = IKARD + 1
C LESE NAECHSTES CMD
IF ( RULER ) WRITE (SYSOUT,5017)
5017 FORMAT (' NEXT COMMAND...'/' ....;....1....;....2....;....3....;'
1 ,'....4....;....5....;....6....;....7....;')
READ(SYSIN,5000) CMD, WPARM
5000 FORMAT (A10, 6F10.0)
c
c-tp write string in upper case
c
c call str$upcase(cmd,cmd)
IF ( CMD .EQ. ' ' ) GOTO 50
if(cmd(1:1) .eq. '*') goto 50 ! '*' means comment in Input file
C ANZAHL DER ARG. ?
K = 1
DO 60 I= 1, 6
IF (WPARM(I)) 58,60,58
58 K = I
60 CONTINUE
C PROTOKOLLIERE
IF ( TEST ) THEN
WRITE (SYSOUT,5001) IKARD, CMD, ( WPARM(I), I=1, K )
5001 FORMAT (' ',10('*')/' ** ',I3,' **',A10,7F15.5)
WRITE (SYSOUT,5004)
5004 FORMAT (' **********')
ENDIF
C SUCHE KOMMANDO
DO 80 I= 1, NNAME
IF ( CMD .EQ. CMDTAB(I) ) GO TO 90
80 CONTINUE
C KOMMANDO NICHT GEF.
WRITE (SYSOUT,5006)
5006 FORMAT ('+',10X,'(COMMAND IGNORED)')
GO TO 50
c-kt erweiterte Tabelle
CC 1 2 3 4 5 6 7 8 9 10
90 GO TO ( 100, 200, 300, 400, 500, 600, 700, 800, 900,1000,
1 1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,
2 2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,
3 3100,3200,3300 ), I
C ----------------------------------------------------------------------
C "BEAM_P" IMPULSVERTEILUNG DES STRAHLS
C ----------------------------------------------------------------------
100 IF ( BEAMP ) THEN
WRITE (SYSOUT,5010)
5010 FORMAT(' MOMENTUM DISTRIBUTION OF BEAM ALREADY DEFINED.')
STOP 12
ENDIF
PMEAN = ABS( WPARM(1) )
RPFWHM = ABS( WPARM(2) )
PFWHM = PMEAN * RPFWHM
PLOW = ABS( WPARM(3) )
PHIGH = ABS( WPARM(4) )
c CALL TBOUND ('R4','P_MEAN ',PMEAN, 1.0, 5000. , 'IN')
c
CALL TBOUND ('R4','P_MEAN ',PMEAN, 0.7, 5000. , 'IN')
CALL TBOUND ('R4','P_FWHM ',PFWHM, 0.0, 5.*PMEAN, 'IN')
CALL TBOUND ('R4','P_LOW ',PLOW , 1.E-1, PMEAN , 'IN')
CALL TBOUND ('R4','P_HIGH ',PHIGH, PMEAN, 5000. , 'IN')
PSIGMA = PFWHM * 0.4246609
BEAMP = .TRUE.
IF ( MASSOK ) THEN
DUMMY = STRAG1( PMEAN, BMASS )
WRITE(SYSOUT,5012) PMEAN, DUMMY
5012 FORMAT(' RELATIVE STRAGGLING FOR ',1PE12.5,' MEV/C = ',1PE12.5)
ENDIF
IF ( TEST ) THEN
WRITE(SYSOUT,5011) PMEAN,PFWHM,RPFWHM,PLOW,PHIGH
5011 FORMAT(' ','BEAM PARAMETER'/
1 6X,'MEAN MOMENTUM IN MEV/C',T60,1PE12.5/
2 6X,'ABSOLUTE FWHM OF THE MOMENTUM DISTRIBUTION',T60,1PE12.5/
3 6X,'RELATIVE FWHM OF THE MOMENTUM DISTRIBUTION',T60,1PE12.5/
4 6X,'LOWEST MOMENTUM IN MEV/C',T60,1PE12.5/
5 6X,'HIGHEST MOMENTUM IN MEV/C',T60,1PE12.5)
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "LAYER" DATEN FUER DIE BESCHAFFENHEIT EINER SCHICHT
C ----------------------------------------------------------------------
200 IF ( LAYON ) THEN
WRITE (SYSOUT,5021)
5021 FORMAT(' PREVIOUS LAYER OR HOLE-COMMAND STILL ACTIVE.')
STOP 12
ENDIF
NLAY = NLAY + 1
NGRP = NGRP + 1
CALL TBOUND ('I4','#LAYER ',NLAY ,1,200,'IN ')
CALL TBOUND ('I4','#GROUP ',NGRP ,1,200,'IN ')
NELEM = IFIX( WPARM(1) + 0.5 )
RO (NLAY) = WPARM (2)
THICK (NLAY) = WPARM (3)
NDIV = MAX0( IFIX( WPARM(4) + 0.5 ), 1 )
RPDIST = ABS( WPARM(5) )
NREP = NDIV - 1
THICK (NLAY) = THICK (NLAY) / FLOAT(NDIV)
READ (SYSIN,5171) COMENT(NGRP)
CALL TBOUND ('I4','#COMPONENT',NELEM,1,20,'IN ')
CALL TBOUND ('R4','DENSITY ',RO(NLAY),1.E-08,100.,'IN')
CALL TBOUND ('R4','THICKNESS ',THICK(NLAY),1.E-08,1.E08,'IN')
CALL TBOUND ('I4','REPETITION',NREP,0,(200-NLAY),'IN ')
CALL TBOUND ('R4','DISTANCE ',RPDIST,0.0, 1.E06,'IN')
RTHICK(NLAY) = RO(NLAY) * THICK(NLAY)
NELEMC = 0
LAYON = .TRUE.
LAYGEO = .FALSE.
TABOK = .FALSE.
NLAYXX = NLAY + NREP
IF ( TEST ) THEN
WRITE (SYSOUT,5022) NLAY, NLAYXX, COMENT(NLAY), NELEM, RO(NLAY),
1 THICK(NLAY), RTHICK(NLAY), NREP, RPDIST
5022 FORMAT (' ','LAYER#',I3,'-',I3,' / ',A20/
1 ' ',5X,'NO OF ELEMENTS',T45,I12/
2 ' ',5X,'DENSITY IN G/CM**3',T45,1PE12.5/
3 ' ',5X,'THICKNESS IN CM',T45,1PE12.5/
4 ' ',5X,'THICKNESS IN G/CM**2',T45,1PE12.5/
5 ' ',5X,'NUMBER OF ADDITIONAL LAYERS',T45,I12/
6 ' ',5X,'DISTANCE IN CM',T45,1PE12.5)
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "GEOMETRY" ALLE GEOMETRIEDATEN EINER SCHICHT
C ----------------------------------------------------------------------
300 IF ( LAYGEO .OR. .NOT.LAYON ) THEN
WRITE (SYSOUT,5031)
5031 FORMAT(' GEOMETRY-DATA ALREADY DEFINED OR NO LAYER-COMMAND',
1 ' ACTIVE.')
STOP 12
ENDIF
LAYGEO = .TRUE.
ZKOO (NLAY) = ABS (WPARM (1))
IGEOT(NLAY) = IFIX ( WPARM(2) + 0.5 )
PGEOT(1,NLAY) = WPARM(3)
CALL TBOUND ('R4','Z-COORD ',ZKOO(NLAY),ZKOMAX,5.E05,'IN ')
CALL TBOUND ('I4','GEO-TYPE ',IGEOT(NLAY),1,2,'IN ')
ZKOMAX = ZKOO(NLAY) + THICK(NLAY)
IF ( IGEOT(NLAY) .EQ. 2 ) THEN
DO 310 LL=2,4
PGEOT(LL,NLAY) = WPARM(LL+2)
310 CONTINUE
CALL TBOUND ('R4','X-SIZE ',(PGEOT(2,NLAY)-PGEOT(1,NLAY)),
1 1.E-05,500.0,'IN ')
CALL TBOUND ('R4','Y-SIZE ',(PGEOT(4,NLAY)-PGEOT(3,NLAY)),
1 1.E-05,500.0,'IN ')
ELSE
PGEOT(2,NLAY) = PGEOT(1,NLAY)**2
CALL TBOUND ('R4','RADIUS ',PGEOT(1,NLAY),1.E-05,500.0,'IN ')
ENDIF
IGRPS(NGRP) = NLAY
IGRPE(NGRP) = NLAY + NREP
IF ( .NOT. SKIP .AND. TEST ) THEN
WRITE (SYSOUT,5032) NLAY, NLAYXX, GEOTYP(IGEOT(NLAY)),
1 IGEOT(NLAY), ZKOO(NLAY)
5032 FORMAT (' ','SHAPE OF LAYER#',I3,'-',I3/
1 ' ',5X,'TYPE',T45,A10,'(',I1,')'/
2 ' ',5X,'Z-COORDINATE',T45,1PG12.5,' CM')
IF ( IGEOT(NLAY) .EQ. 1 ) THEN
WRITE (SYSOUT,5033) PGEOT(1,NLAY)
5033 FORMAT(' ',5X,'RADIUS',T45,1PG12.5,' CM')
ELSE
WRITE (SYSOUT,5034) (PGEOT(LL,NLAY),LL=1,4)
5034 FORMAT(' ',5X,'X-LEFT',T45,1PG12.5,' CM'/
1 ' ',5X,'X-RIGHT',T45,1PG12.5,' CM'/
2 ' ',5X,'Y-BOTTOM',T45,1PG12.5,' CM'/
3 ' ',5X,'Y-TOP',T45,1PG12.5,' CM')
ENDIF
ENDIF
IF ( .NOT. TABOK ) GOTO 50
C LAYER-DATA NOW COMPLETE
350 TABOK = .FALSE.
LAYON = .FALSE.
LAYGEO = .FALSE.
SKIP = .FALSE.
IF ( NREP .LE. 0 ) GOTO 50
C WIEDERHOLE "NREP" - MAL
IS = NLAY + 1
IE = NLAY + NREP
DO 360 IREP=IS,IE
ILAYS(IREP) = ILAYS(NLAY)
ILAYE(IREP) = ILAYE(NLAY)
RO(IREP) = RO(NLAY)
THICK(IREP) = THICK(NLAY)
RTHICK(IREP) = RTHICK(NLAY)
ZKOO(IREP) = ZKOO(IREP-1) + THICK(IREP-1) + RPDIST
IGEOT(IREP) = IGEOT(NLAY)
PGEOT(1,IREP) = PGEOT(1,NLAY)
PGEOT(2,IREP) = PGEOT(2,NLAY)
PGEOT(3,IREP) = PGEOT(3,NLAY)
PGEOT(4,IREP) = PGEOT(4,NLAY)
PGEOT(5,IREP) = PGEOT(5,NLAY)
ZMEAN(IREP) = ZMEAN(NLAY)
AMEAN(IREP) = AMEAN(NLAY)
RMEAN(IREP) = RMEAN(NLAY)
NLIMIT(IREP) = NLIMIT(NLAY)
NLUNIT(IREP) = NLUNIT(NLAY)
IOTYPE(IREP) = IOTYPE(NLAY)
360 CONTINUE
ZDUMMY = ZKOO(IE) - ZKOO(NLAY)
NLAY = IE
ZKOMAX = ZKOO(NLAY) + THICK(NLAY)
WRITE(SYSOUT,5035) NREP,ZDUMMY
5035 FORMAT(' LAYER ',I3,' TIME(S) REPEATED.'/' ',5X,
1 'TOTAL LENGTH IN CM',T45,1PE12.5)
GOTO 50
C ----------------------------------------------------------------------
C "COMPONENT" ALLE DATEN EINER SCHICHTKOMPONENTE
C ----------------------------------------------------------------------
400 IF ( SKIP ) GOTO 50
NELEMC = NELEMC + 1
CALL TBOUND ('I4','#ELEMENTS ',NELEMC,1,NELEM,'IN ')
NZCOMP(NELEMC) = IFIX ( WPARM(1) + 0.5 )
ACOMP(NELEMC) = ABS(WPARM(2))
NATOM (NELEMC) = IFIX ( WPARM(3) + 0.5 )
CALL TBOUND ('I4','Z ',NZCOMP(NELEMC),1,92,'IN ')
CALL TBOUND ('R4','AWT ',ACOMP(NELEMC), 1.0,240.0,'IN ')
CALL TBOUND ('I4','# OF ATOMS',NATOM (NELEMC),1,1000,'IN ')
IF ( TEST ) THEN
WRITE (SYSOUT,5041) NLAY, NLAYXX, NELEMC, NZCOMP(NELEMC),
1 ACOMP(NELEMC), NATOM(NELEMC)
5041 FORMAT(' ','LAYER#',I3,'-',I3,', COMPONENT#',I3/
1 ' ',5X,'Z = ',I3,', A = ',F10.3,', #ATOMS = ',I3)
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "BUILD_TAB" ERSTELLE DIE REICHWEITE-ENERGIE TABELLE
C ----------------------------------------------------------------------
500 IF ( SKIP ) GOTO 590
IF ( LAYON .AND. .NOT. TABOK .AND. ( NELEMC .EQ. NELEM )
1 .AND. MASSOK )
1THEN
c-kt Beginn der Änderung zur Parameterübergabe
PREC = ABS( WPARM(1) )
IF ( PREC .EQ. 0.0 ) THEN
PREC = 1.0E-04
ELSE
CALL TBOUND ('R4','PRECISION ',PREC,1.E-8,1.E-4,'IN')
ENDIF
c-kt Ende
SUMZ = 0.
SUMA = 0.
SUMR = 0.
NSUM = 0
DO 510 LL=1,NELEM
NSUM = NSUM + NATOM(LL)
SUMZ = SUMZ + ( NATOM(LL) *
1 FLOAT(NZCOMP(LL)) * FLOAT(NZCOMP(LL)+1) )
SUMA = SUMA + ( NATOM(LL) * ACOMP(LL) )
SUMR = SUMR + ( NATOM(LL) / RADLEN(NZCOMP(LL)) )
510 CONTINUE
ZMEAN (NLAY) = SQRT( SUMZ / FLOAT(NSUM) )
AMEAN (NLAY) = SUMA / FLOAT(NSUM)
RMEAN (NLAY) = FLOAT(NSUM) / SUMR
NRGSTA = NRGCUR + 1
NRG = NRGMAX - NRGCUR
c-kt Herauskommentieren des vorgefundenen Wertes:
c-kt PREC = 1.0E-04
CALL DGETAB ( BMASS, EKMIN, EKMAX, PREC,
1 NZCOMP, ACOMP, NATOM, NELEM,
2 RGDATA(1,NRGSTA), NRG, IFLAG )
IF ( IFLAG .NE. 0 ) THEN
c-kt Ein Hinweis auf den Grund des Absturzes ist immer angenehm ...
IF ( IFLAG .EQ. 5 ) THEN
c-kt Fortsetzung eines Hinweises, der von DGETAB bei IFLAG=5 gegeben wird:
WRITE (SYSOUT,5555) NLAY
5555 FORMAT(' Nummer ',I3)
ELSE IF ( IFLAG .EQ. 1 ) THEN
c-kt Hilfe bei der Fehlersuche:
WRITE (SYSOUT,5556)
5556 FORMAT(' Zu wenig Einträge in der '//
1 'Energie-Reichweite-Tabelle.'
2 /' Wahrscheinlich muß entweder für diese Schicht mit dem '//
3 'Befehl'
4 /'''BUILD_TAB'' ein Parameter kleiner 1.E-4 übergeben '//
5 'oder der '
6 /' 3. Parameter von ''PARTICLE'' vergrößert werden.')
ELSE
c-kt für die anderen Fälle habe ich keine Informationen bereit:
WRITE (SYSOUT,5557) IFLAG
5557 FORMAT(' ROUTINE ''DGETAB'' mit Fehler ',I1,' verlassen.')
ENDIF
c-kt Ende
STOP 12
ENDIF
ILAYS (NLAY) = NRGSTA
ILAYE (NLAY) = NRG
NRGCUR = NRGCUR + NRG
ELSE
WRITE(SYSOUT,5051)
5051 FORMAT(' PARTICLE, LAYER OR COMPONENT ARE INCOMPLETE.')
STOP 12
ENDIF
IF ( TEST ) THEN
WRITE (SYSOUT,5052) NLAY, NLAYXX, ILAYS(NLAY), ILAYE(NLAY),
1 NRGCUR, ZMEAN(NLAY), AMEAN(NLAY),
2 RMEAN(NLAY)
5052 FORMAT(' ','LAYER#',I3,'-',I3/
1 ' ',5X,'INDEX OF FIRST POINT IN (R,E)-TABLE',T45,I12/
2 ' ',5X,'NUMBER OF POINTS IN (R,E)-TABLE',T45,I12/
3 ' ',5X,'TOTAL NUMBER OF POINTS USED',T45,I12/
4 ' ',5X,'MEAN Z',T45,1PG12.5/
5 ' ',5X,'MEAN ATOMIC WEIGHT',T45,1PG12.5/
6 ' ',5X,'MEAN RADIATION LENGTH (G/CM**2)',T45,1PG12.5)
ENDIF
590 IF ( LAYGEO ) GOTO 350
TABOK = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "HOLE" DATEN FUER EIN STRAHLLOCH ODER EINE STRAHLBLENDE
C ----------------------------------------------------------------------
600 IF ( LAYON ) THEN
WRITE (SYSOUT,5061)
5061 FORMAT(' PREVIOUS LAYER OR HOLE-COMMAND STILL ACTIVE.')
STOP 12
ENDIF
NLAY = NLAY + 1
NGRP = NGRP + 1
CALL TBOUND ('I4','#LAYER ',NLAY ,1,200,'IN ')
CALL TBOUND ('I4','#GROUP ',NGRP ,1,200,'IN ')
THICK (NLAY) = WPARM (1)
CALL TBOUND ('R4','THICKNESS ',THICK(NLAY),1.E-08,1.E08,'IN')
NELEM = 0
RO (NLAY) = 0.0
RTHICK(NLAY) = 0.0
ILAYS(NLAY) = 0
ILAYE(NLAY) = 0
NREP = 0
RPDIST = 0.0
READ (SYSIN, 5171) COMENT(NGRP)
LAYON = .TRUE.
LAYGEO = .FALSE.
TABOK = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "OUTPUT" ALLE PARAMETER FUER DIE AUSGABE
C ----------------------------------------------------------------------
700 IF ( SKIP ) GOTO 50
IF ( NLAY .GE. 1 ) THEN
NLUNIT (NLAY) = IFIX ( WPARM(1) + 0.5 )
IOTYPE (NLAY) = IFIX ( WPARM(2) + 0.5 )
NLIMIT (NLAY) = IABS ( IFIX ( WPARM(3) + 0.5 ) )
ELIMIT (NLAY) = ABS ( WPARM(4) )
IF ( ELIMIT(NLAY) .LE. 0.0 ) ELIMIT(NLAY) = 1.0E37
CALL TBOUND ('I4','LOG UNIT# ',NLUNIT(NLAY) ,1,100,'IN ')
CALL TBOUND ('I4','OUT-TYPE ',IOTYPE(NLAY) ,1,7,'IN ')
IF ( TEST ) THEN
WRITE (SYSOUT,5071) NLAY, NLAYXX, NLUNIT(NLAY),
1 OUTTYP(IOTYPE(NLAY)),
2 IOTYPE(NLAY), NLIMIT(NLAY), ELIMIT(NLAY)
5071 FORMAT(' ','LAYER#',I3,'-',I3/
1 6X,'OUTPUT IS DIRECTED TO UNIT#',T45,I3/
2 6X,'OUTPUT-TYPE',T45,A10,'(',I1,')'/
3 6X,'MAXIMUM NUMBER OF EVENTS',T45,I12/
4 6X,'MAXIMUM ENERGY',T45,E12.5)
ENDIF
ELSE
WRITE (SYSOUT,5072)
5072 FORMAT(' NO LAYER OR HOLE-COMMAND ACTIVE OR GIVEN.')
STOP 12
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "COMPUTE" STARTE DIE EIGENTLICHE RECHNUNG
C ----------------------------------------------------------------------
800 IF ( .NOT.LAYON .AND. BEAMP .AND. SEEDOK .AND. (NLAY .GE. 1)
1 .AND. MASSOK .AND. BEAMXY )
2THEN
NEVCNT = IABS ( IFIX ( WPARM(1) + 0.5 ) )
NSVMAX = IABS ( IFIX ( WPARM(2) + 0.5 ) )
CALL TBOUND ('I4','#EVENTS ',NEVCNT ,1,50000000,'IN ')
CALL TBOUND ('I4','#THROUGH ',NSVMAX ,1,50000000,'IN ')
IZL=IZL+1
c CALL COMPUT (NEVCNT, NSVMAX, IZL, FORM)
CALL COMPUT (NEVCNT, NSVMAX, FORM) ! IZL removed, AH, 17-dec-1996
CALL COMSUM
ELSE
WRITE (SYSOUT,5081)
5081 FORMAT(' ','MASS, BEAM, LAYER OR SEED INCOMPLETE OR MISSING.')
STOP 12
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "EXIT" BEENDE DAS PROGRAMM
C ----------------------------------------------------------------------
c-kt Erweitert um die Ausgabe der Zeitpunkte des Programmstarts und -endes,
c-kt der Dauer der Rechnung und ggf. der Zahl der aus der Rechnung genommenen
c-kt Teilchen:
900 CALL TIME(endezeit)
CALL DATE(endedatum)
dauer = SECNDS(startsec)
IF ( fehler .GT. 0 ) THEN
WRITE (SYSOUT,901) fehler
901 FORMAT(' ',59('='),/T5,'Zahl der aus der Rechnung genommenen'//
1 ' Teilchen:',I6,/' ',59('='))
ENDIF
h = INT(dauer/3600.)
dauer = MOD(dauer,3600.)
m = INT(dauer/60.)
dauer = MOD(dauer,60.)
s = INT (dauer)
WRITE (SYSOUT,902) startzeit,startdatum,endezeit,endedatum
902 FORMAT(' ',//T5,'Programmstart:',3X,A8,2X,A9,/T5,
1 'Programmende:',4X,A8,2X,A9)
WRITE (SYSOUT,903) h,m,s
903 FORMAT(' ',/T5,'Zeitdifferenz:',3X,I2.2,':',I2.2,':',I2.2,/)
c-kt Ende
c-kt Schließen der Standard-Ausgabedatei
c CLOSE(SYSOUT)
c-kt Schließen der Standard-Eingabedatei
c CLOSE(SYSIN)
GOTO 9999
C ----------------------------------------------------------------------
C "RULER" SCHALTET DEN RULER EIN
C ----------------------------------------------------------------------
1000 RULER = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "NORULER" SCHALTET DEN RULER AUS
C ----------------------------------------------------------------------
1100 RULER = .FALSE.
GOTO 50
C ----------------------------------------------------------------------
C "TEST" SCHALTET DEN TEST EIN
C ----------------------------------------------------------------------
1200 TEST = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "NOTEST" SCHALTET DEN TEST AB
C ----------------------------------------------------------------------
1300 TEST = .FALSE.
GOTO 50
C ----------------------------------------------------------------------
C "SEEDS" LIEST DIE STARTWERTE DER ZUFALLSZAHLEN EIN
C ----------------------------------------------------------------------
1400 DO 1410 LL=1,6
SEEDS (LL) = WPARM (LL)
1410 CONTINUE
SEEDOK = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "PARTICLE" MASSE UND ENERGIE-BEREICH DES TEILCHENS
C ----------------------------------------------------------------------
1500 BMASS = ABS( WPARM(1) )
EKMIN = ABS( WPARM(2) )
EKMAX = ABS( WPARM(3) )
CALL TBOUND ('R4','MASS/MEV ',BMASS, 50., 10000., 'IN ')
BMASS2 = BMASS * BMASS
CALL TBOUND ('R4','EK_MIN/MEV',EKMIN, 1.0E-04, 0.1, 'IN ')
CALL TBOUND ('R4','EK_MAX/MEV',EKMAX, 1., 10000.,'IN ')
MASSOK = .TRUE.
IF ( TEST ) THEN
WRITE (SYSOUT,5151) BMASS, EKMIN, EKMAX
5151 FORMAT(' ','PARTICLE PARAMETER'/
1 ' ',5X,'MASS IN MEV/C**2',T45,1PG12.5/
2 ' ',5X,'MINIMUM ENERGY IN MEV',T45,1PG12.5/
3 ' ',5X,'MAXIMUM ENERGY IN MEV',T45,1PG12.5)
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "BEAM_XY" RAEUMLICHE STRUKTUR DES STRAHLS
C ----------------------------------------------------------------------
1600 IF ( BEAMXY ) THEN
WRITE (SYSOUT,5160)
5160 FORMAT(' SPATIAL DISTRIBUTION OF BEAM ALREADY DEFINED.')
STOP 12
ENDIF
XFWHM = ABS( WPARM(1) )
YFWHM = ABS( WPARM(2) )
RMAX = ABS( WPARM(3) )
ibtyp = ifix( wparm(4) + 0.5 ) ! require ibtyp as input, tp
RMAX2 = RMAX * RMAX
CALL TBOUND ('R4','X_FWHM ',XFWHM, 0.0, 1000., 'IN')
CALL TBOUND ('R4','Y-FWHM ',YFWHM, 0.0, 1000., 'IN')
CALL TBOUND ('R4','R_MAX ',RMAX, 0.0, 1000., 'IN')
call tbound ('I4','BEAM-PARA ',ibtyp, 1, 3, 'IN')
c-tp IF ( XFWHM .EQ. 0.0 .OR. YFWHM .EQ. 0.0 ) THEN
c-tp IBTYP = 1
c-tp XBEAM = 0.0
c-tp YBEAM = 0.0
c-tp ELSE
c-tp IBTYP = 2
c-tp XBEAM = XFWHM * 0.4246609
c-tp YBEAM = YFWHM * 0.4246609
c-tp ENDIF
if ( ibtyp .eq. 1 ) then
XBEAM = 0.0
YBEAM = 0.0
else if ( ibtyp .eq. 2 ) then
XBEAM = XFWHM * 0.4246609
YBEAM = YFWHM * 0.4246609
else if (ibtyp .eq. 3 ) then ! new option, TP, 1-jun-95
XBEAM = XFWHM
YBEAM = YFWHM
endif
c
BEAMXY = .TRUE.
c
IF ( TEST ) THEN
WRITE(SYSOUT,5161) BEAMTY(IBTYP),XFWHM,YFWHM,RMAX
5161 FORMAT(' ','BEAM SPATIAL DISTRIBUTION PARAMETER'/
1 6X,'BEAM-SPOT DISTRIBUTION',T47,A10/
2 6X,'X-FWHM OF THE BEAMSPOT IN CM',T45,1PE12.5/
3 6X,'Y-FWHM OF THE BEAMSPOT IN CM',T45,1PE12.5/
4 6X,'MAXIMUM RADIUS OF BEAMSPOT IN CM',T45,1PE12.5)
ENDIF
GOTO 50
C ----------------------------------------------------------------------
C "OLD_LAYER" DATEN EINER SCHICHT SCHON VORHANDEN
C ----------------------------------------------------------------------
1700 IF ( LAYON ) THEN
WRITE (SYSOUT,5021)
STOP 12
ENDIF
NLAY = NLAY + 1
NGRP = NGRP + 1
CALL TBOUND ('I4','#LAYER ',NLAY ,1,200,'IN ')
CALL TBOUND ('I4','#GROUP ',NGRP ,1,200,'IN ')
NLOLDG = IFIX ( WPARM (1) + 0.5 )
CALL TBOUND ('I4','#OLD_GROUP',NLOLDG ,1,(NGRP-1),'IN ')
NLOLD = IGRPS(NLOLDG)
CALL TBOUND ('I4','TAB_START ',ILAYS(NLOLD) ,1,NRGMAX,'IN ')
CALL TBOUND ('I4','TAB_LENGTH',ILAYE(NLOLD) ,1,NRGMAX,'IN ')
RO (NLAY) = WPARM (2)
THICK (NLAY) = WPARM (3)
NDIV = MAX0( IFIX( WPARM(4) + 0.5 ), 1 )
RPDIST = ABS( WPARM(5) )
NREP = NDIV - 1
THICK (NLAY) = THICK (NLAY) / FLOAT(NDIV)
READ (SYSIN,5171) COMENT(NGRP)
5171 FORMAT(A20)
CALL TBOUND ('R4','DENSITY ',RO(NLAY),1.E-08,100.,'IN')
CALL TBOUND ('R4','THICKNESS ',THICK(NLAY),1.E-08,1.E08,'IN')
CALL TBOUND ('I4','REPETITION',NREP,0,(200-NLAY),'IN ')
CALL TBOUND ('R4','DISTANCE ',RPDIST,0.0, 1.E06,'IN')
RTHICK(NLAY) = RO(NLAY) * THICK(NLAY)
NLAYXX = NLAY + NREP
IF ( TEST ) THEN
WRITE (SYSOUT,5170) NLAY, NLAYXX, COMENT(NLAY), NLOLD, RO(NLAY),
1 THICK(NLAY), RTHICK(NLAY), NREP, RPDIST
5170 FORMAT (' ','LAYER#',I3,'-',I3,' / ',A20/
1 ' ',5X,'ATOMIC COMPONENTS COPIED FROM LAYER NO',T45,I12/
2 ' ',5X,'DENSITY IN G/CM**3',T45,1PE12.5/
3 ' ',5X,'THICKNESS IN CM',T45,1PE12.5/
4 ' ',5X,'THICKNESS IN G/CM**2',T45,1PE12.5/
5 ' ',5X,'NUMBER OF ADDITIONAL LAYERS',T45,I12/
6 ' ',5X,'DISTANCE IN CM',T45,1PE12.5)
ENDIF
ILAYS(NLAY) = ILAYS(NLOLD)
ILAYE(NLAY) = ILAYE(NLOLD)
ZMEAN(NLAY) = ZMEAN(NLOLD)
AMEAN(NLAY) = AMEAN(NLOLD)
RMEAN(NLAY) = RMEAN(NLOLD)
TABOK = .TRUE.
LAYGEO = .FALSE.
LAYON = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "STRAG_ON" RANGE STRAGGLING EINSCHALTEN
C ----------------------------------------------------------------------
1800 STRAGG = .TRUE.
RELSTR = ABS( WPARM(1) )
CALL TBOUND ('R4','REL_STRAGG',RELSTR,0.0, 1.0,'IN')
GOTO 50
C ----------------------------------------------------------------------
C "STRAG_OFF" RANGE STRAGGLING AUSSCHALTEN
C ----------------------------------------------------------------------
1900 STRAGG = .FALSE.
GOTO 50
C ----------------------------------------------------------------------
C "SCATT_ON" MULTIPLE SCATTERING EINSCHALTEN
C ----------------------------------------------------------------------
2000 SCATTR = .TRUE.
GOTO 50
C ----------------------------------------------------------------------
C "SCATT_OFF" MULTIPLE SCATTERING AUSSCHALTEN
C ----------------------------------------------------------------------
2100 SCATTR = .FALSE.
GOTO 50
C ----------------------------------------------------------------------
C "SUMMARY" GEBE SUMMARY AUF EIN HILFSFILE AUS
C ----------------------------------------------------------------------
2200 ISUNIT = IFIX( WPARM(1) + 0.5 )
CALL SUMOUT (ISUNIT)
GOTO 50
C ----------------------------------------------------------------------
C "SKIP" UEBERSPRINGE DIE FOLGENDE GRUPPE
C ----------------------------------------------------------------------
2300 IF ( LAYON ) THEN
WRITE (SYSOUT,5021)
STOP 12
ENDIF
SKIP = .TRUE.
NLAY = NLAY + 1
NGRP = NGRP + 1
CALL TBOUND ('I4','#LAYER ',NLAY ,1,200,'IN ')
CALL TBOUND ('I4','#GROUP ',NGRP ,1,200,'IN ')
ILAYS(NLAY) = -1
NREP = 0
READ (SYSIN,5171) COMENT(NGRP)
LAYON = .TRUE.
LAYGEO = .FALSE.
TABOK = .TRUE.
GOTO 50
c-kt Anfang
c ----------------------------------------------------------------------
c "IMPULS" AUSGABEANWEISUNG FÜR DEN AKTUELLEN IMPULS
c ----------------------------------------------------------------------
2400 IF ( .NOT. SKIP ) THEN
IMPOUT(NLAY+1) = IFIX ( WPARM(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',IMPOUT(NLAY+1) ,1,99,'IN ')
ENDIF
GOTO 50
c ----------------------------------------------------------------------
c "WINKEL" AUSGABEANWEISUNG FÜR DEN AKTUELLEN STREUWINKEL
c ----------------------------------------------------------------------
2500 IF ( .NOT. SKIP ) THEN
THEOUT(NLAY+1) = IFIX ( WPARM(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',THEOUT(NLAY+1) ,1,99,'IN ')
ENDIF
GOTO 50
c ----------------------------------------------------------------------
c "STOPPED" AUSGABEANWEISUNG FÜR DIE STOPPVERTEILUNG IN DER SCHICHT
c ----------------------------------------------------------------------
2600 IF ( .NOT. SKIP ) THEN
STPOUT(NGRP) = IFIX ( WPARM(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',STPOUT(NGRP) ,1,99,'IN ')
ENDIF
GOTO 50
c ----------------------------------------------------------------------
c "FOCUS" ANGABEN ÜBER DIE FOKUSSIERUNG DES EINGANGSSTRAHLS
c ----------------------------------------------------------------------
2700 XFOCUS = WPARM(1)
YFOCUS = WPARM(2)
CALL TBOUND ('R4','X-FOCUS/CM ',ABS(XFOCUS), 0., 1000., 'IN ')
CALL TBOUND ('R4','Y-FOCUS/CM ',ABS(YFOCUS), 0., 1000., 'IN ')
GOTO 50
c ----------------------------------------------------------------------
c "WHERE_XY" AUSGABEANWEISUNG FÜR DIE X-Y-KOORDINATEN VOR DER SCHICHT
c ----------------------------------------------------------------------
2800 IF ( .NOT. SKIP ) THEN
XYOUT(NLAY+1) = IFIX ( WPARM(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',XYOUT(NLAY+1) ,1,99,'IN ')
ENDIF
GOTO 50
c ----------------------------------------------------------------------
c "MAGNETIC" Einschalten des Magnetfeldes
c ----------------------------------------------------------------------
2900 MSTART = WPARM(1)
MGAUSS = WPARM(2)
CALL TBOUND ('R4','START B-FELD / CM ',MSTART, 0., 120., 'IN ')
CALL TBOUND ('R4','STÄRKE B-FELD / GAUSS ',
1 MGAUSS, -1.0E+05, 1.0E+05, 'IN ')
GOTO 50
c ----------------------------------------------------------------------
c "ALL_XY" AUSGABEANWEISUNG FÜR DIE X-Y-KOORDINATEN VOR DER SCHICHT
c ----------------------------------------------------------------------
3000 IF ( .NOT. SKIP ) THEN
XYALL(NLAY+1) = IFIX ( WPARM(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',XYALL(NLAY+1) ,1,99,'IN ')
ENDIF
GOTO 50
c
c-kt Ende
c
c ----------------------------------------------------------------------
c "INIT_OUT" AUSGABEANWEISUNG FUER DIE STARTWERTE
c ----------------------------------------------------------------------
3100 IF ( .NOT. SKIP ) then
luninit = ifix ( wparm(1) + 0.5 )
CALL TBOUND ('I4','LOG UNIT# ',LUNINIT ,1,99,'IN ')
endif
goto 50
c ----------------------------------------------------------------------
c "FORMATTED" SET FLAG FOR FORMATTED OUTPUT
c-----------------------------------------------------------------------
3200 form = .true.
write(sysout,'(/,'' OUTPUT will be formatted !!!'',/)')
goto 50
c
c ----------------------------------------------------------------------
c "UFORMATTED" SET FLAG FOR UNFORMATTED OUTPUT
c-----------------------------------------------------------------------
3300 form = .false.
write(sysout,'(/,'' OUTPUT will be un-formatted !!!'',/)')
goto 50
c
c ======================================================================
c H I E R E N D E T D A S H A U P T P R O G R A M M
c ======================================================================
9999 END
C ----------------------------------------------------------------------
C
C NAME: COMPUT
C
C FUNKTION: BERECHNET DIE TRAJEKTORIEN DER TEILCHEN
C
C NMAX I*4 MAXIMALE ANZAHL VON START-TEILCHEN
C NSVMAX I*4 MAXIMALE ANZAHL VON DURCHLAUFENDEN T.
C
C ----------------------------------------------------------------------
c SUBROUTINE COMPUT (NMAX,NSVMAX,IZL,FORM)
SUBROUTINE COMPUT (NMAX,NSVMAX,FORM) ! IZL removed
CHARACTER*1 ICHAR, BACKSC, STOPPD
CHARACTER*20 COMENT
c-kt Hilfsvariable für die Berechnung der Fokussierung
REAL*4 xxf, yyf, temp
c-kt Variablen für MAGNETIC: Beginn und Stärke des Magnetfeldes [cm], [Gauss]
REAL*4 MSTART, MGAUSS
c-kt Hilfsvariable für die Berechnung der Bahn im Magnetfeld:
REAL*4 PXPY, radius, MZDIST, TMFLUG, omega, rotats, PHIXY,
1 XHELIX, YHELIX, PHIHELIXA, PHIHELIXB, PHIHMP,
2 DELTA, DIFA, DIFB, MAXRAD
REAL*8 RGDATA, D1, D2, D3
INTEGER*4 SEEDS
c-kt Zählvariable für die aus der Rechnung genommenen Teilchen:
INTEGER*4 fehler
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY und ALL_XY:
INTEGER*2 IMPOUT, THEOUT, STPOUT, XYOUT, XYALL, LUNINIT
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
INTEGER*2 SYSOUT, SYSIN
INTEGER*4 TEST
LOGICAL SCATTR, STRAGG, STRAG2, FORM
COMMON
1 / EXECFL / SCATTR, STRAGG, RELSTR
COMMON
1 / ATOMIC / RADLEN(92)
COMMON
1 / RETAB / NRGMAX, NRGCUR, RGDATA (5,3000)
2 / LAYER / NLAY, ILAYS(200), ILAYE(200),
3 THICK(200), RTHICK(200), RO(200),
4 ZMEAN(200), AMEAN(200), RMEAN(200),
5 ZKOO(200), IGEOT(200), PGEOT(5,200)
6 / LAYERG / NGRP, IGRPS(200), IGRPE(200)
6 / LAYERC / COMENT(200)
7 / PARTIC / BMASS, BMASS2, EKMIN, EKMAX, IBFLAG, IBUNIT
9 / BEAM / PMEAN, PFWHM, PLOW, PHIGH, XFWHM, YFWHM,
1 RMAX,
2 IBTYP, XBEAM, YBEAM, PSIGMA, RMAX2
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
3 , XFOCUS, YFOCUS
c-kt Beginn und Stärke des Magnetfeldes
4 , MSTART, MGAUSS
9 / KINEMA / ID, X, Y, Z,
1 PX, PY, PZ, PTOTAL,
c-kt Einführen der Variable TTHETA für den Tangens des Streuwinkels:
2 THETA, PHI, STHETA, CTHETA, TTHETA,
3 SPHI, CPHI,
4 TTOTAL, TLAST,
5 ELOSS
COMMON
2 / ZUFALL / R(6)
COMMON
1 / OUTPUT / NLIMIT(200), NCOUNT(200), IOTYPE(200),
8 NLUNIT(200), ELIMIT(200)
3 / COMPON / NZCOMP(20), ACOMP(20), NATOM(20), NELEM, NELEMC
4 / SEEDCO / SEEDS (6)
5 / INFOCO / NPARTI, NTHROU, NSTOP(200), NSTRVO(200),
6 NSTRIN(200), NRUECK(200), NSTRGE(200), NMLRE(200)
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
COMMON
1 / STATUS / SYSOUT, SYSIN, TEST
c 1 / STATUS / SYSOUT
c-kt Zählvariable für die aus der Rechnung genommenen Teilchen:
COMMON
1 / COUNT / fehler
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY, ALL_XY und INIT_OUT:
COMMON
1 / SGLOUT / IMPOUT(201), THEOUT(201), STPOUT(200), XYOUT(201),
2 XYALL(201), LUNINIT
c DATA ZWEIPI /6.283185/, VLIGHT /29.997924/, SQDREI/1.732051/,
c 1 PIHALB /1.570796/
c
c PIHALB is not used, commented out, AH 17-dec-1996
c
DATA ZWEIPI /6.283185/, VLIGHT /29.997924/, SQDREI/1.732051/
DATA BACKSC /'B'/, STOPPD /'S'/
c do i = 1, 1000
c write(sysout,'(5f10.4)') rgdata(1,i),rgdata(2,i),rgdata(3,i),
c 1 rgdata(4,i),rgdata(5,i)
c
c enddo
D2 = BMASS
D3 = D2 * D2
NPARTI = 0
NTHROU = 0
DO 5 I=1,200
NCOUNT(I) = 0
NSTOP(I) = 0
NSTRVO(I) = 0
NSTRIN(I) = 0
NRUECK(I) = 0
NSTRGE(I) = 0
NMLRE (I) = 0
5 CONTINUE
DO 6 I=1,NLAY
IF ( ILAYS(I) .GT. 0 ) GOTO 7
6 CONTINUE
STRAGG = .FALSE.
7 NFIRST = I
WRITE(SYSOUT,2204) SEEDS(1)
2204 FORMAT(' START VALUE OF RANDOM NUMBER = ',I16)
WRITE(SYSOUT,1800)
1800 FORMAT(' CALCULATION STARTED')
IF ( STRAGG ) THEN
WRITE(SYSOUT,1801) RELSTR
1801 FORMAT(' ',5X,'RANGE STRAGGLING ENABLED (',1PE12.5,')')
IF ( RELSTR .EQ. 0.0 ) THEN
STRAG2 = .TRUE.
ELSE
STRAG2 = .FALSE.
c-kt IRELST ist nur für die Ausgabe mit I3 in 'OUTPUT' definiert;
c-kt es scheint mir sinnvoller, dort RELSTR direkt mit (E10.4) auszugeben.
c-kt IRELST = IFIX( RELSTR * 1000. )
ENDIF
ENDIF
IF ( SCATTR ) THEN
WRITE(SYSOUT,1802)
1802 FORMAT(' ',5X,'MULTIPLE SCATTERING ENABLED.')
ENDIF
NSV = 0
ID = 0
DO 100 I=1,NMAX
C WUERFLE STRAHL-TEILCHEN
CALL BEAMIN (NSV)
EBEGIN = SQRT( PTOTAL*PTOTAL + BMASS2 ) - BMASS
IF ( STRAGG ) THEN
C MERKE ANFANGSENERGIE UND -REICHW.
RBEGIN = DRANGE ( EBEGIN, RGDATA(1,ILAYS(NFIRST)),
1 ILAYE(NFIRST) )
c-kt Wenn 'DRANGE' einen negativen Wert übergibt, hatte ein Teilchen einen
c-kt höheren Wert für die kinetische Energie als der größte tabellierte
c-kt Energiewert:
IF (RBEGIN .LT. 0.0) THEN
WRITE (SYSOUT,181) NSV
181 FORMAT(' ',T5,'Startnummer des Teilchens',T50,I22,/T5,
1 'Der 3. Parameter des Befehls PARTICLE muß vergrößert '//
2 'werden.')
STOP 12
ENDIF
c-kt Ende
IF ( STRAG2) THEN
RELSTR = STRAG1( PTOTAL, BMASS )
c-kt IRELST ist nur für die Ausgabe mit I3 in 'OUTPUT' definiert;
c-kt es scheint mir sinnvoller, dort RELSTR direkt mit (E10.4) auszugeben.
c-kt IRELST = IFIX( RELSTR * 1000. )
ENDIF
ENDIF
c-kt Wenn mindestens eine FOCUS-Konstante ungleich Null ist, fokussiere
c-kt den Strahl, sonst setze die üblichen Anfangswerte
IF ( XFOCUS .NE. 0.0 .OR. YFOCUS .NE. 0.0 ) THEN
IF ( XFOCUS .NE. 0.0 ) THEN
xxf = X / XFOCUS
ELSE
xxf = 0.0
ENDIF
IF ( YFOCUS .NE. 0.0 ) THEN
yyf = Y / YFOCUS
ELSE
yyf = 0.0
ENDIF
temp = xxf**2 + yyf**2
CTHETA = 1.0 / SQRT( 1.0 + temp )
STHETA = SQRT( temp / ( 1.0 + temp ) )
TTHETA = SQRT( temp )
THETA = ATAN( TTHETA )
IF ( yyf .NE. 0.0 .OR. xxf .NE. 0.0) THEN
PHI = ATAN2( -yyf, -xxf )
ELSE
PHI = 0.0
ENDIF
CPHI = COS( PHI )
SPHI = SIN( PHI )
PZ = PTOTAL * CTHETA
PX = - PZ * xxf
PY = - PZ * yyf
ELSE
PHI = 0.0
THETA = 0.0
CTHETA = 1.0
STHETA = 0.0
TTHETA = 0.0
CPHI = 1.0
SPHI = 0.0
PZ = PTOTAL
PX = 0.0
PY = 0.0
ENDIF
c-kt Ende
TTOTAL = 0.0
ID = ID + 1
c
c-tp write initial values to file if desired
c
il = 0
if ( luninit .gt. 0 ) then
if ( form ) then
write(luninit,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
1 ebegin,theta,phi,ttotal
else
c
c-tp unformatted output
c
write(luninit) id,il,x,y,z,px,py,pz,ebegin,theta,phi,ttotal
c
endif
endif
c
C ----------------------------------------------------------------------
C TRACK THROUGH LAYERS AND HOLES
C ----------------------------------------------------------------------
DO 200 IL=1,NLAY
IF ( ILAYS(IL) .LT. 0 ) GOTO 200
c-kt (zur Erklärung dieser Zeile: ILAYS ist genau dann negativ, wenn die
c-kt betreffende Schicht mit SKIP aus der Rechnung genommen wurde.)
c-kt Gegebenenfalls Speichern des Impulsbetrages oder des Streuwinkels
IF ( IMPOUT(IL) .NE. 0 ) THEN
WRITE (IMPOUT(IL),201) PTOTAL
201 FORMAT(' ',E13.7)
ENDIF
IF ( THEOUT(IL) .NE. 0 ) THEN
WRITE (THEOUT(IL),202) THETA
202 FORMAT(' ',F13.7)
ENDIF
c-kt Ende
C BERECHNE EINTRITTSPUNKT
ZDIST = ZKOO(IL) - Z
IF ( ZDIST .GT. 0.01 ) THEN
c-kt Wenn der Befehl MAGNETIC gegeben wurde, definiere ein Magnetfeld
c-kt ================================================================
IF ( ( ABS(MGAUSS) .GT. 0.0 ) .AND.
1 ( ZKOO(IL) .GT. MSTART ) ) THEN
c-kt Die Magnetfeldlinien zeigen bei positivem MGAUSS in Strahlrichtung, bei
c-kt negativem entgegen der Strahlrichtung. Demnach ist die Bahn der mu+ eine
c-kt links- bzw. rechtsdrehende Helix.
c-kt Phi wird von der positiven X-Achse aus nach links und rechts positiv bzw.
c-kt negativ gezählt. Falls zwischen PHI und PHINEU die negative X-Achse über-
c-kt quert wurde, fand ein Vorzeichenwechsel statt.
c-kt Wenn das Magnetfeld erst nach dem vorigen Bauteil begonnen hat, dann be-
c-kt rechne die X,Y-Koordinaten und ihre geometr. Summme zu Beginn des Feldes
IF ( Z .LT. MSTART ) THEN
XYPROJ = TTHETA * ( MSTART - Z )
X = CPHI * XYPROJ + X
Y = SPHI * XYPROJ + Y
ENDIF
c-kt Gyrationsradius im Magnetfeld [cm]:
PXPY = SQRT( PX**2 + PY**2 )
radius = PXPY * 3335.668 / ABS( MGAUSS )
c-kt Im Magnetfeld zurückgelegte Strecke [cm]:
MZDIST = ZKOO(IL) - AMAX1( Z, MSTART )
c-kt Flugzeit im Magnetfeld [s]:
TMFLUG = MZDIST * SQRT( PZ*PZ + BMASS2 ) /
1 ( PZ * 0.29979246E+11 )
c-kt Winkelgeschwindigkeit [1/s]:
omega = ABS(MGAUSS) * 0.89874789E+07 / BMASS
c-kt Zahl der bis zum Ort ZKOO(IL) erfolgten Rotationen:
rotats = omega * TMFLUG / 6.2831853
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
WRITE(IL+20,4321) radius, rotats
4321 FORMAT(' ',E13.7,1X,E13.7)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c-kt Richtungswinkel von X,Y zur Helixmittelachse
PHIXY = PHI + SIGN( 1.0, MGAUSS ) * 1.5707963
IF ( PHIXY .GT. 3.1415927 ) PHIXY = PHIXY - 6.2831853
IF ( PHIXY .LT. -3.1415927 ) PHIXY = PHIXY + 6.2831853
c-kt Berechnung der X,Y-Koordinaten der Mittelachse der Helix
XHELIX = X + radius * COS(PHIXY)
YHELIX = Y + radius * SIN(PHIXY)
c-kt Richtungswinkel von der Helixmittelachse zu X,Y
PHIHELIXA = PHIXY + 3.1415927
IF ( PHIHELIXA .GT. 3.1415927 )
1 PHIHELIXA = PHIHELIXA - 6.2831853
c-kt Bei ZKOO(IL) angekommen beträgt der Winkel vom Helixmittelpunkt zum
c-kt Ort des Teilchens PHIHELIXB (Vollständige Rotationen sind belanglos):
PHIHELIXB = PHIHELIXA + SIGN( 1.0, MGAUSS ) *
1 AMOD( rotats, 1.0 ) * 6.28318531
IF ( PHIHELIXB .GT. 3.1415927 )
1 PHIHELIXB = PHIHELIXB - 6.2831853
IF ( PHIHELIXB .LT. -3.1415927 )
1 PHIHELIXB = PHIHELIXB + 6.2831853
c-kt Berechnung der X,Y-Koordinate an ZKOO(IL)
X = XHELIX + radius * COS( PHIHELIXB )
Y = YHELIX + radius * SIN( PHIHELIXB )
c-kt Berechnung von XYPROJ
XYPROJ = SQRT( X**2 + Y**2 )
c-kt Zuweisung von Z:
Z = ZKOO(IL)
c-kt Berechnung des neuen Phasenwinkels Phi an ZKOO(IL)
IF ( Y .NE. 0.0 .OR. X .NE. 0.0) THEN
PHI = ATAN2( Y, X )
ELSE
PHI = 0.0
ENDIF
c-kt Berechnung der Sinus- und Cosinuswerte
CPHI = COS(PHI)
SPHI = SIN(PHI)
c-kt Damit ergeben sich auch PX und PY neu:
PX = PXPY * CPHI
PY = PXPY * SPHI
c-kt Ende der Berechnung des Magnetfeldes.
c-kt =====================================
ELSE
c-kt Variable für TAN(THETA) an die Stelle von STHETA/CTHETA gesetzt:
XYPROJ = TTHETA * ZDIST
X = CPHI * XYPROJ + X
Y = SPHI * XYPROJ + Y
Z = ZKOO(IL)
ENDIF
C BERECHNE FLUGZEIT
BETAZ = PZ / SQRT( PZ*PZ + BMASS2 )
TLAST = ZDIST / (BETAZ * VLIGHT)
TTOTAL = TTOTAL + TLAST
ENDIF ! Dieses ENDIF gehört zu der Z-Distanz<0.01cm-Abfrage
ELOSS = 0.0
c-kt Gegebenenfalls Speichern der X-Y-Koordinaten (aller Teilchen)
IF ( XYALL(IL) .NE. 0 ) THEN
WRITE (XYALL(IL),'(2I6,2F10.5)') ID, IL, X, Y
ENDIF
c-kt Ende
C BERECHNE, OB TEILCHEN TRIFFT
c-kt Berechnung, ob die ermittelten X,Y-Koordinaten in den Grenzen des
c-kt Bauteiles liegen:
GOTO ( 210, 220 ), IGEOT(IL)
C RUNDE SCHICHT
210 IF ( (X*X+Y*Y) .GT. PGEOT(2,IL) ) GOTO 350
GOTO 230
C RECHTECK
220 IF ( X .LT. PGEOT(1,IL) .OR. X .GT. PGEOT(2,IL) ) GOTO 350
IF ( Y .LT. PGEOT(3,IL) .OR. Y .GT. PGEOT(4,IL) ) GOTO 350
230 CONTINUE
c-kt Wenn die Rechnung mit einem Magnetfeld durchgeführt wird, muß der
c-kt Versatz in der X,Y-Ebene durch die Helixbahn berücksichtigt werden.
IF ( ( ABS(MGAUSS) .GT. 0.0 ) .AND.
1 ( ZKOO(IL) .GT. MSTART ) .AND.
2 ( ZDIST .GT. 0.01 ) ) THEN
c-kt Hat das Teilchen auf der Helixbahn das Strahlrohr berührt?
c-kt ==========================================================
c-kt Berechnung des Winkels von der X,Y=0-Achse zur Helixmittelachse. Die
c-kt Helixbahn hat in diesem Winkel den größten Abstand zur X,Y=0-Achse
IF ( YHELIX .NE. 0.0 .OR. XHELIX .NE. 0.0) THEN
PHIHMP = ATAN2( YHELIX, XHELIX )
c-kt Sollte die Helixmittelachse wirklich auf X,Y=0 liegen, ist auch der
c-kt Radius der Teilchenbahn konstant.
ELSE
GOTO 234
ENDIF
c-kt Welcher der beim Flug des Teilchens auf dem Segment der Helixbahn
c-kt überstrichenen Winkel kommt PHIHMP am nächsten und bezeichnet somit
c-kt den am weitesten von der X,Y=0-Achse entfernten Punkt?
c-kt Berechnung der Differenz DELTA dieses Winkels zu PHIHMP
c-kt 1. Es gab mindestens eine vollständige Umdrehung
IF ( rotats .GE. 1.0 ) THEN
DELTA = 0.0
c-kt 2. Es gab keine vollständige Umdrehung
ELSE
c-kt a) Kein Vorzeichenwechsel zwischen PHIHELIXA und PHIHELIXB
IF ( SIGN(1.0,MGAUSS)*PHIHELIXB .GT.
1 SIGN(1.0,MGAUSS)*PHIHELIXA ) THEN
c-kt In dem durchflogenen Helixsegment wurde PHIHMP überstrichen
IF ( SIGN(1.0,MGAUSS)*PHIHELIXB .GE.
1 SIGN(1.0,MGAUSS)*PHIHMP .AND.
2 SIGN(1.0,MGAUSS)*PHIHMP .GE.
3 SIGN(1.0,MGAUSS)*PHIHELIXA ) THEN
DELTA = 0.0
c-kt PHIHMP wurde nicht überstrichen
ELSE
c-kt Welcher der beiden Winkel PHIHELIXA und PHIHELIXB ist PHIHMP näher?
DIFA = ABS( PHIHELIXA - PHIHMP )
IF ( DIFA .GT. 3.14159265 )
1 DIFA = 6.28318531 - DIFA
DIFB = ABS( PHIHELIXB - PHIHMP )
IF ( DIFB .GT. 3.14159265 )
1 DIFB = 6.28318531 - DIFB
IF ( DIFA .LT. DIFB ) THEN
DELTA = DIFA
ELSE
DELTA = DIFB
ENDIF
ENDIF ! PHIHMP überstrichen?
c-kt b) Vorzeichenwechsel zwischen PHIHELIXA und PHIHELIXB
ELSE
c-kt In dem durchflogenen Helixsegment wurde PHIHMP überstrichen
IF ( SIGN(1.0,MGAUSS)*PHIHMP .GE.
1 SIGN(1.0,MGAUSS)*PHIHELIXA .OR.
2 SIGN(1.0,MGAUSS)*PHIHMP .LE.
3 SIGN(1.0,MGAUSS)*PHIHELIXB ) THEN
DELTA = 0.0
c-kt PHIHMP wurde nicht überstrichen
ELSE
c-kt Welcher der beiden Winkel PHIHELIXA und PHIHELIXB ist PHIHMP näher?
DIFA = ABS( PHIHELIXA - PHIHMP )
IF ( DIFA .GT. 3.14159265 )
1 DIFA = 6.28318531 - DIFA
DIFB = ABS( PHIHELIXB - PHIHMP )
IF ( DIFB .GT. 3.14159265 )
1 DIFB = 6.28318531 - DIFB
IF ( DIFA .LT. DIFB ) THEN
DELTA = DIFA
ELSE
DELTA = DIFB
ENDIF
ENDIF ! PHIHMP überstrichen?
ENDIF ! Vorzeichenwechsel von PHIHELIXA auf PHIHELIXB?
ENDIF ! Anzahl Rotationen > 1 ?
c-kt Auswertung: wenn der maximal erreichte Abstand von der X,Y=0-Achse den
c-kt Innenradius des Strahlrohres übersteigt, denn zähle das Teilchen als
c-kt herausgestreut.
IF ( Z .LT. 140.0 ) THEN
MAXRAD = 40.768225
ELSE
MAXRAD = 33.64
ENDIF
IF ( ( SQRT(XHELIX**2+YHELIX**2)
1 + radius * COS( DELTA ) )**2 +
2 ( radius * SIN(DELTA) )**2 .GE. MAXRAD )
2 GOTO 350
ENDIF ! Ist ein Magnetfeld eingeschaltet?
c-kt Ende Berechnung, ob das Teilchen duch das Magnetfeld abgelenkt wird
c-kt ===================================================================
C TEILCHEN KOMMT DURCH
234 CONTINUE
c-kt Gegebenenfalls Speichern der X-Y-Koordinaten (nur von nicht weggestreuten
c-kt Teilchen)
IF ( XYOUT(IL) .NE. 0 ) THEN
WRITE (XYOUT(IL),'(2I6,2F10.5)') ID, IL,X, Y
ENDIF
c-kt Ende
GOTO ( 250, 240, 250, 250, 250, 240, 240 ), IOTYPE(IL)
C DRUCKE 'IN FRONT'-DATEN AUS
240 IF ( EKOUT .GT. ELIMIT(IL) ) GOTO 250
IUNIT = NLUNIT(IL)
c-kt RELSTR wird anstelle von IRELST (=RELSTR*1000) ausgegeben
c WRITE (IUNIT,2100) ID,IL,SCATTR,STRAGG,STRAG2,RELSTR,NLAY,
c 1 EBEGIN,BMASS,X,Y,Z, PX,PY,PZ,THETA,PHI,
c 2 TTOTAL, TLAST, ELOSS
c 2100 FORMAT(' V',I10,2X,I3,2X,3L1,2X,E10.4,1X,I4,1X,1PE12.5 /
c 1 4X,6(1PE12.5) / 4X,6(1PE12.5) )
c
c new output, TP, 1-jun-95
c
e_current = sqrt( px**2 + py**2 + pz**2 + bmass2 ) - bmass
if ( form ) then
write(iunit,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
1 e_current,theta,phi,ttotal
else
c
c-tp unformatted output
c
write(iunit) id,il,x,y,z,px,py,pz,e_current,theta,phi,
1 ttotal
c
endif
NCOUNT(IL) = NCOUNT(IL) + 1
IF ( NCOUNT(IL) .GE. NLIMIT(IL) ) GOTO 420
C SCHICHT ODER LOCH ?
250 IF ( ILAYS(IL) .EQ. 0 ) THEN
c-kt Variable für TAN(THETA) an die Stelle von STHETA/CTHETA gesetzt:
XYPROJ = TTHETA * THICK(IL)
X = CPHI * XYPROJ + X
Y = SPHI * XYPROJ + Y
Z = THICK(IL) + Z
C BERECHNE FLUGZEIT
BETAZ = PZ / SQRT( PZ*PZ + BMASS2 )
TLAST = ZDIST / (BETAZ * VLIGHT)
TTOTAL = TTOTAL + TLAST
ELSE
C SCHICHTDICKE / COS(THETA)
CRANGE = RTHICK (IL) / CTHETA
CTHICK = THICK (IL) / CTHETA
C ENDENERGIE DES TEILCHENS
D1 = PTOTAL
EKIN = DSQRT ( D1 * D1 + D3 ) - D2
ELOSS = ELOSS2 ( EKIN, CRANGE, RGDATA(1,ILAYS(IL)),
1 ILAYE(IL) )
c-kt Wenn ELOSS2 einen negativen Wert übergibt ist einer von drei möglichen
c-kt Fehlern aufgetreten, die hier unterschieden werden:
IF ( ELOSS .LT. 0.0 ) THEN
WRITE (SYSOUT,251) NSV
251 FORMAT(' ',T5,'Startnummer des Teilchens',T50,I22)
IF ( ELOSS .LT. -2.5 ) THEN
c-kt Durch Rundungsfehler wurde keine korrekte Lösung des kubischen
c-kt Spline-Polynoms gefunden:
fehler = fehler + 1
c-kt zu welcher Gruppe gehört die Schicht?
DO 252 LGRP=1,NGRP
IF ( IL .LE. IGRPE(LGRP) ) GOTO 253
252 CONTINUE
253 WRITE (SYSOUT,254) IL, LGRP, COMENT(LGRP)
254 FORMAT(' ',T5,'Nummer der Schicht',T50,I22,/T5,
1 'Nummer der Schichtgruppe',T50,I22,/T5,
2 'Bezeichnung der Schichtgruppe',T52,A20,/T5,
3 'Das Teilchen wurde in der angegebenen Schicht'//
4 ' aus der Rechnung genommen.')
GOTO 100
ELSEIF ( ELOSS .LT. -1.5 ) THEN
c-kt Ein Teilchen hatte einen kleineren Wert für die verbleibende Reichweite
c-kt als der kleinste tabellierte Reichweiten-Wert:
WRITE (SYSOUT,255)
255 FORMAT(' ',T5,'Der 2. Parameter des Befehls PARTICLE '//
1 ' muß verkleinert werden.')
STOP 11
ELSE
c-kt Ein Teilchen hatte einen höheren Wert für die kinetische Energie
c-kt als der größte tabellierte Energiewert:
WRITE (SYSOUT,256)
256 FORMAT(' ',T5,'Der 3. Parameter des Befehls PARTICLE '//
1 ' muß vergrößert werden.')
STOP 12
ENDIF
ENDIF
c-kt Ende
EKOUT = EKIN - ELOSS
IF ( STRAGG ) THEN
C BERECHNE STRAGGLING
RSTRAG = RELSTR * SQRT( ELOSS / EBEGIN )
CALL GGNML ( SEEDS, 1, R(4) )
ERANGE = CRANGE + ( RSTRAG * R(4) * RBEGIN )
IF ( ERANGE .LE. 0.0 ) THEN
EKOUT = EKIN
ELOSS = 0.0
POUT = PTOTAL
NSTRGE(IL) = NSTRGE(IL) + 1
ELSE
ELOSS = ELOSS2 ( EKIN, ERANGE, RGDATA(1,ILAYS(IL)),
1 ILAYE(IL) )
EKOUT = EKIN - ELOSS
IF ( EKOUT .LE. 0.0 ) GOTO 400
POUT = EKOUT * SQRT( 1.0 + (2.0*BMASS/EKOUT) )
ENDIF
ELSE
IF ( EKOUT .LE. 0.0 ) GOTO 400
POUT = EKOUT * SQRT( 1.0 + (2.0*BMASS/EKOUT) )
ENDIF
C AUSTRITTSKOORD. OHNE STREUUNG
c-kt Variable für TAN(THETA) an die Stelle von STHETA/CTHETA gesetzt:
XYPROJ = TTHETA * THICK(IL)
XAUS = CPHI * XYPROJ + X
YAUS = SPHI * XYPROJ + Y
ZAUS = ZKOO(IL) + THICK(IL)
IF ( SCATTR ) THEN
IF ( (ELOSS/EKIN) .GE. 0.10 ) NMLRE(IL) = NMLRE(IL) + 1
C BERECHNE STREUWINKEL THETA
CALL MLR ( ZMEAN(IL), AMEAN(IL), RO(IL), CTHICK,
1 PTOTAL, BMASS, 1.0, THETAS )
C WUERFLE PHI 0-2PI
CALL GGUBS (SEEDS,1,R)
PHIS = R(1) * ZWEIPI
C BERECHNE SIN/COS-WERTE
STHETS = SIN( THETAS )
CTHETS = COS( THETAS )
SPHIS = SIN( PHIS )
CPHIS = COS( PHIS )
C (XAUS,YAUS)-VERSCHIEBUNG
C (SIN(X) = BOGENMASS !)
DXYPRJ = THETAS * CRANGE / SQDREI
X = CPHIS * DXYPRJ + XAUS
Y = SPHIS * DXYPRJ + YAUS
Z = ZAUS
ELSE
X = XAUS
Y = YAUS
Z = ZAUS
STHETS = 0.0
CTHETS = 1.0
SPHIS = 0.0
CPHIS = 1.0
ENDIF
ENDIF
C BERECHNE, OB TEILCHEN TRIFFT
GOTO ( 270, 280 ), IGEOT(IL)
C RUNDE SCHICHT
270 IF ( (X*X+Y*Y) .GT. PGEOT(2,IL) ) GOTO 360
GOTO 290
C RECHTECK
280 IF ( X .LT. PGEOT(1,IL) .OR. X .GT. PGEOT(2,IL) ) GOTO 360
IF ( Y .LT. PGEOT(3,IL) .OR. Y .GT. PGEOT(4,IL) ) GOTO 360
C IMPULSKOMPONENTEN BZGL DES EIGENKOORDINATENSYSTEMS
290 IF ( ILAYS(IL) .EQ. 0 ) GOTO 200
PTOTAL = POUT
PXEI = PTOTAL * CPHIS * STHETS
PYEI = PTOTAL * SPHIS * STHETS
PZEI = PTOTAL * CTHETS
C TRANSFORMIERE INS LABORSYSTEM
PX = PXEI*CPHI*CTHETA - PYEI*SPHI + PZEI*CPHI*STHETA
PY = PXEI*SPHI*CTHETA + PYEI*CPHI + PZEI*SPHI*STHETA
PZ = -PXEI*STHETA + PZEI *CTHETA
c-kt PTOTAL = SQRT(PX*PX + PY*PY + PZ*PZ) !ADDED 22.1.87 HJM, KAW!!!
c-kt Geändert zu höherer Genauigkeit bei der späteren Winkelberechnung
temp = PX*PX + PY*PY + PZ*PZ
PTOTAL = SQRT( temp )
IF ( PZ .LE. 0.0 ) GOTO 370
C BERECHNE NEUE RICHTUNGSWINKEL
IF ( SCATTR ) THEN
CTHETA = PZ / PTOTAL
c-lz 6-dec-91 changed the following sentence to avoid any negative sqrt
c-lz original STHETA = SQRT( 1.0 - CTHETA*CTHETA )
c-kt STHETA = SQRT(amax1(0.0,(1.0 - CTHETA*CTHETA )))
c-kt Geändert zu höherer Genauigkeit bei der Berechnung:
STHETA = SQRT( ( PX*PX + PY*PY ) / temp )
c-kt Berechnung der neu eingeführten Variable für TAN(THETA):
TTHETA = SQRT( ( PX*PX + PY*PY ) / (PZ*PZ) )
THETA = ASIN( STHETA )
c-kt Zugegeben ein unwahrscheinlicher Fall, aber dennoch korrekterweise:
IF ( PY .NE. 0.0 .OR. PX .NE. 0.0) THEN
PHI = ATAN2( PY, PX )
ELSE
PHI = 0.0
ENDIF
CPHI = COS( PHI )
SPHI = SIN( PHI )
ENDIF
GOTO ( 200, 200, 200, 300, 300, 200, 300 ), IOTYPE(IL)
C DRUCKE 'BEHIND'-DATEN AUS
300 IF ( EKOUT .GT. ELIMIT(IL) ) GOTO 200
IUNIT = NLUNIT(IL)
c-kt RELSTR wird anstelle von IRELST (=RELSTR*1000) ausgegeben
c WRITE (IUNIT,2101) ID,IL,SCATTR,STRAGG,STRAG2,RELSTR,NLAY,
c 1 EBEGIN,BMASS,X,Y,Z, PX,PY,PZ,THETA,PHI,
c 2 TTOTAL, TLAST, ELOSS
c 2101 FORMAT(' N',I10,2X,I3,2X,3L1,2X,E10.4,1X,I4,1X,1PE12.5 /
c 1 4X,6(1PE12.5) / 4X,6(1PE12.5) )
c
c new output, TP, 1-jun-95
c
e_current = sqrt( px**2 + py**2 + pz**2 + bmass2 ) - bmass
if ( form ) then
write(iunit,'(I6,2x,I3,10f10.3)') id,il,x,y,z,px,py,pz,
1 e_current,theta,phi,ttotal
else
c
c-tp unformatted output
c
write(iunit) id,il,x,y,z,px,py,pz,e_current,theta,phi,
1 ttotal
c
endif
NCOUNT(IL) = NCOUNT(IL) + 1
IF ( NCOUNT(IL) .GE. NLIMIT(IL) ) GOTO 420
200 CONTINUE
c-kt Gegebenenfalls Speichern des Impulsbetrages oder des Streuwinkels
c-kt nach der letzten Schicht
IF ( IMPOUT(NLAY+1) .NE. 0 ) THEN
WRITE (IMPOUT(NLAY+1),207) PTOTAL
207 FORMAT(' ',E13.7)
ENDIF
IF ( THEOUT(NLAY+1) .NE. 0 ) THEN
WRITE (THEOUT(NLAY+1),208) THETA
208 FORMAT(' ',E13.7)
ENDIF
c-kt Ende
C T. DURCHGEKOMMEN
NTHROU = NTHROU + 1
IF ( NTHROU .GE. NSVMAX ) GOTO 440
GOTO 100
C T. HERAUSGESTREUT VOR DER SCH.
350 NSTRVO(IL) = NSTRVO(IL) + 1
GOTO 100
C T. HERAUSGESTREUT IN DER SCH.
360 NSTRIN(IL) = NSTRIN(IL) + 1
GOTO 100
C GESTOPPT DURCH RUECKSTREUUNG
370 NRUECK(IL) = NRUECK(IL) + 1
ICHAR = BACKSC
GOTO 405
C GESTOPPT IN DER SCHICHT
400 NSTOP(IL) = NSTOP(IL) + 1
ICHAR = STOPPD
405 GOTO ( 100, 100, 410, 100, 410, 410, 410 ), IOTYPE(IL)
C DRUCKE 'STOPPED'-DATEN AUS
410 IUNIT = NLUNIT(IL)
c-kt RELSTR wird anstelle von IRELST (=RELSTR*1000) ausgegeben
WRITE (IUNIT,2102) ICHAR,ID,IL,SCATTR,STRAGG,STRAG2,RELSTR,
1 NLAY,EBEGIN,BMASS,X,Y,Z, PX,PY,PZ,THETA,
2 PHI,TTOTAL, TLAST, ELOSS
2102 FORMAT(' ',A1,I10,2X,I3,2X,3L1,2X,E10.4,1X,I4,1X,1PE12.5 /
1 4X,6(1PE12.5) / 4X,6(1PE12.5) )
NCOUNT(IL) = NCOUNT(IL) + 1
IF ( NCOUNT(IL) .GE. NLIMIT(IL) ) GOTO 420
100 CONTINUE
C ----------------------------------------------------------------------
WRITE (SYSOUT,2103) NMAX
2103 FORMAT(//' MAXIMUM NUMBER OF PARTICLES (',I7,') REACHED.')
GOTO 500
C LIMIT IN EINER SCHICHT ERREICHT
420 WRITE (SYSOUT,2104) NLIMIT(IL), IL
2104 FORMAT(' MAXIMUM NUMBER OF OUTPUT-RECORDS (',I7,') IN LAYER#',I3,
1 ' REACHED.')
GOTO 500
C LIMIT AN START-VERSUCHEN ERREICHT
440 WRITE (SYSOUT,2105) NSVMAX
2105 FORMAT(' MAXIMUM NUMBER OF PASSING PARTICLES (',I7,') REACHED.')
500 NPARTI = ID
WRITE(SYSOUT,2205) SEEDS(1)
2205 FORMAT(' STOP VALUE OF RANDOM NUMBER = ',I16)
RETURN
END
C ----------------------------------------------------------------------
C
C NAME: BEAMIN
C
C FUNKTION: BERECHNET DEN IMPULS- UND DIE ORTSKOORDINATEN
C
C ----------------------------------------------------------------------
SUBROUTINE BEAMIN (NSV)
REAL R(3)
INTEGER*4 SEEDS
COMMON
1 / SEEDCO / SEEDS(6)
2 / BEAM / PMEAN, PFWHM, PLOW, PHIGH, XFWHM, YFWHM,
3 RMAX,
4 IBTYP, XBEAM, YBEAM, PSIGMA, RMAX2
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
5 , XFOCUS, YFOCUS
c-kt Beginn und Stärke des Magnetfeldes
6 , MSTART, MGAUSS
COMMON
1 / KINEMA / ID, X, Y, Z,
2 PX, PY, PZ, PTOTAL,
c-kt Einführen der Variable TTHETA für den Tangens des Streuwinkels:
3 THETA, PHI, STHETA, CTHETA, TTHETA,
4 SPHI, CPHI,
5 TTOTAL, TLAST, ELOSS
C ERSTELLE GAUSS-ZUFALLSZAHLEN
GOTO ( 100, 200, 300 ), IBTYP
C ISOTROPIC DISTRIBUTION
100 CALL GGUBS (SEEDS,2,R)
NSV = NSV + 1
X = ( R(1) - 0.500 ) * RMAX * 2.0
Y = ( R(2) - 0.500 ) * RMAX * 2.0
IF ( (X*X+Y*Y) .GT. RMAX2 ) GOTO 100
GOTO 500
C GAUSSIAN DISTRIBUTION
200 CALL GGNML (SEEDS,2,R)
NSV = NSV + 1
X = XBEAM * R(1)
Y = YBEAM * R(2)
IF ( (X*X+Y*Y) .GT. RMAX2 ) GO TO 200
GOTO 500
c
c isotropic rectangular distribution, TP, 1-jun-1995
c
300 call ggubs(seeds,2,r)
nsv = nsv + 1
x = ( R(1) - 0.500 ) * XBEAM
y = ( R(2) - 0.500 ) * YBEAM
c
500 Z = 0.0
C IMPULSVERTEILUNG
600 CALL GGNML (SEEDS,1,R(3))
IF ( PFWHM .GT. 0.0 ) THEN
PTOTAL = PMEAN + ( PSIGMA * R(3) )
ELSE
PTOTAL = PLOW + ( R(3) * (PHIGH-PLOW) )
ENDIF
IF ( PTOTAL .LT. PLOW .OR. PTOTAL .GT. PHIGH ) GOTO 600
RETURN
END
C ----------------------------------------------------------------------
C
C NAME: COMSUM
C
C FUNKTION: DRUCKT DIE STATISTIK AUS
C
C ----------------------------------------------------------------------
SUBROUTINE COMSUM
CHARACTER*10 BEAMTY(3)
CHARACTER*20 COMENT,COMT
REAL*8 RGDATA
INTEGER*4 SEEDS, TEST
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY und ALL_XY:
INTEGER*2 IMPOUT, THEOUT, STPOUT, XYOUT, XYALL, LUNINIT
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
INTEGER*2 SYSOUT, SYSIN
COMMON
1 / OUTPUT / NLIMIT(200), NCOUNT(200), IOTYPE(200),
2 NLUNIT(200), ELIMIT(200)
3 / COMPON / NZCOMP(20), ACOMP(20), NATOM(20), NELEM, NELEMC
4 / INFOCO / NPARTI, NTHROU, NSTOP(200), NSTRVO(200),
5 NSTRIN(200), NRUECK(200), NSTRGE(200), NMLRE(200)
COMMON
1 / RETAB / NRGMAX, NRGCUR, RGDATA (5,3000)
2 / LAYER / NLAY, ILAYS(200), ILAYE(200),
3 THICK(200), RTHICK(200), RO(200),
4 ZMEAN(200), AMEAN(200), RMEAN(200),
5 ZKOO(200), IGEOT(200), PGEOT(5,200)
6 / LAYERG / NGRP, IGRPS(200), IGRPE(200)
7 / LAYERC / COMENT(200)
8 / PARTIC / BMASS, BMASS2, EKMIN, EKMAX, IBFLAG, IBUNIT
9 / BEAM / PMEAN, PFWHM, PLOW, PHIGH, XFWHM, YFWHM,
1 RMAX,
2 IBTYP, XBEAM, YBEAM, PSIGMA, RMAX2
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
3 , XFOCUS, YFOCUS
c-kt Beginn und Stärke des Magnetfeldes
4 , MSTART, MGAUSS
3 / SEEDCO / SEEDS(6)
c-kt Integer-Arrays für die Unit-Nummern zu den Befehlen IMPULS, WINKEL,
c-kt STOPPED, WHERE_XY und ALL_XY:
COMMON
1 / SGLOUT / IMPOUT(201), THEOUT(201), STPOUT(200), XYOUT(201),
2 XYALL(201), LUNINIT
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
COMMON
1 / STATUS / SYSOUT, SYSIN, TEST
DATA BEAMTY / 'ISOTROPIC ','GAUSSIAN ', 'RECTANGULA' /
C ----------------------------------------------------------------------
WRITE (SYSOUT,2200) BMASS, EKMIN, EKMAX
2200 FORMAT('0','PARTICLE PARAMETER'/
1 ' ',5X,'MASS IN MEV/C**2',T45,1PG12.5/
2 ' ',5X,'MINIMUM ENERGY IN MEV',T45,1PG12.5/
3 ' ',5X,'MAXIMUM ENERGY IN MEV',T45,1PG12.5)
RPFWHM = PFWHM / PMEAN
WRITE(SYSOUT,2201) PMEAN,PFWHM,RPFWHM,PLOW,PHIGH
2201 FORMAT('0','BEAM MOMENTUM DISTRIBUTION PARAMETER'/
1 6X,'MEAN MOMENTUM IN MEV/C',T60,1PE12.5/
2 6X,'ABSOLUTE FWHM OF THE MOMENTUM DISTRIBUTION',T60,1PE12.5/
3 6X,'RELATIVE FWHM OF THE MOMENTUM DISTRIBUTION',T60,1PE12.5/
4 6X,'LOWEST MOMENTUM IN MEV/C',T60,1PE12.5/
5 6X,'HIGHEST MOMENTUM IN MEV/C',T60,1PE12.5)
WRITE(SYSOUT,2202) BEAMTY(IBTYP),XFWHM,YFWHM,RMAX
2202 FORMAT('0','BEAM SPATIAL DISTRIBUTION PARAMETER'/
1 6X,'BEAM-SPOT DISTRIBUTION',T47,A10/
2 6X,'X-FWHM OF THE BEAMSPOT IN CM',T45,1PE12.5/
3 6X,'Y-FWHM OF THE BEAMSPOT IN CM',T45,1PE12.5/
4 6X,'MAXIMUM RADIUS OF BEAMSPOT IN CM',T45,1PE12.5)
C ----------------------------------------------------------------------
MSTOP = 0
MRUECK = 0
MSTRIN = 0
MSTRVO = 0
MCOUNT = 0
MSTRGE = 0
MMLRE = 0
C ----------------------------------------------------------------------
DO 10 I=1,NLAY
MSTOP = MSTOP + NSTOP(I)
MRUECK = MRUECK + NRUECK(I)
MSTRIN = MSTRIN + NSTRIN(I)
MSTRVO = MSTRVO + NSTRVO(I)
MCOUNT = MCOUNT + NCOUNT(I)
MSTRGE = MSTRGE + NSTRGE(I)
MMLRE = MMLRE + NMLRE (I)
10 CONTINUE
C ----------------------------------------------------------------------
WRITE(SYSOUT,1900) NPARTI,NTHROU
1900 FORMAT('0',20X,'S U M M A R Y '//
1 '0','PARTICLES STARTET FROM Z=0',T45,I12/
2 '0','PARTICLES THROUGH ALL AREAS',T45,I12//
3 '0',20X,'LAYER STATISTIC'//
4 '0','GRP/LAY/LUN',1X,'COMMENT',14X,
5 ' STOPPED',1X,'BACK-SCAT',1X,
6 'SCAT_BEF.',1X,'SCAT_INS.',1X,
7 ' STRG_ERR',1X,' MLR_ERR',1X,
8 ' RECORDS'//)
C ----------------------------------------------------------------------
DO 30 LGRP=1,NGRP
KSTOP = 0
KRUECK = 0
KSTRIN = 0
KSTRVO = 0
KCOUNT = 0
KSTRGE = 0
KMLRE = 0
WRITE(SYSOUT,1904)
1904 FORMAT(' ',102('-'))
IFROM = IGRPS(LGRP)
ITO = IGRPE(LGRP)
COMT = COMENT(LGRP)
DO 20 I=IFROM,ITO
WRITE(SYSOUT,1902) LGRP,I,NLUNIT(I),COMT,NSTOP(I),NRUECK(I),
1 NSTRVO(I),NSTRIN(I),
2 NSTRGE(I),NMLRE(I),NCOUNT(I)
1902 FORMAT(' ',I3,'/',I3,'/',I3,1X,A20,7(1X,I9))
c-kt Ausgabe der Stop-Daten des Targets: aufeinanderfolgend in jeder
c-kt Zeile die Anzahl der in der Lage gestoppten Teilchen
IF ( STPOUT(LGRP) .NE. 0 ) THEN
WRITE (STPOUT(LGRP),1909) NSTOP(I)
1909 FORMAT(' ',I6)
ENDIF
c-kt Ende
KSTOP = KSTOP + NSTOP(I)
KRUECK = KRUECK + NRUECK(I)
KSTRIN = KSTRIN + NSTRIN(I)
KSTRVO = KSTRVO + NSTRVO(I)
KCOUNT = KCOUNT + NCOUNT(I)
KSTRGE = KSTRGE + NSTRGE(I)
KMLRE = KMLRE + NMLRE (I)
20 CONTINUE
IF ( IFROM .GE. ITO ) GOTO 30
WRITE(SYSOUT,1904)
WRITE(SYSOUT,1905) LGRP, NLUNIT(I), KSTOP, KRUECK, KSTRVO,
1 KSTRIN, KSTRGE, KMLRE, KCOUNT
1905 FORMAT(' ',I3,'/SUM/',I3,21X,7(1X,I9))
30 CONTINUE
C-----------------------------------------------------------------------
WRITE(SYSOUT,1904)
WRITE(SYSOUT,1903) MSTOP, MRUECK, MSTRVO, MSTRIN, MSTRGE, MMLRE,
1 MCOUNT
1903 FORMAT(' ','TOTAL SUM ',21X,7(1X,I9)//)
RETURN
END
C ----------------------------------------------------------------------
C
C NAME: SUMOUT
C
C FUNKTION: GIBT DIE STATISTIK AUF EINEM HILFSFILE AUS
C (ZUR WEITERVERARBEITUNG MIT SAS)
C
C ----------------------------------------------------------------------
SUBROUTINE SUMOUT(ILUNIT)
CHARACTER*20 COMENT
REAL*8 RGDATA
COMMON
1 / OUTPUT / NLIMIT(200), NCOUNT(200), IOTYPE(200),
2 NLUNIT(200), ELIMIT(200)
3 / COMPON / NZCOMP(20), ACOMP(20), NATOM(20), NELEM, NELEMC
4 / INFOCO / NPARTI, NTHROU, NSTOP(200), NSTRVO(200),
5 NSTRIN(200), NRUECK(200), NSTRGE(200), NMLRE(200)
COMMON
1 / RETAB / NRGMAX, NRGCUR, RGDATA (5,3000)
2 / LAYER / NLAY, ILAYS(200), ILAYE(200),
3 THICK(200), RTHICK(200), RO(200),
4 ZMEAN(200), AMEAN(200), RMEAN(200),
5 ZKOO(200), IGEOT(200), PGEOT(5,200)
6 / LAYERG / NGRP, IGRPS(200), IGRPE(200)
7 / LAYERC / COMENT(200)
8 / PARTIC / BMASS, BMASS2, EKMIN, EKMAX, IBFLAG, IBUNIT
9 / BEAM / PMEAN, PFWHM, PLOW, PHIGH, XFWHM, YFWHM,
1 RMAX,
2 IBTYP, XBEAM, YBEAM, PSIGMA, RMAX2
c-kt Z-Komponente von X- und Y-Strahlfoci für den Befehl FOCUS
3 , XFOCUS, YFOCUS
c-kt Beginn und Stärke des Magnetfeldes
4 , MSTART, MGAUSS
C ----------------------------------------------------------------------
MSTOP = 0
MRUECK = 0
MSTRIN = 0
MSTRVO = 0
MCOUNT = 0
MSTRGE = 0
MMLRE = 0
WRITE(ILUNIT,1900) NPARTI, NTHROU
1900 FORMAT(' ','000/000/000',2(1X,I9))
C ----------------------------------------------------------------------
DO 10 I=1,NLAY
MSTOP = MSTOP + NSTOP(I)
MRUECK = MRUECK + NRUECK(I)
MSTRIN = MSTRIN + NSTRIN(I)
MSTRVO = MSTRVO + NSTRVO(I)
MCOUNT = MCOUNT + NCOUNT(I)
MSTRGE = MSTRGE + NSTRGE(I)
MMLRE = MMLRE + NMLRE (I)
10 CONTINUE
DO 30 LGRP=1,NGRP
IFROM = IGRPS(LGRP)
ITO = IGRPE(LGRP)
DO 20 I=IFROM,ITO
WRITE(ILUNIT,1902) LGRP,I,NLUNIT(I),NSTOP(I),NRUECK(I),
1 NSTRVO(I),NSTRIN(I),
2 NSTRGE(I),NMLRE(I),NCOUNT(I)
1902 FORMAT(' ',I3,'/',I3,'/',I3,7(1X,I9))
20 CONTINUE
30 CONTINUE
C ----------------------------------------------------------------------
WRITE(ILUNIT,1903) MSTOP, MRUECK, MSTRVO, MSTRIN, MSTRGE,
1 MMLRE, MCOUNT
1903 FORMAT(' ','999/000/000',7(1X,I9))
RETURN
END
C ----------------------------------------------------------------------
C
C NAME: TBOUND
C
C FUNKTION: PRUEFT ZAHLEN AUF GUELTIGKEIT
C
C ----------------------------------------------------------------------
SUBROUTINE TBOUND (TYPE, TEXT, VAR, LEFT, RIGHT, INOUT)
CHARACTER*2 TYPE,INOUT
CHARACTER*10 TEXT
REAL*4 VAR, LEFT, RIGHT
INTEGER*4 IVAR, ILEFT, IRIGHT, TEST
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
INTEGER*2 SYSOUT, SYSIN
REAL*4 RVAR, RLEFT, RRIGHT
EQUIVALENCE ( RVAR,IVAR ), ( RLEFT,ILEFT ), ( RRIGHT,IRIGHT )
c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT
COMMON
1 / STATUS / SYSOUT, SYSIN, TEST
C ----------------------------------------------------------------------
RVAR = VAR
RLEFT = LEFT
RRIGHT = RIGHT
C ----------------------------------------------------------------------
IF ( TYPE .EQ. 'I4' ) THEN
IF ( INOUT .EQ. 'IN' ) THEN
IF ( IVAR .LT. ILEFT .OR. IVAR .GT. IRIGHT ) THEN
WRITE(SYSOUT,1800) TEXT,IVAR,ILEFT,IRIGHT
1800 FORMAT(' ',A10,'(',I12,') NOT IN (',I10,'/',I10,')')
STOP 8
ENDIF
ELSE
IF ( IVAR .GE. ILEFT .AND. IVAR .LE. IRIGHT ) THEN
WRITE(SYSOUT,1801) TEXT,IVAR,ILEFT,IRIGHT
1801 FORMAT(' ',A10,'(',I12,') IN (',I10,'/',I10,')')
STOP 8
ENDIF
ENDIF
ELSE
IF ( INOUT .EQ. 'IN' ) THEN
IF ( RVAR .LT. RLEFT .OR. RVAR .GT. RRIGHT ) THEN
WRITE(SYSOUT,1802) TEXT,RVAR,RLEFT,RRIGHT
1802 FORMAT(' ',A10,'(',1PG12.5,') NOT IN (',1PG12.5,
1 '/',1PG12.5,')')
STOP 8
ENDIF
ELSE
IF ( RVAR .GE. RLEFT .AND. RVAR .LE. RRIGHT ) THEN
WRITE(SYSOUT,1803) TEXT,RVAR,RLEFT,RRIGHT
1803 FORMAT(' ',A10,'(',1PG12.5,') IN (',1PG12.5,
1 '/',1PG12.5,')')
STOP 8
ENDIF
ENDIF
ENDIF
RETURN
END
Reference in New Issue View Git Blame Copy Permalink