c********************************************************************
c
c                     i-source.for
c
C##############################################################################
C
      SUBROUTINE RAND(XFE,HDM,EXTRA,IQAL,INST)
      COMMON/SEN1/ FINT(32),FI1B(32),FIS1(7,7,32)
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      COMMON/OUTPUT/OUT(17,8)
      LOGICAL EXTRA,HDM
C
C       RAND LIEFERT MIT XFE EIN MASS DAFUER , WIE STARK DAS SPEKTRUM
C       BEEINTRAECHTIGT WIRD DURCH RANDLAGE DES MAXIMUMS ODER RANDKANAELE
C       MIT GROSSER ZAEHLRATE
C
      IQAL = 0
      EXTRA = .FALSE.
      XMAX = 0.
      XFE1  = 0.
      XFE2  = 0.
      LLC1 = 1
      LLC2 = LC1
      IF ( .NOT. HDM .OR. INST .EQ. 3) GO TO 11
      LLC1 = 9
      LLC2 = LC1-8
   11 CONTINUE
      DO 1 I=LLC1,LLC2
      DO 1 J=1,LC2
      DO 1 K=1,LC3
      IF ( FIS1(K,J,I) .LT. XMAX) GO TO 1
      XMAX = FIS1(K,J,I)
      JMAX = J
      KMAX = K
      IMAX = I
    1 CONTINUE
      IF ( JMAX .EQ. 1 .OR. JMAX .EQ. LC2 ) IQAL = 1
      IF ( KMAX .EQ. 1 .OR. KMAX .EQ. LC3 ) IQAL = 1
      IF ( IMAX .EQ. 1 .OR. IMAX .EQ. LC1 ) IQAL = 1
      IF ( IQAL .EQ. 1 ) GO TO 10
      IF ( INST .EQ. 3 ) GO TO 10
      X1MAX = XMAX*0.01
      X2MAX = XMAX*0.1
      DO 2 J=1,LC2
      DO 2 I=LLC1,LLC2
      IF ( FIS1(1,J,I) .GT. X1MAX) XFE1 = XFE1 + 1.
      IF ( FIS1(LC3,J,I) .GT. X1MAX )  XFE1 = XFE1 + 1.
      IF ( FIS1(1,J,I) .GT. X2MAX ) XFE2 = XFE2 + 1.
      IF ( FIS1(LC3,J,I) .GT. X2MAX ) XFE2 = XFE2 + 1.
    2 CONTINUE
      DO 3 K=1,LC3
      DO 3 I=LLC1,LLC2
      IF ( FIS1(K,1,I)  .GT. X1MAX) XFE1 = XFE1 + 1.
      IF ( FIS1(K,LC2,I) .GT. X1MAX ) XFE1 = XFE1 + 1.
      IF ( FIS1(K,1,I) .GT. X2MAX ) XFE2 = XFE2 + 1.
      IF ( FIS1(K,LC2,I) .GT. X2MAX) XFE2 = XFE2 + 1.
    3 CONTINUE
      IF ( HDM ) GO TO 5
      DO 4 K=1,LC3
      DO 4 J=1,LC2
      IF ( FIS1(K,J,1) .GT. X1MAX) XFE1 = XFE1 + 1.
      IF ( FIS1(K,J,LC1) .GT. X1MAX) XFE1 = XFE1 + 1.
      IF ( FIS1(K,J,1) .GT. X2MAX ) XFE2 = XFE2 + 1.
      IF ( FIS1(K,J,LC1) .GT. X2MAX ) XFE2 = XFE2 + 1.
    4 CONTINUE
    5 CONTINUE
C     WRITE(6,55) XFE1,XFE2
   55 FORMAT(7H  XFE1=E12.4,6H XFE2=E12.4)
C     IM NDM GIBT ES MAXIMAL 230 RANDPUNKTE
C     IM HDM GIBT ES BEI OBIGER METHODE MAXIMAL 448 PUNKTE
C     WENN MEHR ALS 20 ( 5 ) PUNKTE MIT > 10% VOM MAXIMUM->EXTRA = .TRUE
      IF ( HDM ) GO TO 6
      IF ( XFE1 .GT. 50. ) EXTRA = .TRUE.
      IF ( XFE2 .GT. 20. ) EXTRA = .TRUE.
      GO TO 7
    6 IF ( XFE1 .GT. 30. ) EXTRA = .TRUE.
      IF ( XFE2 .GT. 5.  ) EXTRA = .TRUE.
    7 XFE = 1000.*XFE2 + XFE1
C     IF ( JMAX .EQ. 1 .OR. JMAX .EQ. LC2 ) XFE = XFE + 1.E6
C     IF ( KMAX .EQ. 1 .OR. KMAX .EQ. LC3 ) XFE = XFE + 1.E6
      OUT(10,2) = XFE
      IF ( EXTRA ) OUT(11,2) = 1.
   10 RETURN
      END
C
C##############################################################################
C
      SUBROUTINE EXTRAP(HDM,EXTRA,C,MENDIF,MAZDIF,MELDIF)
      COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
      COMMON/KANAL/VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      DOUBLE PRECISION C(33,9,9)
      DIMENSION ENMIT(32),AZMIT(7),ELMIT(7),CMITT(32),DIFF(32),CORRE(32)
     1,CDIFF(32)
      LOGICAL HDM,EXTRA
C
C     DIESE SUBROUTINE EXTRAPOLIERT DIE RANDKANAELE DER GEMESSENEN VER-
C      TEILUNGNS-FUNKTION
C
      DO 5 I=1,LC1
    5 ENMIT(I) = 0.
      DO 6 J=1,LC2
    6 AZMIT(J) = 0.
      DO 7 K=1,LC3
    7 ELMIT(K) = 0.
C     WRITE(6,8) MEN,MAZ,MEL,LC1,LC2,LC3
    8 FORMAT(3X,3H   6I2)
      ZMAX = 0.
      DO 10 I=1,LC1
       DO 10 J=1,LC2
        DO 10 K=1,LC3
      ENMIT(I) = ENMIT(I) + FIS1(K,J,I)
      AZMIT(J) = AZMIT(J) + FIS1(K,J,I)
      ELMIT(K) = ELMIT(K) + FIS1(K,J,I)
      IF ( FIS1(K,J,I) .LT. ZMAX ) GO TO 10
      ZMAX = FIS1(K,J,I)
      IMV = I
      IMA = J
      IME = K
   10 CONTINUE
C      WRITE(6,88) (AZMIT(J),J=1,LC2)
C      WRITE(6,88) (ELMIT(K),K=1,LC3)
C      IF ( HDM ) GO TO 91
C      WRITE(6,90) (ENMIT(I),I=1,LC1)
C      GO TO 93
C   91 CONTINUE
C      WRITE(6,92) (ENMIT(I),I=1,16)
C      WRITE(6,92) (ENMIT(I),I=1,LC1)
C   88 FORMAT(3X,7F7.1)
C   90 FORMAT(3X,9F7.1)
C   92 FORMAT(16F7.1)
C   93 CONTINUE
      CALL RMIMA(ENMIT,LC1,ENMI,ENMA,MIV,MV)
      CALL RMIMA(AZMIT,LC2,AZMI,AZMA,MIA,MA)
      CALL RMIMA(ELMIT,LC3,ELMI,ELMA,MIE,ME)
      MENDIF = MEN
      MAZDIF = MAZ
      MELDIF = MEL
C
C       array expansion
C
      LC2 = LC2 + 2
      LC3 = LC3 + 2
      LC21 = LC2 - 1
      LC31 = LC3 - 1
      LC22 = LC2 - 2
      LC32 = LC3 - 2
      DO 11 I = 1,33
       DO 11 J =1,9
        DO 11 K = 1,9
   11 C(I,J,K) = 0.
      LL = 1
      IF ( HDM ) LL=0
C
C       FALLS HDM = .TRUE.  WIRD KEINE EXTRAPOLATION BEZUEGLICH DER
C       ENERGIE MEHR VORGENONMMEN
C
      LC1 = LC1 + 2*LL
      LC11 = LC1 - LL
      LC12 = LC1 - 2*LL
C
C          EINLESEN VON C(I,J,K) AUS FIS1(K,J,I)
C
      LST = 1 + LL
      L1 = 0
      L2 = 0
      L3 = 0
      IF ( MA .EQ. 1 ) L2 = -1
      IF ( MA .EQ. LC22 ) L2 = +1
      IF ( ME .EQ. 1 ) L3 = -1
      IF ( ME .EQ. LC32 ) L3 = +1
C     WRITE(6,1222) MV,MA,ME
 1222 FORMAT(6H   MV=I2,4H MA=I2,4H ME=I2)
      DO 12 M1 = LST,LC11
       DO 12 M2 = 2,LC21
        DO 12 M3 = 2,LC31
   12 C(M1-L1,M2-L2,M3-L3) = FIS1(M3-1,M2-1,M1-LL)
      IMEN = MEN
      MEN = MEN - 1
      MAZ = MAZ - 1
      MEL = MEL - 1
      IF ( MA .EQ. LC22 ) MAZ = MAZ + 1
      IF ( ME .EQ. LC32 ) MEL = MEL + 1
      IF ( MA .EQ. 1 ) MAZ = MAZ - 1
      IF ( ME .EQ. 1 ) MEL = MEL - 1
      IF ( HDM ) MEN = IMEN
      MENDIF = MENDIF - MEN
      MAZDIF = MAZDIF - MAZ
      MELDIF = MELDIF - MEL
      IF ( HDM ) GO TO 102
C
C      ENERGIEEXTRAPOLATION * SYMMETRISCH
C
      DO 101 M2 = 2,LC21
      DO 101 M3 = 2,LC31
      X2 = C(2,M2,M3)
      XLC11 = C(LC11,M2,M3)
      X1 = X2*X2/(C(3,M2,M3) + 0.1)
      XLC1 = XLC11*XLC11/(C(LC12,M2,M3) + 0.1)
      IF ( X1 .GT. X2 )  X1 = X2*ENMIT(1)/(ENMIT(2) +0.1)
      IF ( XLC1 .GT. XLC11 ) XLC1 = XLC11*ENMIT(LC12)/(ENMIT(LC12-1) +
     1 0.1)
      IF ( X1 .GT. X2 ) X1 = X2
      IF ( XLC1 .GT.  XLC11 ) XLC1 = XLC11
      C(1,M2,M3) = X1
      C(LC1,M2,M3) = XLC1
  101 CONTINUE
  102 CONTINUE
C
C       AZIMUTEXTRAPOLATION
C
      IF ( MA .EQ. 1) GO TO 310
      IF ( MA .EQ. LC22 ) GO TO 330
      DO 201 M1 = 1,LC1
      DO 201 M3 = 2,LC31
      X2 = C(M1,2,M3)
      XLC21 = C(M1,LC21,M3)
      X1 = X2*X2/(C(M1,3,M3) + 0.1)
      XLC2 = XLC21*XLC21/(C(M1,LC22,M3) + 0.1)
      IF ( X1 .GT. X2 ) X1 = X2*AZMIT(1)/(AZMIT(2) + 0.1)
      IF ( XLC2 .GT. XLC21 ) XLC2 = XLC21*AZMIT(LC22)/(AZMIT(LC22-1) +
     1 0.1)
      IF ( X1 .GT. X2 ) X1 = X2
      IF ( XLC2 .GT. XLC21 ) XLC2 = XLC21
      C(M1,1,M3) = X1
      C(M1,LC2,M3) = XLC2
  201 CONTINUE
      GO TO 40
C
C  MAXIMUM AM UNTEREN AZ-KANAL
C
  310 CONTINUE
      DO 32 I=1,LC1
      DO 32 K=2,LC31
       IF(C(I,3,K).LT.C(I,4,K)) GO TO 31
C  I/K MAX TATSAECHLICH AM RAND
      C(I,2,K) = C(I,4,K)
       C(I,1,K) = C(I,2,K)**2/(C(I,3,K)+0.1)
       GO TO 32
C  NORMALE EXTRAPOLATION
   31 CONTINUE
      C(I,2,K) = C(I,3,K)**2/(C(I,4,K)+0.1)
      C(I,1,K) = C(I,2,K)**2/(C(I,3,K)+0.1)
   32 CONTINUE
      GO TO 40
C
C MAXIMUM AM OBEREN AZ-KANAL
C
  330 DO 35 I=1,LC1
      DO 35 K=2,LC31
       IF(C(I,LC22,K).LT.C(I,LC22-1,K)  ) GO TO 34
C  I/K MAX TATSAECHLICH AM RAND
       C(I,LC22+1,K) = C(I,LC22-1,K)
       C(I,LC22+2,K) = C(I,LC22+1,K)**2/(C(I,LC22,K)+0.1)
       GO TO 35
C  NORMALE EXTRAPOLATION
   34  C(I,LC22+1,K) = C(I,LC22,K)**2/(C(I,LC22-1,K)+0.1)
       C(I,LC22+2,K) = C(I,LC22+1,K)**2/(C(I,LC22,K)+0.1)
   35 CONTINUE
C
   40 CONTINUE
C
C           ELEVATIONSEXTRAPOLATION
C
      IF ( ME .EQ. 1 ) GO TO 410
      IF ( ME .EQ. LC32 ) GO TO 430
      DO 301 M1 = 1,LC1
      DO 301 M2 = 1,LC2
      X2 = C(M1,M2,2)
      XLC31 = C(M1,M2,LC31)
      X1 = X2*X2/(C(M1,M2,3) + 0.1)
      XLC3 = XLC31*XLC31/(C(M1,M2,LC32) + 0.1)
      IF ( X1 .GT. X2 ) X1 = X2*ELMIT(1)/(ELMIT(2) + 0.1)
      IF ( XLC3 .GT. XLC31 ) XLC3 = XLC31*ELMIT(LC32)/(ELMIT(LC32-1) +
     1 0.1)
      IF ( X1 .GT. X2 ) X1 = X2
      IF ( XLC3 .GT. XLC31 ) XLC3 = XLC31
      C(M1,M2,1) = X1
      C(M1,M2,LC3) = XLC3
  301 CONTINUE
      GO TO 50
C
C  MAXIMUM AM UNTEREN EL-KANAL
C
  410 CONTINUE
      DO 42 I=1,LC1
      DO 42 J=1,LC2
       IF(C(I,J,3).LT.C(I,J,4)) GO TO 41
C  I/J MAX TATSAECHLICH AM RAND
       C(I,J,2) = C(I,J,4)
       C(I,J,1)=C(I,J,2)**2/(C(I,J,3)+0.1)
       GO TO 42
C  NORMALE EXTRAPOLATION
   41 C(I,J,2) = C(I,J,3)**2/(C(I,J,4)+0.1)
      C(I,J,1) = C(I,J,2)**2/(C(I,J,3)+0.1)
   42 CONTINUE
C
      GO TO 50
C
C  MAXIMUM AM OBEREN EL-KANAL
C
  430 DO 45 I=1,LC1
      DO 45 J=1,LC2
       IF(C(I,J,LC32).LT.C(I,J,LC32-1)) GO TO 44
C  I/J MAX TATSAECHLICH AM RAND
       C(I,J,LC32+1) = C(I,J,LC32-1)
       C(I,J,LC32+2) = C(I,J,LC32+1)**2/(C(I,J,LC32-1)+0.1)
       GO TO 45
   44 CONTINUE
C  NORMALE EXTRAPOLATION
      C(I,J,LC32+1) = C(I,J,LC32)**2/(C(I,J,LC32)+0.1)
      C(I,J,LC32+2) = C(I,J,LC32+1)**2/(C(I,J,LC32)+0.1)
   45 CONTINUE
   50 CONTINUE
C
C        BERECHNUNG DER GESAMTZAEHLRATEN VOR UND NACH DER KORREKTUR UND
C        VERGLEICH MIT DEM INTEGRATIONSZAEHLER
C
      DO 181 I=1,32
      CMITT(I) = 0.
      CORRE(I) = 0.
      CDIFF(I) = 0.
      DIFF(I) = 0.
  181 CONTINUE
      DO 188 I=1,LC1
      DO 189 J=1,LC2
      DO 189 K=1,LC3
  189 CMITT(I) = CMITT(I) + C(I,J,K)
C      WRITE(6,1811) CMITT(I),I,FINT(I+IMEN-1-LL)
  188 CONTINUE
C 1811 FORMAT(5X,7H CMITT=F10.4,3H I=I2,6H FINT=F10.4)
      DO 199 I = 1,LC12
      DIFF(I) = FINT(IMEN-1+I) - ENMIT(I)
      CDIFF(I) = CMITT(I+LL)  - ENMIT(I)
  199 CONTINUE
C      DO 3001 I=1,LC12
C      WRITE(6,3002) FINT(I-1+IMEN),DIFF(I),CDIFF(I),ENMIT(I)
C 3002 FORMAT(6H FINT=F10.2,6H DIFF=F10.2,7H CDIFF=F10.2,7H ENMIT=F10.2)
C 3001 CONTINUE
      IF ( HDM ) GO TO 2071
      IF ( CMITT(1) ) 701,701,702
  701 CORRE(1) = 1.
      GO TO 703
  702 CORRE(1) = FINT(IMEN - 1)/CMITT(1)
  703 IF ( CMITT(LC1) ) 704,704,705
  704 CORRE(LC1) = 1.
      GO  TO 706
  705 CORRE(LC1) = FINT(IMEN + LC1 -2 )/CMITT(LC1)
  706 CONTINUE
      DO 2072 I=2,LC11
      IF (CDIFF(I-1)) 802,801,802
  801 CORRE(I) = 1.
      GO TO 2072
  802 IF ( DIFF(I-1) .LT. 0. )  DIFF(I-1) = 0.
      CORRE(I) = DIFF(I-1)/CDIFF(I-1)
 2072 CONTINUE
      GO TO 2073
 2071 CONTINUE
      DO 2074 I=1,LC1
      IF (CDIFF(I) ) 902,901,902
  901 CORRE(I) = 1.
      GO TO 2074
  902 IF ( DIFF(I) .LT. 0. )  DIFF(I) = 0.
      CORRE(I) = DIFF(I)/CDIFF(I)
 2074 CONTINUE
      IF ( CORRE(1) .EQ. 0. ) CORRE(1) = 1.
      IF ( CORRE(LC1) .EQ. 0. ) CORRE(LC1) = 1.
 2073 CONTINUE
C      DO 2075 I =1,LC1
C      WRITE(6,2076) I,CORRE(I)
C 2076 FORMAT(3H I=I2,7H CORRE=F10.2)
C 2075 CONTINUE
C
C        KORREKTUREN 
C
C             EN correction
C
      DO 2201 M02=2,LC21
      DO 2201 M03=2,LC31
       C(1,M02,M03) = C(1,M02,M03)*CORRE(1)
 2201  C(LC1,M02,M03) = C(LC1,M02,M03)*CORRE(LC1)
C
C             AZ correction
C
      IF(MA.EQ.1) GO TO 203
      IF(MA.EQ.LC22) GO TO 204
      DO 202 M04=1,LC1
      DO 202 M06=2,LC31
       C(M04,1,M06) = C(M04,1,M06)*CORRE(M04)
  202  C(M04,LC2,M06) = C(M04,LC2,M06)*CORRE(M04)
      GO TO 205
  203 DO 213 M04=1,LC1
      DO 213 M06=2,LC31
       C(M04,1,M06) = C(M04,1,M06)*CORRE(M04)
  213  C(M04,2,M06) = C(M04,2,M06)*CORRE(M04)
      GO TO 205
  204 DO 214 M04=1,LC1
      DO 214 M06=2,LC31
       C(M04,LC2,M06) = C(M04,LC2,M06)*CORRE(M04)
  214  C(M04,LC2-1,M06) = C(M04,LC2-1,M06)*CORRE(M04)
  205 CONTINUE
C
C            EL correction
C
      IF(ME.EQ.1) GO TO 207
      IF(ME.EQ.LC32) GO TO 208
      DO 206 M07=1,LC1
      DO 206 M08=1,LC2
       C(M07,M08,1) =   C(M07,M08,1)*CORRE(M07)
  206  C(M07,M08,LC3) = C(M07,M08,LC3)*CORRE(M07)
      GO TO 209
  207 DO 217 M07=1,LC1
      DO 217 M08=1,LC2
       C(M07,M08,1) =   C(M07,M08,1)*CORRE(M07)
  217  C(M07,M08,2) = C(M07,M08,2)*CORRE(M07)
      GO TO 209
  208 DO 218 M07=1,LC1
      DO 218 M08=1,LC2
       C(M07,M08,LC3) =   C(M07,M08,1)*CORRE(M07)
  218  C(M07,M08,LC3-1) = C(M07,M08,LC3-1)*CORRE(M07)
  209 CONTINUE
C
C        check reliability of the extrapolation
C
      DO 600 I=1,LC1
      CMITT(I) = 0.
      DO 601 J=1,LC2
      DO 601 K=1,LC3
  601 CMITT(I) = CMITT(I) + C(I,J,K)
C      WRITE(6,1811) CMITT(I),I,FINT(I+MEN-1)
  600 CONTINUE
C
C       WRITE(6,8) MEN,MAZ,MEL,LC1,LC2,LC3
c      DO 1001 I = 1,LC1
c      DO 1002 J = 1,LC2
c      WRITE(6,1003) ( C(I,J,K) ,K =1,LC3)
c 1003 FORMAT(' ',9F10.2)
c 1002 CONTINUE
c      WRITE(6,1004)
c 1004 FORMAT(' ',/)
c 1001 CONTINUE
C
      RETURN
      END
C
C##############################################################################
C
      SUBROUTINE RMIMA(X,N,RMIN,RMAX,IMIN,IMAX)
C  BESTIMMT MINIMUM UND MAXIMUM WERTE AUS EINEM ARRAY <N> LANG
C  RMIN-RMAX : WERTE
C  IMIN-IMAX : INDEX IHRES AUFRTETENS IM ARRAY
      DIMENSION X(1)
C
      RMIN=X(1)
      IMIN=1
      RMAX=X(1)
      IMAX=1
      IF(N.EQ.1) RETURN
      DO 1 I=2,N
        IF(RMIN.LE.X(I)) GO TO 2
        RMIN=X(I)
        IMIN=I
    2   IF(RMAX.GE.X(I)) GO TO 1
        RMAX=X(I)
        IMAX=I
    1 CONTINUE
      RETURN
      END
C                                                                              
C############################################################################
C                                                                              
      SUBROUTINE GAUSSQ(MV,MA,ME,NPRO1,NPRO2,NPRO3,LL,N,M,L,W)
      COMMON /GQUA/ S(4),NZ1,NZ2,NZ3,LLC,NC,MC,LC
      COMMON/A1/MVC,MAC,MEC,NPRO1C,NPRO2C,NPRO3C
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      COMMON/A3/ A(24),B(24),D(24),EG5,FG5,GG5
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      DOUBLE PRECISION S,Y3,FN3
      DOUBLE PRECISION EG5,FG5,GG5,VC(24),A,B,D
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      EXTERNAL FN3
C
C          INTEGRATIONSROUTINE UND DAMIT ZUSAMMENHAENGENDE UNTERPROGRAMME
C 
      MVC = MV
      MAC = MA
      MEC = ME
      NPRO1C = NPRO1
      NPRO2C = NPRO2
      NPRO3C = NPRO3
      NZ1=1
      NZ2=1
      NZ3=1
      LLC=LL
      NC=N
      MC=M
      LC=L
C      write(6,*) 'GAUSSQ',NPRO3,MVC,MAC,MEC,LL,N,M,L
C        
C      CALL COEFS(DVU,DVO,NPRO1,VC)
C      DO 401 I=1,16
C  401 A(I)=VC(I)
C      EG5=DVO-DVU
C      CALL COEFS(DAZU,DAZO,NPRO2,VC)
C      DO 402 J=1,16
C  402 B(J)=VC(J)
C      FG5=DAZO-DAZU
C      CALL COEFS(DELU,DELO,NPRO3,VC)
C      DO 403 K=1,16
C  403 D(K)=VC(K)
C      GG5=DELO-DELU
C              
      GO TO (1,2,3), NPRO3
    1 CALL DTG8 (DELU,DELO,FN3,Y3)
      GO TO 10
    2 CALL DTG12(DELU,DELO,FN3,Y3)
      GO TO 10
    3 CALL DTG16(DELU,DELO,FN3,Y3)
   10 W = SNGL(Y3)*UMRF
      RETURN
      END
C************************************************:
      DOUBLE PRECISION FUNCTION FN3(X3)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON /GQUA/ S(4),NZ1,NZ2,NZ3,LL,N,M,L
      COMMON/A1/MV,MA,ME,NPRO1,NPRO2,NPRO3
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      DOUBLE PRECISION S,X3,Y2,FN2
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      EXTERNAL FN2
      S(3)=X3
      NZ2=1
      GO TO (1,2,3),NPRO2
    1 CALL DSG8 (DAZU,DAZO,FN2,Y2)
      GO TO 10
    2 CALL DSG12(DAZU,DAZO,FN2,Y2)
      GO TO 10
    3 CALL DSG16(DAZU,DAZO,FN2,Y2)
   10 FN3=Y2*UMRF
      NZ3=NZ3+1
      RETURN
      END
C**************************************
      DOUBLE PRECISION FUNCTION FN2(X2)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON /GQUA/ S(4),NZ1,NZ2,NZ3,LL,N,M,L
      COMMON/A1/MV,MA,ME,NPRO1,NPRO2,NPRO3
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      DOUBLE PRECISION S,X2,Y1,FN1
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      EXTERNAL FN1
      S(2)=X2
      NZ1=1
      GO TO (1,2,3),NPRO1
    1 CALL DRG8 (DVU,DVO,FN1,Y1)
      GO TO 10
    2 CALL DRG12(DVU,DVO,FN1,Y1)
      GO TO 10
    3 CALL DRG16(DVU,DVO,FN1,Y1)
   10 FN2=Y1
      NZ2=NZ2+1
      RETURN
      END
C*********************************************
      DOUBLE PRECISION FUNCTION FN1(X1)
      DIMENSION SP(16,16,16)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON /GQUA/ S(4),NZ1,NZ2,NZ3,LL,N,M,L
      COMMON/A1/MV,MA,ME,NPRO1,NPRO2,NPRO3
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      DOUBLE PRECISION S,X1,CDS2,SDS2,CDS3,SDS3,SP
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      LOGICAL IPL
      IPL = .TRUE.
      CDS2= 1.
      SDS2= 1.
      SDS3= 1.
      S(1) = X1
      LS = LL + N + M + L
C      write(6,*) 'FN1', LS,N,M,L,ipl,umrf
      IF(LS) 100,100,200
C
  100 CALL CINT(S,IPL,IERR)
C      Write(6,*) 'ierr', IERR
      IF ( IERR .EQ. 1 ) GO TO 101
      S(4)=DEXP(S(4))
      GO TO 102
  101 S(4) = 0.D0
C      Write(6,*) 'index',NZ1,NZ2,NZ3
  102 SP(NZ1,NZ2,NZ3) = S(4)
C
C       write(6,*) 'sp', SP(NZ1,NZ2,NZ3) 
C      WRITE(6,1423) (I,S(I),I=1,4 )
C 1423 FORMAT(4(3H S(I1,2H)=E10.4))
C
  200 IF ( N .NE. 0 )   CDS2= DCOS(UMRF*S(2))**N
      IF ( M .NE. 0 )   SDS2= DSIN(UMRF*S(2))**M
      IF ( L .NE. 0 )   SDS3= DSIN(UMRF*S(3))**L
      CDS3= DCOS(UMRF*S(3))**(M + N + 1)
C      write(6,*) 'sin/cos',SDS2,SDS3,CDS2,CDS3
      FN1=SP(NZ1,NZ2,NZ3)*(S(1))**(LS+2)*SDS2*SDS3*CDS2*CDS3
      NZ1=NZ1+1
C      Write(6,*) 'Function FN1', FN1
      RETURN
      END
C###############################################################################
C
      SUBROUTINE GRENZ(F,V,AZ,EL,MV,MA,ME,INST,ICOUNT,IGREN)
c
c*********************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL PI,PIH,UMRF
      DIMENSION F(33,9,9),V(33),AZ(9),EL(9),FELD(33)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON /A7/ MEN,MAZ,MEL,LC1,LC2,LC3
      COMMON/A2/DVO,DVU,DAZO,DAZU,DELO,DELU
      DATA FELD/33*0d0/
      CORR = 2.5D8
      IF ( INST .EQ. 3 ) CORR=1.D0/1.628D-10
      CORR = CORR*FLOAT(ICOUNT)
      WRITE(6,1077) F(MV,MA,ME)
 1077 FORMAT(13H F(MV,MA,ME)=E12.4)

      DO 111 M1=2,LC1
      MV7=LC1-M1+1
      CG = CORR/(V(MV7)**4)
      DO 111 M2=1,LC2
      DO 111 M3=1,LC3
      CGR = CG/DCOS(UMRF*EL(M3))
      IF(F(MV7,M2,M3).GT.CGR) GO TO 112
  111 CONTINUE
      GO TO 820
  112 M1P1=MV7+1
      DVO= V(M1P1)

      DO 113 M1=2,M1P1
      CG = CORR/(V(M1)**4)
      DO 113 M2=1,LC2
      DO 113 M3=1,LC3
      CGR = CG/DCOS(UMRF*EL(M3))
      IF(F(M1,M2,M3).GT.CGR) GO TO 114
  113 CONTINUE
      GO TO 820

  114 M1M1=M1-1
      DVU= V(M1M1)
      DO 115 M2=2,LC2
      MA7=LC2-M2+1
      DO 115 M1=M1M1,M1P1
      CG = CORR/(V(M1)**4)
      DO 115 M3=1,LC3
      CGR = CG/DCOS(EL(M3)*UMRF)
      IF(F(M1,MA7,M3).GT.CGR) GO TO 116
  115 CONTINUE
      GO TO 820

  116 M2P1=MA7+1
      DAZO= AZ(M2P1)
c*************************************************************

      DO 117 M2=2,M2P1
      DO 117 M1=M1M1,M1P1
      CG = CORR/(V(M1)**4)
      DO 117 M3=1,LC3
      CGR = CG/DCOS(UMRF*EL(M3))
      IF(F(M1,M2,M3).GT.CGR) GO TO 118
  117 CONTINUE

      GO TO 820

  118 M2M1=M2-1
      DAZU= AZ(M2M1)
c****************************************************************

      DO 119 M3=2,LC3
      ME7=LC3-M3+1
      CG = CORR/DCOS(UMRF*EL(ME7))
      DO 119 M1=M1M1,M1P1
      CGR = CG/(V(M1)**4)
      DO 119 M2=M2M1,M2P1
      IF(F(M1,M2,ME7).GT.CGR) GO TO 120
  119 CONTINUE
      GO TO 820
  120 M3P1 =ME7+1
      DELO= EL(M3P1)
      DO 121 M3=2,M3P1
      CG = CORR/DCOS(UMRF*EL(M3))
      DO 121 M1=M1M1,M1P1
      CGR = CG/(V(M1)**4)
      DO 121 M2=M2M1,M2P1
      IF(F(M1,M2,M3).GT.CGR) GO TO 122
  121 CONTINUE
      GO TO 820
  122 M3M1=M3-1
      DELU= EL(M3M1)
C      WRITE(6,825) M1M1,M1P1,M2M1,M2P1,M3M1,M3P1
C  825 FORMAT(6I5)
C     IF ( M1P1 .LT. 33 ) GO TO 200
      LLC1 = M1P1 - M1M1 + 1
      LLC2 = M2P1 - M2M1 + 1
      LLC3 = M3P1 - M3M1 + 1
      IF ( LLC1 - 2*(LLC1/2)  .EQ. 0 ) GO TO 701
   71 IF ( LLC2 - 2*(LLC2/2)  .EQ. 0 ) GO TO 702
   72 IF ( LLC3 - 2*(LLC3/2)  .EQ. 0 ) GO TO 703
      GO TO 999
  701 IF ( M1P1 .EQ. LC1)  GO TO 7001
      M1P1 = M1P1 + 1
      LLC1 = LLC1 + 1
      GO TO 71
 7001 M1M1 = M1M1 - 1
      LLC1 = LLC1 + 1
      GO TO 71
  702 IF ( M2P1 .EQ. LC2) GO TO 7002
      M2P1 = M2P1 + 1
      LLC2 = LLC2 + 1
      GO TO 72
 7002 M2M1 = M2M1 - 1
      LLC2 = LLC2 + 1
      GO TO 72
  703 IF ( M3P1 .EQ. LC3) GO TO 7003
      M3P1 = M3P1 + 1
      LLC3  =  LLC3 + 1
      GO TO 999
 7003 M3M1 = M3M1 - 1
      LLC3 = LLC3 + 1
  999 CONTINUE
      WRITE(6,1001) LLC1,LLC2,LLC3
 1001 FORMAT('   Grenze',6H LLC1=I2,6H LLC2=I2,6H LLC3=I2)
      IF ( LLC1 .EQ. LC1 ) GO TO 211
      DO 1002 I=1,LLC1
 1002 FELD(I) = V(M1M1+I-1)
      DO 1003 I=1,LLC1
      V(I) = FELD(I)
 1003 FELD(I) = 0.
  211 IF( LLC2 .EQ. LC2 ) GO TO 222
      DO 1004 J=1,LLC2
 1004 FELD(J) = AZ(M2M1+J-1)
      DO 1005 J=1,LLC2
      AZ(J) = FELD(J)
 1005 FELD(J) = 0.
  222 IF ( LLC3 .EQ. LC3 ) GO TO 233
      DO 1006 K=1,LLC3
 1006 FELD(K) = EL(M3M1+K-1)
      DO 1007 K=1,LLC3
      EL(K) = FELD(K)
 1007 FELD(K) = 0.
      DO 1008 I=1,LC1
      DO 1008 J=1,LC2
      DO 1009 K=1,LLC3
 1009 FELD(K) =  F(I,J,K+M3M1-1)
      DO 1010 K=1,LLC3
      F(I,J,K) = FELD(K)
 1010 FELD(K) = 0.
 1008 CONTINUE
  233 IF( LLC2 .EQ. LC2 ) GO TO 244
      DO 1011 I=1,LC1
      DO 1011 K=1,LLC3
      DO 1012 J=1,LLC2
 1012 FELD(J) = F(I,J+M2M1-1,K)
      DO 1013 J=1,LLC2
      F(I,J,K) = FELD(J)
 1013 FELD(J) = 0.
 1011 CONTINUE
  244 IF (LLC1 .EQ. LC1 ) GO TO 255
      DO 1014 J=1,LLC2
      DO 1014 K=1,LLC3
      DO 1015 I=1,LLC1
 1015 FELD(I) = F(I+M1M1-1,J,K)
      DO 1016 I=1,LLC1
      F(I,J,K) = FELD(I)
 1016 FELD(I) = 0.
 1014 CONTINUE
      LC1 = LLC1
  255 LC2 = LLC2
      LC3 = LLC3
      MV = MV - M1M1 + 1
      MA = MA - M2M1 + 1
      ME = ME - M3M1 + 1
C      WRITE(6,1088) F(MV,MA,ME),MV,MA,ME
C 1088 FORMAT(16H F(MV',MA',ME')=E12.4,3I5)
      IGREN = 0
      GO TO 200
  820 WRITE(6,821)
  821 FORMAT(23H  GRENZE NICHT GEFUNDEN)
      IGREN = 1
  200 RETURN
      END
C******************************************
       SUBROUTINE COEFS(XL,XU,NPRO,VC)
       DOUBLE PRECISION XL,XU,A,B,C,VC(24)
       GO TO (1,2,3)  , NPRO
    1 A=.5D0*(XU+XL)
      B=XU-XL
      C=.48014492824876812D0*B
      VC(1)=A+C
      VC(2)=A-C
      C=.39833323870681337D0*B
      VC(3)=A+C
      VC(4)=A-C
      C=.26276620495816449D0*B
      VC(5)=A+C
      VC(6)=A-C
      C=.9171732124782490D-1*B
      VC(7)=A+C
      VC(8)=A-C
      GO TO 4
    2 A=.5D0*(XU+XL)
      B=XU-XL
      C=.49078031712335963D0*B
      VC(1)=A+C
      VC(2)=A-C
      C=.45205862818523743D0*B
      VC(3)=A+C
      VC(4)=A-C
      C=.38495133709715234D0*B
      VC(5)=A+C
      VC(6)=A-C
      C=.29365897714330872D0*B
      VC(7)=A+C
      VC(8)=A-C
      C=.18391574949909010D0*B
      VC(9)=A+C
      VC(10)=A-C
      C=.62616704255734458D-1*B
      VC(11)=A+C
      VC(12)=A-C
      GO TO 4
    3 A=.5D0*(XU+XL)
      B=XU-XL
      C=.49470046749582497D0*B
      VC(1)=A+C
      VC(2)=A-C
      C=.47228751153661629D0*B
      VC(3)=A+C
      VC(4)=A-C
      C=.43281560119391587D0*B
      VC(5)=A+C
      VC(6)=A-C
      C=.37770220417750152D0*B
      VC(7)=A+C
      VC(8)=A-C
      C=.30893812220132187D0*B
      VC(9)=A+C
      VC(10)=A-C
      C=.22900838882861369D0*B
      VC(11)=A+C
      VC(12)=A-C
      C=.14080177538962946D0*B
      VC(13)=A+C
      VC(14)=A-C
      C=.47506254918818720D-1*B
      VC(15)=A+C
      VC(16)=A-C
    4 RETURN
      END
C*****************************************
      SUBROUTINE DRG8(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.50614268145188130D-1*(FCT(A(1))+FCT(A(2)))
      Y=Y+.11119051722668724D0*(FCT(A(3))+FCT(A(4)))
      Y=Y+.15685332293894364D0*(FCT(A(5))+FCT(A(6)))
      Y=E*(Y+.18134189168918099D0*(FCT(A(7))+FCT(A(8))))
      RETURN
      END
C*****************************************
      SUBROUTINE DSG8(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.50614268145188130D-1*(FCT(B(1))+FCT(B(2)))
      Y=Y+.11119051722668724D0*(FCT(B(3))+FCT(B(4)))
      Y=Y+.15685332293894364D0*(FCT(B(5))+FCT(B(6)))
      Y=F*(Y+.18134189168918099D0*(FCT(B(7))+FCT(B(8))))
      RETURN
      END
C***********************************
      SUBROUTINE DTG8(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.50614268145188130D-1*(FCT(D(1))+FCT(D(2)))
      Y=Y+.11119051722668724D0*(FCT(D(3))+FCT(D(4)))
      Y=Y+.15685332293894364D0*(FCT(D(5))+FCT(D(6)))
      Y=G*(Y+.18134189168918099D0*(FCT(D(7))+FCT(D(8))))
      RETURN
      END
C***********************************
      SUBROUTINE DRG12(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.23587668193255914D-1*(FCT(A(1))+FCT(A(2)))
      Y=Y+.53469662997659215D-1*(FCT(A(3))+FCT(A(4)))
      Y=Y+.8003916427167311D-1*(FCT(A(5))+FCT(A(6)))
      Y=Y+.10158371336153296D0*(FCT(A(7))+FCT(A(8)))
      Y=Y+.11674626826917740D0*(FCT(A(9))+FCT(A(10)))
      Y=E*(Y+.12457352290670139D0*(FCT(A(11))+FCT(A(12))))
      RETURN
      END
C*************************************
      SUBROUTINE DSG12(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.23587668193255914D-1*(FCT(B(1))+FCT(B(2)))
      Y=Y+.53469662997659215D-1*(FCT(B(3))+FCT(B(4)))
      Y=Y+.8003916427167311D-1*(FCT(B(5))+FCT(B(6)))
      Y=Y+.10158371336153296D0*(FCT(B(7))+FCT(B(8)))
      Y=Y+.11674626826917740D0*(FCT(B(9))+FCT(B(10)))
      Y=F*(Y+.12457352290670139D0*(FCT(B(11))+FCT(B(12))))
      RETURN
      END
C****************************************
      SUBROUTINE DTG12(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.23587668193255914D-1*(FCT(D(1))+FCT(D(2)))
      Y=Y+.53469662997659215D-1*(FCT(D(3))+FCT(D(4)))
      Y=Y+.8003916427167311D-1*(FCT(D(5))+FCT(D(6)))
      Y=Y+.10158371336153296D0*(FCT(D(7))+FCT(D(8)))
      Y=Y+.11674626826917740D0*(FCT(D(9))+FCT(D(10)))
      Y=G*(Y+.12457352290670139D0*(FCT(D(11))+FCT(D(12))))
      RETURN
      END
C****************************************************
      SUBROUTINE DRG16(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.13576229705877047D-1*(FCT(A(1))+FCT(A(2)))
      Y=Y+.31126761969323946D-1*(FCT(A(3))+FCT(A(4)))
      Y=Y+.47579255841246392D-1*(FCT(A(5))+FCT(A(6)))
      Y=Y+.62314485627766936D-1*(FCT(A(7))+FCT(A(8)))
      Y=Y+.7479799440828837D-1*(FCT(A(9))+FCT(A(10)))
      Y=Y+.8457825969750127D-1*(FCT(A(11))+FCT(A(12)))
      Y=Y+.9130170752246179D-1*(FCT(A(13))+FCT(A(14)))
      Y=E*(Y+.9472530522753425D-1*(FCT(A(15))+FCT(A(16))))
      RETURN
      END
C*****************************************************
      SUBROUTINE DSG16(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.13576229705877047D-1*(FCT(B(1))+FCT(B(2)))
      Y=Y+.31126761969323946D-1*(FCT(B(3))+FCT(B(4)))
      Y=Y+.47579255841246392D-1*(FCT(B(5))+FCT(B(6)))
      Y=Y+.62314485627766936D-1*(FCT(B(7))+FCT(B(8)))
      Y=Y+.7479799440828837D-1*(FCT(B(9))+FCT(B(10)))
      Y=Y+.8457825969750127D-1*(FCT(B(11))+FCT(B(12)))
      Y=Y+.9130170752246179D-1*(FCT(B(13))+FCT(B(14)))
      Y=F*(Y+.9472530522753425D-1*(FCT(B(15))+FCT(B(16))))
      RETURN
      END
C*****************************************
      SUBROUTINE DTG16(XL,XU,FCT,Y)
       DIMENSION VC(24)
       COMMON/A3/ A(24),B(24),D(24),E,F,G
       DOUBLE PRECISION A,B,D,VC,E,F,G
      DOUBLE PRECISION XL,XU,Y,FCT
      Y=.13576229705877047D-1*(FCT(D(1))+FCT(D(2)))
      Y=Y+.31126761969323946D-1*(FCT(D(3))+FCT(D(4)))
      Y=Y+.47579255841246392D-1*(FCT(D(5))+FCT(D(6)))
      Y=Y+.62314485627766936D-1*(FCT(D(7))+FCT(D(8)))
      Y=Y+.7479799440828837D-1*(FCT(D(9))+FCT(D(10)))
      Y=Y+.8457825969750127D-1*(FCT(D(11))+FCT(D(12)))
      Y=Y+.9130170752246179D-1*(FCT(D(13))+FCT(D(14)))
      Y=G*(Y+.9472530522753425D-1*(FCT(D(15))+FCT(D(16))))
      RETURN
      END
C
C############################################################################
C#########################################################################
C
      SUBROUTINE CINTI (F,XZ,YZ,ZZ,MVZ,MAZZ,MEZ)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL PI,PIH,UMRF
      DIMENSION F(33,9,9),S(4),A(3,3,3),X(33),Y(9),Z(9),VDX(3,3,3),
     1IK(3),JK(3),KK(3),AC2(2,2),BC2(2),AC4V(2,2,2,2),AC4(4,4),
     2BC4V(2,2),BC4(4),AC8V(2,2,2,2,2,2),AC8(8,8),BC8(8),FV2(3),FV3(3),
     3FV4(2,2,2),AQA(3,3,3,16,4,4),AQM(3,3,3),XZ(33),YZ(9),ZZ(9)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      LOGICAL IPL
      write(6,*) ' entry of cinti'
C
C       ROUTINEN ZUR BESTIMMUNG DER VERTEILUNGSFUNKTION 
C
      DO 101 I=1,LC1
  101 X(I) = XZ(I)
      DO 102 J=1,LC2
  102 Y(J) = YZ(J)
      DO 103 K=1,LC3
  103 Z(K) = ZZ(K)
      MV = MVZ
      MA = MAZZ
      ME = MEZ
      IK(1)=1
      IK(2)=0
      IK(3)=0
      JK(1)=0
      JK(2)=1
      JK(3)=0
      KK(1)=0
      KK(2)=0
      KK(3)=1
C 
c      WRITE(6,1113) (((F(I,J,K),I,J,K,I=1,LC1),J=1,LC2),K=1,LC3)
c 1113 FORMAT(4H  F=D20.8,4H  I=I5,4H  J=I5,4H  K=I5)
c      WRITE(6,1114) LC1,LC2,LC3
c 1114 FORMAT(6H  LC1=I5,6H  LC2=I5,6H  LC3=I5)
c      WRITE(6,1115)   (   X(I),I,I=1,LC1)
c 1115 FORMAT(7H  DGES=D20.8,4H  I=I5)
c      WRITE(6,1116) (  Y(J),J,J=1,LC2)
c 1116 FORMAT(6H  DAZ=D20.8,4H  J=I5)
c      WRITE(6,1117) (  Z(K),K,K=1,LC3)
c 1117 FORMAT(6H  DEL=D20.8,4H  K=I5)
C 
      DO 601 IQ=1,LH1
      I=2*IQ+1
      DO 601 JQ=1,LH2
      J=2*JQ+1
      DO 601 KQ=1,LH3
      K=2*KQ+1
      A(1,1,1) = F(I-2,J-2,K-2)
      VDX(1,1,1) =0.
      VDX(2,1,1)= X(I-1)-X(I-2)
      VDX(3,1,1) = X(I)-X(I-2)
      VDX(1,2,1) = Y(J-1)-Y(J-2)
      VDX(1,3,1) = Y(J)-Y(J-2)
      VDX(1,1,2) = Z(K-1) - Z(K-2)
      VDX(1,1,3) = Z(K) - Z(K-2)
C 
      DO 1 I1=1,3
      DO 2 L1=1,2
      IA=L1*IK(I1)
      JA=L1*JK(I1)
      KA=L1*KK(I1)
      DO 3 L2=1,2
    3 AC2(L1,L2)=(VDX(1+IA,1+JA,1+KA))**L2
    2 BC2(L1)=(F(I-2+IA,J-2+JA,K-2+KA)-F(I-2,J-2,K-2))
C
C     write(6,*) 'call DGELG 1'
      CALL DGELG(BC2,AC2,2,1,1.0D-15,IER)
C     write(6,*) 'ier 1', IER
C
      DO 4 L1=1,2
      IB=L1*IK(I1)+1
      JB=L1*JK(I1)+1
      KB=L1*KK(I1)+1
      IG = IB+JB+KB-3
    4 A(IB,JB,KB)=BC2(IG)
    1 CONTINUE
C
      DO 5 I1=1,3
      L1=1
      DO6 L12=1,2
      DO6 L11=1,2
      JO=L11*IK(I1)
      KO=L11*JK(I1)
      IO=L11*KK(I1)
      IT=L12*IK(I1)
      JT=L12*JK(I1)
      KT=L12*KK(I1)
      L2=1
      DO 7 L22=1,2
      DO7 L21=1,2
      AC4V(L11,L12,L21,L22)=(VDX(1+IO,1+JO,1+KO))**L21*(VDX(1+IT,1+JT,
     11+KT))**L22
      AC4(L1,L2)=AC4V(L11,L12,L21,L22)
    7 L2=L2+1
      FV1=F(I-2,J-2,K-2)
      FV2(1)=F(I-2+IO,J-2+JO,K-2+KO)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2+IT,J-2+JT,K-2+KT)-F(I-2,J-2,K-2)
      BC4V(L11,L12)=F(I-2+IO+IT,J-2+JO+JT,K-2+KO+KT)-FV2(1)-FV2(2)-FV1
      BC4(L1)=BC4V(L11,L12)
    6 L1=L1+1
C 
C     Write(6,*) 'call DGELG 2'
      CALL DGELG(BC4,AC4,4,1,1.0D-15,IER)
C    
C     write(6,*) 'ier 2', IER
      LZ=1
      DO 8 L12=1,2
      DO 8 L11=1,2
      ID= L11*KK(I1)+L12*IK(I1)+1
      JD= L11*IK(I1)+L12*JK(I1)+1
      KD= L11*JK(I1)+L12*KK(I1)+1
      A(ID,JD,KD)=BC4(LZ)
    8 LZ = LZ+1
    5 CONTINUE
      L = 1
      DO 9 L1=1,2
      DO 9 L2=1,2
      DO 9 L3=1,2
      FV2(1)=F(I-2+L1,J-2,K-2)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2,J-2+L2,K-2)-F(I-2,J-2,K-2)
      FV2(3)=F(I-2,J-2,K-2+L3)-F(I-2,J-2,K-2)
      FV3(1)=F(I-2+L1,J-2+L2,K-2)-FV2(1)-FV2(2)-FV1
      FV3(2)=F(I-2+L1,J-2,K-2+L3)-FV2(1)-FV2(3)-FV1
      FV3(3)=F(I-2,J-2+L2,K-2+L3)-FV2(2)-FV2(3)-FV1
      FV4(L1,L2,L3)=F(I-2+L1,J-2+L2,K-2+L3)-FV3(1)-FV3(2)-FV3(3)-FV2(1)-
     1FV2(2)-FV2(3)-FV1
      BC8(L)=FV4(L1,L2,L3)
    9 L=L+1
      L1=1
      DO 10 L11=1,2
      DO 10 L12=1,2
      DO 10 L13=1,2
      L2=1
      DO 11 L21=1,2
      DO 11 L22=1,2
      DO 11 L23=1,2
      AC8V(L11,L12,L13,L21,L22,L23)=(VDX(1+L11,1,1))**L21*(VDX(1,1+L12,1
     1))**L22*(VDX(1,1,1+L13))**L23
      AC8(L1,L2) = AC8V(L11,L12,L13,L21,L22,L23)
   11 L2=L2+1
   10 L1=L1+1
C 
C     Write(6,*) 'call DGELG 3' 
      CALL DGELG(BC8,AC8,8,1,1.0D-15,IER)
C
C      write(6,*) 'ier 3', IER
      LZ=1
      DO 13 I13=2,3
      DO 13 J13=2,3
      DO 13 K13=2,3
      A(I13,J13,K13)=BC8(LZ)
   13 LZ=LZ+1
      DO 501 IQ1=1,3
      DO 501 JQ1=1,3
      DO 501 KQ1=1,3
      AQA(IQ1,JQ1,KQ1,IQ,JQ,KQ)=A(IQ1,JQ1,KQ1)
C      Write(6,*) 'cinti 501', AQA(IQ1,JQ1,KQ1,IQ,JQ,KQ)
  501 CONTINUE
C
  601 CONTINUE
      IF(MV.EQ.1.OR.MV.EQ.LC1) GO TO 270
      IF(MA.EQ.1.OR.MA.EQ.LC2) GO TO 270
      IF(ME.EQ.1.OR.ME.EQ.LC3) GO TO 270
      I = MV+1
      J = MA+1
      K = ME+1
      A(1,1,1) = F(I-2,J-2,K-2)
      VDX(1,1,1) =0.
      VDX(2,1,1)= X(I-1)-X(I-2)
      VDX(3,1,1) = X(I)-X(I-2)
      VDX(1,2,1) = Y(J-1)-Y(J-2)
      VDX(1,3,1) = Y(J)-Y(J-2)
      VDX(1,1,2) = Z(K-1) - Z(K-2)
      VDX(1,1,3) = Z(K) - Z(K-2)
      DO 701 I1=1,3
      DO 702 L1=1,2
      IA=L1*IK(I1)
      JA=L1*JK(I1)
      KA=L1*KK(I1)
      DO703 L2=1,2
  703 AC2(L1,L2)=(VDX(1+IA,1+JA,1+KA))**L2
  702 BC2(L1)=(F(I-2+IA,J-2+JA,K-2+KA)-F(I-2,J-2,K-2))
C 
C     Write(6,*) 'call DGELG 4' 
      CALL DGELG(BC2,AC2,2,1,1.0D-15,IER)
C
C      write(6,*) 'ier 4', IER
      DO 704 L1=1,2
      IB=L1*IK(I1)+1
      JB=L1*JK(I1)+1
      KB=L1*KK(I1)+1
      IG = IB+JB+KB-3
  704 A(IB,JB,KB)=BC2(IG)
  701 CONTINUE
      DO 705 I1=1,3
      L1=1
      DO706 L12=1,2
      DO706 L11=1,2
      JO=L11*IK(I1)
      KO=L11*JK(I1)
      IO=L11*KK(I1)
      IT=L12*IK(I1)
      JT=L12*JK(I1)
      KT=L12*KK(I1)
      L2=1
      DO 707 L22=1,2
      DO707 L21=1,2
      AC4V(L11,L12,L21,L22)=(VDX(1+IO,1+JO,1+KO))**L21*(VDX(1+IT,1+JT,
     11+KT))**L22
      AC4(L1,L2)=AC4V(L11,L12,L21,L22)
  707 L2=L2+1
      FV1=F(I-2,J-2,K-2)
      FV2(1)=F(I-2+IO,J-2+JO,K-2+KO)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2+IT,J-2+JT,K-2+KT)-F(I-2,J-2,K-2)
      BC4V(L11,L12)=F(I-2+IO+IT,J-2+JO+JT,K-2+KO+KT)-FV2(1)-FV2(2)-FV1
      BC4(L1)=BC4V(L11,L12)
  706 L1=L1+1
C
C     Write(6,*) 'call DGELG 5'
      CALL DGELG(BC4,AC4,4,1,1.0D-15,IER)
C
C     write(6,*) 'ier 5', IER
      LZ=1
      DO 708 L12=1,2
      DO 708 L11=1,2
      ID= L11*KK(I1)+L12*IK(I1)+1
      JD= L11*IK(I1)+L12*JK(I1)+1
      KD= L11*JK(I1)+L12*KK(I1)+1
      A(ID,JD,KD)=BC4(LZ)
  708 LZ = LZ+1
  705 CONTINUE
      L = 1
      DO 709 L1=1,2
      DO 709 L2=1,2
      DO 709 L3=1,2
      FV2(1)=F(I-2+L1,J-2,K-2)-F(I-2,J-2,K-2)
      FV2(2)=F(I-2,J-2+L2,K-2)-F(I-2,J-2,K-2)
      FV2(3)=F(I-2,J-2,K-2+L3)-F(I-2,J-2,K-2)
      FV3(1)=F(I-2+L1,J-2+L2,K-2)-FV2(1)-FV2(2)-FV1
      FV3(2)=F(I-2+L1,J-2,K-2+L3)-FV2(1)-FV2(3)-FV1
      FV3(3)=F(I-2,J-2+L2,K-2+L3)-FV2(2)-FV2(3)-FV1
      FV4(L1,L2,L3)=F(I-2+L1,J-2+L2,K-2+L3)-FV3(1)-FV3(2)-FV3(3)-FV2(1)-
     1FV2(2)-FV2(3)-FV1
      BC8(L)=FV4(L1,L2,L3)
  709 L=L+1
      L1=1
      DO710 L11=1,2
      DO710 L12=1,2
      DO710 L13=1,2
      L2=1
      DO711 L21=1,2
      DO711 L22=1,2
      DO711 L23=1,2
      AC8V(L11,L12,L13,L21,L22,L23)=(VDX(1+L11,1,1))**L21*(VDX(1,1+L12,1
     1))**L22*(VDX(1,1,1+L13))**L23
      AC8(L1,L2) = AC8V(L11,L12,L13,L21,L22,L23)
  711 L2=L2+1
  710 L1=L1+1
C
C     Write(6,*) 'call DGELG 6'
      CALL DGELG(BC8,AC8,8,1,1.0D-15,IER)
C
C      write(6,*) 'ier 6', IER
      LZ=1
      DO713 I=2,3
      DO713 J=2,3
      DO713 K=2,3
      A(I,J,K)=BC8(LZ)
  713 LZ=LZ+1
      DO 505 I=1,3
      DO 505 J=1,3
      DO 505 K=1,3
      AQM(I,J,K)=A(I,J,K)
c      Write(6,*) 'cinti 505', AQM(I,J,K)
  505 CONTINUE
C      write(6,*) 'return'
  270 RETURN
C
C************ EINGANG IN DIE SUBROUTINE CINT *******************
      ENTRY CINT(S,IPL,IERR)
  130 XS=S(1)
      YS=S(2)
      ZS=S(3)
      IF( IPL ) GO TO 1270
      IF(MV.EQ.1.OR.MV.EQ.LC1) GO TO 1270
      IF(MA.EQ.1.OR.MA.EQ.LC2) GO TO 1270
      IF(ME.EQ.1.OR.ME.EQ.LC3) GO TO 1270
      IF(XS.LT.X(MV-1).OR.YS.LT.Y(MA-1).OR.ZS.LT.Z(ME-1)) GO TO 1270
      IF(XS.GT.X(MV+1).OR.YS.GT.Y(MA+1).OR.ZS.GT.Z(ME+1))GO TO 1270
      I = MV+1
      J = MA+1
      K = ME+1
      XO=X(I-2)
      YO=Y(J-2)
      ZO=Z(K-2)
      DO 1501 I=1,3
      DO 1501 J=1,3
      DO 1501 K=1,3
 1501 A(I,J,K)=AQM(I,J,K)
      GO TO 193
 1270 IF(XS.LT.X(1).OR.YS.LT.Y(1).OR.ZS.LT.Z(1)) GO TO 170
      DO 140 I=3,LC1,2
      IF(XS.LE.X(I)) GO TO 150
  140 CONTINUE
      GO TO 170
  150 XO=X(I-2)
      DO 160 J=3,LC2,2
      IF(YS.LE.Y(J)) GO TO 151
  160 CONTINUE
      GO TO 170
  151 YO=Y(J-2)
      DO 180 K=3,LC3,2
      IF(ZS.LE.Z(K)) GO TO 200
  180 CONTINUE
  170 IERR = 1
      RETURN
C 
  200 ZO = Z(K-2)
      IQ=(I-1)/2
      JQ=(J-1)/2
      KQ=(K-1)/2
      DO 502 I=1,3
      DO 502 J=1,3
      DO 502 K=1,3
  502 A(I,J,K)=AQA(I,J,K,IQ,JQ,KQ)
  193 SU = 0.
      DO 14 I=2,3
      DO 14 J=2,3
      DO 14 K=2,3
   14 SU = SU + A(I,J,K)*(S(1)-XO)**(I-1)*(S(2)-YO)**(J-1)*(S(3)-ZO)**
     1(K-1)
      DO15 J=2,3
      DO15 K=2,3
   15 SU=SU+A(1,J,K)*(S(2)-YO)**(J-1)*(S(3)-ZO)**(K-1)
      DO16 I=2,3
      DO 16 K=2,3
   16 SU=SU+A(I,1,K)*(S(1)-XO)**(I-1)*(S(3)-ZO)**(K-1)
      DO17 I=2,3
      DO17 J=2,3
   17 SU=SU+A(I,J,1)*(S(1)-XO)**(I-1)*(S(2)-YO)**(J-1)
      DO 18 I=2,3
   18 SU=SU+A(I,1,1)*(S(1)-XO)**(I-1)
      DO 19 J=2,3
   19 SU=SU+A(1,J,1)*(S(2)-YO)**(J-1)
      DO 20 K=2,3
   20 SU=SU+A(1,1,K)*(S(3)-ZO)**(K-1)
C      write(6,*) ' coefficient A(1,1,1),SU in CINT', A(1,1,1), SU
      SU = SU+A(1,1,1)
      S(4)=SU
      IERR=0
      RETURN
      END
C
C
C########################################################################
C
       SUBROUTINE RALFA(FALF,ISHIF,HDM,ALFA,EXTRA,F,MENDIF,MAZDIF,
     1 MELDIF,VER,MAXALF,MZAHL,KZAHL)
C
C       DIESE SUBROUTINE TRENNT MIT HILFE EINER INTERPOLATION DIE ALFA-
C       ZAEHLRATEN VON DEN PROTONENZAEHLRATEN UND SPEICHERT SIE IM FELD
C       FALF(9,9,9)  ; ENTSPRECHEND  WERDEN DIE PROTONENZAEHLRATEN KORRI
C       GIERT UND NEU IN FIS1(7,7,32) GESPEICHERT ; DIE  GROESSEN  NEN,
C       NAZ,NEL GEBEN DEN ECKPUNKT DES ALFA-MESSBEREICHS AN UND LA1,LA2,
C       LA3 SEINE LAENGE ( COMMONBLOCK AA ) ; WENN EXTRA = .TRUE.,D.H.
C       WENN RANDEXTRAPOLATION , DANN WIRD F ANSTELLE VON FIS1 BENUTZT.
C       
       COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
       COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
       COMMON /KANAL/ VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
       COMMON/AA/NEN,NAZ,NEL,LA1,LA2,LA3
       COMMON/MAXIMUM/MAXINT,RATINT,MAX3D(3),RAT3D,MAXI1B,RATI1B
       DIMENSION PAZ(3),PEN(3),PEL(3),FALF(9,9,9),XF(6),YF(6),SIGMAY(3)
       DOUBLE PRECISION COEF(3),X(3),Y(3),F(33,9,9),FH,YY,DJ
       LOGICAL HDM,ALFA,EXTRA
       DIMENSION FMAX(5),JMIN(5),FMIN(5),JMAX(5),I2M(5),I3M(5),MAXALF(3)
       INTEGER  PAZ,PEN,PEL
       MN=MEN+LC1-1
       MZ=MAZ+LC2-1
       ML=MEL+LC3-1
       VER = 0.
       K=0
       NTERMS = 3
       IFEHL = 0
       IF ( KZAHL .LT. MZAHL ) GO TO 5555
       LMAX = MEN - 1 + LC1
       LMIN = MEN - 1 + LC1
C
C****** EINLESEN VON FIS1 IN F ****************************
C
 5555 IF ( EXTRA ) GO TO 190
      DO 191 I=1,LC1
      DO 191 J=1,LC2
      DO 191 K=1,LC3
  191 F(I,J,K) = FIS1(K,J,I)
C 
  190 CONTINUE   
C     Write(6,*) (((F(i,j,k),i=5,7),j=5,7),k=7,11)
C
C-----------------------------------------------------------------------
C----- SUCHE EXTREMA IN 3D UND INTEGRATIONSZAEHLERDATEN ---------------
C-----------------------------------------------------------------------
       FLA=0.
       J=MAXINT+9
       IF(J.GT.32) J=32
C       WRITE(6,3007) J
 3007  FORMAT(4H  J=I2)
       DO 10 I=1,6
       J=J-1
       IF(FINT(J+1).LT.15) GO TO 10
       FGL=FINT(J+1)*0.9
       JS=J+1
       JM=JS
       IF(FINT(J).LT.FGL) GO TO 11
       JM=JS+1
       IF(FINT(J).LT.FLA) GO TO 11
       FLA=FGL
10     CONTINUE
       WRITE(6,602)
602    FORMAT(' NO ALFA-MAXIMUM FOUND IN INTEGRATED DATA 1A ')
       IFEHL = 1
       KZAHL = KZAHL + 1
       GO TO 7000
11     CONTINUE
       DO 12 I=1,5
       J=JM-I
       IF ( J.LE. 1 ) GO TO 12
       FGL=FINT(J)
       KM = J
       IF(FINT(J-1).GT.FGL) GO TO 13
12    CONTINUE
      WRITE(6,603)
603   FORMAT(' NO ALFA-MINIMUM FOUND IN INTEGRATED DATA 1A')
       IFEHL = 1
       KZAHL = KZAHL + 1
      GO TO 7000
13    CONTINUE
      KZAHL = 0
C
C******EXAKTE BESTIMMUNG DER  EXTREMA AUS DEN INTEGRATIONSZAEHLERN ***
C     INSTRUMENT 1A
      I1 = JM - 1
      I2 = JM + 2
      IF(I2.GT.31) I2=31
C      WRITE(6,3008) I1,I2,KM
 3008 FORMAT(5H  I1=I3,3HI2=I3,3HKM=I3)
      DO 14 I=I1,I2
      FGL=FINT(I)
      LMAX=I
C      WRITE(6,5463) I
 5463 FORMAT(3H I=I2)
      IF(FINT(I+1).LT.FGL) GO TO 15
14    CONTINUE
      WRITE(6,604) I2,JM
604   FORMAT(' THIS SHOULD NOT HAPPEN  I2=',I3,'JM=',I3)
      ALFA = .FALSE.
      GO TO 900
15    DO 114 L = 1,4
      J = LMAX -L
C      WRITE(6,5464) J
 5464 FORMAT(3H J=I2)
      IF ( J .LE. 1 ) GO TO 114
      FGL = FINT(J)
      LMIN = J
      IF ( FINT(J-1) .GT. FGL ) GO TO 1144
  114 CONTINUE
C     WRITE(6,6043)
 6043 FORMAT( 'THIS SHOULD NOT HAPPEN' )
      ALFA = .FALSE.
      GO TO 900
 1144 CONTINUE
      WRITE(6,6044) LMAX,LMIN,FINT(LMAX),FINT(LMIN),MAXINT,RATINT
 6044 FORMAT(' Alpha separation'/6H LMAX=I2,6H LMIN=I2,5H MAX=F8.2,
     1 5H MIN=F8.2,8H MAXINT=I2,8H RATINT=F8.2)
      IF ( .NOT. HDM .AND. LMIN  .GE. (MEN+LC1) ) GO TO 6615
C
C     INSTRUMENT 1B
      DO 3014 I=I1,I2
      FGL=FI1B(I)
      LBMAX=I
C      WRITE(6,3011) I
 3011 FORMAT(3H I=I2)
      IF(FI1B(I+1).LT.FGL) GO TO 3015
 3014 CONTINUE
      WRITE(6,3604)I2,JM
 3604 FORMAT(' THIS SHOULD NOT HAPPEN  I2=',I3,'JM=',I3)
      IF ( IFEHL .EQ. 0 ) GO TO 3015
      ALFA = .FALSE.
      GO TO 900
C
 3015 DO 3114 L= 1,4
      J = LBMAX-L
C      WRITE(6,3064) J
 3064 FORMAT(3H J=I2)
      IF ( J .LE. 1 ) GO TO 3044
      FGL = FI1B(J)
      LBMIN = J
      IF ( FI1B(J-1) .GT. FGL ) GO TO 3044
 3114 CONTINUE
      WRITE(6,3043)
 3043 FORMAT( 'THIS SHOULD NOT HAPPEN' )
      IF ( IFEHL .EQ. 0 ) GO TO 3044
      ALFA = .FALSE.
      GO TO 900
 3044 CONTINUE
      WRITE(6,3444) LBMAX,LBMIN,FI1B(LBMAX),FI1B(LBMIN)
 3444 FORMAT(7H LBMAX=I2,6HLBMIN=I2,5H MAX=F8.2,5H MIN=F8.2)
      BMAXIN = FI1B(LBMAX)
      DO 3055 I=2,LBMAX
      K = LBMAX - I + 1
      IF ( FI1B(K) .LT. BMAXIN ) GO TO 3055
      BMAXIN = FI1B(K)
 3055 CONTINUE
      VER1A = 0.
      VER1B = 0.
      IF (FINT(LMAX) .GT. 0.) VER1A = RATINT/FINT(LMAX)
      VHI = FI1B(LBMAX) - 25.
      IF (VHI .GT. 0.) VER1B = (BMAXIN-25.)/VHI
      IF (VER1B .GT. 0.) VER = VER1A/VER1B
C      WRITE(6,3057) VER1A,VER1B,VER
 3057 FORMAT(7H VER1A=E10.4,6HVER1B=E10.4,6H  VER=E10.4)
C     IF ( VER .LT. 1.2  .AND. VER  .GT. 0.0 ) GO TO 62
 7000 CONTINUE
C
C***** SUCHE NACH ALFA-MAXIMUM UND MINIMUM IN RICHTUNG DES PROTONEN-
C***** MAXIMUMS; FALLS KEINES GEFUNDEN , SUCHE IN DEN NACHBARKANAELEN *
      FLA = 0.
      I1 = MAX3D(3) + MENDIF
      I2 = MAX3D(2) + MAZDIF
      I3 = MAX3D(1) + MELDIF
C      WRITE(6,6677) (  MAX3D(I),I ,I=1,3)
 6677 FORMAT(6I4)
      I2M(1) = I2
      I3M(1) = I3
      I2M(2) = I2 + 1
      I3M(2) = I3
      I2M(3) = I2
      I3M(3) = I3 + 1
      I2M(4) = I2
      I3M(4) = I3 - 1
      I2M(5) = I2 - 1
      I3M(5) = I3
      DO 708 KK = 1,5
      I2 = I2M(KK)
      I3 = I3M(KK)
C      WRITE(6,703) I1,I2,I3
  703 FORMAT(4H    3I4)
      J = I1 + 9
      IF ( J.GT. LC1) J=LC1
      DO 701 I = 1,6
      J = J - 1
      IF ( J .LT. I1 + 3 )  GO TO  709
      IF ( F(J+1,I2,I3)  .LE. 5.D0) GO TO 701
      FGL = F(J+1,I2,I3)*0.8
      LS = J + 1
      LM = LS
      IF ( F(J,I2,I3)  .LT. FGL*1.D0) GO TO 710
      LM = LS + 1
      IF ( F(J,I2,I3) .LT. FLA*1.D0) GO TO 710
      FLA = FGL
  701 CONTINUE
      FLA = 0.
  708 CONTINUE
  709 CONTINUE
      WRITE(6,6022)
 6022 FORMAT('NO ALFA MAXIMUM FOUND IN 3D DATA ')
C     IF ( VER .GT. 1.8 .AND. IFEHL .EQ. 0 ) GO TO 719
      IF ( IFEHL .EQ. 0) GO TO 719
      ALFA = .FALSE.
      GO TO 900
  710 DO 702 I  = 1,5
      J = LM-I
      IF ( J .LE. 1 ) GO TO 702
      FGL = F(J,I2,I3)
      IF ( F(J-1,I2,I3) .GT. FGL*1.D0) GO TO 720
  702 CONTINUE
      WRITE(6,6033)
 6033 FORMAT('  NO MINIMUM FOUND IN 3D DATA ')
      IF ( IFEHL .EQ. 0 ) GO TO 719
      ALFA = .FALSE.
      GO TO 900
  719 N3DMAX = LMAX - MEN + 1
      N3DMIN = LMIN - MEN + 1
      WRITE(6,7199)
 7199 FORMAT(' TRENNE ALFAS BEIM MINIMUM DES INTEGRATIONSZAEHLERS AB')
      GO TO 722
  720 CONTINUE
      N3DMAX = LM
      N3DMIN = J
      IF ( IFEHL .EQ. 1 .AND.(N3DMAX -N3DMIN) .GT. 2) GO TO 6613
  721 MAXDIF = N3DMAX - LMAX + MEN -1
      MINDIF = N3DMIN - LMIN + MEN -1
      WRITE(6,6055) N3DMAX,N3DMIN,MAXDIF,MINDIF
 6055 FORMAT(8H  3DMAX=I3,8H  3DMIN=I3,8H MAXDIF=I3,8H MINDIF=I3)
      IF ( MAXDIF .NE. 0 )  WRITE(6,6066)
      IF ( MINDIF .NE. 0 ) WRITE(6,6077)
 6066 FORMAT(' ALFA-MAXIMUM AUS INTEGRATIONS- UND 3D-DATEN STIMMT  NICHT
     1 UEBEREIN ')
 6077 FORMAT(' MINIMUM ZWISCHEN ALFA UND PROTONEN PEAK AUS INTEGRATIONS
     1 UND 3D-DATEN ERMITTELT STIMMT NICHT UEBEREIN ')
      IF ( ABS(1.*MINDIF) .GT. 2..AND. IFEHL .EQ. 1) GO TO 725
      IF ( MINDIF .GE. 1 .OR. MINDIF .LT. -1) GO TO 719
C*****CHECK IF SUFFICIENT 3D DATA EXISTS*********************
  722 IF(N3DMAX .LT. LC1 .AND. N3DMIN .GT. 0 ) GO TO 16
      WRITE(6,605) LC1,N3DMAX
605   FORMAT(' INSUFFICIENT 3D DATA  LC1,N3DMAX=',2I3)
  725 ALFA =.FALSE.
      GO TO 900
16    CONTINUE
C
C*****SEARCH 3D DATA FOR ALPHA MAXIMUM ********************************
      I1=N3DMAX
      FH= F(I1,I2,I3)
      N1=I1-1
      IF(N1.LT.1) N1=1
      N2=I1+1
      IF(N2.GT.LC1) N2=LC1
      L1=I2-1
      IF(L1.LT.1) L1=1
      L2=I2+1
      IF(L2.GT.LC2) L2=LC2
      M1=I3-1
      IF(M1.LT.1) M1=1
      M2=I3+1
      IF(M2.GT.LC3) M2=LC3
      K=0
17    DO 20 I1=N1,N2
      DO 20 I2=L1,L2
      DO 20 I3=M1,M2
      IF( F(I1,I2,I3)-FH )  20,21,22
22    FH= F(I1,I2,I3)
      K=0
21    K=K+1
      IF(K.GE.4) GO TO 23
C***** PEN , PAZ , PEL ENTHAELT DIE MAXIMUMS ADRESSE ******************
      PAZ(K)=I2
      PEL(K)=I3
      PEN(K)=I1
20    CONTINUE
      MAXALF(1) = I1
      MAXALF(2) = I2
      MAXALF(3) = I3
      GO TO 24
23    K=K-1
      WRITE(6,606) K
606   FORMAT(' THIS OUGHT NOT TO HAPPEN  K=',I3)
      ALFA = .FALSE.
      GO TO 900
C
24    CONTINUE
      IF(FH.GT. 5.D0) GO TO 25
      WRITE(6,615)  FH
615   FORMAT('  MAX 3D ALFA FLUX KLEINER ALS  5;FLUX=',D10.4)
      ALFA =.FALSE.
      GO TO 900
25    FM3D=FH
      FMZ=0.3*FM3D
C      WRITE(6,616) FH,PEN(1),PAZ(1),PEL(1)
  616 FORMAT(5X,'ALFA-MAXIMUM=',D10.2,4H EN=I2,4H AZ=I2,4H EL=I2)
C***** GEGENWAERTIG WIRD DEN ALFAS NUR DAS ENERGIE-INTERVALL LIMITIERT**
      L2 = 1
      M2 = LC2
      L3 = 1
      M3 = LC3
      L1 = N3DMIN
      M1 = L1 + 8
      IF ( M1 .GT. LC1 ) M1 = LC1
      WRITE(6,607) L1,M1,L2,M2,L3,M3
607   FORMAT(' LIMITS ARE (EN,AZ,EL)',6I5)
C***** STARTKANAL NUMMERN  DES ALFA-MESSAUSSCHNITTES BEZOGEN AUF DAS
C***** GROSSE RASTER ***************
      NEN=L1+MEN-1
      NAZ=L2+MAZ-1
      NEL=L3+MEL-1
      LA1=M1-L1+1
      LA2=M2-L2+1
      LA3=M3-L3+1
      IF ( LA1 .LT. 3 ) GO TO 613
C**** AUSSCHREIBEN DER ZAEHLRATEN *************************************
C          L5=L1-5
C          IF(L5.LT.1) L5=1
C          DO 94 I3=  L3,M3
C          DO 94 I2=  L2,M2
C  94      WRITE(6,620)( F(I1,I2,I3)  ,I1=L5,M1)
C
C      WRITE(6,600)      LA1,LA2,LA3,NEN,NAZ,NEL
C 600 FORMAT(3H   6I4)
C*****PRESET FALF ARRAY TO ZERO**********************************
      DO 4000 I1=1,9
       DO 4000 I2=1,9
        DO 4000 I3 =1,9
4000     FALF(I1,I2,I3)= 0.0
       ALFA = .TRUE.
C
C-----------------------------------------------------------------------
C*****JEDES ENERGIESPEKTRUM WIRD BEI FESTEM WINKELKANAL AUF EXTREMA
C***** UNTERSUCHT UND DIE COUNTS AUF FALF UND F VERTEILT .
C-----------------------------------------------------------------------
      DO 700 I3=L3,M3
      DO 700 I2=L2,M2
      IMAX = N3DMAX
      IMIN = N3DMIN
C***** ES WIRD EINE MINDESTZAEHLRATE VON 1 ANGENOMMEN *****************
      DJ=1.D0
C     CHECK IF DATA >DJ EXISTS
      LX= M1 - 1
      DO 40 I1=L1,LX
      IF( F(I1,I2,I3) .GE. DJ) GO TO 41
40    CONTINUE
C     NO DATA >1
      GO TO 700
41    CONTINUE
      IF ( I1 .GE. L1 + 1 )  GO TO  598
      CMIN = F(L1,I2,I3)
      IF ( F(L1+1,I2,I3) .NE. 0. ) GO TO 43
      IF ( F(L1-1,I2,I3) .GT. CMIN ) IMIN = IMIN + 1
      GO TO 598
   43 CONTINUE
      IF ( F(L1+2,I2,I3) .NE. 0. ) GO TO 44
      IF ( F(L1-1,I2,I3) .GE. CMIN ) GO TO 390
      GO TO 598
   44 CONTINUE
      IF ( F(L1-1,I2,I3) .GE. CMIN ) GO TO 45
      GO TO 390
   45 IF ( F(L1+1,I2,I3) .LT. CMIN )  GO TO 46
      GO TO 390
   46 IMIN = IMIN + 1
      IF ( F(L1+2,I2,I3) .LT. F(L1+1,I2,I3) ) IMIN = IMIN + 1
C
C*****ABZIEHEN DER EXTRAPOLIERTEN ZAEHLRATEN DER PROTONEN ************
  390 ISTRT = IMIN - 4
C      WRITE(6,9066) IMIN,ISTRT,ISHIF
 9066 FORMAT(6H IMIN=I2,7H ISTRT=I2,7H ISHIF=I2)
      DO 840  I = 1,3
      JJ = ISTRT + I + 1 + MEN - 1
      LL = ISTRT + I
      X(I) = VOLT( JJ,ISHIF)*1.D0
      YY = F(LL,I2,I3)
C      WRITE(6,9077) YY
 9077 FORMAT(4H YY=D12.3)
      IF ( YY .LE. 1.D0) YY=1.D0
C     Y(I) = DBLE(YY)
      Y(I) = DLOG(YY)
      COEF(I) = 0.D0
  840 SIGMAY(I) = 1./SQRT(SNGL(YY))
C     DO 841 I=1,3
C     WRITE(6,906) X(I),Y(I),SIGMAY(I),I
C 906 FORMAT(2D12.3,E12.3,I5)
C 841 CONTINUE
      NPTS = 3
      CALL POLFIT(X,Y,SIGMAY,NPTS,NTERMS,+1,COEF,CHISQR)
      NPTS = NPTS + 3
      DO 842 I= 4,NPTS
      JJ = ISTRT + 1 +I + MEN - 1
      XF(I) = VOLT(JJ,ISHIF)
      YY = (COEF(1) + COEF(2)*XF(I) + COEF(3)*(XF(I)**2))
      YF(I) = DEXP( YY  )
C      WRITE(6,8422) XF(I),YF(I),I,COEF(I-3)
 8422 FORMAT( 2E12.3 ,I5,D12.3)
  842 CONTINUE
C***** ERZEUGUNG DES ALFA-ZAEHLRATENFELDES ****************************
      IANF = IMIN
      IEND = IMIN + 2
      JJ = 3
C      WRITE(6,8455) IANF,IEND
 8455 FORMAT(' ianf,iend',2(I4))
      DO 843 I = IANF,IEND
      JJ = JJ + 1
      FF = F(I,I2,I3)    - YF(JJ)
      IF ( FF .LT. 0. ) FF = F(I,I2,I3)
      K = I -L1 + 1
      FALF(K,I2,I3) = FF
      F(I,I2,I3)  =  F(I,I2,I3)  -  FF*1.D0
C      WRITE(6,8433)   F(I,I2,I3) ,FALF(K,I2,I3),FF,I,K,I2,I3
 8433 FORMAT( 1D12.3,2E12.3 ,3X,4I2)
  843 CONTINUE
      GO TO 7043
  598 CONTINUE
      IEND = IMIN - 1
C***** NACH ABZUG DER FITZAEHLRATE DER PROTONEN VON DER ALFA-ZAEHLRATE
C***** WIRD FUER I GROESSER ALS IEND DIE PROTONENZAEHLRATE NULL GESETZT
 7043 IEND = IEND + 1
      DO 844 I = IEND , M1
      K = I - L1 + 1
      IF ( I .LT. L1 ) GO TO 844
      FALF(K,I2,I3) = F(I,I2,I3)
C      WRITE(6,8404) FALF(K,I2,I3),K,I2,I3
 8404 FORMAT(10X,5HFALF=F8.2,3I5)
  844 F(I,I2,I3) =0.D0
      IF ( M1 .EQ. LC1) GO TO 846
      MM1 = M1 +1
      DO 845 I=MM1,LC1
  845 F(I,I2,I3) = 0.D0
  846 CONTINUE
      GO TO 700
  613 WRITE(6,614)
  614 FORMAT('  BAD DATA : FELDLAENGE LA1 < 3 ')
      ALFA = .FALSE.
      GO TO 900
 6613 WRITE(6,6614)
 6614 FORMAT(' BAD DATA : MAX-MIN-DIFFERENZ IN 2D GROESSER ALS 2 UND
     1 IFEHL = 1')
      ALFA = .FALSE.
      GO TO 900
 6615 WRITE(6,6616)
 6616 FORMAT(' MINIMUM ZWISCHEN ALFAS UND PROTONEN AUS INSTRUMENT 1A
     1 LIEGT AUSSERHALB DES MESSARRAYS')
      GO TO 900
700   CONTINUE
C-----------------------------------------------------------------------
C****** SCHLEIFENENDE **************************************************
C-----------------------------------------------------------------------
C      FELDLAENGE IN DER ENERGIE SOLL UNGERADE SEIN
       IF ( LA1 .LT. 9 .AND. (LA1-2*(LA1/2)) .EQ. 0 ) LA1 = LA1 +1
       IF ( LA3 .EQ. 6 .OR. LA3 .EQ. 8 ) LA3 = LA3 + 1
      IF ( EXTRA ) GO TO 902
      DO 93 I=1,LC1
      DO 93 J=1,LC2
      DO 93 K=1,LC3
   93 FIS1(K,J,I)= F(I,J,K)
  902 CONTINUE
C           DO 92 I3=1,LA3
C           DO 92 I2=1,LA2
C   92      WRITE(6,620)(FALF(I1,I2,I3),I1=1,LA1)
C     GO TO 900
   62 WRITE(6,6226) VER
 6226 FORMAT(' VERHAELTNIS DER PEAKS VER=',E10.4)
C     ALFA = .FALSE.
C
900    CONTINUE
       WRITE(6,4563) ALFA
 4563  FORMAT(6H ALFA=L3)
C      WRITE(6,630) FINT,FH ,FMZ
C      WRITE(6,631) MAX3D
C      WRITE(6,631) MN,MZ,ML,MEN,MAZ,MEL,LC1,LC2,LC3
C      WRITE(6,631) M1,M2,M3,L1,L2,L3
C      WRITE(6,631) (PEN(I),PAZ(I),PEL(I),I=1,K),LMAX,JM
C      WRITE(6,631) LA1,LA2,LA3,NEN,NAZ,NEL
620    FORMAT(15F7.2)
621         FORMAT(10E10.2)
630    FORMAT(16F7.1/16F7.1/2F7.1)
631    FORMAT(10I5)
640          FORMAT(9(D10.3,I3))
641          FORMAT(3I5,E11.3, I5,E11.3)
690    FORMAT(10X,6I6)
      RETURN
      END
C
C***********************************************************************
C      SUBROUTINE POLFIT
C      LEAST-SQUARE FIT TO DATA WITH POLYNOMIAL CURVE
C      Y = A(1) +A(2)*X + A(3)*X***2 + .....
       SUBROUTINE POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)
       DIMENSION X(1),Y(1),SIGMAY(1),A(1)
       DIMENSION SUMX(19),SUMY(10),ARRAY(10,10)
       DOUBLE PRECISION X,Y,SUMX,SUMY,XTERM,YTERM,ARRAY,CHISQ
       DOUBLE PRECISION DETERM,DELTA,A
C
   11 NMAX = 2*NTERMS - 1
      DO 13 N=1,NMAX
   13   SUMX(N) = 0.
      DO 15 J=1,NTERMS
   15   SUMY(J) = 0.
        CHISQ = 0.
21     DO 50 I=1,NPTS
          XI = X(I)
          YI = Y(I)
31     IF ( MODE ) 32,37,39
  32  IF (YI) 35,37,33
  33  WEIGHT = 1./YI
      GO TO 41
  35  WEIGHT = 1./(-YI)
      GO TO 41
  37  WEIGHT = 1.
      GO TO 41
  39  WEIGHT = 1./ SIGMAY(I)**2
41    XTERM = WEIGHT
      DO 44 N=1,NMAX
      SUMX(N) = SUMX(N) + XTERM
  44  XTERM = XTERM*XI
  45  YTERM = WEIGHT*YI
      DO 48 N=1,NTERMS
       SUMY(N) = SUMY(N) + YTERM
  48  YTERM = YTERM*XI
  49  CHISQ = CHISQ + WEIGHT*YI**2
   50 CONTINUE
C
51    DO 54 J=1,NTERMS
       DO 54 K=1,NTERMS
      N = J +K -1
  54  ARRAY(J,K) = SUMX(N)
      DELTA = DETERM(ARRAY,NTERMS)
      IF ( DELTA ) 61,57,61
  57  CHISQR = 0.
      DO 59 J=1,NTERMS
  59  A(J) = 0.
      GO TO 80
  61  DO 70 L =1,NTERMS
       DO 66 J=1,NTERMS
        DO 65 K =1,NTERMS
      N = J+ K -1
  65  ARRAY(J,K) = SUMX(N)
  66  ARRAY(J,L) = SUMY(J)
  70  A(L) = DETERM ( ARRAY,NTERMS )/DELTA
C
      IF (NPTS .LE. NTERMS ) GO TO 90
   71 DO 75 J=1,NTERMS
       CHISQ = CHISQ - 2.*A(J)*SUMY(J)
       DO 75 K =1,NTERMS
        N = K + J -1
   75  CHISQ = CHISQ + A(J)*A(K)*SUMX(N)
   76 FREE = NPTS - NTERMS
   77 CHISQR = CHISQ  / FREE
80    RETURN
90    CHISQR = 1.E5
      RETURN
      END
C
C**************  FUNCTION DETERM  **************************************
      REAL FUNCTION DETERM*8  (ARRAY,NORDER)
      DOUBLE PRECISION ARRAY,SAVE
      DIMENSION ARRAY(10,10)
10    DETERM = 1.D0
11    DO 50 K = 1,NORDER
      IF ( ARRAY(K,K) )  41,21,41
21    DO 23 J=K,NORDER
      IF ( ARRAY(K,J) )  31,23,31
   23 CONTINUE
      DETERM = 0.D0
      GO TO 60
   31 DO 34 I = K,NORDER
       SAVE = ARRAY(I,J)
       ARRAY(I,J) = ARRAY(I,K)
   34 ARRAY(I,K) = SAVE
      DETERM = - DETERM
C
   41  DETERM = DETERM*ARRAY(K,K)
      IF ( K-NORDER)  43,50,50
   43 K1 = K  + 1
      DO 46 I=K1,NORDER
      DO 46 J=K1,NORDER
   46  ARRAY(I,J) = ARRAY(I,J) - ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
50    CONTINUE
60    RETURN
      END
C
C#######################################################################
C
      SUBROUTINE VERT(F,V,AZ,EL,MV,MA,ME,VXS,VYS,VZS,HDM,BAS,EXTRA,INST,
     1 MENDIF,MAZDIF,MELDIF,ALFA,TEST,ISHIF,Hhelios)
      COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
      COMMON/A7/MEN,MAZ,MEL,LC1,LC2,LC3
      COMMON/KANAL/VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
      COMMON/KANAL3/VOLT3(18),AZIMU3(18,2),ELEVA3(11)
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/MAXIMUM/MAXINT,RATINT,MAX3D(3),RAT3D,MAXI1B,RATI1B
      DOUBLE PRECISION V(33),AZ(9),EL(9),F(33,9,9),DF,DFN,
     1 XMAX,ZMAX,FMAX,DEE,ELZ(9),FZZ(9)
      LOGICAL HDM,EXTRA,ALFA,TEST
      DIMENSION CHACOR(11,2),CHACO3(11,2)
      INTEGER Hhelios
C
C******** BERECHNUNG DER PROTONENVERTEILUNGSFUNKTION *******************
C
C********EICH- UND KORREKTURFAKTOREN FUER INSTRUMENT 1A ***************
C
      DATA CHACOR/4.,4.626,4.956,4.687,4.527,4.393,4.588,4.731,4.591,4.9
     150,4.,4.,4.532,4.609,4.768,4.748,4.876,4.673,4.622,4.540,4.039,4./
C     DELEL = 0.08727
C     DELAZ = 0.09774
C
C********EICH- UND KORREKTURFAKTOREN FUER INSTRUMENT 3 ****************
      DATA CHACO3/1.,0.9298,0.9756,1.0718,1.0574,0.9613,1.0477,1.1077,
     10.94177,0.90586,1.,1.,0.9862,1.0039,1.0177,0.9933,0.9967,0.9972,
     21.0313,0.9737,1.,1./
      DELZ  = 0.01
      GMOA = 4.E-7
      GMOA3 = 1.628E-8
      CORR = 1./(GMOA*DELZ)
      CORR3 = 1./(GMOA3*DELZ)
C      WRITE(6,1250) CORR ,HDM
 1250 FORMAT(7H  CORR=E12.4,6H  HDM=L3)
      MKK = 1
      MJJ = 1
      MII = 1
      IF ( TEST ) GO TO 3
C     MAXIMUMSINDIZES
      MKK = MAX3D(1) + MELDIF
      MJJ = MAX3D(2) + MAZDIF
      MII = MAX3D(3) + MENDIF
C      WRITE(6,1255) MKK,MJJ,MII,Hhelios
 1255 FORMAT(' MAXINDIZES=',3I3,8H HELIOS=I3)
      IF ( .NOT. TEST ) GO TO 3
      IF ( INST .EQ. 3 ) GO TO 3
C     BEGRENZUNG DER FELDLAENGE IM HDM
      IF ( .NOT. HDM .OR. ALFA ) GO TO 3
      FINMAX = FINT(MAXINT)
      XFRAC = FINMAX*0.05
      INDEX = MAXINT + 5
      IF ( FINT(INDEX) .LT. XFRAC  .AND. FINT(INDEX+1) .LT. XFRAC )
     1 GO TO 1
      GO TO 3
    1 INDEX = MII + 5
      IF ( EXTRA) GO TO 111
      F3DMAX = FIS1(MKK,MJJ,MII)
      XFRAC = F3DMAX*0.05
      IF (FIS1(MKK,MJJ,INDEX).LT.XFRAC .AND. FIS1(MKK,MJJ,INDEX+1).LT.XF
     1RAC)
     1 GO TO 2
      GO TO 3
  111 F3DMAX = F(MII,MJJ,MKK)
      XFRAC = F3DMAX*0.1
      IF ( F(INDEX,MJJ,MKK) .LT. XFRAC .AND. F(INDEX+1,MJJ,MKK) .LT. XFR
     1AC)
     1 GO TO  2
      GO TO 3
    2 LC1 = INDEX - 1
      IF ( LC1 -2*(LC1/2) .NE. 0 ) LC1 = LC1 + 1
    3 CONTINUE
      IF ( TEST ) GO TO 335
      IF ( .NOT. HDM ) GO TO 335
      IF ( INST .EQ. 3) GO TO 335
      IF ( EXTRA ) GO TO 3155
      DO 500 I=1,8
      DO 500 J=1,LC2
      DO 500 K=1,LC3
      IF ( FIS1(K,J,I) .GT. 2.5) GO TO 335
      FIS1(K,J,I) = 0.
  500 CONTINUE
      GO TO 335
 3155 DO 501 I=1,8
      DO 501 J=1,LC2
      DO 501 K=1,LC3
      IF ( F(I,J,K) .GT. 2.5) GO TO 335
      F(I,J,K) = 0.
  501 CONTINUE
C
  335 LL = 0
      IF ( HDM .OR. LC1 .EQ.8 .OR. LC1 .EQ. 6 )  LL = 1
      IANF = MEN + 1
      IEND = IANF + LC1 - 1
      JANF = MAZ + 1
      JEND = JANF + LC2 - 1
      KANF = MEL + 1
      KEND = KANF + LC3 - 1
C      WRITE (6,677) IANF,IEND,JANF,JEND,KANF,KEND
C      Write(6,*) 'Proton array coordinates'
C      WRITE (6,677) MEN,MAZ,MEL,LC1,LC2,LC3
C  677 FORMAT( 6(I3) )
      IF ( INST .EQ. 3) GO TO 800
      DO 11 I=IANF,IEND
C      write(6,*)' 22@@ ', ISHIF, VOLT(I,ISHIF)
   11 V(I+LL+1-IANF) = VOLT(I,ISHIF)*1.D5
      DO 12 J =JANF,JEND
   12 AZ(J+1-JANF) = AZIMUT(J,ISHIF)*1.D0
      DO 13 K = KANF,KEND
   13 EL(K+1-KANF) =  ELEVAT(K)*1.D0
      IF ( EXTRA ) GO TO 1000
      GO TO 900
  800 CONTINUE
C      WRITE(6,8787) ( AZIMU3(J,ISHIF),J=1,9)
C      WRITE(6,8787) ( AZIMU3(J,ISHIF),J=10,18)
C 8787 FORMAT(9(2X,E10.4))
      DO 811 I=IANF,IEND
  811 V(I+LL+1-IANF) = VOLT3(I)*1.D5
      DO 812 J=JANF,JEND
  812 AZ(J+1-JANF) = AZIMU3(J,ISHIF)*1.D0
      DO 813 K=KANF,KEND
  813 EL(K+1-KANF) =  ELEVA3(K)*1.D0
      DO 16 I = 1,LC1
      DO 16 J = 1,LC2
      DO 16 K = 1,LC3
      DFN =   V(I+LL)**4*DCOS(EL(K)*UMRF)
   16 F(I+LL,J,K) = ( FIS1(K,J,I)+BAS)*1./CHACO3(K,Hhelios)*CORR3/DFN
C      WRITE(6,8788) ( AZ(J),J=1,9)
C 8788 FORMAT(9(2X,E10.4))
      GO TO 1002
  900 CONTINUE
      DO 14 I = 1,LC1
      DO 14 J = 1,LC2
      DO 14 K = 1,LC3
      DFN =   V(I+LL)**4*DCOS(EL(K)*UMRF)
   14 F(I+LL,J,K) = ( FIS1(K,J,I)+BAS)*4./CHACOR(K,Hhelios)*CORR/DFN
      IF ( .NOT. HDM ) GO TO 100
      GO TO 1001
 1000 DO 15 J=1,LC2
      DO 15 K=1,LC3
      DO 15 I=1,LC1
      MI = LC1 - I + 1
      DFN = V(MI+LL)**4*DCOS(EL(K)*UMRF)
   15 F(MI+LL,J,K) =(F(MI,J,K)+BAS)*4./CHACOR(K,Hhelios)*CORR/DFN
 1002 IF ( LC1 .EQ. 8 .OR. LC1 .EQ. 6) GO TO 1001
      IF ( .NOT. HDM ) GO TO 100
C
C***** RANDERGAENZUNG FUER LC3 = 6 (LC3=8) UND  FUER  DEN 33. ( 17. FUER
C     INSTRUMENT 3 ) KANAL***************************
 1001 V(1) = VOLT(IANF-1,ISHIF)*1.D5
      IF ( INST .EQ. 3) V(1) = VOLT3(IANF-1)*1.D5
      DO 21 J = 1,LC2
      DO 21 K = 1,LC3
      ZMAX = 0.D0
C      
      DO 211 I = 1,LC1
      IF ( F(I+1,J,K) .LT. ZMAX ) GO TO 211
      ZMAX = F(I+1,J,K)
      IMAX = I + 1
  211 CONTINUE
      IF ( ZMAX .EQ. 0.D0) GO TO 2111
      DEE = 1
      IF ( IMAX .EQ. 2 ) GO TO 2113
      DEE = ((V(1)-V(IMAX))/(V(2)-V(IMAX)))**2
 2113 F(1,J,K) = ZMAX*( (F(2,J,K)/ZMAX)**DEE )
C      WRITE(6,212) J,K,F(1,J,K),J,K,F(2,J,K),DEE
  212 FORMAT(5H F(1,I2,1H,I2,2H)=D12.4,5H F(2,I2,1H,I2,2H)=D12.4,5H DEE=
     1D12.4)
      GO TO 21
 2111 F(1,J,K) = 0.D0
   21 CONTINUE
      LC1 = LC1 + 1
      IF ( LC3 .EQ. 5 .OR. LC3  .EQ. 7 .OR. LC3 .EQ. 9  )  GO TO 100
      DO 31 I = 1,LC1
      DO 31 J = 1,LC2
      YMAX = 0.D0
      DO 311 K = 1,LC3
      IF ( F(I,J,K) .LT. YMAX ) GO TO 311
      YMAX = F(I,J,K)
      KK = K
  311 CONTINUE
      IF ( KK - LC3 ) 3001,3002,3002
 3001 IF ( YMAX .EQ. 0.D0 ) GO TO 3002
      F(I,J,LC3+1) = F(I,J,LC3)**2/F(I,J,LC3-1)
      GO TO 31
 3002 F(I,J,LC3+1) = YMAX
   31 CONTINUE
      EL(LC3+1) =ELEVAT(KEND+1)*1.D0
      IF ( INST .EQ. 3) EL(LC3+1) = ELEVA3(KEND+1)*1.D0
      LC3 = LC3 +1
C*********************************************************************
  100 CONTINUE
      FMAX = F(MII,MJJ,MKK)
C      WRITE(6,8790) FMAX
 8790 FORMAT(6H FMAX=E10.4)
C     UMRECHNUNG DER  WINKELZAEHLUNG  VON - NACH +
      IF ( Hhelios .EQ. 2 ) GO TO 47
      DO 40 I = 1,LC1
      DO 40 J = 1,LC2
      DO 41 K = 1,LC3
   41 FZZ(K)= F(I,J,LC3-K+1)
      DO 42 K = 1,LC3
   42 F(I,J,K) = FZZ(K)
   40 CONTINUE
      DO 50 K=1,LC3
   50 ELZ(K)= EL(LC3-K+1)
      DO 51 K=1,LC3
   51 EL(K) = ELZ(K)
      GO TO 46
   47 CONTINUE
      DO 43 I=1,LC1
      DO 43 K=1,LC3
      DO 44 J=1,LC2
   44 FZZ(J) = F(I,LC2-J+1,K)
      DO 45 J=1,LC2
   45 F(I,J,K) = FZZ(J)
   43 CONTINUE
      DO 52 J=1,LC2
   52 ELZ(J) = AZ(LC2-J+1)
      DO 53 J=1,LC2
   53 AZ(J) = ELZ(J)
C     MAXIMUM -BESTIMMUNG   DER VERTEILUNGSFUNKTION
   46 I1 = 1
      I2 = LC1
      IF ( HDM .AND. INST .NE. 3 ) I1=9
      DO 20 I = I1,I2
      DO 20 J = 1,LC2
      DO 20 K = 1,LC3
      IF ( F(I,J,K) .LT. FMAX ) GO TO 20
      FMAX = F(I,J,K)
      MV = I
      MA = J
      ME = K
   20 CONTINUE
C      WRITE(6,8789) ( AZ(J),J=1,9)
C 8789 FORMAT(9(2X,E10.4))
      WRITE(6,1875) FMAX
 1875 FORMAT(7H  FMAX=D12.4)
C*****BESTIMMUNG DER CARTESISCHEN GESCHWINDIGKEITSKOMPONENTEN DES
C     MAXIMUMS
      AZB = AZ(MA)*UMRF
      ELB = EL(ME)*UMRF
      COSA = COS(AZB)
      COSE = COS(ELB)
      SINA = SIN(AZB)
      SINE = SIN(ELB)
      VXS = V(MV)*COSA*COSE
      VYS = V(MV)*SINA*COSE
      VZS = V(MV)*SINE
      WRITE(6,9002) VXS,VYS,VZS,MV,MA,ME,LC1,LC2,LC3,MEN,MAZ,MEL
 9002 FORMAT(6H  VXS=E15.4,6H  VYS=E15.4,6H  VZS=E15.4/4H MV=I2,
     14H MA=I2,4H ME=I2,6H  LC1=I2,4HLC2=I2,4HLC3=I2/6H  MEN=I2,
     24HMAZ=I2,4HMEL=I2)
  888 RETURN
      END
C
C######################################################################
C
      SUBROUTINE  ZTEST(TN,VXT,VYT,VZT,TXT,TYT,TZT,MEN,MAZ,MEL,HDM,
     1LC1,LC2,LC3,ISHIF,Hhelios)
      COMMON /KONST/ PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/SEN1/FINT(32),FI1B(32),FIS1(7,7,32)
      COMMON/KANAL/VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
      DIMENSION CHACOR(11,2)
      LOGICAL HDM
      INTEGER Hhelios
C
C
C***** SUBROUTINE ZUR BESTIMMUNG DER TESTZAEHLRATEN
C
C********EICH- UND KORREKTURFAKTOREN FUER DIE VERTEILUNGSFUNKTION *****
      DATA CHACOR/4.,4.626,4.956,4.687,4.527,4.393,4.588,4.731,4.591,4.9
     150,4.,4.,4.532,4.609,4.768,4.748,4.876,4.673,4.622,4.540,4.039,4./
C     WRITE(6,123) TN,TXT,TYT,TZT,VXT,VYT,VZT,LC1,LC2,LC3,MEN,MAZ,MEL
  123 FORMAT(3H   4(E12.4)/3H   3(E10.4)/3H   6(I3))
      CORR = 2.5E+8
      XK=(1.67/(2.*1.38))*1.0E-08
C     WRITE(6,124) XK
  124 FORMAT(6H   XK=E10.4)
      DO 3 I=1,LC1
      XV=VOLT(MEN+I,ISHIF)*1.E5
C     WRITE(6,125) XV,I
  125 FORMAT(6H   XV=E12.4,3H I=I2)
      DO 3 J=1,LC2
      XAZ=AZIMUT(MAZ+J,ISHIF)
C     WRITE(6,126) XAZ,J
  126 FORMAT(6H  XAZ=E10.4,3H J=I2)
      DO 3 K=1,LC3
      XEL=ELEVAT(MEL+K)
C     WRITE(6,127) XEL,K
  127 FORMAT(6H  XEL=E10.4,3H K=I2)
      FIS1(K,J,I)=TN*(XK/PI)**1.5*SQRT(1./(TXT*TYT*TZT))*
     1EXP((-XK)*((XV*COS(UMRF*XEL)*COS(UMRF*XAZ)-
     2VXT)**2/TXT+(XV*COS(UMRF*XEL)*SIN(UMRF*XAZ)-
     3VYT)**2/TYT+(XV*SIN(UMRF*XEL)-VZT)**2/TZT))
      FIS1(K,J,I)=FIS1(K,J,I)*XV**4*COS(UMRF*XEL)/(CORR*4.)*CHACOR(K,
     1Hhelios)
C     WRITE(6,128) FIS1(K,J,I)
  128 FORMAT(2H  E12.4)
    3 CONTINUE
      RETURN
      END
C
C#########################################################################
C
      SUBROUTINE MOMENT(MV,MA,ME,RHO,BV,BAZ,BEL,T1,T2,T3,T1AZ,T1EL,T2AZ,
     1T2EL,T3AZ,T3EL,NPRO1,NPRO2,NPRO3,VX,VY,VZ,QX,QY,QZ,ALFA,
     2NALFA,R,PHIQ,EPSQ,MOFE)
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      COMMON/A3/ A(24),B(24),D(24),EG5,FG5,GG5
      COMMON/A2/ DVO,DVU,DAZO,DAZU,DELO,DELU
      COMMON/OUTPUT/OUT(17,8)
      DOUBLE PRECISION A,B,D,VC(24),EG5,FG5,GG5
      DOUBLE PRECISION DVO,DVU,DAZO,DAZU,DELO,DELU
      DIMENSION WERT(21),TXYZ(20),XMASS(2),R(9)
      LOGICAL ALFA
C     Definitions for calling the eigenvalue subroutine
      INTEGER IFAIL,INTGER(3)
      DOUBLE PRECISION AA(3,3),RR(3),RI(3),VR(3,3),VI(3,3),RXY
C
      DATA VC /24*0d0/
C
C
C***** BERECHNUNG DER MOMENTE  *****************************************
C
      BK = 1.381E-16
      PROM = 1.673E-24
      ALFM = 4.*PROM
      XMASS(1) = PROM
      XMASS(2) = ALFM
C
      KMASS = 1
      IF ( ALFA .AND.(NALFA.EQ.1)) KMASS=2
      FAKT = XMASS(KMASS)/BK
      MOFE = 0
C         
      CALL COEFS(DVU,DVO,NPRO1,VC)
      DO 401 I=1,16
  401 A(I)=VC(I)
      EG5=DVO-DVU
      CALL COEFS(DAZU,DAZO,NPRO2,VC)
      DO 402 J=1,16
  402 B(J)=VC(J)
      FG5=DAZO-DAZU
      CALL COEFS(DELU,DELO,NPRO3,VC)
      DO 403 K=1,16
  403 D(K)=VC(K)
      GG5=DELO-DELU
C         
      LL = 0
      N = 0
      M = 0
      L = 0
       DO 16 NL=1,21,2
       LL = NL - 1
       IF (NL - 5) 13,18,19
   18  LL = 0
       N = 1
       GO TO 13
   19  IF (NL - 7) 100,10,11
   10  N = 0
       M = 1
       LL = 0
       GO TO 13
   11  IF (NL - 9) 100,12,21
   12  LL = 0
       L = 1
       M = 0
       GO TO 13
   21  IF(NL - 11) 100,22,23
   22  LL = 0
       L = 0
       N = 2
       GO TO 13
   23  IF(NL - 13) 100,24,25
   24  LL = 0
       N = 0
       M = 2
       GO TO 13
   25  IF(NL - 15) 100,26,27
   26  LL = 0
       M = 0
       L = 2
       GO TO 13
   27  IF(NL -17) 100,28,29
   28  LL = 0
       L = 0
       N = 1
       M = 1
       GO TO 13
   29  IF(NL - 19) 100,30,31
   30  LL = 0
       M = 0
       L = 1
       N = 1
       GO TO 13
   31  IF(NL - 21) 100,32,100
   32  LL = 0
       M = 1
       N = 0
       L = 1
   13 CALL GAUSSQ(MV,MA,ME,NPRO1,NPRO2,NPRO3,LL,N,M,L,W)
      WERT(NL) = W
C      WRITE(6,1366) NL,WERT(NL)
 1366 FORMAT(4H NL=I3,4H  W=E10.4)
   16  CONTINUE
       IF ( WERT(1) .LT. 0.001) GO TO 107
       DO 20 JJJ=3,21,2
   20  WERT(JJJ) = WERT(JJJ)/WERT(1)
       WRITE (6,14) (WERT(NL1),NL1=1,15,2)
   140 FORMAT(6H RHO =,1PE10.3,7H V**2 =,1PE10.3/5H VX =,1PE10.3,
     1    5H VY = ,1PE10.3,5H VZ =,1PE10.3/8H VX**2 =,1PE10.3,
     2  8H VY**2 =,1PE10.3/8H VZ**2 =,1PE10.3)
       DO 40 I10=11,15,2
   40  TXYZ(I10) = (WERT(I10) - WERT(I10-6)**2)*FAKT
       WRITE (6,41) (TXYZ(I11),I11=11,15,2)
   41  FORMAT(5H TX =,1PE10.3,5H TY =,1PE10.3,5H TZ =,1PE10.3)
       AA(1,1) = WERT(11) - WERT(5)**2
       AA(2,2) = WERT(13) - WERT(7)**2
       AA(3,3) = WERT(15) - WERT(9)**2
       AA(1,2) = WERT(17) - WERT(5)*WERT(7)
       AA(1,3) = WERT(19) - WERT(5)*WERT(9)
       AA(2,3) = WERT(21) - WERT(7)*WERT(9)
       AA(3,2) = AA(2,3)
       AA(3,1) = AA(1,3)
C
       TXX = AA(1,1)*FAKT
       TYY = AA(2,2)*FAKT
       TZZ = AA(3,3)*FAKT
       TXY = AA(1,2)*FAKT
       TXZ = AA(1,3)*FAKT
       TYZ = AA(2,3)*FAKT
       TZX = TXZ
       TZY = TYZ
       TYX = TXY
C      BERECHNUNG DER WAERMESTROMDICHTE
       FAC = WERT(1)*XMASS(KMASS)/2.
       EKIN = WERT(11) + WERT(13) + WERT(15)
       LL = 2
       N = 1
       M = 0
       L = 0
      CALL GAUSSQ(MV,MA,ME,NPRO1,NPRO2,NPRO3,LL,N,M,L,W)
       W = W/WERT(1)
      QX = FAC*( W - EKIN*WERT(5) -
     & 2.*(AA(1,1)*WERT(5) + AA(2,1)*WERT(7) + AA(3,1)*WERT(9)))
       N = 0
       M = 1
      CALL GAUSSQ(MV,MA,ME,NPRO1,NPRO2,NPRO3,LL,N,M,L,W)
       W = W/WERT(1)
      QY = FAC*( W - EKIN*WERT(7) - 
     & 2.*(AA(1,2)*WERT(5) + AA(2,2)*WERT(7) + AA(3,2)*WERT(9)))
       M = 0
       L = 1
      CALL GAUSSQ(MV,MA,ME,NPRO1,NPRO2,NPRO3,LL,N,M,L,W)
       W = W/WERT(1)
      QZ = FAC*(W - EKIN*WERT(9) - 
     & 2.*(AA(1,3)*WERT(5) + AA(2,3)*WERT(7) + AA(3,3)*WERT(9)))
C
       PHIQ = 1./UMRF*ATAN2(QY,QX)
       EPSQ = 1./UMRF*ATAN2(QZ,SQRT(QX*QX+QY*QY))
       Q = SQRT( QX*QX + QY*QY + QZ*QZ )
C
      WRITE (6,42) AA(1,1),AA(1,2),AA(1,3),AA(2,2),AA(2,3),AA(3,3)
   42 FORMAT(15H MATRIX AA(I,J)/6(2XE10.3))
C
      IFAIL = 0
      CALL F02AGF(AA,3,3,RR,RI,VR,3,VI,3,INTGER,IFAIL)
      WRITE(6,*) ' eigenvalues', IFAIL 
C       
      WRITE (6,43)  RR(1),VR(1,1),VR(2,1),VR(3,1),RR(2),VR(1,2),VR(2,2),
     & VR(3,2), RR(3), VR(1,3),VR(2,3),VR(3,3)
   43 FORMAT(3(12H EIGENWERT =,E10.3,4X,14H EIGENVEKTOR =,3E10.3/))
C   
       VX = WERT(5)
       VY = WERT(7)
       VZ = WERT(9)
       RHO = WERT(1)
C
       BV = SQRT(VX**2 + VY**2 + VZ**2)/1.E5
       BXY = SQRT(VX**2 + VY**2)
       BAZ = 1./UMRF*ATAN2(VY,VX)
       BEL = 1./UMRF*ATAN2(VZ,BXY)
C
      T1  = RR(1)*FAKT
      T2  = RR(2)*FAKT
      T3  = RR(3)*FAKT
C
      T1AZ = 1./UMRF*DATAN2(VR(2,1),VR(1,1))
      RXY = DSQRT(VR(1,1)**2 + VR(2,1)**2)
      T1EL = 1./UMRF*DATAN2(VR(3,1),RXY)
C
      T2AZ = 1./UMRF*DATAN2(VR(2,2),VR(1,2))
      RXY = DSQRT(VR(1,2)**2 + VR(2,2)**2)
      T2EL = 1./UMRF*DATAN2(VR(3,2),RXY)
C
      T3AZ = 1./UMRF*DATAN2(VR(2,3),VR(1,3))
      RXY = DSQRT(VR(1,3)**2 + VR(2,3)**2)
      T3EL = 1./UMRF*DATAN2(VR(3,3),RXY)
C
      WRITE(6,9004) RHO,BV,BAZ,BEL,T1,T2,T3,T1AZ,T1EL,T2AZ,T2EL,T3AZ,T3EL
 9004 FORMAT(6H  RHO=E10.4,5H  BV=E10.4,6H  BAZ=E10.4,
     &6H  BEL=E10.4/6H   T1=E15.4,6H   T2=E15.4,6H   T3=E15.4/
     &8H   T1AZ=E15.4,8H   T1EL=E15.4/
     &8H   T2AZ=E15.4,8H   T2EL=E15.4/8H   T3AZ=E15.4,8H   T3EL=E15.4)
C              
      WRITE(6,9005) QX,QY,QZ,PHIQ,EPSQ,Q
 9005 FORMAT(6H   QX=E12.4,5H  QY=E12.4,5H  QZ=E12.4/6H PHIQ=E10.4,
     16H EPSQ=E10.4,4H  Q=E12.4)
C     Eigen vectors of temperature matrix
      R(1) = VR(1,1)
      R(2) = VR(2,1)
      R(3) = VR(3,1) 
C
      R(4) = VR(1,2)
      R(5) = VR(2,2)
      R(6) = VR(3,2) 
C
      R(7) = VR(1,3)
      R(8) = VR(2,3)
      R(9) = VR(3,3) 
C      
      IF ( ALFA .AND. NALFA .EQ. 1) GO TO 800
C     EINLESEN DER PROTONEN-MOMENTE    INS OUT-ARRAY
      OUT(1,3) = RHO
      OUT(2,3) = BV
      OUT(3,3) = BAZ
      OUT(4,3) = BEL
      OUT(11,3)= T1
      OUT(12,3)= T2
      OUT(13,3)= T3
      OUT(14,3)= T1AZ
      OUT(15,3)= T1EL
      OUT(16,3)= T2AZ
      OUT(17,3)= T2EL
      OUT(1,4) = T3AZ
      OUT(2,4) = T3EL
      OUT(3,4) = Q
      OUT(4,4) = PHIQ
      OUT(5,4) = EPSQ
      IF ( .NOT. ALFA ) GO TO 801
      OUT(14,6) = -VX
      OUT(15,6) = -VY
      OUT(16,6) = -VZ
  801 CONTINUE
      GO TO 100
  800 CONTINUE
C     EINLESEN DER MOMENTE DER ALFA-TEILCHEN
      OUT(1,5) = RHO
      OUT(2,5) = BV
      OUT(3,5) = BAZ
      OUT(4,5) = BEL
      OUT(11,5)= T1
      OUT(12,5)= T2
      OUT(13,5)= T3
      OUT(14,5)= T1AZ
      OUT(15,5)= T1EL
      OUT(16,5)= T2AZ
      OUT(17,5)= T2EL
      OUT(1,6) = T3AZ
      OUT(2,6) = T3EL
      OUT(3,6) = Q
      OUT(4,6) = PHIQ
      OUT(5,6) = EPSQ
      IF ( .NOT. ALFA ) GO TO 100
      OUT(14,6) = OUT(14,6) + VX
      OUT(15,6) = OUT(15,6) + VY
      OUT(16,6) = OUT(16,6) + VZ
      OUT(17,6) =(SQRT(OUT(14,6)**2 + OUT(15,6)**2 +OUT(16,6)**2))*1.E-5
  100 RETURN
  107 MOFE = 1
      RETURN
      END
C
C########################################################################
C
      SUBROUTINE ALVER(F,V,AZ,EL,FALF,MVALF,MAALF,MEALF,AVXS,AVYS,AVZS,
     1BAS,ISHIF,Hhelios)
      COMMON/AA/NEN,NAZ,NEL,LA1,LA2,LA3
      COMMON/KANAL/VOLT(34,2),AZIMUT(18,2),ELEVAT(11)
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      DOUBLE PRECISION V(33),AZ(9),EL(9),F(33,9,9),DFN,
     1 FMAX,ELZ(9),FZZ(9)
      DIMENSION CHACOR(11,2),FALF(9,9,9)
      INTEGER Hhelios
C
C******************** ALFAVERTEILUNGSFUNKTION **************************
C
C********EICH- UND KORREKTURFAKTOREN FUER DIE VERTEILUNGSFUNKTION *****
      DATA CHACOR/4.,4.626,4.956,4.687,4.527,4.393,4.588,4.731,4.591,4.9
     150,4.,4.,4.532,4.609,4.768,4.748,4.876,4.673,4.622,4.540,4.039,4./
C     DELEL = 0.08727
C     DELAZ = 0.09774
      DELZ = 0.01
      GMOA = 4.E-7
      CORR = 1./(GMOA*DELZ)
      WURZ = 1./SQRT(2.)
C
      IANF = NEN + 1
      IEND = IANF + LA1 - 1
      JANF = NAZ + 1
      JEND = JANF + LA2 - 1
      KANF = NEL + 1
      KEND = KANF + LA3 - 1
C
      write(6,*)  'alpha array dimensions'
      WRITE (6,677) IANF,IEND,JANF,JEND,KANF,KEND
 677  FORMAT( 6(I3) )
      write(6,*)  'length', LA1,LA2,LA3
C      write(6,*) Hhelios,ISHIF
C
       DO 11 I=IANF,IEND
   11  V(I+1-IANF) = VOLT(I,ISHIF)*1.D5*WURZ
       DO 12 J =JANF,JEND
   12  AZ(J+1-JANF) = AZIMUT(J,ISHIF)*1.D0
       DO 13 K = KANF,KEND
   13  EL(K+1-KANF) =  ELEVAT(K)*1.D0
        DO 14 I = 1,LA1
        DO 14 J = 1,LA2
        DO 14 K = 1,LA3
        DFN =   V(I)**4*DCOS(EL(K)*UMRF)
C        write(6,*) 'DFN=', DFN
        F(I,J,K) = (FALF(I,J,K) + BAS)*4.D0/CHACOR(K,Hhelios)*CORR/DFN
C        write(6,*) 'F=', F(I,J,K)
   14   CONTINUE
C*********************************************************************
  100 CONTINUE
      IF ( Hhelios .EQ. 2) GO TO 47
C     UMRECHNUNG DER  WINKELZAEHLUNG IN DER ELEVATION VON - NACH +
      DO 40 I = 1,LA1
      DO 40 J = 1,LA2
      DO 41 K = 1,LA3
   41 FZZ(K)= F(I,J,LA3-K+1)
      DO 42 K = 1,LA3
   42 F(I,J,K) = FZZ(K)
   40 CONTINUE
      DO 50 K=1,LA3
   50 ELZ(K)= EL(LA3-K+1)
      DO 51 K=1,LA3
   51 EL(K) = ELZ(K)
      GO TO 46
   47 CONTINUE
      DO 43 I=1,LA1
      DO 43 K=1,LA3
      DO 44 J=1,LA2
   44 FZZ(J) = F(I,LA2-J+1,K)
      DO 45 J=1,LA2
   45 F(I,J,K) = FZZ(J)
   43 CONTINUE
      DO 52 J=1,LA2
   52 ELZ(J) = AZ(LA2-J+1) 
      DO 53 J=1,LA2
   53 AZ(J) = ELZ(J)
C     MAXIMUM -BESTIMMUNG   DER VERTTEILUNGSFUNKTION
   46 FMAX = 0.D0
      DO 20 I = 1,LA1
      DO 20 J =1,LA2
      DO 20 K = 1,LA3
      IF ( F(I,J,K) .LT. FMAX ) GO TO 20
      FMAX = F(I,J,K)
      MVALF = I
      MAALF = J
      MEALF = K
   20 CONTINUE
      WRITE(6,1875) FMAX ,I,J,K
 1875 FORMAT(11H  ALFAFMAX=D12.4,3I5)
C*****BESTIMMUNG DER CARTESISCHEN GESCHWINDIGKEITSKOMPONENTEN DES
C     MAXIMUMS
      AZB = AZ(MAALF)*UMRF
      ELB = EL(MEALF)*UMRF
      COSA = COS(AZB)
      COSE = COS(ELB)
      SINA = SIN(AZB)
      SINE = SIN(ELB)
      AVXS = V(MVALF)*COSA*COSE
      AVYS = V(MVALF)*SINA*COSE
      AVZS = V(MVALF)*SINE
      write(6,*)'alpha particle maximum parameters'
      WRITE(6,9002) AVXS,AVYS,AVZS,MVALF,MAALF,MEALF,LA1,LA2,LA3
 9002 FORMAT(6H AVXS=E15.4,6H AVYS=E15.4,6H AVZS=E15.4/4H MV=I2,
     14H MA=I2,4H ME=I2,4HLA1=I2,4HLA2=I2,4HLA3=I2)
  888 RETURN
      END
C
C     MAXIMUMSSUCHROUTINE
      SUBROUTINE MAXSU(VXS,VYS,VZS,FMAX,VL,IPL)
      LOGICAL IPL
      COMMON/KONST/PI,UMRF,LH1,LH2,LH3,PIH
      DIMENSION Z(3)
      DOUBLE PRECISION S(4)
      NZAHL = 0
      FMAX = EXP(FMAX)
      VL1 = VL
      VXS1 = VXS
      VYS1 = VYS
      VZS1 = VZS
      VXS2 = VXS1
      VYS2 = VYS1
      VZS2 = VZS1
      FMAX1 = FMAX
      FMAX2 = FMAX1
  100 CONTINUE
      NZAHL = NZAHL + 1
C     WRITE(6,101) NZAHL
  101 FORMAT(7H NZAHL=I1)
      Z(1) = -VL1
      Z(2) = 0.
      Z(3) = +VL1
      DO 1  I=1,3
      DO 1  J=1,3
      DO 1  K=1,3
      VX = VXS1 + Z(I)
      VY = VYS1 + Z(J)
      VZ = VZS1 + Z(K)
      R = VX**2 + VY**2 + VZ**2
      R1 = SQRT(VX**2 + VY**2)
      S(1) = 1.D0*SQRT(R)
      S(2) = 1.D0*ATAN2(VY,VX)/UMRF
      S(3) = 1.D0*ATAN2(VZ,R1)/UMRF
      CALL CINT(S,IPL,IERR)
      FFMAX = 0.
      IF ( IERR .EQ. 0. ) FFMAX = SNGL(DEXP(S(4)))
C     IF ( IERR .EQ. 1 ) FFMAX= 0.
C     WRITE(6,102) VX,VY,VZ,FFMAX
  102 FORMAT(5X,3E10.3,6HFFMAX=E10.3)
      IF ( FFMAX .LE.  FMAX2 ) GO TO 1
      FMAX2 = FFMAX
      VXS2 = VX
      VYS2 = VY
      VZS2 = VZ
    1 CONTINUE
      IF ( FMAX2 .EQ. FMAX1 ) GO TO 120
      EPS = ABS( (FMAX2-FMAX1)/FMAX1 )
C     WENN DIE ABWEICHUNG DES NEUEN MAXIMUMS VOM ALTEN <0.2% SPRING ANS
C     ENDE
      IF ( EPS .LT. 0.02 ) GO TO 200
      GO TO 121
  120 IF ( FMAX2 .NE. FMAX ) GO TO 200
  121 VL1 = VL1/2.
      FMAX1 = FMAX2
      VXS1 = VXS2
      VYS1 = VYS2
      VZS1 = VZS2
      IF ( NZAHL .GT. 3 ) GO TO 200
      GO TO 100
  200 VXS = VXS2
      VYS = VYS2
      VZS = VZS2
      FMAX = ALOG(FMAX2)
  300 RETURN
      END
C
      SUBROUTINE KORORB(PARA,WEL,WAZ,AWEL,AWAZ,Hhelios,
     & ora4,ora8,ora9,ora15,ora17)
C***********************************************************************
C
C  KORREKTUR DER PARAMETER MIT O/A-DATEN
C
C[B             STAND VOM 19.11.76
C             UNTERSCHIED HELIOS1/2 FUER DIE WINKEL SOWIE KORREKTUR
C             WEGEN DES ENDLICHEN SONNENRADIUS EINGEBAUT.  RS, 7.3.77
C***********************************************************************
C
      DIMENSION PARA(3)
      INTEGER Hhelios
      REAL*4 ora4,ora8,ora9,ora15,ora17,
     & ORAL8,ORAL9,PAR1,WEL,WAZ,D
C
      IF (PARA(1).EQ.0.)   RETURN
      AWEL = 0.
      AWAZ = 0.
      KWAZ = 1
      PAR1 = PARA(1)
      IF (ABS(WAZ).LT.50.) GOTO 300
      WAZA = 0.
      KWAZ = 0
 300  WAZA = WAZ*3.14159265/180.
C    
      ORAL9 = ORA9*149590000./24./60./60.
      ORAL8 = ORA8*149590000./24./60./60.
C
      write(6,*) waza,oral8,oral9,par1,ora15
      IF(WAZ.EQ.0.) GO TO 400
C
      PARA(1) = SQRT((PAR1*SIN(WAZA) + ORAL9)**2
     *             + (PAR1*COS(WAZA) + ORAL8)**2)
      PARA(3) = PARA(3)/ORA17*60.
C
      WAZB   =  ATAN2((PAR1*SIN(WAZA) + ORAL9),(PAR1*COS(WAZA) + ORAL8))
C
      IF (Hhelios .EQ. 1) AWAZ = WAZB/3.14159265*180. - 0.25646/ORA4 
      IF (Hhelios .EQ. 2) AWAZ = WAZB/3.14159265*180. + 0.25646/ORA4
C
      IF ( KWAZ.EQ.0 ) AWAZ = 99.
      IF ( ABS(WEL) .GT. 50 ) GOTO 100
  350 IF (Hhelios .EQ. 1) BWEL = ( WEL + ORA15 - 90.)*3.141959265/180.
      IF (Hhelios .EQ. 2) BWEL = ( WEL - ORA15 + 90.)*3.141959265/180.
C
      write(6,*) BWEL, PAR1,PARA(1),ora15
C
      AWEL = ASIN( BWEL*PAR1/PARA(1) )/3.14159265*180.
      IF(PAR1 .LE. 0.) PARA(1) = 0.
      GOTO 200
 100  AWEL = 99.
 200  RETURN
C
C   KORR. WENN KEINE WINKEL BERECHNET ODER VORHANDEN
C
  400 D = ORAL9/PAR1
      IF( D .GT.  1.0) D =  1.0
      IF( D .LT. -1.0) D = -1.0
      WAZA = ASIN(D)
      PARA(1) = PAR1*COS(WAZA) + ORAL8
      PARA(3) = PARA(3)/ORA17*60.
      IF(WEL.NE.0.) GO TO 350
      RETURN
      END
C
