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 Eintraege in der '//
     1        'Energie-Reichweite-Tabelle.'
     2      /' Wahrscheinlich muss entweder fuer diese Schicht mit dem '//
     3      ,'Befehl',
     4      /'''BUILD_TAB'' ein Parameter kleiner 1.E-4 uebergeben '//
     5       'oder der '
     6      /' 3. Parameter von ''PARTICLE'' vergroessert 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 muss vergroessert '//
     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) I, NSTOP(I)
 1909       FORMAT(' ',I3, ' ',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