CCY   PROGRAM GAMMA (TAPE 1, TAPE 2, TAPE 3)
CAT   PROGRAM GAMMA
      DIMENSION EGRENZ(10), NITER(9)
      DIMENSION M(9), H(9), NSTR(9)
      DIMENSION SIJ(9,9), SA(9,9), D(9,9), RHO(9), XNE(9) 
      DIMENSION FLUSS(9,500), QUELLE(9) 
      DIMENSION NABS(100), H1(100) ,H2(100)
      DIMENSION STROM(9,100), STRAHL(9,100), DOSIS(100) 
C  DEFINITION DER I/O DATEIEN
CAT   OPEN (1, FILE="BEREIN.DAT")
CAT   OPEN (2, FILE="BERMAT.DAT")
CAT   OPEN (3, FILE="BERAUS.DAT", POSITION="APPEND")
C  AUSDRUCK DER AUFGABENSTELLUNG
      CALL COPY
C  EINLESEN DER MATERIALDATEN 
      READ (1,100) X
       NENERG=INT(X+.5) 
      READ (1,100) (EGRENZ(N),N=1,NENERG+1) 
      READ (1,100) X
       NZONEN=INT(X+.5) 
      DO 12 NZ=1,NZONEN 
       READ (1,100) X1,H(NZ),X2 
       M(NZ)=INT(X1+.5) 
       NSTR(NZ)=INT(X2+.5)
 12   CONTINUE
      READ (1,100) X1, X2, X3 
       NVAR=INT(X1+.5)
       NEND=INT(X2+.5)
       NSP =INT(X3+SIGN(.5,X3)) 
      READ (1,100) XKAPPA 
      READ (1,100) X
       NI =INT(X+.5)
      READ (1,100) R1 
      READ (1,100) (QUELLE(NE),NE=1,NENERG) 
      READ (1,100) PARA

      CALL MAKOMP (NENERG, EGRENZ, NMATE, SIJ, SA, D, XNE, RHO, PARA)
      IF (NVAR.LT.1.OR.NVAR.GE.NZONEN) THEN 
       DO 31 NE=1,NENERG 
 31     NITER(NE)=NI
       CALL SCHIRM (NENERG, NZONEN, XKAPPA, NITER, SIJ, SA, D, XNE, 
     1              NSTREI, NMATE, M, H, NSTR, R1, QUELLE, FLUSS) 
        WRITE (3,200) NZONEN, NSTREI
        WRITE (3,300) (NSTR(NZ),NZ=1,NZONEN)
        WRITE (3,400) NENERG
CCY     WRITE (3,500) (EGRENZ(NE),NE=1       ,NENERG+1)
CAT     WRITE (3,500) (EGRENZ(NE),NE=NENERG-2,NENERG+1)
         NULL= 0
CCY     WRITE (3,300) NULL, (NITER(NE),NE=1       ,NENERG) 
CAT     WRITE (3,300) NULL, (NITER(NE),NE=NENERG-2,NENERG) 
         DO 32 NS=1,NSTREI
CCY32   WRITE (3,600) NS,(FLUSS(NE,NS),NE=1       ,NENERG) 
CAT32   WRITE (3,600) NS,(FLUSS(NE,NS),NE=NENERG-2,NENERG) 
      ELSE
        NANF=NSTR(NVAR) 
        NBER= 0 
        WRITE (3,400) NENERG
CCY     WRITE (3,500) (EGRENZ(NE),NE=1       ,NENERG+1)
CAT     WRITE (3,500) (EGRENZ(NE),NE=NENERG-2,NENERG+1)
       DO 4 NLOOP=NANF,NEND,NSP 
        NBER = NBER+1 
        NSTR(NVAR)= NLOOP 
         DO 41 NE=1,NENERG 
 41     NITER(NE)=NI
       CALL SCHIRM (NENERG, NZONEN, XKAPPA, NITER, SIJ, SA, D, XNE, 
     1              NSTREI, NMATE, M, H, NSTR, R1, QUELLE, FLUSS) 
         NABS(NBER)=     NSTR(NVAR)
         H1(NBER)= FLOAT(NSTR(NVAR  )) *H(NVAR) 
         HGEWEB =  FLOAT(NSTR(NZONEN)) *H(NZONEN)
          H2(NBER) = R1
          DO 42 NZ=1,NZONEN-1
 42      H2(NBER)= H2(NBER) + FLOAT(NSTR(NZ)) * H(NZ)
          DOSIS(NBER)= 0. 
          NWAND= NSTREI -NSTR(NZONEN)
          DO 43 NE=1,NENERG 
         EM= (EGRENZ(NE+1) +EGRENZ(NE)) /2.
         STROM(NE,NBER) =FLUSS(NE,NWAND)
          PSIEIN= FLUSS(NE,NWAND)-FLUSS(NE,NWAND+1)
          PSIAUS= FLUSS(NE,NSTREI-1)-FLUSS(NE,NSTREI)
          PSIEIN= PSIEIN *NSTR(NZONEN) /(NSTR(NZONEN)-2)
         STRAHL(NE,NBER) =(PSIEIN-PSIAUS)/H(NZONEN)*D(M(NZONEN),NE)
     1                   /HGEWEB /RHO(M(NZONEN)) *1.602E-13 *3600.  
 43      DOSIS(NBER) = DOSIS(NBER) + STRAHL(NE,NBER)
 4     CONTINUE 
       WRITE (3,700)
        DO 5 N=1,NBER
CCY 5  WRITE (3,600) NABS(N),(STROM(NE,N),NE=1       ,NENERG)
CAT 5  WRITE (3,600) NABS(N),(STROM(NE,N),NE=NENERG-2,NENERG)
       WRITE (3,100)
       WRITE (3,800)
CCY    WRITE (3,900) (NE,NE=1       ,NENERG)
CAT    WRITE (3,900) (NE,NE=NENERG-2,NENERG)
        DO 6 N=1,NBER
CCY 6  WRITE (3,500) H1(N),H2(N),DOSIS(N),(STRAHL(NE,N),NE=1,NENERG)
CAT 6  WRITE (3,500) H1(N),H2(N),DOSIS(N),
CAT  1               (STRAHL(NE,N),NE=NENERG-2,NENERG)
       ENDIF 

100   FORMAT (1X,F12.3) 
200   FORMAT (1X,10HZONENZAHL:,I5,15H; STREIFENZAHL:,I5,7H; NSTR:)
300   FORMAT (12(1X,I10)) 
400   FORMAT (1X,19HZAHL DER ENERGIEGR:,I3,16H; GRENZENERGIEN:) 
500   FORMAT (12(1X,1PE10.3)) 
600   FORMAT (1X,I10,11(1X,1PE10.3))
700   FORMAT (1X,10HSTR(ABSCH),3X,30H        GRUPPENFLUESSE        ,
     1        2X,20H [M**-2*S**-1]      )
800   FORMAT (1X,10H   SCHIRM ,1X,10HENTFERNUNG,1X,10HGESAMTDOS.,
     1        3X,30H               GRUPPENDOSEN   )
900   FORMAT (2(1X,10H    [M]   ),1X,10H  [SV/H]  ,9(1X,I10))
      STOP
      END 

C  UNTERPROGRAMM ZUR ERRECHNUNG DER MATERIALDATEN 
      SUBROUTINE MAKOMP (NENERG,EGRENZ,NMATE,SIJ,SA,D,XNE,RHO,PARA)
      DIMENSION EGRENZ(10),ZORD(19),AMOL(19),DICHTE(9,19) 
      DIMENSION SIJ(9,9), SA(9,9), D(9,9), RHO(9), XNE(9) 
      DOUBLE PRECISION S, SS, DI
      DOUBLE PRECISION O, U, A, B, DE, XM, XG, XI 
      DOUBLE PRECISION EI, E, E1, E2, E3, E4, ZLN 
      DO 1 NK=1,19
 1     READ (2,100) ZORD(NK), AMOL(NK)
      READ (2,100) X
       NMATE= INT(X+.5) 
      DO 2 NM=1,NMATE 
 2     READ (2,100) (DICHTE(NM,NK),NK=1,19) 
      IF (PARA.GE.1.) THEN
        DO 21 NE=1,NENERG 
 21      READ (2,200) (SIJ(NE,J),J =1,NENERG)
        DO 22 NM=1,NMATE
         READ (2,200) (SA(NM,NE),NE=1,NENERG) 
         READ (2,200) (D (NM,NE),NE=1,NENERG) 
         READ (2,200) RHO (NM)
 22      READ (2,200) XNE (NM)
      ELSE
        EL  = .5110034
        R   = .2494668E-28
        AVOG= 6.022045E26 
      DO 3 NE=1,NENERG
      DO 4 J=1,NENERG 
        O= EGRENZ(J+1)
        U= EGRENZ(J)
        A= EGRENZ(NE) 
        B= EGRENZ(NE+1) 
        DE= B-A 
        S= 0. 
        XM= 10. 
        XG= 10. 
       IF (2.*O.LT.EL) XM= EL*O/ (EL-2.*O)
       IF (2.*U.LT.EL) XG= EL*U/ (EL-2.*U)
       IF (A.GE.XM .OR. A.LT.(O-DE/1000.)) GOTO 12
       IF (B.GT.XM) B=XM
       IF (XG.LT.A) GOTO 11 
       IF (XG.GT.B) XG=B
        S=S +DLOG(O/U) *DLOG(XG/A)
        S=S +(-2.*EL*DLOG(O/U)+EL*EL*(1./U-1./O)) *(1./A-1./XG) 
        S=S +((O*O-U*U)/2. +2.*EL*(O-U) -2.*EL*EL*DLOG(O/U))
     1       /2. *(1./A/A -1./XG/XG)
        S=S +EL*EL*(O-U)/3. *(1./A**3 -1./XG**3)
 11    IF (XG.GE.B) GOTO 12 
       IF (XG.LT.A) XG=A
        S=S +(1./XG**3 -1./B**3) *EL*EL*O/3.
        S=S +(1./XG/XG -1./B/B) *(EL*EL/2. +EL*O +O*O/4.) 
        S=S +(1./XG -1./B) *(2.*EL -EL*EL/O)
        S=S +EL/2. *(1./(EL+2.*XG) -1./(EL+2.*B)) 
        S=S +(DLOG( B/EL)+2.*EL/ B+EL*EL/ B/ B) *DLOG(O/ B+2.*O/EL) 
        S=S -(DLOG(XG/EL)+2.*EL/XG+EL*EL/XG/XG) *DLOG(O/XG+2.*O/EL) 
        S=S +(DLOG((EL+2.*B)/(EL+2.*XG)) -DLOG(B/XG)) /2. 
        S=S +ALOG(2.) *DLOG((2.*B+EL)/(2.*XG+EL)) 
        S=S +((DLOG(B/EL)**2) -(DLOG(XG/EL)**2)) /2.
        S=S -((DLOG((2.*B+EL)/EL))**2 -(DLOG((2.*XG+EL)/EL))**2) /2.
       DO 13 I=1,25 
        XI= FLOAT(I)
 13     S=S +1./XI/XI *((EL/(2*XG+EL))**I -(EL/(2.*B+EL))**I) 
 12    SIJ(NE,J)= S*R*EL/DE 
 4    CONTINUE
 3    CONTINUE
      DO 7 NM=1,NMATE 
       RHO (NM) =0.
       XNE (NM) =0. 
      DO 5 NK=1,19
       RHO(NM)=RHO(NM) +DICHTE(NM,NK)
 5     XNE(NM)=XNE(NM) +DICHTE(NM,NK) /AMOL(NK) *AVOG *ZORD(NK) 
       IMAX= 50
       FIM= FLOAT(IMAX) 
       DICHTE(NM,1)= DICHTE(NM,1)/2.
      DO 7 NE=1,NENERG
       SA(NM,NE)=0. 
       D (NM,NE)=0. 
       DE= (EGRENZ(NE+1)-EGRENZ(NE)) /FIM 
       EI = EGRENZ(NE) -DE/2. 
      DO 7 I=1,IMAX 
       EI= EI +DE 
       E= EI
       IF (EI.LT..1) E= 10.*EI
       SS= 0. 
      DO 6 NK=1,19
       Z= ZORD(NK)
       ZLN= DLOG(DBLE(Z)) 
       S=   (1.6268E-9 -2.683E-12*Z) /(1. +4.173E-2*Z) /E 
       S=S +(1.5274E-9 -0.511E-12*Z) /(1. +1.027E-2*Z) /E/E 
       S=S +(1.1330E-9 -2.177E-12*Z) /(1. +2.013E-2*Z) /E**3.5
       S=S -(0.0912E-9             )                   /E**4
       S=S *Z**5 *(1. +.01481*ZLN*ZLN -.788E-3*ZLN*ZLN*ZLN) 
       S=S *1.E-28 *DICHTE(NM,NK) *AVOG /AMOL(NK) 
       SS=SS +S 
 6    CONTINUE
       IF (EI.LT..1) SS= SS*1000. 
       E1= E/EL 
       E2= E1*E1
       E3= E2*E1
       E4= E3*E1
       DI=   (2.*E2-2.*E1-6.)/E3 +2./(2.*E1+1.)**2
       DI=DI +(-E2+4.*E1+3.) /E4 *DLOG(2.*E1+1.)
       DI=DI *R *XNE(NM)
       DI= 1./(SS+DI)/3.
      SA(NM,NE)= SA(NM,NE) +SS/FIM
      D (NM,NE)= D (NM,NE) +DI/FIM
 7    CONTINUE
      DO 31 NE=1,NENERG 
 31    WRITE (2,200) (SIJ(NE,J),J=1,NENERG) 
      DO 32 NM=1,NMATE
       WRITE (2,200) (SA(NM,NE),NE=1,NENERG)
       WRITE (2,200) (D (NM,NE),NE=1,NENERG)
       WRITE (2,200) RHO(NM)
       WRITE (2,200) XNE(NM)
 32   CONTINUE
      ENDIF
      RETURN
100   FORMAT (1X,F12.3) 
200   FORMAT (1X,E20.14)
      END 

      SUBROUTINE SCHIRM (NENERG,NZONEN,XKAPPA,NITER,SIJ,SA,D,XNE, 
     1                   NSTREI, NMATE, M, H, NSTR, R1, QUELLE, FLUSS)
      DIMENSION SIJ(9,9), SA(9,9), ST(9), D(9,9), XNE(9)
      DIMENSION M(9), H(9), NSTR(9)
      DIMENSION NITER(9), QUELLE(9) 
      DIMENSION S(9,500), FLUSS(9,500)
      DIMENSION A(3,500), Q(500)
      REAL      B(3,500),FI(500)
       NSTREI= 1
       DO 11 NZ=1,NZONEN
 11     NSTREI= NSTREI +NSTR(NZ)
       DO 12 NE=1,7 
       IF (R1.LT.H(1)) THEN 
         S(NE,1)= 0.5 *QUELLE(NE) /XNE(M(1))
        DO 13 NS=2,NSTR(1)
 13      S(NE,NS)= QUELLE(NE) /XNE(M(1))
        NS = NSTR(1) + 1
         S(NE,NS)= QUELLE(NE) *H(1) /(H(1)*XNE(M(1))+H(2)*XNE(M(2)))
        DO 14 NS=NSTR(1)+2,NSTREI 
 14      S(NE,NS)=0.
       ELSE 
         S(NE,1)= QUELLE(NE) *(R1/H(1) -XKAPPA/2.) /XNE(M(1)) 
        DO 15 NS=2,NSTREI 
 15      S(NE,NS)=0.
       ENDIF
 12    CONTINUE 
C  GRUPPENWEISE BERECHNUNG DES FLUSSVERLAUFES 
       DO 1 NE=NENERG,1,-1 
       DO 21 NM=1,NMATE 
       ST(NM)=SA(NM,NE) 
       DO 21 J=1,NENERG 
       ST(NM)=ST(NM)+SIJ(NE,J)*XNE(NM)
 21   CONTINUE
C  MATRIX, QUELLVEKTOR
      NS=1
       R=R1 
      IF (R1.LT.H(1)) THEN
       R= 0.
       A(1,NS)= 0.
       A(2,NS)= D(M(1),NE)/H(1)/H(1) *(2.+2.*XKAPPA) +.5*ST(M(1)) 
       A(3,NS)=-D(M(1),NE)/H(1)/H(1) *(2.+2.*XKAPPA)
      ELSE
       A(1,NS)=0. 
       A(2,NS)=D(M(1),NE)/H(1)* (1./H(1)+XKAPPA/2./(R+H(1)/2.)) 
     1        + ST(M(1))
       A(3,NS)=-D(M(1),NE)/H(1)* (1./H(1)+XKAPPA/2./(R+H(1)/2.))
      ENDIF 
       Q(NS)=S(NE,NS)*XNE(M(1)) 
      DO 2 NZ=1,NZONEN
      DO 3 NST=2,NSTR(NZ) 
      NS=NS+1
       R=R +H(NZ) 
       A(1,NS)=-D(M(NZ),NE)/H(NZ)*(1./H(NZ)-XKAPPA/2./(R-H(NZ)/2.)) 
       A(2,NS)=D(M(NZ),NE)/H(NZ)
     1        *(2./H(NZ) +XKAPPA/2.*(1./(R+H(NZ)/2.)-1./(R-H(NZ)/2.)))
     2        +ST(M(NZ))
       A(3,NS)=-D(M(NZ),NE)/H(NZ)*(1./H(NZ)+XKAPPA/2./(R+H(NZ)/2.)) 
       Q(NS)=S(NE,NS)*XNE(M(NZ))
 3    CONTINUE
      NS=NS+1 
       R= R+H(NZ) 
      IF (NS.EQ.NSTREI) THEN
       A(1,NS)=-D(M(NZ),NE)/H(NZ)*(1./H(NZ)-XKAPPA/2./(R-H(NZ)/2.)) 
       A(2,NS)=D(M(NZ),NE)/H(NZ) *(1./H(NZ)-XKAPPA/2./(R-H(NZ)/2.)) 
     1        +(1./H(NZ) +XKAPPA/2./R) +ST(M(NZ)) 
       A(3,NS)=0. 
       Q(NS)=S(NE,NS) *XNE(M(NZ)) 
      ELSE
       A(1,NS)=-D(M(NZ),NE)/H(NZ)*(2./(H(NZ)+H(NZ+1)) 
     1        -XKAPPA/2./(R-H(NZ)/2.))
       A(2,NS)=D(M(NZ+1),NE)/H(NZ+1) 
     1        *(2./(H(NZ)+H(NZ+1)) +XKAPPA/2./(R+H(NZ+1)/2.)) 
     2        +D(M(NZ),NE)/H(NZ)
     3        *(2./(H(NZ)+H(NZ+1)) -XKAPPA/2./(R-H(NZ)/2.)) 
     4        +(ST(M(NZ+1))*H(NZ+1)+ST(M(NZ))*H(NZ))
     5        /(H(NZ+1)+H(NZ))
       A(3,NS)=-D(M(NZ+1),NE)/H(NZ+1) 
     1        *(2./(H(NZ)+H(NZ+1)) +XKAPPA/2./(R+H(NZ+1)/2.)) 
       Q(NS)=S(NE,NS) /(H(NZ)+H(NZ+1))
     1      *(XNE(M(NZ))*H(NZ) +XNE(M(NZ+1))*H(NZ+1)) 
      ENDIF 
 2    CONTINUE

C  KONTROLLE DER MATRIX 
CCC   CALL MATRIX (NE,NZONEN,NSTR,M,A)

C  LOESUNG MIT GAUSSELIMINATION 
       DO 22 NS=1,NSTREI 
      FI(NS)=Q(NS) 
       DO 22 I=1,3
      B(I,NS)=A(I,NS) 
 22    CONTINUE 
      CALL GAUSS (NSTREI, B, FI)

C  ITERATIONSAUFRUF 
      NIT=NITER(NE) 
      CALL ITER (NSTREI, A, B, Q, FI, NIT)

C  FLUSS ABLEGEN UND QUELLEN ERMITTELN
      DO 4 NS=1,NSTREI
       FLUSS(NE,NS)=FI(NS)
      DO 4 J=1,NE
       S(J,NS)=S(J,NS)+FI(NS)*SIJ(NE,J) 
 4    CONTINUE
       NITER(NE)=NIT
 1    CONTINUE
      RETURN
      END 

C  UNTERPROGRAMM ZUR GAUSSELIMINATION 
      SUBROUTINE GAUSS (NSTREI, B ,FI)
      REAL B(3,500),FI(500) 
        DO 10 NS=1,NSTREI-1 
       B(1,NS+1)=B(1,NS+1)/B(2,NS)
       B(2,NS+1)=B(2,NS+1) - B(1,NS+1)*B(3,NS)
       FI(NS+1)=FI(NS+1) - B(1,NS+1)*FI(NS) 
 10     CONTINUE
       FI(NSTREI)=FI(NSTREI)/B(2,NSTREI)
        DO 11 NS=NSTREI-1,1,-1
       FI(NS)= (FI(NS) - B(3,NS)*FI(NS+1)) / B(2,NS)
 11     CONTINUE
      RETURN
      END 


C  UNTERPROGRAMM ZUR NACHITERATION
      SUBROUTINE ITER (NSTREI, A, B, Q, FI, NIT)
      DIMENSION A(3,500),Q(500) 
      REAL B(3,500),FI(500),RES(500)
      DO 1 NITER=1,NIT
       NS=1 
        RES(NS)=Q(NS) -A(2,NS)*FI(NS) -A(3,NS)*FI(NS+1) 
       DO 11 NS=2,NSTREI-1
        RES(NS)=Q(NS)-A(1,NS)*FI(NS-1)-A(2,NS)*FI(NS)-A(3,NS)*FI(NS+1)
 11    CONTINUE 
       NS=NSTREI
        RES(NS)= Q(NS) -A(1,NS)*FI(NS-1) -A(2,NS)*FI(NS)
       DO 12 NS=2,NSTREI
        RES(NS)=RES(NS) -B(1,NS)*RES(NS-1)
 12    CONTINUE 
        RES(NSTREI)= RES(NSTREI)/B(2,NSTREI)
       DO 13 NS=NSTREI-1,1,-1 
        RES(NS)=RES(NS) - B(3,NS)*RES(NS+1) 
        RES(NS)= RES(NS)/B(2,NS)
 13    CONTINUE 
       DO 14 NS=1,NSTREI
        FI(NS)=FI(NS)+RES(NS) 
 14    CONTINUE 
 1    CONTINUE
        DELTA=0.
       DO 2 NS=1,NSTREI
        DELTA = DELTA + RES(NS)/FI(NS)*RES(NS)/FI(NS)
 2     CONTINUE
       NIT= INT (-.5+10.*ALOG10 (SQRT (DELTA/NSTREI)))
      RETURN
      END 


C  UNTERPROGRAMM ZUR KONTROLLDOKUMENTATION
      SUBROUTINE MATRIX (NE,NZONEN,NSTR,M,A)
      DIMENSION NSTR(9), M(9)
      DIMENSION A(3,500)
        NULL= 0 
        NS=1 
       WRITE (3,900) NE,NULL,NULL,NS,(A(I,NS),I=1,3) 
       DO 11 NZ=1,NZONEN
       DO 12 NST=2,NSTR(NZ) 
        NS=NS+1 
 12    WRITE (3,900) NE,NZ,M(NZ),NS,(A(I,NS),I=1,3)
        NS=NS+1 
 11    WRITE (3,900) NE,NULL,NULL,NS,(A(I,NS),I=1,3) 
900   FORMAT(4(1X,I5),3(1X,E16.10)) 
      RETURN
      END

C  UNTERPROGRAMM ZUM AUSDRUCK DER AUFGABENSTELLUNG
      SUBROUTINE COPY
      CHARACTER*65 TITEL
 1    READ (1,100,END=2) TITEL
      WRITE (3,200) TITEL
      GOTO 1
 2    CONTINUE
      REWIND 1
      RETURN
100   FORMAT (A65)
200   FORMAT (1X,A65)
      END