C                                                                       00010000
C                THIS PROGRAM CALCULATES                                00020000
C THE GAMMA-RAY ANGULAR DISTRIBUTION AND THE DOUBLY-DIFFERENTIAL        00030000
C CROSS SECTION FOR THE SCATTERING OF AN INELASTIC PARTICLE AND         00040000
C               AN EMISSION OF A GAMMA RAY                              00050000
C                                                                       00060002
C         >>>>>>> FOR PROJECTILE EXCITATION <<<<<<<                     00070002
C         >>>>>>> OR FOR TARGET  EXCITATION <<<<<<<                     00080002
C                                                                       00090002
C THE NUCLEAR LEVELS ARE NOT ONLY POPULATED BY COULOMB EXCITATION       00100000
C BUT ALSO BY CONVERSION AND GAMMA TRANSITIONS CASCADING DOWN           00110000
C HIGHER STATES                                                         00120000
C                                                                       00130000
C THE FORMULAS ARE TAKEN FROM A.WINTHER AND J.DE BOER,                  00140000
C "COULOMB EXCITATION",P. 310-314                                       00150000
C AND D.SCHWALM ET AL., NUCL.PHYSICS A192 (1972) 449-495                00160000
C                                                                       00170000
C THE PROGRAM IS COMPUTED BY                                            00180000
C                              HANS-JUERGEN WOLLERSHEIM                 00190000
C                              NUCLEAR PHYSICS GROUP 2                  00200000
C                     GESELLSCHAFT FUER SCHWERIONENFORSCHUNG MBH        00210000
C                              POSTFACH 110541                          00220000
C                              6100 DARMSTADT 11                        00230000
C                              WEST GERMANY                             00240000
C                                                                       00250000
C                              GSI-DARMSTADT 8.12.88                    00260002
C                                                                       00270000
C **********************************************************************00280000
C >>>>>>>>>>  ONLY E1-,E2-,E3- AND M1-TRANSITIONS ARE POSSIBLE <<<<<<<<<00290000
C **********************************************************************00300000
C             MAX. NUMBER OF STATES = 35                                00310000
C             MAX. NUMBER OF POSSIBLE TRANSITIONS = 199                 00320000
C **********************************************************************00330000
C                                                                       00340003
C                   COULEX - CALCULATION :                              00350000
C                   --------------------                                00360003
C SCATTERING ANGLE IN THE CM-SYSTEM : THETA (1) > THETA (2)             00370000
C                                                                       00380000
C THE NEW COULOMB EXCITATION PROGRAM FROM A.LELL HAS TO BE USED, THERE  00390000
C THE MATRIX-ELEMENTS FOR ELECTRIC MULTIPOLES ARE NOT MULTIPLIED WITH   00400000
C                    >>>> (I ** L)<<<<                                  00410000
C                 ("COULOMB EXCITATION" P.310)                          00420000
C IN THE FOLLOWING GAMMA-DECAY PROGRAM THE E2- AND E3-MATRIXELEMENTS    00430000
C ARE MULTIPLIED WITH (I ** 2).                                         00440002
C IF ONE USES ANOTHER COULOMB EXCITATION PROGRAM IN WHICH THE E2- AND   00450000
C E3-MATRIXELEMENTS ARE MULTIPLIED WITH (I ** L),ONE HAS TO CHANGE THE  00460000
C SIGN OF THESE MATRIXELEMENTS                                          00470000
C                                                                       00480000
C **********************************************************************00490000
C                                                                       00500000
C                        SUBROUTINES :                                  00510003
C                        -----------                                    00520003
C    ANGPAR , FK , FCT , RACAHI , CLEBI , FIT , SPLINE , DQTFG , DQG8   00530000
C SUBROUTINE FIT AND SPLINE ARE PUBLISHED BY                            00540000
C                                    HAGER+SELTZER NUCL.DATA A4,1-235   00550000
C SUBROUTINE DQTFG AND DQG8 ARE FROM THE MATH.LIB. OF THE IBM (GSI)     00560000
C                                                                       00570000
C **********************************************************************00580003
C                                                                       00590003
C                **** INPUT OF THE DATA ****                            00600000
C                     -----------------                                 00610003
C DATA CARD #1                                                          00620000
C ------------                                                          00630003
C I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I24,I25           00640002
C                                                       FORMAT(15(1XI1))00650002
C >>>> OUTPUT AND INPUT CONTROLS <<<<                                   00660000
C                                                                       00670000
C ALL OUTPUT CONTROLS ARE FIXED-POINT NUMBERS TAKING THE VALUES 1 OR 0  00680000
C ACCORDING TO WHETHER THE OUTPUT IN  QUESTION IS DESIRED OR NOT        00690000
C EXCEPTION : I22 = 0,1,2,3,4,5 ; I24 = 0,1,>1                          00700000
C                                                                       00710000
C I11=OUTPUT OF THE CONVERSION COEFFICIENTS (E2,M1,E1,E3)               00720000
C I12=POSSIBLE DECAYS OF A STATE                                        00730000
C I13=HN (LIFETIME OF THE STATE)       "COULOMB EXCITATION" P.311 (48)  00740000
C I14=GK(N,M)                          "COULOMB EXCITATION" P.311 (47)  00750000
C I15=FK(N,M)                          "COULOMB EXCITATION" P.311 (46)  00760000
C I16=ANG(K,KK)+FEEDING                "COULOMB EXCITATION" P.311 (45)  00770000
C I17=SPIN(N),SPIN(K),W(N,K)           "COULOMB EXCITATION" P.311 (44)  00780000
C I18=SPIN(N),SPIN(K),DS(N,K)          "COULOMB EXCITATION" P.314 (58)  00790000
C I19=GAMMA-RAY ANGULAR DISTRIBUTION TT=0,180,5 VGAMMA=0.0 AND PI=3.14  00800000
C I20=GAMMA-RAY ANGULAR DISTRIBUTION PLOTT-OUTPUT                       00810000
C I21=EXCITATION PROBABILITIES, CM-CROSS SECTIONS, AND                  00820000
C     ANGULAR DISTRIBUTION TENSORS                                      00830000
C I22=1 INPUT OF FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR  00840000
C I22=2 F.SOLID-ANGLE CORR. + INPUT OF THE DEORIENTATION CORRECTION     00850000
C I22=3 F.SOLID-ANGLE CORR. + DEORIENTATION CORRECTION (SCHERER-BLUME)  00860000
C I22=4 F.SOLID-ANGLE CORR. + DEORIENTATION CORRECTION (ABRAGAM-POUND)  00870000
C I22=5 F.SOLID-ANGLE CORR. + G2=0.96304-0.08579*ALOG(TAU)              00880000
C I23=INPUT OF M1-MATRIXELEMENTS AND M1-CONVERSION COEFFICIENTS         00890000
C I24=CALCULATION OF THE GAMMA-RAY ANGULAR DISTRIBUTION IN A COORDINATE 00900000
C     SYSTEM FOR WHICH THE EXCITED NUCLEUS IS AT REST                   00910000
C >>> I24=1 <<<                                                         00920000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  00930000
C     OF THE GAMMA - DETECTOR IN THE LABORATORY COORDINATE SYSTEM       00940000
C >>> I24>1 <<<                                                         00950000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  00960000
C     OF THE GAMMA - DETECTOR IN THE REST-COORDINATE SYSTEM             00970000
C I25=CALCULATION FOR PROJECTILE EXCITATION                             00980003
C                                                                       00990000
C DATA CARD #2 TO #8 INPUT FOR HAGER+SELTZER INTERPOLATION PROGRAM      01000000
C ------------------                                                    01010003
C CONVERSION COEFFICIENTS ARE CALCULATED FOR E1-,E2-,E3-,M1-TRANSITIONS 01020000
C                                                                       01030000
C NCC = NUMBER OF VALUES GIVEN IN EACH CASE (K,L,M) IN THE TABLE        01040000
C CCE1 = LOWEST TABULATED ENERGY OF SUBSHELL TO BE INTERPOLATED, EXCEPT 01050000
C        FOR L OR M SHELLS IN WHICH CASE E1 = -1.0                      01060000
C CCE2 = LOWEST TABULATED ENERGY OF THE K,L2,M5 SUBSHELLS FOR THE       01070000
C        INTERPOLATION IN THE K, L, M SHELLS,RESPECTIVELY               01080000
C CC MIN = MIN. GAMMA-ENERGY GIVEN IN THE TABLE OF THE                  01090000
C          CONVERSION COEFFICIENTS                                      01100000
C CC MAX = MAX. GAMMA-ENERGY GIVEN IN THE TABLE OF THE                  01110000
C          CONVERSION COEFFICIENTS                                      01120000
C                                                                       01130000
C EXAMPLE :                                                             01140000
C                                                                       01150000
C FOR K-SHELL: CCE1,CCE1,CCE1,CCKMAX (K-SHELL)      (CCE1=CCE2=CCKMIN)  01160000
C FOR L-SHELL: -1. ,CCE2 (L2),CCLMIN (TOTAL L-),CCLMAX (TOTAL L-SHELL)  01170000
C FOR M-SHELL: -1. ,CCE2 (M5),CCMMIN (TOTAL M-),CCMMAX (TOTAL M-SHELL)  01180000
C                                                                       01190000
C ALPHA = LIST OF CONVERSION COEFFICIENTS FROM THE TABLE                01200000
C                                                                       01210000
C        >>>> IF E1- AND E3-TRANSITIONS ARE POSSIBLE <<<<               01220000
C    LIST OF CONVERSION COEFFICIENTS IN THE FOLLOWING SEQUENCE :        01230000
C           A. ELECTRIC (1. E1-,2. E2-,3. E3-) AND                      01240000
C           B. MAGNETIC (4. M1-) TRANSITIONS                            01250000
C  INPUT OF CONVERSION COEFFICIENTS FOR M1-TRANSITIONS ONLY IN          01260000
C             CONNECTION WITH >>> I23 = 1 <<<                           01270000
C  USE ALWAYS DATA CARD #2 TO #8 FOR THE DIFFERENT TRANSITIONS          01280000
C                                                                       01290000
C DATA CARD #2                                                          01300000
C ------------                                                          01310003
C NCCK,NCCL,NCCM                                    FORMAT(3I5)         01320000
C DATA CARD #3                                                          01330000
C ------------                                                          01340003
C CCE1K,CCE2K,CCKMIN,CCKMAX   (MULTIPOLE)           FORMAT(4E10.3)      01350000
C DATA CARD #4                                                          01360000
C ------------                                                          01370003
C CCE1L,CCE2L,CCLMIN,CCLMAX   (MULTIPOLE)               "    "          01380000
C DATA CARD #5                                                          01390000
C ------------                                                          01400003
C CCE1M,CCE2M,CCMMIN,CCMMAX   (MULTIPOLE)               "    "          01410000
C DATA CARD #6 + FOLLOWING DATA CARDS                                   01420000
C ------------                                                          01430003
C (ALPHAK(I),I=1,NCCK)        (MULTIPOLE)           READ(5,*)           01440000
C DATA CARD #7 + FOLLOWING DATA CARDS                                   01450000
C ------------                                                          01460003
C (ALPHAL(I),I=1,NCCL)        (MULTIPOLE)               "               01470000
C DATA CARD #8 + FOLLOWING DATA CARDS                                   01480000
C ------------                                                          01490003
C (ALPHAM(I),I=1,NCCM)        (MULTIPOLE)               "               01500000
C                                                                       01510000
C DATA CARD #9 + FOLLOWING DATA CARDS                                   01520000
C ------------                                                          01530003
C IXYZ,JXYZ,MM1(IXYZ,JXYZ)                          FORMAT(2I5,F10.3)   01540000
C                                                                       01550000
C INPUT OF M1-MATRIXELEMENTS ONLY IN CONNECTION WITH >>> I23 = 1 <<<    01560000
C                    AND DATA CARD #9                                   01570000
C                                                                       01580000
C EACH CARD CONTAINS THE INFORMATION ABOUT THE NUMBER OF THE NUCLEAR    01590000
C STATE : IXYZ (INITIAL-), JXYZ (FINAL-STATE), AND THE M1-MATRIX -      01600000
C ELEMENT                                                               01610000
C THE NUMBER OF THE NUCLEAR STATE IS IDENTICAL WITH THE ONE IN THE      01620000
C COULEX-CALCULATION                                                    01630000
C IF ( IXYZ.EQ.0 ) END OF INPUT FOR M1-MATRIXELEMENTS                   01640000
C                                                                       01650000
C DATA CARD #10                                                         01660000
C -------------                                                         01670003
C TT,VGAMMA,VI1,VI2,K1LAB                           FORMAT(4F10.3,1I2)  01680000
C                                                                       01690000
C TT     =  POLAR ANGLE OF THE GAMMA - DETECTOR      (IN DEGREE)        01700000
C VGAMMA =  AZIMUTHAL ANGLE OF THE GAMMA - DETECTOR  (IN DEGREE)        01710000
C VI1 AND VI2 ARE THE AZIMUTHAL ANGLES OF THE PARTICLE DETECTOR         01720000
C VI1    =  LOWER BOUNDARY OF THE PARTICLE DETECTOR  (IN DEGREE)        01730025
C VI2    =  UPPER BOUNDARY OF THE PARTICLE DETECTOR  (IN DEGREE)        01740025
C K1LAB  =  SPECIFY THE NUMBER OF THE STATE,FOR WHICH THE TRANSFORMATION01750000
C           FROM THE CM-SYSTEM TO THE LAB-SYSTEM IS DONE                01760000
C IF >>> I24=0 OR 1 <<<                                                 01770000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  01780000
C     OF THE GAMMA - DETECTOR IN THE LABORATORY COORDINATE SYSTEM       01790000
C IF >>> I24>1 <<<                                                      01800000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  01810000
C     OF THE GAMMA - DETECTOR IN THE REST-COORDINATE SYSTEM             01820000
C                                                                       01830000
C                                                                       01840000
C DATA CARD #11                                                         01850000
C -------------                                                         01860003
C (Q(I),I=1,5,2)                                    FORMAT(3F10.3)      01870000
C                                                                       01880000
C INPUT OF THE FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR    01890000
C ONLY IN CONNECTION WITH >>> I22 = 1 <<< AND DATA CARD #11             01900000
C                                                                       01910000
C FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR                 01920000
C ASSUMPTION : Q0,Q2,Q4 ARE INDEPENDENT FROM THE GAMMA-ENERGY           01930000
C D.C.CAMP AND A.L.VAN LEHN NUCLEAR INSTRUMENTS AND METHODS 76 (1969),  01940000
C P.192-240                                                             01950000
C J.BARRETTE,G.LAMOUREUX,S.MONARO NUCLEAR INSTRUMENTS AND METHODS 93    01960000
C (1971), P.1-11                                                        01970000
C                                                                       01980000
C DATA CARD #11A                                                        01990000
C --------------                                                        02000003
C IXYZ,QGS(3,IXYZ),QGS(5,IXYZ)                     FORMAT(1I5,2F10.3)   02010000
C                                                                       02020000
C DATA CARD #11B                                                        02030000
C --------------                                                        02040003
C AJ,GZA,TAUC,H0(IN)                                                    02050000
C                                                                       02060000
C IF >>> I22 = 2 <<< DATA CARD #11 + DATA CARD #11A                     02070000
C INPUT OF DEORIENTATION CORRECTION FOR THE GAMMA-RAY ANGULAR           02080000
C                   DISTRIBUTION                                        02090000
C LARGE MAGNETIC HYPERFINE FIELDS ACT IN HIGHLY IONIZED AND EXCITED     02100000
C HEAVY IONS.THE ATTENUATION FACTORS G2 AND G4 CAN BE CALCULATED        02110000
C USING THE SCHERER-BLUME THEORY OR ABRAGAM-POUND THEORY                02120000
C                                                                       02130000
C EACH CARD CONTAINS THE INFORMATION ABOUT THE NUMBER OF THE NUCLEAR    02140000
C STATE : IXYZ , AND THE ATTENUATION - COEFFICIENTS :                   02150000
C            G2=QGS(3,IXYZ) G4=QGS(5,IXYZ)                              02160000
C                                                                       02170000
C THE NUMBER OF THE NUCLEAR STATE IS IDENTICAL WITH THE ONE IN THE      02180000
C COULEX-CALCULATION                                                    02190000
C IF ( IXYZ.EQ.0 ) END OF INPUT FOR ATTENUATION - COEFFICIENTS          02200000
C                                                                       02210000
C IF >>> I22 = 3 <<< DATA CARD #11 + DATA CARD #11B                     02220000
C DEORIENTATION CORRECTION FOR THE GAMMA-RAY ANGULAR DISTRIBUTION       02230000
C            USING THE SCHERER-BLUME THEORY                             02240000
C                                                                       02250000
C  AJ = ATOMIC SPIN                                                     02260000
C  GZA = G-FACTOR                                                       02270000
C  TAUC = REORIENTATION TIME IN VACUUM (PSEC)                           02280000
C  H0 = MAGNETIC HYPERFINE FIELD (MEGA GAUSS)                           02290000
C                                                                       02300000
C IF >>> I22 = 4 <<< DATA CARD #11 + DATA CARD #11B                     02310000
C DEORIENTATION CORRECTION FOR THE GAMMA-RAY ANGULAR DISTRIBUTION       02320000
C            USING THE ABRAGAM-POUND THEORY                             02330000
C                                                                       02340000
C  AJ = ATOMIC SPIN                                                     02350000
C  GZA = G-FACTOR                                                       02360000
C  TAUC = REORIENTATION TIME IN VACUUM (PSEC)                           02370000
C  H0 = MAGNETIC HYPERFINE FIELD (MEGA GAUSS)                           02380000
C                                                                       02390000
C DATA CARD #12                                                         02400000
C -------------                                                         02410003
C MZAHL,NORM                                        FORMAT(2I5)         02420000
C                                                                       02430000
C  MZAHL = NUMBER OF THETA - INTEGRATIONS FOR THE PROJECTILE ANGULAR    02440000
C          DISTRIBUTION                                                 02450000
C  NORM  = POS.VALUE: NUMBER OF THE GAMMA-TRANSITION ON WHICH ALL OTHER 02460003
C                     GAMMA-YIELDS ARE NORMALIZED                       02470003
C        = NEG.VALUE: NORMALIZED ON THE RUTHERFORD CROSS SECTION        02480003
C                                                                       02490003
C  EXAMPLE :                                                            02500000
C         NORM=0 FOR THE LAST POSSIBLE TRANSITION IN THE GAMMA CASCADE  02510000
C         NORM=1 FOR THE ONE BEFORE IN THE GAMMA CASCADE  E.T.C.        02520000
C                    NORM=0 FOR THE 2-->0 TRANSITION                    02530000
C  THE SEQUENCE OF THE GAMMA TRANSITIONS IS GIVEN FOR OUTPUT CONTROL    02540000
C  I17=1 OR I18=1                                                       02550000
C                                                                       02560000
C DATA CARD #13 + FOLLOWING DATA CARDS                                  02570000
C -------------                                                         02580003
C XA,XE                                             FORMAT(2F10.2)      02590000
C                                                                       02600000
C  XA,XE = SCATTERING ANGLE OF THE PROJECTILS IN THE CM - SYSTEM        02610000
C          INTEGRATION FROM XE TO XA (XA > XE)                          02620000
C                                                                       02630000
C  THE DIFFERENTIAL CROSS SECTION IS TRANSFORMED INTO THE LAB - SYSTEM  02640000
C  AND THE INTEGRATION IS DONE FOR THE SCATTERING ANGLE OF THE          02650000
C  PROJECTILE OR TARGET NUCLEUS IN THE LAB SYSTEM                       02660004
C  EXAMPLE:                                                             02670004
C  THETA(LAB-RECOIL NUCLEUS) = ( 180. - THETA(CM-PROJECTILE) ) / 2.     02680004
C                                                                       02690000
C **********************************************************************02700004
      REAL ME,MM1,MEE2                                                  02710000
      INTEGER OUAMP,OUPROW,OUXI,OUPSI                                   02720000
      INTEGER OUANG0,OUANG1,OUANG2,OUANG3                               02730000
      COMMON/C1/ E(35),SPIN(35),EE(35),MEE2(35,35),MM1(35,35),NMAX      02740000
      COMMON/C2/ ALPHA1(35,35),ALPHA2(35,35),ANG(35,5,5)                02750000
      COMMON/C12/ ALPHE1(35,35),ALPHE3(35,35),DELTA1(199)               02760000
      COMMON/C3/ I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I2402770000
     1,I25                                                              02780005
      COMMON/C4/ Z1,Z2,A1,A2,EP,THETAP,QG(5,35),QGS(5,35),IPAR(35)      02790000
      COMMON/C14/ AJ,GZA,TAUC,H0(35)                                    02800000
C **********************************************************************02810000
C   >>>> FOR THE HAGER AND SELTZER INTERPOLATION PROGRAM  <<<<          02820000
C **********************************************************************02830000
      COMMON/C5/ ALPHAK(23,4),ALPHAL(23,4),ALPHAM(23,4),NCCK,NCCL,NCCM  02840000
      COMMON/C6/ CCE1K(4),CCE2K(4),CCE1L(4),CCE2L(4),CCE1M(4),CCE2M(4)  02850000
      COMMON/C8/ ITHETA                                                 02860000
      COMMON/C9/  CCKMIN(4),CCLMIN(4),CCMMIN(4)                         02870000
      COMMON/C19/ CCKMAX(4),CCLMAX(4),CCMMAX(4)                         02880000
C **********************************************************************02890000
      COMMON/RUTH/SIGLAB                                                02900007
      COMMON/C10/ VGAMMA,Q(5),ME(35,35,4),IITRAN(199),IFTRAN(199)       02910000
      DIMENSION AA(35),W(35,35),AYY(199,199),BYY(199,199)               02920000
      DIMENSION DS(199),AXX(199),YE(199),XX(199)                        02930000
      DIMENSION YY(199),ZZ(199),AAXX(199),XXX(199),SINLAB(199)          02940000
      DIMENSION P(35),DSIG(35),CROSSN(199)                              02950000
      DIMENSION RUTHER(199),BRUTH(199)                                  02960000
C **********************************************************************02970000
      open (5,file='input')
      open (6,file='output')
      open (17,file='coulex')
C                                                                       02980006
C >>>> OUTPUT AND INPUT CONTROLS <<<<                                   02990006
C                                                                       03000006
      READ(5,52) I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I2403010006
     1,I25                                                              03020006
  52  FORMAT(15(1X,I1))                                                 03030006
C **********************************************************************03040006
C                                                                       03050005
C               >>>> INPUT FROM COULEX-PROGRAM <<<<                     03060000
C                                                                       03070000
C          THE FIRST PART OF THIS PROGRAM - INPUT OF DATA -             03080000
C              IS TAKEN FROM A DISK >>>> READ(17) <<<<                  03090000
C                                                                       03100005
C **********************************************************************03110000
C                                                                       03120000
      READ(17,8000) Z1,A1,Z2,A2,EP,QPAR                                 03130000
 8000 FORMAT(6F10.3)
C                                                                       03140000
C Z1 = CHARGE NUMBER OF THE PROJECTILE                                  03150000
C A1 = MASS OF THE PROJECTILE IN A.M.U.                                 03160000
C Z2 = CHARGE NUMBER OF THE TARGET NUCLEUS                              03170000
C A2 = MASS OF THE TARGET NUCLEUS IN A.M.U.                             03180000
C EP = LABORATORY ENERGY OF THE PROJECTILE IN MEV                       03190000
C QPAR = E1-POLARISATION                                                03200000
C                                                                       03210000
      IF(I25.EQ.0) GOTO 1000                                            03220005
      Z1P=Z2                                                            03230000
      A1P=A2                                                            03240000
      EP=EP*A2/A1                                                       03250000
      Z2=Z1                                                             03260000
      A2=A1                                                             03270000
      Z1=Z1P                                                            03280000
      A1=A1P                                                            03290000
C                                                                       03300005
 1000 READ(17,8001) NMAX                                                03310005
 8001 FORMAT(I10)
C                                                                       03320000
C NMAX = NUMBER OF NUCLEAR STATES                                       03330000
C                                                                       03340000
      DO 7001 IWOLLE=1,NMAX                                             03350000
      READ(17,8002) SPIN(IWOLLE),E(IWOLLE),IPAR(IWOLLE)                 03360000
 8002 FORMAT(2F10.3,I10)
 7001 CONTINUE                                                          03370000
C                                                                       03380000
C SPIN(N) = SPIN QUANTUM NUMBER OF THE N-TH NUCLEAR STATE               03390000
C E(N) = EXCITATION ENERGY OF THE N-TH NUCLEAR STATE IN MEV             03400000
C IPAR(N) = PARITY OF THE N-TH NUCLEAR STATE                            03410000
C        IPAR = 1 POSITIVE PARITY, IPAR=-1 NEGATIVE PARITY              03420000
C                                                                       03430000
      READ(17,8003) INTERV,NCM,EMMAX,XIMAX,ACCUR,NTIME                  03440000
 8003 FORMAT(2I10,3F10.5,I10)
C                                                                       03450000
C OUTPUT CONTROLS FOR THE COULEX CALCULATION                            03460000
C                                                                       03470000
      READ(17,8004) OUXI,OUPSI,OUAMP,OUPROW,OUANG0,OUANG1,OUANG2,OUANG3 03480000
 8004 FORMAT(8I10)
C                                                                       03490000
C OUTPUT CONTROLS FOR THE COULEX CALCULATION                            03500000
C                                                                       03510000
      DO 7003 IME14=1,4                                                 03520000
      DO 7002 IWOLLE=1,NMAX                                             03530000
      READ(17,8005) (ME(IWOLLE,IW,IME14),IW=1,NMAX)                     03540000
 8005 FORMAT(10F10.3)
C                                                                       03550006
C INPUT OF THE E1-,E2-,E3-, AND E4-MATRIX                               03560006
C                                                                       03570006
      IF(IME14.NE.4) GOTO 7002                                          03580000
      DO 7004 IW=1,NMAX                                                 03590000
      IF(ABS(ME(IWOLLE,IW,4)).GT.0.0) IME4=1                            03600000
 7004 CONTINUE                                                          03610000
 7002 CONTINUE                                                          03620000
 7003 CONTINUE                                                          03630000
C                                                                       03640000
      READ(17,8006) THETA                                               03650000
 8006 FORMAT(F10.3)
C                                                                       03660000
C THETA = SCATTERING ANGLE OF THE PROJECTILE (IN DEGREE) IN THE         03670000
C         CENTER-OF-MASS SYSTEM                                         03680000
C                                                                       03690000
      WRITE(6,430)                                                      03700000
  430 FORMAT(54H1DE BOER - WINTHER MULTIPLE COULOMB EXCITATION PROGRAM, 03710000
     12X,27HFOR LAMBDA  1 , 2 , 3 AND 4)                                03720000
      WRITE(6,431) Z1,A1,EP                                             03730000
  431 FORMAT(32H0PROJECTILE CHARGE NUMBER  Z1 = ,F6.2,14H,  MASS  A1 = ,03740000
     1F8.3,23HAMU,  LAB ENERGY  EP = ,F8.3,3HMEV)                       03750000
      WRITE(6,434) THETA                                                03760000
  434 FORMAT(52H0SCATTERING ANGLE IN CENTER OF MASS SYSTEM  THETA = ,F7.03770000
     12,8H DEGREES)                                                     03780000
      WRITE(6,437) Z2,A2                                                03790000
  437 FORMAT(28H0TARGET CHARGE NUMBER  Z2 = ,F6.2,14H,  MASS  A2 = ,F8.303800000
     1,3HAMU/)                                                          03810000
      WRITE(6,438)                                                      03820000
  438 FORMAT(11X,'NO.',5X,'  SPIN   ',3X,' ENERGY    PARITY ( POS.=1 /  03830000
     1NEG.=-1 )'/)                                                      03840000
      NE1E3=0                                                           03850000
C                                                                       03860000
C     NE1E3 = 1 IF NEGATIVE PARITY STATES ARE INCLUDED IN THE COULEX-   03870000
C     CALCULATION ===> OUTPUT OF E1- AND E3-MATRIX                      03880000
C                                                                       03890000
      DO 7175 I175=1,NMAX                                               03900000
      IF(IPAR(I175).LT.0) NE1E3=1                                       03910000
 7175 WRITE(6,7176) I175,SPIN(I175),E(I175),IPAR(I175)                  03920000
 7176 FORMAT(10X,I3,7X,F5.1,4X,F9.4,2X,I5)                              03930000
      WRITE(6,7177) QPAR                                                03940000
 7177 FORMAT(1X/1X,'E1-POLARISATION,QPAR = ',F8.4/)                     03950000
      WRITE(6,442)                                                      03960000
  442 FORMAT(84H0PERFORMANCE CONTROLS  NMAX,   INTERV,   NCM,   EMMAX,  03970000
     1 XIMAX,     ACCUR,     NTIME)                                     03980000
      WRITE(6,443)NMAX,INTERV,NCM,EMMAX,XIMAX,ACCUR,NTIME               03990000
  443 FORMAT(1H0,19X,I6,2X,2I8,2F9.2,F13.7,I9)                          04000000
      WRITE(6,444)                                                      04010000
  444 FORMAT(81H0OUTPUT CONTROLS       OUXI, OUPSI, OUAMP, OUPROW, OUANG04020000
     10, OUANG1, OUANG2, OUANG3)                                        04030000
      WRITE(6,445) OUXI,OUPSI,OUAMP,OUPROW,OUANG0,OUANG1,OUANG2,OUANG3  04040000
  445 FORMAT(1H0,17X,3I7,5I8)                                           04050000
      WRITE(6,43)                                                       04060000
   43 FORMAT(67H0MULTIPOLE MATRIX ELEMENTS IN BARNS**(LAM/2) , READ FROM04070000
     1 DATA CARDS)                                                      04080000
C                                                                       04090000
C THE SIGN OF THE MATRIXELEMENT FOR ELECTRIC MULTIPOLES HAS TO BE       04100000
C CHANGED, IF ONE USES THE NEW COULOMB EXCITATION PROGRAM FROM A.LELL   04110000
C THE MATRIX-ELEMENTS FOR ELECTRIC MULTIPOLES ARE NOT MULTIPLIED WITH   04120000
C (I ** L) IN THIS COULOMB EXCITATION PROGRAMM                          04130000
C ("COULOMB EXCITATION" P.310)                                          04140000
C THE CHANGE OF THE SIGN IS ONLY DONE FOR THE E2- AND E3-MATRIXELEMENTS 04150000
C                                                                       04160000
      DO 50 N=1,NMAX                                                    04170000
      DO 51 M=1,NMAX                                                    04180000
      MEE2(N,M)=-ME(N,M,2)                                              04190000
      ME(N,M,3)=-ME(N,M,3)                                              04200000
  51  CONTINUE                                                          04210000
  50  CONTINUE                                                          04220000
C                                                                       04230000
C **********************************************************************04240000
C            >>>> END OF INPUT DATA FROM COULEX <<<<                    04250000
C **********************************************************************04260000
C                                                                       04270000
C >>>> OUTPUT AND INPUT CONTROLS <<<<                                   04280000
C      -------------------------                                        04290009
      WRITE(6,53)                                                       04300000
  53  FORMAT(1H0,1X,'I11 I12 I13 I14 I15 I16 I17 I18 I19 I20 I21 I22 I2304310000
     1 I24 I25'/)                                                       04320006
      WRITE(6,54) I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I204330000
     14,I25                                                             04340006
  54  FORMAT(15(2X,I2)/)                                                04350006
C                                                                       04360000
C      >>>> INPUT FOR HAGER+SELTZER INTERPOLATION PROGRAM <<<<          04370009
C  CONVERSION COEFFICIENTS ARE CALCULATED FOR E1-,E2-,E3-TRANSITIONS    04380000
C                                                                       04390000
      IF(NE1E3.EQ.0) GOTO 500                                           04400000
C                                                                       04410000
      MULT=1                                                            04420000
      GOTO 501                                                          04430000
  500 MULT=2                                                            04440000
  501 READ(5,100) NCCK,NCCL,NCCM                                        04450000
  100 FORMAT(3I5)                                                       04460000
      READ(5,101) CCE1K(MULT),CCE2K(MULT),CCKMIN(MULT),CCKMAX(MULT)     04470000
      READ(5,101) CCE1L(MULT),CCE2L(MULT),CCLMIN(MULT),CCLMAX(MULT)     04480000
      READ(5,101) CCE1M(MULT),CCE2M(MULT),CCMMIN(MULT),CCMMAX(MULT)     04490000
  101 FORMAT(4E10.3)                                                    04500000
      READ(5,*) (ALPHAK(I,MULT),I=1,NCCK)                               04510000
      READ(5,*) (ALPHAL(I,MULT),I=1,NCCL)                               04520000
      READ(5,*) (ALPHAM(I,MULT),I=1,NCCM)                               04530000
      IF(NE1E3.EQ.0) GOTO 502                                           04540000
      IF(MULT.EQ.1) GOTO 503                                            04550000
      IF(MULT.EQ.2) GOTO 504                                            04560000
      GOTO 502                                                          04570000
  503 MULT=2                                                            04580000
      GOTO 501                                                          04590000
  504 MULT=3                                                            04600000
      GOTO 501                                                          04610000
  502 IF(I23.EQ.0.OR.MULT.EQ.4) GOTO 511                                04620000
      MULT=4                                                            04630000
      GOTO 501                                                          04640000
  511 CONTINUE                                                          04650000
C                                                                       04660000
C M1-MATRIXELEMENTS : M1(SPIN(INITIAL) --> SPIN(FINAL)) = 0.0           04670000
C                                                                       04680000
      DO 551 I=1,NMAX                                                   04690000
      DO 552 J=1,NMAX                                                   04700000
      MM1(I,J)=0.0                                                      04710000
 552  CONTINUE                                                          04720000
 551  CONTINUE                                                          04730000
C                                                                       04740000
      IF(I23.EQ.0) GOTO 5432                                            04750000
C                                                                       04760009
C      >>>> INPUT OF M1-MATRIXELEMENTS <<<<                             04770009
C IXYZ (INITIAL-),JXYZ (FINAL-STATE),XYZ (M1-MATRIX ELEMENT)            04780009
C IF ( IXYZ.EQ.0 ) END OF INPUT FOR M1-MATRIXELEMENTS                   04790009
C                                                                       04800009
  55  READ(5,56) IXYZ,JXYZ,XYZ                                          04810000
  56  FORMAT(2I5,F10.3)                                                 04820000
      IF(IXYZ.EQ.0) GOTO 5432                                           04830000
      MM1(IXYZ,JXYZ)=XYZ                                                04840000
      GOTO 55                                                           04850000
 5432 CONTINUE                                                          04860000
C                                                                       04870000
C OUTPUT >>> E1-MATRIX <<<                                              04880000
C                                                                       04890000
      IF(NE1E3.NE.1) GOTO 5434                                          04900000
      WRITE(6,5435)                                                     04910000
 5435 FORMAT(1H0,1X,'ME1(N,M)'/)                                        04920000
      DO 5436 I=1,NMAX                                                  04930000
 5436 WRITE(6,59) (ME(I,J,1),J=1,NMAX)                                  04940000
  59  FORMAT(1X,10F8.3)                                                 04950000
C                                                                       04960000
C OUTPUT >>> E2-MATRIX <<<                                              04970000
C                                                                       04980000
 5434 WRITE(6,57)                                                       04990000
  57  FORMAT(1H0,1X,'ME2(N,M)'/)                                        05000000
      DO 58 I=1,NMAX                                                    05010000
  58  WRITE(6,59) (MEE2(I,J),J=1,NMAX)                                  05020000
C                                                                       05030000
C OUTPUT >>> E3-MATRIX <<<                                              05040000
C                                                                       05050000
      IF(NE1E3.NE.1) GOTO 5437                                          05060000
      WRITE(6,5438)                                                     05070000
 5438 FORMAT(1H0,1X,'ME3(N,M)'/)                                        05080000
      DO 5439 I=1,NMAX                                                  05090000
 5439 WRITE(6,59) (ME(I,J,3),J=1,NMAX)                                  05100000
 5437 CONTINUE                                                          05110000
C                                                                       05120000
C OUTPUT >>> E4-MATRIX <<<                                              05130000
C                                                                       05140000
      IF(IME4.NE.1) GOTO 5443                                           05150000
      WRITE(6,5441)                                                     05160000
 5441 FORMAT(1H0,1X,'ME4(N,M)'/)                                        05170000
      DO 5442 I=1,NMAX                                                  05180000
 5442 WRITE(6,59) (ME(I,J,4),J=1,NMAX)                                  05190000
 5443 CONTINUE                                                          05200000
C                                                                       05210000
C OUTPUT >>> M1-MATRIX <<<                                              05220000
C                                                                       05230000
      IF(I23.EQ.0) GOTO 5433                                            05240000
      WRITE(6,60)                                                       05250000
  60  FORMAT(1H0,1X,'MM1(N,M)'/)                                        05260000
      DO 61 I=1,NMAX                                                    05270000
  61  WRITE(6,59) (MM1(I,J),J=1,NMAX)                                   05280000
 5433 CONTINUE                                                          05290000
C                                                                       05300000
C TOTAL M1-CONVERSION COEFFICIENT:ALPHA(SPIN(INITIAL)-->SPIN(FINAL))=0.005310000
C                                                                       05320000
      DO 661 I=1,NMAX                                                   05330000
      DO 662 J=1,NMAX                                                   05340000
      ALPHA1(I,J)=0.0                                                   05350000
 662  CONTINUE                                                          05360000
 661  CONTINUE                                                          05370000
      MMAX=0                                                            05380000
C                                                                       05390000
C THE COMPLETE ANGULAR DISTRIBUTION IS CALCULATED ONLY FOR              05400009
C THE FIRST SCATTERING ANGLE : MMAX = 1                                 05410009
C     -----                                                             05420009
      IF(I19.EQ.1) MMAX=1                                               05430000
C                                                                       05440000
      READ(5,67) TT,VGAMMA,VI1,VI2,K1LAB                                05450000
  67  FORMAT(4F10.3,2I2)                                                05460000
C                                                                       05470010
C TT     =  POLAR ANGLE OF THE GAMMA - DETECTOR      (IN DEGREE)        05480010
C VGAMMA =  AZIMUTHAL ANGLE OF THE GAMMA - DETECTOR  (IN DEGREE)        05490010
C VI1 AND VI2 ARE THE AZIMUTHAL ANGLES OF THE PARTICLE DETECTOR         05500010
C VI1    =  LOWER BOUNDARY OF THE PARTICLE DETECTOR  (IN DEGREE)        05510025
C VI2    =  UPPER BOUNDARY OF THE PARTICLE DETECTOR  (IN DEGREE)        05520025
C K1LAB  =  SPECIFY THE NUMBER OF THE STATE,FOR WHICH THE TRANSFORMATION05530010
C           FROM THE CM-SYSTEM TO THE LAB-SYSTEM IS DONE                05540010
C                                                                       05550010
      IF(I24.GT.1) WRITE(6,72)                                          05560000
  72  FORMAT(1X/1X,'ANGLE OF THE GE(LI)-DETECTOR IN THE REST-SYSTEM :') 05570000
      WRITE(6,68) TT,VGAMMA,VI1,VI2,NMAX,K1LAB                          05580000
  68  FORMAT(1X/1X,'THETA = ',F8.2,1X,'VGAMMA = ',F8.2//                05590000
     11X,'VI1 = ',F8.2,1X,'VI2 = ',F8.2//                               05600000
     21X,'NUMBER OF TENSORS + 1 (FOR EVERY SCATTERING ANGLE) = ',I2//   05610000
     31X,'TRANSFORMATION FROM THE CM- TO THE LAB-SYSTEM ',              05620000
     4'(INDEX OF THE STATE) = ',I2/)                                    05630000
C                                                                       05640000
      VGAMMA=VGAMMA*3.1415926/180.                                      05650000
C ----------------------------------------------------------------------05660012
C                                                                       05670000
C INPUT OF THE FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR    05680000
C                                                                       05690000
      IF(I22.EQ.0) GOTO 1807                                            05700000
C                                                                       05710000
      DO 123 IXYZ=1,35                                                  05720000
      DO 124 JXYZ=1,35                                                  05730000
      QG(1,IXYZ)=1.                                                     05740000
      QG(3,IXYZ)=1.                                                     05750000
      QG(5,IXYZ)=1.                                                     05760000
      QGS(1,IXYZ)=1.                                                    05770000
      QGS(3,IXYZ)=1.                                                    05780000
      QGS(5,IXYZ)=1.                                                    05790000
  124 CONTINUE                                                          05800000
  123 CONTINUE                                                          05810000
C                                                                       05820000
      READ(5,1808) (Q(IH),IH=1,5,2)                                     05830000
 1808 FORMAT(3F10.3)                                                    05840000
C                                                                       05850012
C FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR                 05860012
C ASSUMPTION : Q0,Q2,Q4 ARE INDEPENDENT FROM THE GAMMA-ENERGY           05870012
C                                                                       05880012
      WRITE(6,1809) (Q(IH),IH=1,5,2)                                    05890000
 1809 FORMAT(1H0,1X,'FINITE SOLIDE-ANGLE CORRECTION FOR THE GE(LI)-DETEC05900000
     1TOR'/1X,' Q0 = ',F10.4,' Q2 = ',F10.4,' Q4 = ',F10.4/)            05910000
C                                                                       05920000
C >>>> INPUT FOR THE DEORIENTATION CORRECTION <<<<                      05930011
C                                                                       05940000
      IF(I22.EQ.2) GOTO 121                                             05950000
      IF(I22.EQ.3.OR.I22.EQ.4) GOTO 126                                 05960000
C                                                                       05970000
      GOTO 1807                                                         05980000
  121 READ(5,122) IXYZ,XYZ1,XYZ2                                        05990000
  122 FORMAT(1I5,2F10.3)                                                06000000
C                                                                       06010012
C IXYZ (NUMBER OF NUCLEAR STATE),XYZ1 (G2 - ATTENUATION - COEFFICIENT), 06020012
C XYZ2 (G4 - ATTENUATION - COEFFICIENT)                                 06030012
C IF ( IXYZ.EQ.0 ) END OF INPUT FOR ATTENUATION - COEFFICIENTS          06040012
C                                                                       06050012
      IF(IXYZ.EQ.0) GOTO 1807                                           06060000
      QGS(1,IXYZ)=1.                                                    06070000
      QGS(3,IXYZ)=XYZ1                                                  06080000
      QGS(5,IXYZ)=XYZ2                                                  06090000
      WRITE(6,125) IXYZ,QGS(1,IXYZ),QGS(3,IXYZ),QGS(5,IXYZ)             06100000
  125 FORMAT(1X,'ATTENUATION COEFFICIENTS #',I2,' G0 = ',F6.3,' G2 = '  06110000
     1,F6.3,' G4 = ',F6.3)                                              06120000
      GOTO 121                                                          06130000
  126 DO 132 IN=2,NMAX                                                  06140000
      READ(5,127) AJ,GZA,TAUC,H0(IN)                                    06150000
  127 FORMAT(2F10.3,2E10.3)                                             06160000
C                                                                       06170011
C  AJ = ATOMIC SPIN                                                     06180011
C  GZA = G-FACTOR                                                       06190011
C  TAUC = REORIENTATION TIME IN VACUUM (PSEC)                           06200011
C  H0 = MAGNETIC HYPERFINE FIELD (MEGA GAUSS)                           06210011
C                                                                       06220011
      ALAMD2=2.*((AJ+1.)/AJ)*((4789.44*GZA*H0(IN))**2)*TAUC*1.E-12      06230000
      IF(I22.EQ.3) WRITE(6,128) AJ,GZA,TAUC,H0(IN)                      06240000
  128 FORMAT(1X,' CALCULATION OF THE ATTENUATION COEFFICIENTS > GK < U  06250000
     1SING THE SCHERER-BLUME THEORY'//1X,'ATOMIC SPIN = ',F3.1,2X,'G-FAC06260000
     2TOR = ',F5.2,2X,'REORIENTATION TIME = ',E10.3,' (PSEC) ',2X,'MAGN.06270000
     3FIELD = ',E10.3,' (MEGA GAUSS)'/)                                 06280000
      IF(I22.EQ.4) WRITE(6,130) AJ,GZA,TAUC,H0(IN),ALAMD2               06290000
  130 FORMAT(1X,' CALCULATION OF THE ATTENUATION COEFFICIENTS > GK < U  06300000
     1SING THE ABRAGAM-POUND THEORY'//1X,'ATOMIC SPIN = ',F3.1,2X,'G-FAC06310000
     2TOR = ',F5.2,2X,'REORIENTATION TIME = ',E10.3,' (PSEC) ',2X,'MAGN.06320000
     3FIELD = ',E10.3,' (MEGA GAUSS)'/2X,' LAMDA2 = ',E10.3,            06330000
     4' (1/PSEC)'/)                                                     06340000
      H0(IN)=H0(IN)*1.E+6                                               06350000
      TAUC=TAUC*1.E-12                                                  06360000
  132 CONTINUE                                                          06370000
 1807 CONTINUE                                                          06380000
C ----------------------------------------------------------------------06390012
C                                                                       06400000
C THETAP = THETA (SCATTERING ANGLE IN THE CM-SYSTEM)                    06410000
C                                                                       06420000
      THETAP=THETA                                                      06430000
      ITHETA=0                                                          06440000
      I=0                                                               06450000
  80  ITHETA=ITHETA+1                                                   06460000
      I=I+1                                                             06470000
      DO 83 J=2,NMAX                                                    06480000
      DO 84 J1=1,5,2                                                    06490000
      DO 85 J2=1,5                                                      06500000
      ANG(J,J1,J2)=0.0                                                  06510000
  85  CONTINUE                                                          06520000
  84  CONTINUE                                                          06530000
  83  CONTINUE                                                          06540000
C                                                                       06550000
C INPUT OF THE ANGULAR DISTRIBUTION TENSORS FROM THE DISK >>READ(17)<<  06560000
C                                                                       06570000
      IF(I21.EQ.0) GOTO 1805                                            06580000
      WRITE(6,1010)                                                     06590000
 1010 FORMAT(2X,'ANG(0,0)   ANG(2,2)   ANG(2,1)   ANG(2,0)   ANG(4,4) ' 06600000
     1,'  ANG(4,3)   ANG(4,2)   ANG(4,1)   ANG(4,0)'/)                  06610000
 1805 DO 6 J=2,NMAX                                                     06620000
      READ(17,8007) ((ANG(J,J1,J2),J2=1,J1),J1=1,5,2)                   06630000
 8007 FORMAT(1X,9E12.4)
C                                                                       06640000
      P(J)=ANG(J,1,1)                                                   06650000
C                                                                       06660000
      IF(I21.EQ.0) GOTO 1806                                            06670000
      WRITE(6,1030) ANG(J,1,1),ANG(J,3,3),ANG(J,3,2),ANG(J,3,1),        06680000
     1ANG(J,5,5),ANG(J,5,4),ANG(J,5,3),ANG(J,5,2),ANG(J,5,1)            06690000
 1030 FORMAT(9(1X,E10.4))                                               06700000
 1806 CONTINUE                                                          06710000
   6  CONTINUE                                                          06720000
C                                                                       06730000
C THETAP = SCATTERING ANGLE OF THE PROJECTILES IN THE CM - SYSTEM       06740000
C          FOR WHICH THE ANGULAR DISTRIBUTION TENSORS ARE CALCULATED    06750000
C                                                                       06760000
      WRITE(6,1301) THETAP                                              06770000
 1301 FORMAT(1H0,'THETAP = ',F10.2/)                                    06780000
C                                                                       06790000
      IF(I21.EQ.0) GOTO 1800                                            06800000
      PSUM=0.0                                                          06810000
      DO 1801 J=2,NMAX                                                  06820000
      PSUM=PSUM+P(J)                                                    06830000
C                                                                       06840000
C CALCULATION OF THE RUTHERFORD CROSS SECTION  (IN BARNS/STERADIAN)     06850000
C                                                                       06860000
      VF1=((1.+A1/A2)*Z1*Z2/EP)**2                                      06870000
      VF2=SQRT(EP/(EP-(1.+A1/A2)*E(J)))                                 06880000
      DSR=1.2958E-3*VF1*VF2*(1./SIN(THETAP*3.1415926/(180.*2.))**4)     06890000
C                                                                       06900000
 1801 DSIG(J)=P(J)*DSR                                                  06910000
      P(1)=1.-PSUM                                                      06920000
      VF1=((1.+A1/A2)*Z1*Z2/EP)**2                                      06930000
      VF2=SQRT(EP/(EP-(1.+A1/A2)*E(1)))                                 06940000
      DSR=1.2958E-3*VF1*VF2*(1./SIN(THETAP*3.1415926/(180.*2.))**4)     06950000
C                                                                       06960000
C CALCULATION OF THE CM - CROSS SECTION (IN BARNS/STERADIAN)            06970000
C                                                                       06980000
      DSIG(1)=P(1)*DSR                                                  06990000
      WRITE(6,1802)                                                     07000000
 1802 FORMAT(1H0,1X,'LEVEL INDEX',6X,'EXCITATION',6X,'CM CROSS SECTION'/07010000
     117X,'PROBABILITIES',5X,'BARNS/STERADIAN'//                        07020000
     26X,'N',15X,'P(N)',13X,'DSIG(N)')                                  07030000
      DO 1803 J=1,NMAX                                                  07040000
 1803 WRITE(6,1804) J,P(J),DSIG(J)                                      07050000
 1804 FORMAT(5X,I2,11X,E11.4,7X,E11.4)                                  07060000
 1800 CONTINUE                                                          07070000
C                                                                       07080000
C SUBROUTINE FOR CALCULATION OF THE DOUBLE DIFFERENTIAL CROSS SECTION   07090000
C INCLUDING THE INTEGRATION OVER THE FINAL PARTICLE DETECTOR            07100000
C (AZIMUTHAL-ANGLE)                                                     07110000
C                                                                       07120000
      CALL ANGPAR(TT,VI1,VI2,DS,LAUF,SINLAB)                            07130000
C                                                                       07140000
      LAUF=LAUF-1                                                       07150000
      THET=THETAP*3.1415926/180.                                        07160000
      DO 90 K1=1,LAUF                                                   07170000
C                                                                       07180006
C AXX = THETA (CM. ANGLE OF PROJECTILE)                                 07190006
C                                                                       07200006
      AXX(I)=THETAP                                                     07210000
C                                                                       07220006
      IF(I25.EQ.0) GOTO 1001                                            07230006
      ETA=((180.-THETAP)/2.)*3.1415926/180.                             07240006
      AAXX(I)=ATAN(SIN(2.*ETA)/((A1/A2)-COS(2.*ETA)))                   07250006
      if(aaxx(i).lt.0) aaxx(i)=3.1415926+aaxx(i)
C                                                                       07260006
C ETA = ETA  (LAB. ANGLE OF TARGET NUCLEUS)                             07270006
C AAXX = PSI (LAB. ANGLE OF PROJECTILE)                                 07280006
C                                                                       07290006
      GOTO 90                                                           07300006
 1001 AAXX(I)=((180.-THETAP)/2.)*3.1415926/180.                         07310006
C                                                                       07320000
C AAXX = ETA (LAB. ANGLE OF TARGET NUCLEUS)                             07330006
C                                                                       07340006
  90  AYY(I,K1)=DS(K1)*SINLAB(K1LAB)                                    07350000
C-----------------------------------------------------------------------07360012
C                                                                       07370008
C RUTHERFORD CROSS SECTION IN LAB.SYSTEM * SIN(THETA-LAB)               07380008
C                                                                       07390008
      VFR=((1.+A1/A2)*Z1*Z2/EP)**2                                      07400006
      DSR=1.2958E-3*VFR*(1./SIN(THETAP*3.1415926/(180.*2.))**4)         07410006
      RUTHER(I)=DSR*SIGLAB*SINLAB(K1LAB)                                07420006
C                                                                       07430012
C-----------------------------------------------------------------------07440012
C                                                                       07450008
      IF(MMAX.EQ.1) GOTO 82                                             07460000
C                                                                       07470000
C INPUT OF THE NEXT SCATTERING ANGLE IN THE CM-SYSTEM                   07480008
C              ----                                                     07490008
      READ(17,8006) THETAP                                              07500000
      IF(THETAP.LT.0.0001) GOTO 82                                      07510000
      GOTO 80                                                           07520000
C **********************************************************************07530008
  82  CONTINUE                                                          07540000
      MMAX=ITHETA                                                       07550000
      WRITE(6,81) MMAX                                                  07560000
  81  FORMAT(1H0,1X,'MAX. NUMBER OF SCATTERING ANGLES -THETA- =',I3/)   07570000
C                                                                       07580000
C OUTPUT >>> E2-CONVERSION COEFFICIENTS <<<                             07590000
C                                                                       07600000
      IF(I11.EQ.0) GOTO 103                                             07610000
      WRITE(6,62)                                                       07620000
  62  FORMAT(1H0,1X,'CONVERSION-COEFFICIENTS - E2 -'/)                  07630000
      DO 63 IME=1,NMAX                                                  07640000
  63  WRITE(6,64) (ALPHA2(IME,JME),JME=1,NMAX)                          07650000
  64  FORMAT(1X,10F8.3)                                                 07660000
C                                                                       07670000
C OUTPUT >>> M1-CONVERSION COEFFICIENTS <<<                             07680000
C                                                                       07690000
      IF(I23.EQ.0) GOTO 5440                                            07700000
      WRITE(6,65)                                                       07710000
  65  FORMAT(1H0,1X,'CONVERSION COEFFICIENTS - M1 -'/)                  07720000
      DO 104 IMM=1,NMAX                                                 07730000
 104  WRITE(6,64) (ALPHA1(IMM,JMM),JMM=1,NMAX)                          07740000
 5440 IF(NE1E3.EQ.0) GOTO 103                                           07750000
C                                                                       07760000
C OUTPUT >>> E1-CONVERSION COEFFICIENTS <<<                             07770000
C                                                                       07780000
      WRITE(6,508)                                                      07790000
 508  FORMAT(1H0,1X,'CONVERSION-COEFFICIENTS - E1 -'/)                  07800000
      DO 505 IME=1,NMAX                                                 07810000
 505  WRITE(6,64) (ALPHE1(IME,JME),JME=1,NMAX)                          07820000
C                                                                       07830000
C OUTPUT >>> E3-CONVERSION COEFFICIENTS <<<                             07840000
C                                                                       07850000
      WRITE(6,506)                                                      07860000
 506  FORMAT(1H0,1X,'CONVERSION-COEFFICIENTS - E3 -'/)                  07870000
      DO 507 IME=1,NMAX                                                 07880000
 507  WRITE(6,64) (ALPHE3(IME,JME),JME=1,NMAX)                          07890000
 103  CONTINUE                                                          07900000
C                                                                       07910000
C FOR A COMPLETE ANGULAR DISTRIBUTION :  I19 = 1                        07920000
C OR A CALCULATION INCLUDING ONLY ONE SCATTERING ANGLE                  07930000
C THE INTEGRATION OVER THE PARTICLE-DETECTOR IS NOT POSSIBLE            07940000
C                                                                       07950000
      IF(I19.EQ.1) GOTO 1900                                            07960000
      IF(MMAX.EQ.1) GOTO 1900                                           07970000
C                                                                       07980000
      READ(5,5100) MZAHL,NORM                                           07990000
 5100 FORMAT(2I5)                                                       08000000
C                                                                       08010012
C  MZAHL = NUMBER OF THETA - INTEGRATIONS FOR THE PROJECTILE ANGULAR    08020012
C          DISTRIBUTION                                                 08030012
C  NORM  = POS.VALUE: NUMBER OF THE GAMMA-TRANSITION ON WHICH ALL OTHER 08040012
C                     GAMMA-YIELDS ARE NORMALIZED                       08050012
C        = NEG.VALUE: NORMALIZED ON THE RUTHERFORD CROSS SECTION        08060012
C                                                                       08070012
      IF(NORM.GE.0) WRITE(6,5110) MZAHL,NORM                            08080012
 5110 FORMAT(1H0,1X,'NUMBER OF THETA-INTEGRATIONS = ',I3/               08090000
     11X,'NUMBER OF THE TRANSITION FOR THE NORMALIZATION = ',I3/)       08100000
      IF(NORM.LT.0) WRITE(6,5111) MZAHL,NORM                            08110012
 5111 FORMAT(1H0,1X,'NUMBER OF THETA-INTEGRATIONS = ',I3/               08120012
     11X,'NORMALIZED TO THE RUTHERFORD CROSS SECTION = ',I3/)           08130012
C ----------------------------------------------------------------------08140012
      DO 7654 JZAHL=1,MZAHL                                             08150000
      READ(5,1302) XA,XE                                                08160000
 1302 FORMAT(2F10.2)                                                    08170000
C                                                                       08180012
C  XA,XE = SCATTERING ANGLE OF THE PROJECTILS IN THE CM - SYSTEM        08190012
C          INTEGRATION FROM XE TO XA (XA > XE)                          08200012
C                                                                       08210012
      WRITE(6,1303) XA,XE                                               08220000
      WRITE(20,1303) XA,XE                                              08230027
 1303 FORMAT(1X,'INTEGRATION FROM ',F10.2,' TO ',F10.2/)                08240000
      DO 91 K1=1,MMAX                                                   08250000
      K2=K1+1                                                           08260000
      IF(AXX(K1).GT.XA.AND.AXX(K2).LT.XA) GOTO 92                       08270000
      IF(AXX(K1).EQ.XA) GOTO 92                                         08280000
      GOTO 93                                                           08290000
C                                                                       08300000
C LINEAR INTERPOLATION AT THE BOUNDARIES OF THE PARTICLE DETECTOR       08310000
C                                                                       08320000
  92  DO 94 K3=1,LAUF                                                   08330000
  94  BYY(1,K3)=(AYY(K2,K3)-AYY(K1,K3))*(XA-AXX(K1))/                   08340000
     1(AXX(K2)-AXX(K1))+AYY(K1,K3)                                      08350000
      BRUTH(1)=((RUTHER(K2)-RUTHER(K1))*(XA-AXX(K1))/                   08360000
     1(AXX(K2)-AXX(K1)))+RUTHER(K1)                                     08370000
C     WRITE(6,705) XA,(BYY(1,K3),K3=1,LAUF)                             08380000
 705  FORMAT(1X,11E10.4)                                                08390000
      GOTO 91                                                           08400000
  93  IF(AXX(K1).GT.XE.AND.AXX(K2).LT.XE) GOTO 95                       08410000
      IF(AXX(K1).EQ.XE) GOTO 95                                         08420000
      GOTO 91                                                           08430000
C                                                                       08440000
C LINEAR INTERPOLATION AT THE BOUNDARIES OF THE PARTICLE DETECTOR       08450000
C                                                                       08460000
  95  DO 96 K3=1,LAUF                                                   08470000
  96  YE(K3)=(AYY(K2,K3)-AYY(K1,K3))*(XE-AXX(K1))/                      08480000
     1(AXX(K2)-AXX(K1))+AYY(K1,K3)                                      08490000
      ERUTH=((RUTHER(K2)-RUTHER(K1))*(XE-AXX(K1))/                      08500000
     1(AXX(K2)-AXX(K1)))+RUTHER(K1)                                     08510000
C     WRITE(6,705) XE,(YE(K3),K3=1,LAUF)                                08520000
  91  CONTINUE                                                          08530000
C                                                                       08540012
C XXX = THETA (CM. ANGLE OF PROJECTILE)                                 08550012
C                                                                       08560012
      XXX(1)=XA                                                         08570000
C                                                                       08580012
      IF(I25.EQ.0) GOTO 1002                                            08590012
C                                                                       08600012
C ETA = ETA  (LAB. ANGLE OF TARGET NUCLEUS)                             08610012
C XX = PSI (LAB. ANGLE OF PROJECTILE)                                   08620012
C                                                                       08630012
      ETA=((180.-XA)/2.)*3.1415926/180.                                 08640000
      XX(1)=ATAN(SIN(2.*ETA)/((A1/A2)-COS(2.*ETA)))                     08650000
      if(xx(1).lt.0) xx(1)=3.1415926+xx(1)
      GOTO 1003                                                         08660012
C                                                                       08670012
 1002 XX(1)=((180.-XA)/2.)*3.1415926/180.                               08680012
C                                                                       08690012
C XX = ETA  (LAB. ANGLE OF TARGET NUCLEUS)                              08700012
C                                                                       08710012
 1003 I=2                                                               08720012
      DO 97 K1=1,MMAX                                                   08730000
      IF(AXX(K1).GT.XA) GOTO 97                                         08740000
      IF(AXX(K1).LT.XE) GOTO 98                                         08750000
      XXX(I)=AXX(K1)                                                    08760000
      XX(I)=AAXX(K1)                                                    08770000
      BRUTH(I)=RUTHER(K1)                                               08780000
      DO 99 K2=1,LAUF                                                   08790000
  99  BYY(I,K2)=AYY(K1,K2)                                              08800000
      I=I+1                                                             08810000
  97  CONTINUE                                                          08820000
  98  XXX(I)=XE                                                         08830000
      BRUTH(I)=ERUTH                                                    08840000
C                                                                       08850000
      IF(I25.EQ.0) GOTO 1004                                            08860012
C                                                                       08870012
C ETA = ETA  (LAB. ANGLE OF TARGET NUCLEUS)                             08880000
C XX = PSI (LAB. ANGLE OF PROJECTILE)                                   08890000
C                                                                       08900000
      ETA=((180.-XE)/2.)*3.1415926/180.                                 08910000
      XX(I)=ATAN(SIN(2.*ETA)/((A1/A2)-COS(2.*ETA)))                     08920000
      if(xx(i).lt.0) xx(i)=3.1415926+xx(i)
      GOTO 1005                                                         08930012
C                                                                       08940000
 1004 XX(I)=((180.-XE)/2.)*3.1415926/180.                               08950012
C                                                                       08960012
C XX = ETA  (LAB. ANGLE OF TARGET NUCLEUS)                              08970012
C                                                                       08980012
 1005 DO 700 K2=1,LAUF                                                  08990012
 700  BYY(I,K2)=YE(K2)                                                  09000000
C ----------------------------------------------------------------------09010012
      IF(NORM.LT.0) GOTO 1006                                           09020012
      NNORM=LAUF-NORM                                                   09030000
      DO 5101 I1=1,I                                                    09040000
 5101 YY(I1)=BYY(I1,NNORM)                                              09050000
C                                                                       09060000
C SUBROUTINE FOR THE INTEGRATION OVER THE PARTICLE-DETECTOR             09070000
C (THETAP)                                                              09080000
C                                                                       09090000
      CALL DQTFG(XX,YY,ZZ,I)                                            09100000
      ANORM=ZZ(I)                                                       09110000
      GOTO 1007                                                         09120012
C                                                                       09130000
C SUBROUTINE FOR THE RUTHERFORD INTEGRATION OVER THE PARTICLE-DETECTOR  09140000
C (THETAP)                                                              09150000
C                                                                       09160000
 1006 CALL DQTFG(XX,BRUTH,ZZ,I)                                         09170012
      ANORM=ZZ(I)*(VI1-VI2)/57.2958                                     09180026
      WRITE(6,1008) ANORM                                               09190012
 1008 FORMAT(1X,'RUTHERFORD CROSS-SECTION=',E15.8)                      09200012
C                                                                       09210024
C RUTHERFORD CROSS-SECTION * 4. * 3.1415926 TO GET PROBABILITIES        09220024
C                                                                       09230024
      ANORM=ANORM/12.5664                                               09240024
C                                                                       09250000
 1007 DO 701 K2=1,LAUF                                                  09260012
      DO 702 I1=1,I                                                     09270000
      YY(I1)=BYY(I1,K2)                                                 09280000
c      write(6,1771) xx(i1),yy(i1)
C 1771 format(1x,2e10.3)
 702  CONTINUE                                                          09290000
C                                                                       09300000
C SUBROUTINE FOR THE INTEGRATION OVER THE PARTICLE-DETECTOR             09310000
C (THETAP)                                                              09320000
C                                                                       09330000
      CALL DQTFG(XX,YY,ZZ,I)                                            09340000
      CROSSN(K2)=ZZ(I)/ANORM                                            09350000
      IF(JZAHL.GT.1) GOTO 704                                           09360000
      WRITE(6,706) DELTA1(K2),IITRAN(K2),IFTRAN(K2),ZZ(I),CROSSN(K2)    09370000
 706  FORMAT(1X,'TRAN.-ENERGY (KEV) = ',F6.1,                           09380000
     13X,'CROSS-SECTION (#',I2,'-->#',I2,') = ',E15.8,3X,F10.5)         09390012
C     WRITE(20) IITRAN(K2),IFTRAN(K2),ZZ(I),CROSSN(K2)                  09400027
      WRITE(20,9980) DELTA1(K2),CROSSN(K2)                              09410027
 9980 FORMAT(1X,F6.1,E15.8)                                             09420027
      GOTO 701                                                          09430000
C                                                                       09440000
 704  WRITE(6,703) IITRAN(K2),IFTRAN(K2),ZZ(I),CROSSN(K2)               09450000
C     WRITE(20) IITRAN(K2),IFTRAN(K2),ZZ(I),CROSSN(K2)                  09460027
      WRITE(20,9980) DELTA1(K2),CROSSN(K2)                              09470027
 703  FORMAT(1X,'CROSS-SECTION (#',I2,'-->#',I2,') = ',E15.8,5X,F10.5)  09480012
 701  CONTINUE                                                          09490000
C     WRITE(20) (CROSSN(IW12),IW12=1,LAUF)                              09500000
7654  CONTINUE                                                          09510000
1900  STOP                                                              09520000
      END                                                               09530000
      SUBROUTINE ANGPAR(TT,VI1,VI2,DS,LAUF,SINLAB)                      09540000
C     -----------------                                                 09550013
C SUBROUTINE FOR CALCULATION OF THE DOUBLE DIFFERENTIAL CROSS SECTION   09560013
C INCLUDING THE INTEGRATION OVER THE FINAL PARTICLE DETECTOR            09570013
C (AZIMUTHAL-ANGLE)                                                     09580013
C                                                                       09590012
      EXTERNAL FCT                                                      09600000
      COMMON/C1/ E(35),SPIN(35),EE(35),ME2(35,35),MM1(35,35),NMAX       09610000
      COMMON/C2/ ALPHA1(35,35),ALPHA2(35,35),ANG(35,5,5)                09620000
      COMMON/C12/ ALPHE1(35,35),ALPHE3(35,35),DELTA1(199)               09630000
      COMMON/C3/ I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I2409640000
     1,I25                                                              09650020
      COMMON/C4/ Z1,Z2,A1,A2,EP,THETAP,QG(5,35),QGS(5,35),IPAR(35)      09660000
      COMMON/C14/ AJ,GZA,TAUC,H0(35)                                    09670000
C **********************************************************************09680000
C   >>>>  FOR THE HAGER AND SELTZER INTERPOLATION PROGRAM  <<<<         09690000
C **********************************************************************09700000
      COMMON/C5/ ALPHAK(23,4),ALPHAL(23,4),ALPHAM(23,4),NCCK,NCCL,NCCM  09710000
      COMMON/C6/ CCE1K(4),CCE2K(4),CCE1L(4),CCE2L(4),CCE1M(4),CCE2M(4)  09720000
      COMMON/C8/ ITHETA                                                 09730000
      COMMON/C9/ CCKMIN(4),CCLMIN(4),CCMMIN(4)                          09740000
      COMMON/C19/ CCKMAX(4),CCLMAX(4),CCMMAX(4)                         09750000
      DIMENSION CCEE(23),ALPHA(23)                                      09760000
C **********************************************************************09770000
      COMMON/RUTH/SIGLAB                                                09780007
      COMMON/C7/ T,RESTLA,KLAUF,KKLAUF,NHILF,IKHILF                     09790000
      COMMON/C10/ VGAMMA,Q(5),ME(35,35,4),IITRAN(199),IFTRAN(199)       09800000
      DIMENSION E1(35),IK(35),Y(5,5),G(35,35,5),F(5,35,35),W(35,35)     09810000
      DIMENSION NN(35),DS(199),SINLAB(199)                              09820000
      DIMENSION TTP(80),VGP(80),DSP(80),TAU(35)                         09830000
      REAL ME,ME2,MM1                                                   09840000
      REAL*4 DSNORM(80),TTPP(80)                                        09850000
C                                                                       09860000
C ONLY FOR THE PLOT OUTPUT NECESSARY                                    09870000
C                                                                       09880000
c      INTEGER*2 TAB(100)/                                              09890000
c     A ' 0',' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9',              09900000
c     1 '10','11','12','13','14','15','16','17','18','19',              09910000
c     2 '20','21','22','23','24','25','26','27','28','29',              09920000
c     3 '30','31','32','33','34','35','36','37','38','39',              09930000
c     4 '40','41','42','43','44','45','46','47','48','49',              09940000
c     5 '50','51','52','53','54','55','56','57','58','59',              09950000
c     6 '60','61','62','63','64','65','66','67','68','69',              09960000
c     7 '70','71','72','73','74','75','76','77','78','79',              09970000
c     8 '80','81','82','83','84','85','86','87','88','89',              09980000
c     9 '90','91','92','93','94','95','96','97','98','99'/              09990000
c      INTEGER*2 TAB1(100)/                                             10000000
c     A '0 ','1 ','2 ','3 ','4 ','5 ','6 ','7 ','8 ','9 ',              10010000
c     1 '10','11','12','13','14','15','16','17','18','19',              10020000
c     2 '20','21','22','23','24','25','26','27','28','29',              10030000
c     3 '30','31','32','33','34','35','36','37','38','39',              10040000
c     4 '40','41','42','43','44','45','46','47','48','49',              10050000
c     5 '50','51','52','53','54','55','56','57','58','59',              10060000
c     6 '60','61','62','63','64','65','66','67','68','69',              10070000
c     7 '70','71','72','73','74','75','76','77','78','79',              10080000
c     8 '80','81','82','83','84','85','86','87','88','89',              10090000
c     9 '90','91','92','93','94','95','96','97','98','99'/              10100000
C                                                                       10110000
C IPLOTT IS A SCALER FOR THE PLOTT OUTPUT                               10120000
C                                                                       10130000
c      IPLOTT=0                                                         10140000
C                                                                       10150000
C LAUF IS A SCALER FOR THE CROSS SECTION                                10160000
C                                                                       10170000
      LAUF=1                                                            10180000
C                                                                       10190000
C TT = ANGLE OF THE GAMMA DETECTOR                                      10200000
C THETAP = ANGLE OF THE SCATTERED PARTICLE                              10210000
C VI2 > VI1 BOUNDARY OF THE PARTICLE COUNTER (AZIMUTHAL ANGLE)          10220025
C                                                                       10230000
      PI=3.1415926                                                      10240000
      T=PI*TT/180.                                                      10250000
      THET=PI*THETAP/180.                                               10260000
      V1=PI*VI1/180.                                                    10270000
      V2=PI*VI2/180.                                                    10280000
C                                                                       10290000
C ARRANGEMENT OF THE ENERGY LEVELS  E(IN) > E(IN+1)                     10300000
C                                                                       10310000
      DO 1 J=1,NMAX                                                     10320000
    1 E1(J)=E(J)                                                        10330000
      N=1                                                               10340000
    4 NN(N)=0                                                           10350000
      EE(N)=-0.0001                                                     10360000
      DO 2 J=1,NMAX                                                     10370000
      IF(EE(N)-E1(J)) 3,2,2                                             10380000
    3 NN(N)=J                                                           10390000
      EE(N)=E1(J)                                                       10400000
    2 CONTINUE                                                          10410000
      E1(NN(N))=-0.0001                                                 10420000
      N=N+1                                                             10430000
      IF(N-NMAX) 4,4,5                                                  10440000
C                                                                       10450000
C THE POSSIBLE DECAYS OF THE NUCLEAR STATE                              10460000
C                                                                       10470000
    5 JJJ1=NMAX-1                                                       10480000
C                                                                       10490000
C                                                                       10500000
C                                                                       10510000
      DO 6 JJJ=1,JJJ1                                                   10520000
C                                                                       10530000
C                                                                       10540000
C                                                                       10550000
      J1=0                                                              10560000
      J2=JJJ+1                                                          10570000
      DO 7 J=J2,NMAX                                                    10580000
C                                                                       10590000
      IF((IPAR(NN(JJJ))+IPAR(NN(J))).NE.0) GOTO 402                     10600000
C                                                                       10610000
      IF (ABS(SPIN(NN(JJJ))-SPIN(NN(J))).GT.3.0) GOTO 7                 10620000
      GOTO 403                                                          10630000
  402 IF (ABS(SPIN(NN(JJJ))-SPIN(NN(J))).GT.2.0) GOTO 7                 10640000
  403 J1=J1+1                                                           10650000
      IK(J1)=NN(J)                                                      10660000
    7 CONTINUE                                                          10670000
C                                                                       10680000
      IF(I12.EQ.0) GOTO 201                                             10690000
      WRITE(6,1003)                                                     10700000
 1003 FORMAT(1H0,1X,'THE EXCITED STATE AND ITS POSSIBLE DECAYS :    ',/ 10710000
     11X,'SPIN(I)   SPIN(F)  SPIN(F)  SPIN(F)  SPIN(F)  SPIN(F)')       10720000
      WRITE(6,1004) SPIN(NN(JJJ)),(SPIN(IK(J)),J=1,J1)                  10730000
 1004 FORMAT(2X,10(F5.2,4X))                                            10740000
      WRITE(6,1901)                                                     10750000
 1901 FORMAT(1H0)                                                       10760000
  201 CONTINUE                                                          10770000
C                                                                       10780000
C CALCULATION OF THE QUANTITY :  H(N)                                   10790000
C                                                                       10800000
      N=JJJ                                                             10810000
      IF(ITHETA.GT.1) GOTO 530                                          10820000
      HN=0.0                                                            10830000
      DO 8 J=1,J1                                                       10840000
C **********************************************************************10850000
C        >>>> HAGER AND SELTZER INTERPOLATION PROGRAM <<<<              10860000
C **********************************************************************10870000
      IF(ITHETA.GT.1) GOTO 524                                          10880000
      MULT=0                                                            10890000
C                                                                       10900000
C NN(N) = INITIAL-STATE                                                 10910000
C IK(J) = FINAL-STATE                                                   10920000
C LEVEL ENERGY (INITIAL-) > LEVEL ENERGY (FINAL-STATE)                  10930000
C                                                                       10940000
      CCEE(1)=(E(NN(N))-E(IK(J)))*1000.                                 10950000
      IF((IPAR(NN(N))+IPAR(IK(J))).NE.0) GOTO 517                       10960000
      IF(ABS(SPIN(NN(N))-SPIN(IK(J))).GT.1.) GOTO 519                   10970000
      MULT=1                                                            10980000
      GOTO 518                                                          10990000
  517 MULT=2                                                            11000000
  518 IF(CCEE(1).LT.CCKMAX(MULT).AND.CCEE(1).GT.CCKMIN(MULT)) GOTO 501  11010000
      IF(CCEE(1).GT.CCKMAX(MULT)) GOTO 502                              11020000
      WRITE(6,503)                                                      11030000
  503 FORMAT(1X,'GAMMA-ENERGY LESS THAN CCKMIN >>> AKTOT=0.0 <<<')      11040000
      AKTOT=0.0                                                         11050000
      GOTO 504                                                          11060000
  502 AKTOT=0.0                                                         11070000
      GOTO 504                                                          11080000
  501 CONTINUE                                                          11090000
      DO 505 IW=1,NCCK                                                  11100000
  505 ALPHA(IW)=ALPHAK(IW,MULT)                                         11110000
      CCE1=CCE1K(MULT)                                                  11120000
      CCE2=CCE2K(MULT)                                                  11130000
C                                                                       11140000
      CALL FIT(CCE1,CCE2,1,MULT,NCCK,ALPHA,1,CCEE)                      11150000
C                                                                       11160000
      AKTOT=ALPHA(1)                                                    11170000
  504 CONTINUE                                                          11180000
      IF(CCEE(1).LT.CCLMAX(MULT).AND.CCEE(1).GT.CCLMIN(MULT)) GOTO 506  11190000
      IF(CCEE(1).GT.CCLMAX(MULT)) GOTO 507                              11200000
      WRITE(6,508)                                                      11210000
  508 FORMAT(1X,'GAMMA-ENERGY LESS THAN CCLMIN >>> ALTOT=0.0 <<<')      11220000
      ALTOT=0.0                                                         11230000
      GOTO 509                                                          11240000
  507 ALTOT=0.0                                                         11250000
      GOTO 509                                                          11260000
  506 CONTINUE                                                          11270000
      DO 510 IW=1,NCCL                                                  11280000
  510 ALPHA(IW)=ALPHAL(IW,MULT)                                         11290000
      CCE1=CCE1L(MULT)                                                  11300000
      CCE2=CCE2L(MULT)                                                  11310000
C                                                                       11320000
      CALL FIT(CCE1,CCE2,2,MULT,NCCL,ALPHA,1,CCEE)                      11330000
C                                                                       11340000
      ALTOT=ALPHA(1)                                                    11350000
  509 CONTINUE                                                          11360000
      IF(CCEE(1).LT.CCMMAX(MULT).AND.CCEE(1).GT.CCMMIN(MULT)) GOTO 511  11370000
      IF(CCEE(1).GT.CCMMAX(MULT)) GOTO 512                              11380000
      WRITE(6,513)                                                      11390000
  513 FORMAT(1X,'GAMMA-ENERGY LESS THAN CCMMIN >>> AMTOT=0.0 <<<')      11400000
      AMTOT=0.0                                                         11410000
      GOTO 514                                                          11420000
  512 AMTOT=0.0                                                         11430000
      GOTO 514                                                          11440000
  511 CONTINUE                                                          11450000
      DO 515 IW=1,NCCM                                                  11460000
  515 ALPHA(IW)=ALPHAM(IW,MULT)                                         11470000
      CCE1=CCE1M(MULT)                                                  11480000
      CCE2=CCE2M(MULT)                                                  11490000
C                                                                       11500000
      CALL FIT(CCE1,CCE2,3,MULT,NCCM,ALPHA,1,CCEE)                      11510000
C                                                                       11520000
      AMTOT=ALPHA(1)                                                    11530000
  514 CONTINUE                                                          11540000
      IF(MULT.EQ.1) GOTO 520                                            11550000
      IF(MULT.EQ.2) GOTO 521                                            11560000
      IF(MULT.EQ.3) GOTO 522                                            11570000
      IF(MULT.EQ.4) GOTO 526                                            11580000
      GOTO 500                                                          11590000
  520 ALPHE1(NN(N),IK(J))=AKTOT+ALTOT+AMTOT                             11600000
      ALPHE1(IK(J),NN(N))=ALPHE1(NN(N),IK(J))                           11610000
      GOTO 523                                                          11620000
  521 ALPHA2(NN(N),IK(J))=AKTOT+ALTOT+AMTOT                             11630000
      ALPHA2(IK(J),NN(N))=ALPHA2(NN(N),IK(J))                           11640000
      GOTO 523                                                          11650000
  522 ALPHE3(NN(N),IK(J))=AKTOT+ALTOT+AMTOT                             11660000
      ALPHE3(IK(J),NN(N))=ALPHE3(NN(N),IK(J))                           11670000
      GOTO 523                                                          11680000
  526 ALPHA1(NN(N),IK(J))=AKTOT+ALTOT+AMTOT                             11690000
      ALPHA1(IK(J),NN(N))=ALPHA1(NN(N),IK(J))                           11700000
      WRITE(6,527) AKTOT,ALTOT,AMTOT                                    11710000
  527 FORMAT(1X,'CCK = ',F10.3,' CCL = ',F10.3,' CCM = ',F10.3,' -- M1 -11720000
     1-'/)                                                              11730000
      GOTO 519                                                          11740000
  523 WRITE(6,516) AKTOT,ALTOT,AMTOT,MULT                               11750000
  516 FORMAT(1X,'CCK = ',F10.3,' CCL = ',F10.3,' CCM = ',F10.3,' -- E', 11760000
     1I1,' --'/)                                                        11770000
  519 IF((IPAR(NN(N))+IPAR(IK(J))).NE.0) GOTO 500                       11780000
      IF(MULT.EQ.3) GOTO 524                                            11790000
      MULT=3                                                            11800000
      GOTO 518                                                          11810000
  500 IF(I23.EQ.0.OR.MULT.EQ.4) GOTO 524                                11820000
      IF(ABS(SPIN(NN(N))-SPIN(IK(J))).GT.1.) GOTO 524                   11830000
      MULT=4                                                            11840000
      GOTO 518                                                          11850000
  524 CONTINUE                                                          11860000
C **********************************************************************11870000
C        >>>> HAGER AND SELTZER INTERPOLATION PROGRAM <<<<              11880000
C **********************************************************************11890000
C                                                                       11900000
C IF THE INITIAL STATE OR THE FINAL STATE HAS NEGATIVE PARITY THE SUM   11910000
C OVER L INCLUDES >>> E1 AND E3 <<< MULTIPOLES                          11920000
C                                                                       11930000
C IF THE INITIAL STATE AND THE FINAL STATE HAS THE SAME PARITY THE SUM  11940000
C OVER L INCLUDES >>> E2 AND M1 <<< MULTIPOLES                          11950000
C                                                                       11960000
      IF((IPAR(NN(N))+IPAR(IK(J))).NE.0) GOTO 400                       11970000
C                                                                       11980000
      AA=1.59034E17*(E(NN(N))-E(IK(J)))**3*ME(NN(N),IK(J),1)**2         11990000
     1  /(2.*SPIN(NN(N))+1.)*(1.+ALPHE1(NN(N),IK(J)))                   12000000
     2  +5.70924E8*(E(NN(N))-E(IK(J)))**7*ME(NN(N),IK(J),3)**2          12010000
     3  /(2.*SPIN(NN(N))+1.)*(1.+ALPHE3(NN(N),IK(J)))                   12020000
      GOTO 401                                                          12030000
  400 AA=1.22537E13*(E(NN(N))-E(IK(J)))**5*ME2(NN(N),IK(J))**2          12040000
     1  /(2.*SPIN(NN(N))+1.)*(1.+ALPHA2(NN(N),IK(J)))                   12050000
     2  +1.75836E13*(E(NN(N))-E(IK(J)))**3*MM1(NN(N),IK(J))**2          12060000
     3  /(2.*SPIN(NN(N))+1.)*(1.+ALPHA1(NN(N),IK(J)))                   12070000
  401 HN=HN+AA                                                          12080000
    8 CONTINUE                                                          12090000
      HN=1./HN                                                          12100000
      TAU(NN(N))=HN                                                     12110000
C                                                                       12120000
      IF(I13.EQ.0) GOTO 202                                             12130000
      WRITE(6,1005) NN(N),HN                                            12140000
 1005 FORMAT(1X,'HN = LIFETIME (SEC) OF THE STATE #',I2,' = ',E15.7/)   12150000
  202 CONTINUE                                                          12160000
  530 CONTINUE                                                          12170000
C                                                                       12180000
C DECAY CALCULATION OF THE NUCLEAR STATE WITH SPIN(NN) AND ITS          12190000
C INTENSITY DISTRIBUTION ON NUCLEAR STATES WITH SPIN(IK) AND            12200000
C LOWER EXCITATION ENERGY                                               12210000
C                                                                       12220000
C                                                                       12230000
C                                                                       12240000
      DO 9 J=1,J1                                                       12250000
C                                                                       12260000
C                                                                       12270000
C                                                                       12280000
      IF(ITHETA.GT.1) GOTO 531                                          12290000
C                                                                       12300000
      DO 10 K=1,5,2                                                     12310000
C                                                                       12320000
C CALCULATION OF THE QUANTITY :  GK (IM,IN)                             12330000
C                                                                       12340000
      VORFAC=HN*SQRT((2.*SPIN(IK(J))+1.)*(2.*SPIN(NN(N))+1.))           12350000
C                                                                       12360000
      IF((IPAR(NN(N))+IPAR(IK(J))).NE.0) GOTO 405                       12370000
C                                                                       12380000
      S1=(-1.)**(IFIX(3.*SPIN(IK(J))+3.*SPIN(NN(N))+1.00001))           12390000
      S2=(-1.)**(IFIX(3.*SPIN(IK(J))+3.*SPIN(NN(N))+3.00001))           12400000
      AME1=SQRT(1.59034E17*(E(NN(N))-E(IK(J)))**3)*ME(NN(N),IK(J),1)    12410000
     1  /SQRT(2.*SPIN(NN(N))+1.)                                        12420000
      AME3=SQRT(5.70924E8*(E(NN(N))-E(IK(J)))**7)*ME(NN(N),IK(J),3)     12430000
     1  /SQRT(2.*SPIN(NN(N))+1.)                                        12440000
      G(IK(J),NN(N),K)=VORFAC*(S1*RACAHI(SPIN(IK(J)),SPIN(IK(J)),       12450000
     1 SPIN(NN(N)),SPIN(NN(N)),FLOAT(K-1),1.)*AME1**2                   12460000
     2 *(1.+ALPHE1(NN(N),IK(J)))                                        12470000
     3 +S2*RACAHI(SPIN(IK(J)),SPIN(IK(J)),SPIN(NN(N)),SPIN(NN(N)),      12480000
     4 FLOAT(K-1),3.)*(AME3**2)*(1.+ALPHE3(NN(N),IK(J))))               12490000
       IF((IPAR(NN(N))+IPAR(IK(J))).EQ.0) GOTO 10                       12500000
C                                                                       12510000
  405 S1=(-1.)**(IFIX(3.*SPIN(IK(J))+3.*SPIN(NN(N))+0.00001))           12520000
      S2=(-1.)**(IFIX(3.*SPIN(IK(J))+3.*SPIN(NN(N))+1.00001))           12530000
      AM2=SQRT(1.22537E13*(E(NN(N))-E(IK(J)))**5)*ME2(NN(N),IK(J))      12540000
     1   /SQRT(2.*SPIN(NN(N))+1.)                                       12550000
C     AM1=SQRT(1.75836E13*(E(NN(N))-E(IK(J)))**3)*MM1(NN(N),IK(J))      12560000
C    1   /SQRT(2.*SPIN(NN(N))+1.)                                       12570000
      AM1=-SQRT(1.75836E13*(E(NN(N))-E(IK(J)))**3)*MM1(NN(N),IK(J))     12580000
     1   /SQRT(2.*SPIN(NN(N))+1.)                                       12590000
      G(IK(J),NN(N),K)=VORFAC*(S1*RACAHI(SPIN(IK(J)),SPIN(IK(J)),       12600000
     1 SPIN(NN(N)),SPIN(NN(N)),FLOAT(K-1),2.)*AM2**2                    12610000
     2 *(1.+ALPHA2(NN(N),IK(J)))                                        12620000
     3 +S2*RACAHI(SPIN(IK(J)),SPIN(IK(J)),SPIN(NN(N)),SPIN(NN(N)),      12630000
     4 FLOAT(K-1),1.)*(AM1**2)*(1.+ALPHA1(NN(N),IK(J))))                12640000
   10 CONTINUE                                                          12650000
C                                                                       12660000
C DELTA1 = TRANSITION-ENERGY                                            12670000
C                                                                       12680000
      DELTA1(LAUF)=(E(NN(N))-E(IK(J)))*1000.                            12690000
C                                                                       12700000
      IF(I14.EQ.0) GOTO 203                                             12710000
      WRITE(6,1006)                                                     12720000
 1006 FORMAT(1X,'  G0(IF,II)   G2(IF,II)   G4(IF,II)'/)                 12730000
      WRITE(6,1007) (G(IK(J),NN(N),K),K=1,5,2)                          12740000
 1007 FORMAT(1X,3(E11.4,1X)//)                                          12750000
  203 CONTINUE                                                          12760000
  531 CONTINUE                                                          12770000
C                                                                       12780000
      IF(I22.LT.2) GOTO 126                                             12790000
C ----------------------------------------------------------------------12800013
C CALCULATION OF THE ATTENUATION COEFFICIENTS >> QG << USING THE        12810000
C     SCHERER-BLUME THEORY OR ABRAGAM-POUND THEORY                      12820000
C     (NOT VELOCITY-DEPENDED)                                           12830000
C                                                                       12840000
      VCM=0.04634*SQRT(EP/A1)/(1.+A2/A1)                                12850000
      DELTA=E(NN(N))                                                    12860000
      AKA=(1.+A1/A2)*DELTA/EP                                           12870000
      AKAPA1=(A2/A1)*SQRT(1.-AKA)                                       12880013
      AKAPA2=SQRT(1.-AKA)                                               12890013
C                                                                       12900013
      IF(I25.EQ.0) GOTO 2000                                            12910013
C                                                                       12920013
C VV2 = VELOCITY OF THE PROJECTILE                                      12930013
C                                                                       12940013
      VV2=VCM*SQRT(1.+(AKAPA1**2)+2.*AKAPA1*COS(THET))                  12950013
      GOTO 2001                                                         12960013
C                                                                       12970013
 2000 VV2=VCM*SQRT(1.+(AKAPA2**2)-2.*AKAPA2*COS(THET))                  12980013
C                                                                       12990013
C VV2 = VELOCITY OF THE TARGET NUCLEUS                                  13000013
C                                                                       13010013
C                                                                       13020000
C VELOCITY DEPENDENCE OF DEORIENTATION COEFFICIENT G2,G4 :              13030000
C                                                                       13040000
C G2 = 0.96304 - 0.08579 * LN(TAU) , TAU {PSEC}                         13050000
C G4 = 3 * G2 / (10 - 7 * G2)                                           13060000
C                                                                       13070000
 2001 G2PAR=0.96304-0.08579*ALOG(TAU(NN(N))*1.E+12)                     13080013
C                                                                       13090013
      IF(G2PAR.GT.1.0) G2PAR=1.                                         13100000
      IF(G2PAR.LT.0.0) G2PAR=0.0                                        13110000
      IF(I22.EQ.5) QG(1,NN(N))=1.                                       13120000
      IF(I22.EQ.5) QGS(3,NN(N))=G2PAR                                   13130000
C                                                                       13140000
      XVV2=1.                                                           13150000
      AVV2=(1.-QGS(3,NN(N)))/(QGS(3,NN(N))*(0.062**XVV2))               13160000
      G2VV2=1./(1.+AVV2*(VV2**XVV2))                                    13170000
      G4VV2=3.*G2VV2/(10.-7.*G2VV2)                                     13180000
C                                                                       13190000
      IF(I22.EQ.2.OR.I22.EQ.5) QG(3,NN(N))=G2VV2                        13200000
      IF(I22.EQ.2.OR.I22.EQ.5) QG(5,NN(N))=G4VV2                        13210000
      IF(I22.EQ.2.OR.I22.EQ.5) GOTO 131                                 13220000
C                                                                       13230000
C     OMEGA=4789.4413*GZA*H0                                            13240000
      OMEGA=4789.4413*GZA*H0(NN(N))*(VV2**0.6)                          13250000
      OMEGAJ=SQRT((AJ+1.)/AJ)*OMEGA                                     13260000
      IF(I22.EQ.4) GOTO 130                                             13270000
      AKK2=1./(1.+OMEGAJ*TAUC*(((2.*SPIN(NN(N))+1.)**2/(4.*AJ*(AJ+1.)))+13280000
     117.*(2.*AJ-1)*(2.*AJ+3.)/(20.*AJ*(AJ+1.))))                       13290000
      AKK4=1./(1.+OMEGAJ*TAUC*(((2.*SPIN(NN(N))+1.)**2/(4.*AJ*(AJ+1.)))+13300000
     159.*(2.*AJ-1)*(2.*AJ+3.)/(20.*AJ*(AJ+1.))))                       13310000
      QG(1,NN(N))=1.                                                    13320000
      QG(3,NN(N))=1./(1.+2.*(OMEGAJ**2)*TAUC*TAU(NN(N))*AKK2)           13330000
      QG(5,NN(N))=1./(1.+(20./3.)*(OMEGAJ**2)*TAUC*TAU(NN(N))*AKK4)     13340000
      GOTO 131                                                          13350000
  130 QG(1,NN(N))=1.                                                    13360000
      QG(3,NN(N))=1./(1.+2.*(OMEGAJ**2)*TAUC*TAU(NN(N)))                13370000
      QG(5,NN(N))=1./(1.+(20./3.)*(OMEGAJ**2)*TAUC*TAU(NN(N)))          13380000
  131 WRITE(6,125) NN(N),QG(1,NN(N)),QG(3,NN(N)),QG(5,NN(N))            13390000
  125 FORMAT(1X,'ATTENUATION COEFFICIENTS #',I2,' G0 = ',F6.3,' G2 = '  13400000
     1,F6.3,' G4 = ',F6.3)                                              13410000
C ----------------------------------------------------------------------13420013
C                                                                       13430000
C CALCULATION OF THE QUANTITY :  A(K,KK)  K=0,2,4   -K <KK< +K          13440000
C                                                                       13450000
  126 DO 11 K=1,5,2                                                     13460000
      DO 12 KK=1,K                                                      13470000
C                                                                       13480000
C DEORIENTATION CORRECTION OF THE GAMMA-RAY ANGULAR DISTRIBUTION        13490000
C                                                                       13500000
      IF(I22.LT.2) GOTO 127                                             13510000
      ANG(IK(J),K,KK)=ANG(IK(J),K,KK)+ANG(NN(N),K,KK)*G(IK(J),NN(N),K)* 13520000
     1 QG(K,NN(N))                                                      13530000
      GOTO 12                                                           13540000
  127 ANG(IK(J),K,KK)=ANG(IK(J),K,KK)+ANG(NN(N),K,KK)*G(IK(J),NN(N),K)  13550000
   12 CONTINUE                                                          13560000
   11 CONTINUE                                                          13570000
      IF(I16.EQ.0) GOTO 205                                             13580000
      WRITE(6,1010)                                                     13590000
 1010 FORMAT(2X,'ANG(0,0)   ANG(2,2)   ANG(2,1)   ANG(2,0)   ANG(4,4) ' 13600000
     1 ,'  ANG(4,3)   ANG(4,2)   ANG(4,1)   ANG(4,0)'/)                 13610000
      WRITE(6,1011) ANG(IK(J),1,1),ANG(IK(J),3,3),ANG(IK(J),3,2),       13620000
     1 ANG(IK(J),3,1),ANG(IK(J),5,5),ANG(IK(J),5,4),ANG(IK(J),5,3),     13630000
     2 ANG(IK(J),5,2),ANG(IK(J),5,1)                                    13640000
 1011 FORMAT(10(1X,E10.4))                                              13650000
  205 CONTINUE                                                          13660000
C                                                                       13670000
C CALCULATION OF THE QUANTITY :  FK (IM,IN)                             13680000
C                                                                       13690000
      IF(ITHETA.GT.1) GOTO 532                                          13700000
      DO 13 K=1,5,2                                                     13710000
C                                                                       13720000
      IF((IPAR(NN(N))+IPAR(IK(J))).NE.0) GOTO 406                       13730000
C                                                                       13740000
      F(K,IK(J),NN(N))=HN*(                                             13750000
     1 AME1**2*FK(SPIN(IK(J)),1.,1.,SPIN(NN(N)),FLOAT(K-1))             13760000
     2 +AME1*AME3*FK(SPIN(IK(J)),1.,3.,SPIN(NN(N)),FLOAT(K-1))          13770000
     3 +AME1*AME3*FK(SPIN(IK(J)),3.,1.,SPIN(NN(N)),FLOAT(K-1))          13780000
     4 +AME3**2*FK(SPIN(IK(J)),3.,3.,SPIN(NN(N)),FLOAT(K-1)))           13790000
      GOTO 408                                                          13800000
  406 F(K,IK(J),NN(N))=HN*(                                             13810000
     1 AM2**2*FK(SPIN(IK(J)),2.,2.,SPIN(NN(N)),FLOAT(K-1))              13820000
     2 +AM2*AM1*FK(SPIN(IK(J)),2.,1.,SPIN(NN(N)),FLOAT(K-1))            13830000
     3 +AM1*AM2*FK(SPIN(IK(J)),1.,2.,SPIN(NN(N)),FLOAT(K-1))            13840000
     4 +AM1**2*FK(SPIN(IK(J)),1.,1.,SPIN(NN(N)),FLOAT(K-1)))            13850000
  408 CONTINUE                                                          13860000
C                                                                       13870000
C FINITE SOLID-ANGLE CORRECTION FOR THE GE(LI)-DETECTOR                 13880000
C                                                                       13890000
      IF(I22.EQ.0) GOTO 13                                              13900000
      F(K,IK(J),NN(N))=F(K,IK(J),NN(N))*Q(K)                            13910000
   13 CONTINUE                                                          13920000
C                                                                       13930000
      IF(I15.EQ.0) GOTO 204                                             13940000
      WRITE(6,1008)                                                     13950000
 1008 FORMAT(1H0,1X,' F0(IF,II)   F2(IF,II)   F4(IF,II)'/)              13960000
      WRITE(6,1009) (F(K,IK(J),NN(N)),K=1,5,2)                          13970000
 1009 FORMAT(1X,3(E11.4,1X)//)                                          13980000
  204 CONTINUE                                                          13990000
  532 CONTINUE                                                          14000000
C                                                                       14010000
C CALCULATION OF THE GAMMA-RAY ANGULAR DISTRIBUTION :  W (IN->IM)       14020000
C INCLUDING THE INTEGRATION OVER THE PARTICLE COUNTER                   14030000
C                                                                       14040000
C                                                                       14050000
C GAMMA-RAY ANGULAR DISTRIBUTION                                        14060000
C                                                                       14070000
      IF(I19.EQ.0) GOTO 1970                                            14080000
      VGAMMA=0.0                                                        14090000
      I18=0                                                             14100000
      ITTP=1                                                            14110000
      ITMAX=185                                                         14120000
      WRITE(6,1910) SPIN(NN(N)),SPIN(IK(J))                             14130000
      WRITE(20,1910) SPIN(NN(N)),SPIN(IK(J))                            14140000
 1910 FORMAT(1H0,1X,F5.2,' -->',F5.2,' TRANSITION :'/)                  14150000
      WRITE(6,1920)                                                     14160000
 1920 FORMAT(1H0,1X,'  TG     VG       DSIGMA       DS'/)               14170000
 1950 IT1=5                                                             14180000
C                                                                       14190000
C DO-LOOP  IT1 = 5,ITMAX,5      NO. 1900                                14200000
C                                                                       14210000
 1900 CONTINUE                                                          14220000
      ITT=IT1                                                           14230000
      IF(VGAMMA.LT.3.0) GOTO 1980                                       14240000
      ITT=180-IT1+5                                                     14250000
 1980 T=PI*FLOAT(ITT-5)/180.                                            14260000
      TT=FLOAT(ITT-5)                                                   14270000
 1970 CONTINUE                                                          14280000
C                                                                       14290000
C                                                                       14300000
C                                                                       14310000
      W(NN(N),IK(J))=0.0                                                14320000
      VORFAC=1./SQRT(4.*3.1415926)                                      14330000
      DO 14 K=1,5,2                                                     14340000
      DO 15 KK=1,K                                                      14350000
      KLAUF=K                                                           14360000
      KKLAUF=KK                                                         14370000
      NHILF=NN(N)                                                       14380000
      IKHILF=IK(J)                                                      14390000
C                                                                       14400000
C SUBROUTINE FOR THE INTEGRATION OVER THE FINAL PARTICLE COUNTER        14410000
C                                                                       14420000
      CALL DQG8(V2,V1,FCT,YYY)                                          14430000
C                                                                       14440000
      IF(I22.LT.2) GOTO 128                                             14450000
C                                                                       14460000
C DEORIENTATION CORRECTION OF THE GAMMA-RAY ANGULAR DISTRIBUTION        14470000
C                                                                       14480000
      W(NN(N),IK(J))=W(NN(N),IK(J))+ANG(NN(N),K,KK)*F(K,IK(J),NN(N))    14490000
     1 *YYY*QG(K,NN(N))                                                 14500000
      GOTO 15                                                           14510000
  128 W(NN(N),IK(J))=W(NN(N),IK(J))+ANG(NN(N),K,KK)*F(K,IK(J),NN(N))    14520000
     1 *YYY                                                             14530000
c      write(6,9999) w(nn(n),ik(j)),ang(nn(n),k,kk),f(k,ik(j),nn(n)),yyy
 9999 format(4e10.3)
   15 CONTINUE                                                          14540000
   14 CONTINUE                                                          14550000
C                                                                       14560000
      W(NN(N),IK(J))=VORFAC*W(NN(N),IK(J))                              14570000
C                                                                       14580000
      IF(I17.EQ.0) GOTO 206                                             14590000
      WRITE(6,1012) SPIN(NN(N)),SPIN(IK(J)),W(NN(N),IK(J))              14600000
 1012 FORMAT(1X,'DW(',F5.2,'-->',F5.2,') = ',E15.7/21X,15(1H*)//)       14610000
  206 CONTINUE                                                          14620000
C                                                                       14630000
C **********************************************************************14640000
C CALCULATION OF THE DOUBLE DIFFERENTIAL CROSS SECTION : DS(LAUF)       14650000
C "COULOMB EXCITATION", PAGE 314                                        14660000
C **********************************************************************14670000
C                                                                       14680000
C 1. CALCULATION OF THE RUTHERFORD CROSS SECTION   (BARNS/STERADIAN)    14690000
C                                                                       14700000
      VF1=((1.+A1/A2)*Z1*Z2/EP)**2                                      14710000
      VF2=SQRT(EP/(EP-(1.+A1/A2)*E(NN(N))))                             14720000
      DSR=1.2958E-3*VF1*VF2*(1./SIN(THET/2.)**4)                        14730000
C                                                                       14740000
C 2. TRANSFORMATION OF THE RUTHERFORD CROSS SECTION FROM THE            14750000
C    CM- INTO THE LAB-SYSTEM (FOR PROJECTILE OR TARGET NUCLEUS)         14760014
C                                                                       14770000
      QQQ=-E(NN(N))                                                     14780000
      ETT=EP+QQQ                                                        14790000
      IF(A1.LT.A2) GOTO 2100                                            14800014
C                                                                       14810000
C A1 > A2                                                               14820000
C                                                                       14830000
      AAA=A1*A1*(EP/ETT)/((A1+A2)**2)                                   14840000
      BBB=A1*A2*(EP/ETT)/((A1+A2)**2)                                   14850000
      CCC=A2*A2*(1.+(A1*QQQ/(A2*ETT)))/((A1+A2)**2)                     14860000
      DDD=A2*A1*(1.+(A1*QQQ/(A2*ETT)))/((A1+A2)**2)                     14870000
C                                                                       14880014
      IF(I25.EQ.0) GOTO 2002                                            14890014
C                                                                       14900014
C PROJECTILE EXCITATION                                                 14910014
C                                                                       14920014
      E4ETT=AAA+CCC+2.*SQRT(AAA*CCC)*COS(THET)                          14930014
      THSIN=SQRT(A2*DDD/(A1*E4ETT))*SIN(PI-THET)                        14940014
      SINLAB(LAUF)=THSIN                                                14950014
      SIGLAB=SQRT(AAA*CCC*((CCC/AAA)-THSIN**2))/E4ETT                   14960014
      SIGLAB=1./SIGLAB                                                  14970018
      GOTO 2200                                                         14980018
C                                                                       14990014
C TARGET EXCITATION                                                     15000014
C                                                                       15010014
 2002 E3ETT=BBB+DDD+2.*SQRT(AAA*CCC)*COS(PI-THET)                       15020014
      THSIN=SQRT(DDD/E3ETT)*SIN(PI-THET)                                15030014
      SINLAB(LAUF)=THSIN                                                15040014
      SIGLAB=SQRT(AAA*CCC*((DDD/BBB)-THSIN**2))/E3ETT                   15050014
      SIGLAB=1./SIGLAB                                                  15060018
      GOTO 2200                                                         15070014
C                                                                       15080000
C A1 < A2                                                               15090000
C                                                                       15100000
 2100 AAA=A1*A2*(EP/ETT)/((A1+A2)**2)                                   15110014
      BBB=A1*A1*(EP/ETT)/((A1+A2)**2)                                   15120000
      CCC=A2*A1*(1.+(A1*QQQ/(A2*ETT)))/((A1+A2)**2)                     15130000
      DDD=A2*A2*(1.+(A1*QQQ/(A2*ETT)))/((A1+A2)**2)                     15140000
C                                                                       15150014
      IF(I25.EQ.0) GOTO 2004                                            15160014
C                                                                       15170014
C PROJECTILE EXCITATION                                                 15180014
C                                                                       15190014
      E3ETT=BBB+DDD+2.*SQRT(AAA*CCC)*COS(THET)                          15200014
      THSIN=SQRT(DDD/E3ETT)*SIN(THET)                                   15210014
      SINLAB(LAUF)=THSIN                                                15220014
      SIGLAB=SQRT(AAA*CCC*((DDD/BBB)-THSIN**2))/E3ETT                   15230014
      SIGLAB=1./SIGLAB                                                  15240018
      GOTO 2200                                                         15250018
C                                                                       15260014
C TARGET EXCITATION                                                     15270014
C                                                                       15280014
 2004 E4ETT=AAA+CCC+2.*SQRT(AAA*CCC)*COS(PI-THET)                       15290014
      THSIN=SQRT(A1*DDD/(A2*E4ETT))*SIN(THET)                           15300014
      SINLAB(LAUF)=THSIN                                                15310014
      SIGLAB=SQRT(AAA*CCC*((CCC/AAA)-THSIN**2))/E4ETT                   15320014
      SIGLAB=1./SIGLAB                                                  15330018
C                                                                       15340000
C                                                                       15350000
 2200 DS(LAUF)=W(NN(N),IK(J))*DSR*SIGLAB                                15360014
C                                                                       15370000
C                                                                       15380000
C                                                                       15390000
C >>> GAMMA-RAY ANGULAR DISTRIBUTION<<<                                 15400000
C                                                                       15410000
      IF(I19.EQ.0) GOTO 1960                                            15420000
      IF(VGAMMA.LT.3.0) GOTO 1903                                       15430000
      TTPP(ITTP)=(360.-TT)*PI/180.                                      15440000
      GOTO 1902                                                         15450000
 1903 TTPP(ITTP)=TT*PI/180.                                             15460000
 1902 CONTINUE                                                          15470000
      TTP(ITTP)=TT                                                      15480000
      VGP(ITTP)=VGAMMA                                                  15490000
      DSP(ITTP)=DS(LAUF)                                                15500000
      ITTP=ITTP+1                                                       15510000
      IT1=IT1+5                                                         15520000
      IF(IT1.GT.ITMAX) GOTO 1999                                        15530000
      GOTO 1900                                                         15540000
C                                                                       15550000
C END OF DO-LOOP IT1 = 5,ITMAX,5      NO. 1900                          15560000
C                                                                       15570000
 1999 IF(VGAMMA.GT.3.0) GOTO 1940                                       15580000
      VGAMMA=3.1415926                                                  15590000
      ITMAX=180                                                         15600000
      GOTO 1950                                                         15610000
 1940 CONTINUE                                                          15620000
      ITTPN=ITTP-1                                                      15630000
      DSPN=DSP(1)                                                       15640000
      DO 1941 IWOLL=2,ITTPN                                             15650000
      IF(DSP(IWOLL).GT.DSPN) DSPN=DSP(IWOLL)                            15660000
 1941 CONTINUE                                                          15670000
      DSPNN=3.5/DSPN                                                    15680000
      DO 1942 IWOLL=1,ITTPN                                             15690000
      DSNORM(IWOLL)=DSP(IWOLL)*DSPNN                                    15700000
      WRITE(6,1930) TTP(IWOLL),VGP(IWOLL),DSP(IWOLL),DSNORM(IWOLL)      15710000
      WRITE(20,1930) TTP(IWOLL),VGP(IWOLL),DSP(IWOLL),DSNORM(IWOLL)     15720000
 1930 FORMAT(2(1X,F6.2),1X,E15.7,F6.2)                                  15730000
 1942 CONTINUE                                                          15740000
C ----------------------------------------------------------------------15750014
C                                                                       15760000
C >>> GAMMA-RAY ANGULAR DISTRIBUTION PLOTT-OUTPUT <<<                   15770000
C                                                                       15780000
      IF(I20.EQ.0) GOTO 1960                                            15790000
C                                                                       15800000
c      IF(IPLOTT.GT.0) GOTO 1915                                        15810000
c      CALL GRAPHC('P')                                                 15820000
c      CALL PLOTS(14.,72.)                                              15830000
c 1915 IF(IPLOTT.EQ.1) GOTO 1911                                        15840000
c      IF(IPLOTT.EQ.2) GOTO 1912                                        15850000
c      IF(IPLOTT.EQ.3) GOTO 1916                                        15860000
c      IF(IPLOTT.EQ.4) GOTO 1917                                        15870000
c      IPLOTT=1                                                         15880000
c      CALL PLOT(7.,7.,-3)                                              15890000
c      GOTO 1913                                                        15900000
c 1911 IPLOTT=2                                                         15910000
c      CALL PLOT(0.,14.,-3)                                             15920000
c      GOTO 1913                                                        15930000
c 1912 IPLOTT=3                                                         15940000
c      CALL PLOT(0.,14.,-3)                                             15950000
c      GOTO 1913                                                        15960000
c 1916 IPLOTT=4                                                         15970000
c      CALL PLOT(0.,14.,-3)                                             15980000
c      GOTO 1913                                                        15990000
c 1917 IPLOTT=0                                                         16000000
c      CALL PLOT(0.,14.,-3)                                             16010000
c 1913 CALL POLAR(DSNORM,TTPP,73,2,0,2,-3.5,1.)                         16020000
c      CALL SYMBOL(-6.5,0.,0.75,3,0.,-1)                                16030000
c      CALL SYMBOL(0.,-6.5,0.75,3,0.,-1)                                16040000
c      CALL SYMBOL(0.,6.5,0.75,3,0.,-1)                                 16050000
c      CALL SYMBOL(6.5,0.,0.75,3,0.,-1)                                 16060000
c      INDEX=IFIX(SPIN(NN(N)))+1                                        16070000
c      CALL SYMBOL(-6.25,4.75,0.5,TAB(INDEX),0.,2)                      16080000
c      CALL SYMBOL(-4.5,5.,0.5,24,0.,-1)                                16090000
c      INDEX=IFIX(SPIN(IK(J)))+1                                        16100000
c      CALL SYMBOL(-3.75,4.75,0.5,TAB1(INDEX),0.,2)                     16110000
c      IF(IPLOTT.GT.0) GOTO 1914                                        16120000
c      CALL PLTEND                                                      16130000
 1914 CONTINUE                                                          16140000
 1960 CONTINUE                                                          16150000
C ----------------------------------------------------------------------16160014
C                                                                       16170000
C                                                                       16180000
C                                                                       16190000
C                                                                       16200000
      IF(I18.EQ.0) GOTO 207                                             16210000
      WRITE(6,1013) SPIN(NN(N)),SPIN(IK(J)),DS(LAUF)                    16220000
 1013 FORMAT(1X,'DS(',F5.2,'-->',F5.2,') = ',E15.7/21X,15(1H*)//)       16230000
  207 CONTINUE                                                          16240000
      IITRAN(LAUF)=NN(N)                                                16250000
      IFTRAN(LAUF)=IK(J)                                                16260000
      LAUF=LAUF+1                                                       16270000
      IF(LAUF.GT.199) GOTO 533                                          16280000
C                                                                       16290000
    9 CONTINUE                                                          16300000
C                                                                       16310000
    6 CONTINUE                                                          16320000
C                                                                       16330000
      GOTO 535                                                          16340000
  533 WRITE(6,534)                                                      16350000
  534 FORMAT(1H0,72(1H*)/1X,2(1H*),14X,'MAX. NUMBER OF POSSIBLE TRANSITI16360000
     1ONS > 199',14X,2(1H*)/1X,72(1H*))                                 16370000
      STOP                                                              16380000
c  535 IF(I20.EQ.1) CALL PLTEND                                         16390000
  535 continue
      RETURN                                                            16400000
      END                                                               16410000
      FUNCTION FK(AJF,AL1,AL2,AJI,AK)                                   16420000
C                                                                       16430000
C  FK(JF,L1,L2,JI)                                                      16440000
C                                                                       16450000
      FK=((-1.)**(IFIX(AJF-AJI-1.00000)))*SQRT((2.*AL1+1.)*             16460000
     1 (2.*AL2+1.)*(2.*AJI+1.))*CLEBI(AL1,AL2,AK,1.,-1.,0.)             16470000
     2 *RACAHI(AJI,AJI,AL1,AL2,AK,AJF)                                  16480000
      RETURN                                                            16490000
      END                                                               16500000
      FUNCTION FCT(X)                                                   16510000
C                                                                       16520014
C                                                                       16530014
C                                                                       16540014
      REAL ME                                                           16550000
      COMMON/C1/ E(35),SPIN(35),EE(35),ME2(35,35),MM1(35,35),NMAX       16560000
      COMMON/C2/ ALPHA1(35,35),ALPHA2(35,35),ANG(35,5,5)                16570000
      COMMON/C12/ ALPHE1(35,35),ALPHE3(35,35),DELTA1(199)               16580000
      COMMON/C3/ I11,I12,I13,I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I2416590000
     1,I25                                                              16600014
      COMMON/C4/ Z1,Z2,A1,A2,EP,THETAP,QG(5,35),QGS(5,35),IPAR(35)      16610000
      COMMON/C14/ AJ,GZA,TAUC,H0(35)                                    16620000
      COMMON/C7/ T,RESTLA,KLAUF,KKLAUF,NHILF,IKHILF                     16630000
      COMMON/C10/ VGAMMA,Q(5),ME(35,35,4),IITRAN(199),IFTRAN(199)       16640000
C **********************************************************************16650000
C GAMMA-DETECTOR IN THE XZ-PLANE                                        16660000
C                                                                       16670000
C EXPERIMENT GSI : THETA-GAMMA = 150. , VGAMMA = 0.                     16680000
C                  THETA-GAMMA = 30.  , VGAMMA = 180.                   16690000
C **********************************************************************16700000
C                                                                       16710000
C 1. TRANSFORMATION OF THE ANGLES FROM THE LAB- INTO THE REST-SYSTEM    16720000
C "NUCL. PHYS. A192 (1972), 449"                                        16730000
C                                                                       16740000
      PI=3.1415926                                                      16750000
      THET=THETAP*PI/180.                                               16760000
      VCM=0.04634*SQRT(EP/A1)/(1.+A2/A1)                                16770000
      DELTA=E(NHILF)                                                    16780000
      AKA=(1.+A1/A2)*DELTA/EP                                           16790000
      AKAPA2=SQRT(1.-AKA)                                               16800000
      AKAPA1=(A2/A1)*SQRT(1.-AKA)                                       16810000
C                                                                       16820014
      IF(I25.EQ.0) GOTO 2000                                            16830014
C                                                                       16840014
C PROJECTILE EXCITATION                                                 16850014
C                                                                       16860014
      VV2=VCM*SQRT(1.+(AKAPA1**2)+2.*AKAPA1*COS(THET))                  16870025
      T2COS=(VCM/VV2)*(1.+AKAPA1*COS(THET))                             16880025
      T2SIN=SQRT(1.-T2COS*T2COS)                                        16890014
      TG2COS=COS(T)*T2COS+T2SIN*SIN(T)*COS(VGAMMA-X)                    16900014
C                                                                       16910014
C CALCULATION OF ANGLE FOR THE GAMMA-DETECTOR IN THE REST-SYSTEM        16920014
C                                                                       16930014
      TGCOS=(COS(T)*SQRT(1.-VV2*VV2)-T2COS*VV2+T2COS*TG2COS*(1.-SQRT(   16940025
     1 1.-VV2*VV2)))/(1.-VV2*TG2COS)                                    16950025
      TGSIN=SQRT(1.-TGCOS*TGCOS)                                        16960014
      VICOS=(+1.)*(TG2COS-T2COS*TGCOS-VV2*(1.-T2COS*TGCOS*TG2COS))      16970025
     1 /(TGSIN*T2SIN*(1.-VV2*TG2COS))                                   16980025
      GOTO 2001                                                         16990014
C                                                                       17000014
C TARGET EXCITATION                                                     17010014
C                                                                       17020014
 2000 VV2=VCM*SQRT(1.+(AKAPA2**2)-2.*AKAPA2*COS(THET))                  17030025
      T2COS=(VCM/VV2)*(1.-AKAPA2*COS(THET))                             17040025
      T2SIN=SQRT(1.-T2COS*T2COS)                                        17050014
      TG2COS=COS(T)*T2COS-T2SIN*SIN(T)*COS(VGAMMA-X)                    17060014
C                                                                       17070000
C CALCULATION OF ANGLE FOR THE GAMMA-DETECTOR IN THE REST-SYSTEM        17080000
C                                                                       17090000
      TGCOS=(COS(T)*SQRT(1.-VV2*VV2)-T2COS*VV2+T2COS*TG2COS*(1.-SQRT(   17100025
     1 1.-VV2*VV2)))/(1.-VV2*TG2COS)                                    17110025
      TGSIN=SQRT(1.-TGCOS*TGCOS)                                        17120000
      VICOS=(-1.)*(TG2COS-T2COS*TGCOS-VV2*(1.-T2COS*TGCOS*TG2COS))      17130025
     1 /(TGSIN*T2SIN*(1.-VV2*TG2COS))                                   17140025
C                                                                       17150000
C 2. TRANSFORMATION OF THE CROSS SECTION FROM THE REST- INTO            17160000
C THE LAB-SYSTEM                                                        17170000
C                                                                       17180000
 2001 RESTLA=ABS((1.-VV2*VV2)/((1.-VV2*TG2COS)**2))                     17190025
C                                                                       17200000
C RESTLA = FUNCTION OF PARTICLE ANGLE (AZIMUTHAL-ANGLE = X )            17210000
C                                                                       17220000
C CALCULATION OF THE SPHERICAL HARMONIC FUNCTION :  Y(K,KK)             17230000
C K=0,2,4  -K <KK< +K                                                   17240000
C INCLUDING THE TRANSFORMATION OF THE CROSS-SECTION FROM THE REST-      17250000
C INTO THE LAB-SYSTEM                                                   17260000
C                                                                       17270000
      IF(I24.GT.0) RESTLA=1.                                            17280000
      IF(I24.GT.1) GOTO 12                                              17290000
      GOTO 11                                                           17300000
   12 TGCOS=COS(T)                                                      17310000
      VICOS=COS(VGAMMA-X)                                               17320000
   11 CONTINUE                                                          17330000
C                                                                       17340000
C I24=CALCULATION OF THE GAMMA-RAY ANGULAR DISTRIBUTION IN A COORDINATE 17350000
C     SYSTEM FOR WHICH THE EXCITED NUCLEUS IS AT REST                   17360000
C >>> I24=1 <<<                                                         17370000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  17380000
C     OF THE GAMMA - DETECTOR IN THE LABORATORY COORDINATE SYSTEM       17390000
C >>> I24>1 <<<                                                         17400000
C     INPUT OF THE POLAR ANGLE ( TT ) AND THE AZIMUTHAL ANGLE (VGAMMA)  17410000
C     OF THE GAMMA - DETECTOR IN THE REST-COORDINATE SYSTEM             17420000
C                                                                       17430000
      XSAVE=X                                                           17440000
      VI2COS=2.*VICOS*VICOS-1.                                          17450000
      VI3COS=4.*(VICOS**3)-3.*VICOS                                     17460000
      VI4COS=8.*(VICOS**4)-8.*(VICOS**2)+1.                             17470000
      K=KLAUF                                                           17480000
      KK=KKLAUF                                                         17490000
      IF(KK.EQ.1.AND.K.EQ.1) GOTO 1                                     17500000
      IF(KK.EQ.1.AND.K.EQ.3) GOTO 2                                     17510000
      IF(KK.EQ.2.AND.K.EQ.3) GOTO 3                                     17520000
      IF(KK.EQ.3.AND.K.EQ.3) GOTO 4                                     17530000
      IF(KK.EQ.1.AND.K.EQ.5) GOTO 5                                     17540000
      IF(KK.EQ.2.AND.K.EQ.5) GOTO 6                                     17550000
      IF(KK.EQ.3.AND.K.EQ.5) GOTO 7                                     17560000
      IF(KK.EQ.4.AND.K.EQ.5) GOTO 8                                     17570000
      IF(KK.EQ.5.AND.K.EQ.5) GOTO 9                                     17580000
C                                                                       17590000
C Y(0,0)                                                                17600000
C                                                                       17610000
    1 FCT=1./SQRT(4.*PI)                                                17620000
      FCT=FCT*RESTLA                                                    17630000
      GOTO 10                                                           17640000
C                                                                       17650000
C Y(2,0)                                                                17660000
C                                                                       17670000
    2 FCT=SQRT(5./(16.*PI))*(3.*TGCOS**2-1.)                            17680000
      FCT=FCT*RESTLA                                                    17690000
      GOTO 10                                                           17700000
C                                                                       17710000
C Y(2,-1) + Y(2,1)                                                      17720000
C                                                                       17730000
    3 FCT=-SQRT(15./(2.*PI))*TGCOS*TGSIN*VICOS                          17740000
      FCT=FCT*RESTLA                                                    17750000
      GOTO 10                                                           17760000
C                                                                       17770000
C Y(2,-2) + Y(2,2)                                                      17780000
C                                                                       17790000
    4 FCT=SQRT(15./(8.*PI))*TGSIN**2*VI2COS                             17800000
      FCT=FCT*RESTLA                                                    17810000
      GOTO 10                                                           17820000
C                                                                       17830000
C Y(4,0)                                                                17840000
C                                                                       17850000
    5 FCT=SQRT(9./(256.*PI))*(35.*TGCOS**4-30.*TGCOS**2+3.)             17860000
      FCT=FCT*RESTLA                                                    17870000
      GOTO 10                                                           17880000
C                                                                       17890000
C Y(4,-1) + Y(4,1)                                                      17900000
C                                                                       17910000
    6 FCT=-SQRT(45./(16.*PI))*(7.*TGCOS**2-3.)*TGSIN*TGCOS*             17920000
     1VICOS                                                             17930000
      FCT=FCT*RESTLA                                                    17940000
      GOTO 10                                                           17950000
C                                                                       17960000
C Y(4,-2) + Y(4,2)                                                      17970000
C                                                                       17980000
    7 FCT=SQRT(45./(32.*PI))*(7.*TGCOS**2-1.)*TGSIN**2*                 17990000
     1VI2COS                                                            18000000
      FCT=FCT*RESTLA                                                    18010000
      GOTO 10                                                           18020000
C                                                                       18030000
C Y(4,-3) + Y(4,3)                                                      18040000
C                                                                       18050000
    8 FCT=-SQRT(315./(16.*PI))*TGSIN**3*TGCOS*VI3COS                    18060000
      FCT=FCT*RESTLA                                                    18070000
      GOTO 10                                                           18080000
C                                                                       18090000
C Y(4,-4) + Y(4,4)                                                      18100000
C                                                                       18110000
    9 FCT=SQRT(315./(128.*PI))*TGSIN**4*VI4COS                          18120000
      FCT=FCT*RESTLA                                                    18130000
   10 X=XSAVE                                                           18140000
      RETURN                                                            18150000
      END                                                               18160000
      FUNCTION CLEBI(AJ1,AJ2,AJ3,AL1,AL2,AL3)                           18170000
C                                                                       18180000
C CLEB(J1,J2,J3,L1;L2,L3)                                               18190000
C                                                                       18200000
      DIMENSION FCT(179),FACT(179)                                      18210000
      NA=IFIX(2.*AJ1)                                                   18220000
      NB=IFIX(2.*AJ2)                                                   18230000
      NC=IFIX(2.*AJ3)                                                   18240000
      ND=IFIX(2.*AL1)                                                   18250000
      NE=IFIX(2.*AL2)                                                   18260000
      NFF=IFIX(2.*AL3)                                                  18270000
C                                                                       18280000
C     INITIALISE FACTORIAL ARRAY FCT                                    18290000
C                                                                       18300000
      FCTZ=1.                                                           18310000
      FACTZ=1.                                                          18320000
      FCT(1)=ALOG(1.0)                                                  18330000
      FACT(1)=ALOG(1.0)                                                 18340000
      DO 101 I=1,178                                                    18350000
      FACT(I+1)=FACT(I)+ALOG(SQRT(FLOAT(I)))                            18360000
  101 FCT(I+1)=ALOG(FLOAT(I))-ALOG(10.)+FCT(I)                          18370000
C                                                                       18380000
C                                                                       18390000
      NF=-NFF                                                           18400000
      CLEBI=0.                                                          18410000
      IF(ND+NE+NF .NE. 0) RETURN                                        18420000
      L3=NB+NE                                                          18430000
      L8=NA+ND                                                          18440000
      L10=NC+NF                                                         18450000
      IF(MOD(L3,2) .NE. 0) RETURN                                       18460000
      IF(MOD(L8,2) .NE. 0) RETURN                                       18470000
      IF(MOD(L10,2) .NE. 0) RETURN                                      18480000
      L3=L3/2                                                           18490000
      L8=L8/2                                                           18500000
      L10=L10/2                                                         18510000
      L0=L3+L8+L10                                                      18520000
      IF(L0 .GT.177) GO TO 2                                            18530000
      IF(IABS(ND) .GT. NA  .OR.                                         18540000
     X   IABS(NE) .GT. NB  .OR.                                         18550000
     X     IABS(NF) .GT. NC) RETURN                                     18560000
      L1=L0-NC                                                          18570000
      L2=L8-ND                                                          18580000
      L9=L3-NE                                                          18590000
      L11=L10-NF                                                        18600000
      L6=L0-NB                                                          18610000
      L7=L0-NA                                                          18620000
      L4=L6-L2                                                          18630000
      L5=L7-L3                                                          18640000
      KMIN=MAX0(0,-L4,-L5)                                              18650000
      KMAX=MIN0(L1,L2,L3)                                               18660000
      IF(KMAX .LT. KMIN)RETURN                                          18670000
C     CO1=FACT(L1+1)*FACT(L6+1)*FACT(L7+1)*FACT(L8+1)                   18680000
C    X    *FACT(L2+1)*FACT(L3+1)                                        18690000
C    X*FACT(L9+1)*FACT(L10+1)*FACT(L11+1)/FACT(L0+2)                    18700000
C    X   *FACT(NC+2)/FACT(NC+1)                                         18710000
      CO1=FACT(L1+1)+FACT(L6+1)+FACT(L7+1)+FACT(L8+1)                   18720000
     X    +FACT(L2+1)+FACT(L3+1)                                        18730000
     X+FACT(L9+1)+FACT(L10+1)+FACT(L11+1)-FACT(L0+2)                    18740000
     X   +FACT(NC+2)-FACT(NC+1)                                         18750000
      SGN=1.                                                            18760000
      DO 1 K=KMIN,KMAX                                                  18770000
      L12=L1-K+1                                                        18780000
      L13=L2-K+1                                                        18790000
      L14=L3-K+1                                                        18800000
      L15=L4+K+1                                                        18810000
      L16=L5+K+1                                                        18820000
C     DEN=FACT(K+1)*FACT(L12)*FACT(L13)*FACT(L14)*FACT(L15)*FACT(L16)   18830000
      DEN=FACT(K+1)+FACT(L12)+FACT(L13)+FACT(L14)+FACT(L15)+FACT(L16)   18840000
C     CLEBI=CLEBI+SGN*CO1/DEN**2                                        18850000
      CLEBI=CLEBI+SGN*EXP(CO1-2.*DEN)                                   18860000
      SGN=-SGN                                                          18870000
    1 CONTINUE                                                          18880000
      IF(MOD(KMIN,2) .NE. 0) CLEBI=-CLEBI                               18890000
      RETURN                                                            18900000
    2 CLEBI=0.                                                          18910000
      STOP 1                                                            18920000
      END                                                               18930000
      FUNCTION RACAHI(AJ1,AJ2,AL2,AL1,AJ3,AL3)                          18940000
C                                                                       18950000
C    W(J1,J2,L2,L1;J3,L3)                                               18960000
C                                                                       18970000
      DIMENSION FCT(179),FACT(179)                                      18980000
      NA=IFIX(2.*AJ1)                                                   18990000
      NB=IFIX(2.*AJ2)                                                   19000000
      NCC=IFIX(2.*AL2)                                                  19010000
      ND=IFIX(2.*AL1)                                                   19020000
      NEE=IFIX(2.*AJ3)                                                  19030000
      NF=IFIX(2.*AL3)                                                   19040000
C                                                                       19050000
C     INITIALISE FACTORIAL ARRAY FCT                                    19060000
C                                                                       19070000
      FCTZ=1.                                                           19080000
      FACTZ=1.                                                          19090000
      FCT(1)=ALOG(1.0)                                                  19100000
      FACT(1)=ALOG(1.0)                                                 19110000
      DO 101 I=1,178                                                    19120000
      FACT(I+1)=FACT(I)+ALOG(SQRT(FLOAT(I)))                            19130000
  101 FCT(I+1)=ALOG(FLOAT(I))-ALOG(10.)+FCT(I)                          19140000
C                                                                       19150000
C                                                                       19160000
      NC=NEE                                                            19170000
      NE=NCC                                                            19180000
      RACAHI=0.                                                         19190000
      N1=NA+NB+NC                                                       19200000
      N2=ND+NE+NC                                                       19210000
      N3=NA+NE+NF                                                       19220000
      N4=ND+NB+NF                                                       19230000
      IF(MOD(N1,2) .NE. 0) RETURN                                       19240000
      IF(MOD(N2,2) .NE. 0) RETURN                                       19250000
      IF(MOD(N4,2) .NE. 0) RETURN                                       19260000
      N1=N1/2                                                           19270000
      N2=N2/2                                                           19280000
      N3=N3/2                                                           19290000
      N4=N4/2                                                           19300000
      N1A=N1-NC                                                         19310000
      N2A=N2-NC                                                         19320000
      N3A=N3-NF                                                         19330000
      N4A=N4-NF                                                         19340000
      L6=NE-N2A-N3A                                                     19350000
      L7=NA-N1A-N3A                                                     19360000
      L1=N3A+N4A+NF+1                                                   19370000
      KMIN=MAX0(0,-L6,-L7)                                              19380000
      KMAX=MIN0(L1,N1A,N2A,N3A,N4A)                                     19390000
      IF(KMAX .LT. KMIN) RETURN                                         19400000
      IF(L1 .GT. 178) GO TO 2                                           19410000
      N1B=N1-NB                                                         19420000
      N1C=N1-NA                                                         19430000
      N2B=N2-NE                                                         19440000
      N2C=N2-ND                                                         19450000
      N3B=N3-NE                                                         19460000
      N3C=N3-NA                                                         19470000
      N4B=N4-NB                                                         19480000
      N4C=N4-ND                                                         19490000
      SGN=1.                                                            19500000
      DO 3 K=KMIN,KMAX                                                  19510000
      LL1=L1-K                                                          19520000
      LL2=N1A-K                                                         19530000
      LL3=N2A-K                                                         19540000
      LL4=N3A-K                                                         19550000
      LL5=N4A-K                                                         19560000
      LL6=L6+K                                                          19570000
      LL7=L7+K                                                          19580000
C     TERM=FACT(LL1+1)/(FACT(K+1)*FACT(LL2+1)                           19590000
C    X   *FACT(LL3+1)*FACT(LL4+1)                                       19600000
C    X   *FACT(LL5+1)*FACT(LL6+1)                                       19610000
C    X  *FACT(LL7+1))                                                   19620000
      TERM=FACT(LL1+1)-(FACT(K+1)+FACT(LL2+1)                           19630000
     X   +FACT(LL3+1)+FACT(LL4+1)                                       19640000
     X   +FACT(LL5+1)+FACT(LL6+1)                                       19650000
     X  +FACT(LL7+1))                                                   19660000
      WOLLE2=                                                           19670000
     X FACT(N1A+1)+FACT(N1B+1)+FACT(N1C+1)-FACT(N1+2)                   19680000
     X +FACT(N2A+1)+FACT(N2B+1)+FACT(N2C+1)-FACT(N2+2)                  19690000
     X +FACT(N3A+1)+FACT(N3B+1)+FACT(N3C+1)-FACT(N3+2)                  19700000
     X +FACT(N4A+1)+FACT(N4B+1)+FACT(N4C+1)-FACT(N4+2)                  19710000
C     RACAHI=RACAHI+SGN*TERM**2                                         19720000
      RACAHI=RACAHI+SGN*EXP(WOLLE2+2.*TERM)                             19730000
      SGN=-SGN                                                          19740000
    3 CONTINUE                                                          19750000
C     RACAHI=RACAHI*                                                    19760000
C    X FACT(N1A+1)*FACT(N1B+1)*FACT(N1C+1)/FACT(N1+2)                   19770000
C    X *FACT(N2A+1)*FACT(N2B+1)*FACT(N2C+1)/FACT(N2+2)                  19780000
C    X *FACT(N3A+1)*FACT(N3B+1)*FACT(N3C+1)/FACT(N3+2)                  19790000
C    X *FACT(N4A+1)*FACT(N4B+1)*FACT(N4C+1)/FACT(N4+2)                  19800000
      IF(MOD(KMIN,2) .NE. 0) RACAHI=-RACAHI                             19810000
      RETURN                                                            19820000
    2 RACAHI=0.                                                         19830000
      STOP 2                                                            19840000
      END                                                               19850000
      SUBROUTINE FIT(E1,E2,NS,L,N,ALPHA,NP,EE)                          19860000
      COMMON/C11/A(23),B(4),C(23),D(4,23),E(23),F(23),H(23),S(23),T(23),19870000
     1X(23),Y(23),ST(23)                                                19880000
C     COMMON/C22/E1,E2,NS,L,N,ALPHA(23),NP,EE(23)                       19890000
      DIMENSION ALPHA(23),EE(23)                                        19900000
      DIMENSION ERG(23,3),EN(23)                                        19910000
      DATA ERG/1.0,1.7,3.0,6.2,10.,14.,20.,28.,40.,53.,70.,83.,100.,123.19920000
     1,150.,215.,300.,390.,500.,730.,1000.,1250.,1500.,1.,1.7,3.0,6.2,1019930000
     2.,14.,20.,28.,40.,53.,70.,83.,100.,123.,150.,215.,300.,390.,500.,719940000
     330.,1000.,1500.,0.,1.,2.,4.,8.,15.,25.,40.,52.,70.,103.,150.,280.,19950000
     4500.,10*0./                                                       19960000
      ES=E1                                                             19970000
      IF(NS.EQ.1) E2=E1                                                 19980000
      E1=-E1+1.                                                         19990000
      E2=-E2+1.                                                         20000000
      NE2=-E2+.501                                                      20010000
      NN=N                                                              20020000
      IF(ES.LT.0.) NN=N+5                                               20030000
      DO 15 I=6,NN                                                      20040000
      M=ERG(I,NS)+.01                                                   20050000
      J=M+NE2                                                           20060000
      IF(M.LT.150) GOTO 13                                              20070000
      IF(M.GE.390) GOTO 11                                              20080000
      J=5*((10*J+25)/50)                                                20090000
      GOTO 14                                                           20100000
   11 IF(M.GE.1000) GOTO 12                                             20110000
      J=10*((J+5)/10)                                                   20120000
      GOTO 14                                                           20130000
   12 J=50*((J+25)/50)                                                  20140000
      GOTO 14                                                           20150000
   13 J=(10*J+5)/10                                                     20160000
   14 EN(I)=J                                                           20170000
   15 C(I)=EN(I)+E1                                                     20180000
      IF(ES.LT.0.) GOTO 30                                              20190000
      DO 16 I=1,5                                                       20200000
      C(I)=ERG(I,NS)                                                    20210000
   16 EN(I)=ERG(I,NS)-E1                                                20220000
      DO 17 I=1,N                                                       20230000
   17 C(I)=ALOG(SQRT((1021.952+C(I))*C(I)))                             20240000
      DO 18 I=1,NP                                                      20250000
      W=EE(I)+E1                                                        20260000
   18 X(I)=ALOG(SQRT((1021.952+W)*W))                                   20270000
      M=2*L+1                                                           20280000
      DO 19 I=1,N                                                       20290000
   19 A(I)=ALOG(ALPHA(I)*EN(I)**M)                                      20300000
      CALL SPLINE(N,NP)                                                 20310000
      DO 20 I=1,NP                                                      20320000
   20 ALPHA(I)=EXP(Y(I))/EE(I)**M                                       20330000
      GOTO 33                                                           20340000
   30 DO 31 I=1,N                                                       20350000
      C(I)=ALOG(EN(I+5))                                                20360000
   31 A(I)=ALOG(ALPHA(I))                                               20370000
      DO 35 I=1,NP                                                      20380000
   35 X(I)=ALOG(EE(I))                                                  20390000
      CALL SPLINE(N,NP)                                                 20400000
      DO 32 I=1,NP                                                      20410000
   32 ALPHA(I)=EXP(Y(I))                                                20420000
   33 CONTINUE                                                          20430000
      L1=0                                                              20440000
      IF(ES.LT.0.) L1=5                                                 20450000
      M1=N+L1                                                           20460000
      DO 37 I=1,NP                                                      20470000
      IF(EE(I).LT.EN(L1+1)) ALPHA(I)=-1.                                20480000
   37 IF(EE(I).GT.EN(M1)) ALPHA(I)=-1.                                  20490000
      RETURN                                                            20500000
      END                                                               20510000
      SUBROUTINE SPLINE(N,NP)                                           20520000
      COMMON/C11/A(23),B(4),C(23),D(4,23),E(23),F(23),H(23),S(23),T(23),20530000
     1X(23),Y(23),ST(23)                                                20540000
      M=N-1                                                             20550000
      DO 7 I=1,M                                                        20560000
      E(I)=C(I+1)-C(I)                                                  20570000
    7 F(I)=A(I+1)-A(I)                                                  20580000
      MS=N-2                                                            20590000
      DO 8 I=1,MS                                                       20600000
      D(1,I)=E(I+1)                                                     20610000
      D(2,I)=2.*(E(I)+E(I+1))                                           20620000
      D(3,I)=E(I)                                                       20630000
    8 D(4,I)=3.*(E(I+1)*F(I)/E(I)+E(I)*F(I+1)/E(I+1))                   20640000
      S(1)=-1.                                                          20650000
      T(1)=2.*F(1)/E(1)                                                 20660000
      DO 9 I=2,M                                                        20670000
      L=I-1                                                             20680000
      S(I)=-(D(1,L)+D(2,L)*S(L))/(D(3,L)*S(L))                          20690000
    9 T(I)=(-T(L)*(D(2,L)+D(3,L)*S(I))+D(4,L))/D(3,L)                   20700000
      AA=2.*F(M)/E(M)                                                   20710000
      H(N)=(T(M)+AA*S(M))/(S(M)+1.)                                     20720000
      DO 10 I=1,M                                                       20730000
      K=N-I                                                             20740000
   10 H(K)=(H(K+1)-T(K))/S(K)                                           20750000
      DO 12 I=1,NP                                                      20760000
      DO 13 K=1,M                                                       20770000
      IF(X(I).LT.C(K+1)) GOTO 14                                        20780000
   13 CONTINUE                                                          20790000
   14 K1=K                                                              20800000
      Z=(C(K1+1)-C(K1))/2.                                              20810000
      B(3)=(H(K1+1)-H(K1))/(4.*Z)                                       20820000
      B(4)=.25*((H(K1+1)+H(K1))*Z-A(K1+1)+A(K1))/Z**3                   20830000
      B(1)=(A(K1+1)+A(K1)-2.*B(3)*Z**2)/2.                              20840000
      B(2)=(H(K1+1)+H(K1)-6.*B(4)*Z**2)/2.                              20850000
      ZP=X(I)-(C(K1+1)+C(K1))/2.                                        20860000
      Y(I)=B(1)+ZP*(B(2)+ZP*(B(3)+ZP*B(4)))                             20870000
   12 CONTINUE                                                          20880000
      RETURN                                                            20890000
      END                                                               20900000
      SUBROUTINE DQTFG(X,Y,Z,NDIM)                                      20910000
C                                                                       20920000
C     ..................................................................20930000
C                                                                       20940000
C     PURPOSE                                                           20950000
C        TO COMPUTE THE VECTOR OF INTEGRAL VALUES FOR A GIVEN           20960000
C        GENERAL TABLE OF ARGUMENT AND FUNCTION VALUES.                 20970000
C                                                                       20980000
C     USAGE                                                             20990000
C        CALL DQTFG(X,Y,Z,NDIM)                                         21000000
C                                                                       21010000
C     DESCRIPTION OF PARAMETERS                                         21020000
C        X      - DOUBLE PRECISION INPUT VECTOR OF ARGUMENT VALUES.     21030000
C        Y      - DOUBLE PRECISION INPUT VECTOR OF FUNCTION VALUES.     21040000
C        Z      - RESULTING DOUBLE PRECISION VECTOR OF INTEGRAL         21050000
C                 VALUES. Z MAY BE IDENTICAL WITH X OR Y.               21060000
C        NDIM   - THE DIMENSION OF VECTORS X,Y,Z.                       21070000
C                                                                       21080000
C     REMARKS                                                           21090000
C        NO ACTION IN CASE NDIM LESS THAN 1.                            21100000
C                                                                       21110000
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     21120000
C        NONE                                                           21130000
C                                                                       21140000
C     METHOD                                                            21150000
C        BEGINNING WITH Z(1)=0. EVALUATION OF VECTOR Z IS DONE BY       21160000
C        MEANS OF TRAPEZOIDAL RULE (SECOND ORDER FORMULA).              21170000
C        FOR REFERENCE, SEE                                             21180000
C        F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,            21190000
C        MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.75.             21200000
C                                                                       21210000
C     ..................................................................21220000
C                                                                       21230000
      DIMENSION X(1),Y(1),Z(1)                                          21240000
c     DOUBLE PRECISION X,Y,Z,SUM1,SUM2                                  21250000
C                                                                       21260000
      SUM2=0.E0                                                         21270000
      IF(NDIM-1)4,3,1                                                   21280000
C                                                                       21290000
C     INTEGRATION LOOP                                                  21300000
    1 DO 2 I=2,NDIM                                                     21310000
      SUM1=SUM2                                                         21320000
      SUM2=SUM2+.5E0*(X(I)-X(I-1))*(Y(I)+Y(I-1))                        21330000
    2 Z(I-1)=SUM1                                                       21340000
    3 Z(NDIM)=SUM2                                                      21350000
    4 RETURN                                                            21360000
      END                                                               21370000
      SUBROUTINE DQG8(XL,XU,FCT,Y)                                      21380000
C                                                                       21390000
C     ..................................................................21400000
C                                                                       21410000
C     SUBROUTINE DQG8                                                   21420000
C                                                                       21430000
C     PURPOSE                                                           21440000
C        TO COMPUTE INTEGRAL(FCT(X), SUMNED OVER X FROM XL TO XU)       21450000
C                                                                       21460000
C     USAGE                                                             21470000
C        CALL DQG8 (XL,XU,FCT,Y)                                        21480000
C        PARAMETER FCT REQUIRES AN EXTERAL STATEMENT                    21490000
C                                                                       21500000
C     DESCRIPTION O PARAMETERS                                          21510000
C        XL     - DOUBLE PRECISION LOWER BOUND OF THE INTERVAL.         21520000
C        XU     - DOUBLE PRECISION UPPER BOUND OF THE INTERVAL.         21530000
C        FCT    - THE NAME OF AN EXTERNAL DOUBLE PRECISION FUNCTION     21540000
C                 SUBPROGRAM USED.                                      21550000
C        Y      - THE RESULTING DOUBLE PRECISION INTEGRAL VALUE.        21560000
C                                                                       21570000
C     REMARKS                                                           21580000
C        NONE                                                           21590000
C                                                                       21600000
C     SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                     21610000
C        THE EXTERNAL DOUBLE PRECISION FUNCTION SUBPROGRAM FCT(X)       21620000
C        MUST BE FURNISHED BY THE USER.                                 21630000
C                                                                       21640000
C     METHOD                                                            21650000
C        EVALUATION IS DONE BY MEANS OF 8-POINT GAUSS QUADRATURE        21660000
C        FORMULA, WHICH INTEGRATES POLYNOMIALS UP TO DEGREE 15          21670000
C        EXACTLY. FOR REFERENCE, SEE                                    21680000
C        V.I.KRYLOV, APPROXIMATE CALCULATION OF INTEGRALS,              21690000
C        MACMILLAN, NEW YORK/LONDON, 1962, PP.100-111 AND 337-340.      21700000
C                                                                       21710000
C     ..................................................................21720000
C                                                                       21730000
c     DOUBLE PRECISION XL,XU,Y,A,B,C,FCT                                21740000
C                                                                       21750000
      A=.5e0*(XU+XL)                                                    21760000
      B=XU-XL                                                           21770000
      C=.48014492824876812e0*B                                          21780000
      Y=.50614268145188130e-1*(FCT(A+C)+FCT(A-C))                       21790000
      C=.39833323870681337e0*B                                          21800000
      Y=Y+.11119051722668724e0*(FCT(A+C)+FCT(A-C))                      21810000
      C=.26276620495816449e0*B                                          21820000
      Y=Y+.15685332293894364e0*(FCT(A+C)+FCT(A-C))                      21830000
      C=.9171732124782490e-1*B                                          21840000
      Y=B*(Y+.18134189168918099e0*(FCT(A+C)+FCT(A-C)))                  21850000
      RETURN                                                            21860000
      END                                                               21870000
