      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      DIMENSION TITLE(10)                                                       
      DIMENSION IDUM(52),DUM(52,12)                                             
      DATA END/'END DATA'/                                                      
C                                                                               
C     ESSENTIAL CONSTANTS                                                       
C                                                                               
      AVO=6.0225D23                                                             
      PI=3.14159265D0                                                           
      FF=PI/1.8D2                                                               
      AMU=9.31502D2                                                             
      RMP=1.00727647D0                                                           
      RME=0.5486D-3                                                             
      E=3.7946D-7                                                               
      C1=4.72D-4                                                                
      C2=4.95D-1                                                                
      C=2.9979D10                                                               
      H=6.58195D-22                                                             
      HSQ=H**2                                                                  
      CSQ=C**2                                                                  
      E4=E**4                                                                   
      CCON=E4/(CSQ*HSQ)                                                         
      BRO=1.0D14/C 
      PERM=1.0D0                                                             
      CALL EP (RME,RMP,EPSP)                                                    
C                                                                               
C                                                                               
C                                                                               
C     INPUT DATA                                                                
C                                                                               
  999 READ(5,1)TITLE                                                            
    1 FORMAT(10A8)                                                              
      IF(TITLE(1).EQ.END) GO TO 1000                                            
      TARG=0.0D0                                                                
      TARGB=0.0D0                                                               
      RRM2=0.0D0                                                                
      READ(5,2)KINEM,KRCOIL,KSPRED,KSMEAR,KRSOVR,KDTDX,KTHRSH                   
C                                                                               
C                                                                               
C                                                                               
C       PLEASE NOTE. IMPORTANT                                                  
C     DO NOT INPUT KRSOVR.GT.0 IF KINEM.GT.0   OR VICE VERSA                    
C     DO NOT INPUT KDTDX .GT.0 IF KINEM.GT.0   OR VICE VERSA                    
C     DO NOT INPUT KTHRSH.GT.0 IF KINEM.GT.0   OR VICE VERSA                    
C     NOT THAT KRCOIL CAN ONLY BE USED WITH KINEM                               
C     NOTE THAT KSPRED CAN BE USED WITH KRSOVR AND KDTDX                        
C                                                                               
C                                                                               
C                                                                               
    2 FORMAT(7I2)                                                               
      IF(KINEM.EQ.2) GO TO 1204                                                 
      READ(5,3)RM1,Z1,TINC,RM2,Z2,T2L                                           
    3 FORMAT(7F10.6)                                                            
      IF(KDTDX.GT.0) GO TO 1203                                                 
      READ(5,3)RM3,Z3,QE3,RM4,Z4,QE4                                            
      IF(KRSOVR.GT.1) READ(5,3)RM1B,Z1B,TINCB,RM2B,Z2B                          
      IF(KRSOVR.GT.0) READ(5,3)RM3B,Z3B,QE3B,RM4B,Z4B,QE4B                      
 1203 IF(KRSOVR.EQ.1.OR.KDTDX.GT.0) READ(5,3)TSTOP,DELT                         
      IF(KINEM.EQ.1)                READ(5,3)ANGS,ANGF,ANGI,      CMORLB        
      IF(KSPRED.EQ.1.OR.KDTDX.GT.0)  READ(5,3)TARG,TANGL,DENS,ANMOL,      
     1EXPOT,RRM2,ZZ2                                                            
      IMP=0                                                                     
      IF(RRM2.GT.0.0D0) IMP=1                                                   
      IF(KSPRED.EQ.1.AND.KRSOVR.GT.1) READ(5,4)TARGB,TANGLB,DENSB,       
     1ANMOLB,EXPOTB                                                             
      IF(KSMEAR.EQ.1) READ(5,4)BEAM,DIVERG,SPOT,HORIZ,RADIUS,DETEC              
    4 FORMAT(6F10.6)                                                            
      IF(KSPRED.EQ.0.AND.KDTDX.EQ.0) GO TO 10                                   
       IF(EXPOT.LT.0.0D0) READ(5,6)DTDX                                         
      IF(KRSOVR.LT.2) GO TO 10                                                  
      IF(EXPOTB.LT.0.0D0.AND.TARGB.GE.0.0D0) READ(5,6)DTDXB                     
    6 FORMAT(F10.4)                                                             
   10 CONTINUE                                                                  
      GO TO 1205                                                                
 1204 READ(5,3)RM3,Z3,QE3,RM4,Z4,QE4                                            
 1205 CONTINUE                                                                  
      IEX=IDINT(EXPOT)                                                          
      IEXB=IDINT(EXPOTB)                                                        
      IZ1=IDINT(Z1)                                                             
      IZ2=IDINT(Z2)                                                             
      IZ3=IDINT(Z3)                                                             
      IZ4=IDINT(Z4)                                                             
      ITB=IDINT(TARGB)                                                          
      IF(KRSOVR.EQ.0) GO TO 21205                                               
      IZ1B=IDINT(Z1B)                                                           
      IZ2B=IDINT(Z2B)                                                           
      IZ3B=IDINT(Z3B)                                                           
      IZ4B=IDINT(Z4B)                                                           
21205 CONTINUE                                                                  
C                                                                               
C     INPUT DATA PRINTED OUT                                                    
      WRITE(6,1121)                                                             
C                                                                               
 1121 FORMAT('1'//'0 -----------------------------'/                            
     1            ' |',29X,'|'/' |',29X,'|'/                                    
     1' |  K   K  III  N    N  K   K  |'/                                       
     2' |  K  K    I   NN   N  K  K   |'/                                       
     3' |  KK      I   N N  N  KK     |'/                                       
     4' |  K K     I   N  N N  K K    |'/                                       
     5' |  K  K    I   N   NN  K  K   |'/                                       
     6' |  K   K  III  N    N  K   K  |'/                                       
     6' |',29X,'|'/' |',29X,'|'/'  -----------------------------'//)            
C                                                                               
C                                                                               
C     THIS PROGRAM WAS WRITTEN BY N.M. CLARKE, KING'S COLLEGE,LONDON            
C                                                                               
C                                                                               
      WRITE(6,1221)                                                             
 1221 FORMAT('0INPUT DATA'/' __________'/)                                      
      WRITE(6,21)TITLE                                                          
   21 FORMAT('0',10A8//)                                                        
      WRITE(6,22)KINEM,KRCOIL,KSPRED,KSMEAR,KRSOVR,KDTDX,KTHRSH                 
   22 FORMAT('0 KINEM =',I3/'  KRCOIL=',I3/'  KSPRED=',I3/'  KSMEAR=',I3/       
     1/'  KRSOVR=',I3/'  KDTDX =',I3/'  KTHRSH=',I3//)                          
      IF(KINEM.EQ.2) GO TO 1325                                                 
      WRITE(6,23 )                                                              
   23 FORMAT('0',10X,'INCIDENT PARTICLE',19X,'TARGET PARTICLE'/' ',4X,'M        
     1ASS',7X,'CHARGE',6X,'ENERGY',7X,'MASS',7X,'CHARGE',6X,'ENERGY')           
      WRITE(6,24)RM1,Z1,TINC,RM2,Z2,T2L                                         
   24 FORMAT('  ',6(F10.6,2X)//)                                                
      IF(KDTDX.GT.0) WRITE(6,1229)TSTOP,DELT                                    
 1229 FORMAT('0 FINAL ENERGY',3X,'ENERGY STEP'/'  ',F10.4,5X,F10.4)             
 1325 CONTINUE                                                                  
      IF(KDTDX.GT.0) GO TO 1222                                                 
      WRITE(6,25 )                                                              
   25 FORMAT('0',11X,'PARTICLE THREE',22X,'PARTICLE FOUR'/' ',4X,'MASS',7X,'CHAR
     17X,'CHARGE',4X,'EXCITATION',4X,'MASS',7X,'CHARGE',4X,'EXCITATION')        
      WRITE(6,24)RM3,Z3,QE3,RM4,Z4,QE4                                          
      IF(KINEM.EQ.2) GO TO 1328                                                 
 1222 CONTINUE                                                                  
      IF(KRSOVR-1) 1224,1226,1226                                               
 1226 WRITE(6,1223 )                                                            
 1223 FORMAT('0 SECOND SET OF PARTICLES'/'  _______________________')           
      IF(KRSOVR-1) 1224,1227,1228                                               
 1228 WRITE(6,23)                                                               
      WRITE(6,1324 )RM1B,Z1B,TINCB,RM2B,Z2B                                     
 1324 FORMAT('  ',5(F10.6,2X)//)                                                
 1227 WRITE(6,25)                                                               
      WRITE(6,24)RM3B,Z3B,QE3B,RM4B,Z4B,QE4B                                    
      IF(KRSOVR.EQ.1) WRITE(6,1229)TSTOP,DELT                                   
 1224 CONTINUE                                                                  
      IF(KINEM.EQ.1)                WRITE(6,26)ANGS,ANGF,ANGI,CMORLB            
   26 FORMAT('0  FIRST ANGLE    LAST ANGLE    INCREMENT   CM. OR LAB.'/         
     1' ',4(4X,F7.3,3X)//)                                                      
      IF(KSPRED.EQ.1.OR.KDTDX.GT.0)WRITE(6,27)TARG,TANGL,DENS,ANMOL,       
     1 EXPOT
   27 FORMAT('0  TARGET THICKNESS  = ',F8.3 ,'   MILLIGRAMS/SQ.CM.'/            
     1'   TARGET ANGLE      = ',F8.3,'   DEG.'/                                 
     1'   TARGET DENSITY    = ',F8.3 ,'   GRAMS/CC.'/                           
     2'   NO. OF ATOMS/MOL. = ',F8.3 /'   EXCITATION POTL.  = ',F8.3 ,' ELEC    
     3  ELECTRON VOLTS'//)                                                      
      IF(IMP.GT.0) WRITE(6,3927)RRM2,ZZ2                                        
 3927 FORMAT('   TARGET MASS       =',F8.3/'   TARGET CHARGE     =',F8.3        
     1/'  KINEMATICS ARE FOR IMPURITY IN ABOVE TARGET'//)                       
 1327 FORMAT('0 SECOND TARGET'/)                                                
      IF(KSPRED.EQ.1.AND.KRSOVR.GT.1) WRITE(6, 1327)                            
      IF(KSPRED.EQ.1.AND.KRSOVR.GT.1) WRITE(6, 27)TARGB,TANGLB,DENSB,      
     1ANMOLB,EXPOTB                                                             
      IF (KSMEAR.EQ.1) WRITE(6,28)BEAM,DIVERG,SPOT,HORIZ,RADIUS,DETEC           
   28 FORMAT('0  FWHM. OF BEAM             = ',F8.4 ,'  MEV  '/                 
     1'1  DIVERGENCE OF BEAM        = ',F8.4 ,'  DEGREES'/                      
     2'   BEAM SPOT WIDTH           =',                                '        
     1 ',F8.4 , '  CMS.'/'   HORIZONTAL APERTURE WIDTH = ',F8.4 ,'  CMS'/       
     3/'   RADIUS OF APERTURE        = ',F8.4 ,'  CMS.'/                        
     4'   DETECTOR RESOLUTION FWHM  = ',F8.4 ,'  MEV')                          
      IF(KDTDX.GT.0.AND.IEX.EQ.-1)      WRITE(6,29)DTDX                         
      IF(KSPRED.GT.0.AND.IEX.EQ.-2)      WRITE(6,1029)DTDX                      
   29 FORMAT('0  STOPPING POWER, DTDX = ',F12.4,'  KEV/MICRON OR MEV/MM. ')     
     1 ')                                                                       
 1029 FORMAT('0  STOPPING POWER, DTDX = ',F12.4,'  KEV/MG./CMSQ OR MEV/G        
     1M/CMSQ')                                                                  
      IF(KRSOVR.GT.1.AND.IEXB.EQ.-1)     WRITE(6, 1329)DTDXB                    
      IF(KRSOVR.GT.1.AND.IEXB.EQ.-2)     WRITE(6,11329)DTDXB                    
 1329 FORMAT('0  SECOND TARGET'/'   STOPPING POWER, DTDXB= ',F12.4,' KEV/MICRON 
     1/MICRON OR MEV/MM. ')                                                     
11329 FORMAT('0  SECOND TARGET'/'   STOPPING POWER, DTDXB= ',F12.4,' KEV        
     1/MG/CMSQ OR MEV/GM/CMSQ')                                                 
 1328 CONTINUE                                                                  
C                                                                               
C                                                                               
C     PRINT  ,THE OPTIONS SPECIFIED BY THE CONTROL PARAMETERS                   
C                                                                               
C                                                                               
      WRITE(6,30 )                                                              
   30 FORMAT('0'//' OPTIONS SPECIFIED'/' _________________'//)                  
      IF(KINEM.EQ.1) WRITE(6,31 )                                               
      IF(KINEM.EQ.2) WRITE(6,1331)                                              
 1331 FORMAT('0-NEW VALUES FOR PARTICLES THREE AND FOUR,SAME PROJECTILE AND TARG
     1AND TARGET')                                                              
   31 FORMAT('0-RELATIVISTIC TWO BODY KINEMATICS')                              
      IF(KRCOIL.GT.0) WRITE(6,32  )                                             
   32 FORMAT('0-INCLUDING KINEMATICS FOR RECOIL PARTICLE FOUR')                 
      IF(KRCOIL.EQ.3) WRITE(6,34 )                                              
   34 FORMAT('0-ALSO KINEMATICS')                                               
      IF(KRCOIL.GE.2) WRITE(6,33)                                               
   33 FORMAT('0-AT IDENTICAL ANGLES AS PARTICLE THREE')                         
      IF(KRSOVR.EQ.1) WRITE(6,37 )                                              
   37 FORMAT('0-KINEMATIC CROSSOVER (ENERGY) CALCULATION')                      
      IF(KRSOVR.EQ.2) WRITE(6,38 )                                              
   38 FORMAT('0-CALCULATION OF KINEMATIC CROSSOVER ANGLE (ENERGY) FOR TW        
     1O DIFFERENT REACTIONS')                                                   
      IF(KRSOVR.EQ.3) WRITE(6,39 )                                              
   39 FORMAT('0-CALCULATION OF KINEMATIC CROSSOVER ANGLE (MAGNETIC RIGIDTY      
     ITY) FOR TWO DIFFERENT REACTIONS')                                         
      IF(KDTDX.GT.0) WRITE(6,41)                                                
   41 FORMAT('0-CALCULATION OF STOPPING POWERS,RANGES,ENERGY SPREAD AND EFFECT
     1EFFECTIVE CHARGES')                                                       
      IF(KTHRSH.GT.0) WRITE(6,42 )                                              
   42 FORMAT('0-CALCULATION OF REACTION THRESHOLDS')                            
      IF(KSPRED.GT.0) WRITE(6,35 )                                              
   35 FORMAT('0-WITH CALCULATION OF ENERGY LOSS AND SPREAD IN TARGET')        
      IF(KSMEAR.GT.0) WRITE(6,36)                                               
   36 FORMAT('0-AND WITH CALCULATION OF ENERGY SPREAD FROM INCIDENT BEAM      
     1,DETECTOR GEOMETRY AND RESOLUTION')                                       
      IF(KRSOVR.GE.2.AND.TARGB.LT.0.0D0) WRITE(6,1235  )                        
 1235 FORMAT('0-SECOND TARGET IDENTICAL TO FIRST TARGET. SECOND REACTION        
     1DUE TO SMALL IMPURITY IN FIRST TARGET')                                   
 9211 FORMAT('1',10A8/)                                                         
C                                                                               
C                                                                               
C     PRINT  ,DERIVED QUANTITIES                                                
C                                                                               
      WRITE(6,11210  )                                                          
11210 FORMAT('0'//' DERIVED QUANTITIES'/' __________________'//)                
C                                                                               
C                                                                               
C     DERIVE BASIC QUANTITIES FOR KINEMATICS                                    
C                                                                               
C                                                                               
C                                                                               
C                                                                               
C     NEXT SECTION DERIVES QUANTITIES  ASSOCIATED WITH CALCULATIONS             
C     OF ENERGY LOSS AND ENERGY SPREAD                                          
C                                                                               
 1340 IF(KSPRED.EQ.0.AND.KDTDX.EQ.0) GO TO 260                                  
      IF(IMP.EQ.0) GO TO 3712                                                   
      ORM2=RM2                                                                  
       OZ2=Z2                                                                   
      RM2=RRM2                                                                  
      Z2=ZZ2                                                                    
 3712 CONTINUE                                                                  
   11 IF(IEX  )12,14,13                                                         
   12 IF(IZ1.EQ.0)    GO TO 1314                                                
      CALL BFROMT(AMU,RM1,TINC,BTA,BSQ)                                         
      CALL EP (RME,RM1,EPS1)                                                    
      CALL AKAY(PERM,AVO,DENS,E4,PI,ANMOL,RME,RM1,RM2,Z1,Z2,AK1)                     
      IF(IEX.EQ.-1)        DTDX=DTDX*1.0D1                                      
      IF(IEX.EQ.-2)        DTDX =DTDX *DENS                                     
      CALL EFCHRG(AK1,RM1,TINC,Z1,AKE,ZTASQ)                                    
      U=0.0D0                                                                   
      EXPOT=0.0D0                                                               
      CALL CKP(Z2,EXPOT,RZ2,Z42,Z52,XZ2,UZ2)                                    
      CALL CKZ(BSQ,U    ,Z2,Z2,Z42,Z52,XZ2,UZ2,RZ2,CKBYZ)                       
      CALL EXFRDT(AKE,AMU,BSQ,CKBYZ,DTDX,EPS1,RM1,EXPOT,U)                      
      WRITE(6,220)EXPOT,U                                                       
  220 FORMAT('0  EXCITATION POTL.   = ',F12.6,'  ELECTRON VOLTS'/'   REDUC      
     1UCED ENERGY   = ',F12.6)                                                  
      GO TO 16                                                                  
 1314 AK1=0.0D0                                                                 
      U1=0.0D0                                                                  
      ZE=0.0D0                                                                  
      DTDX=0.0D0                                                                
      TT1=0.0D0                                                                 
      GO TO 16                                                                  
   14 CALL EXFRZ(Z2,EXPOT)                                                      
      WRITE(6,221 )EXPOT                                                        
 221  FORMAT('0  EXCITATION POTL.   = ',F12.6,' ELECTRON VOLTS')                
   13 IF(IZ1.EQ.0)    GO TO 1314                                                
      CALL EP (RME,RM1,EPS1)                                                    
      CALL AKAY(PERM,AVO,DENS,E4,PI,ANMOL,RME,RM1,RM2,Z1,Z2, AK1)                    
   16 IF(KDTDX.GT.0) GO TO 1341                                                 
      IF(IZ3.EQ.0)    GO TO 1343                                                
      CALL EP (RME,RM3,EPS3)                                                    
      CALL AKAY(PERM,AVO,DENS,E4,PI,ANMOL,RME,RM3,RM2,Z3,Z2, AK3)                    
 1341 TDX=TARG*1.0D-3/DENS                                                      
      CALL CKP(Z2,EXPOT,RZ2,Z42,Z52,XZ2,UZ2)                                    
      WRITE(6,19222)AK1                                                         
19222 FORMAT('0  AK1 = ',D12.6)                                                 
      IF(KDTDX) 1330,1330,1310                                                  
 1310 KLOKA=0                                                                   
      KLOKB=0                                                                   
      DO 1312 I=1,52                                                            
      DO 1312 J=1,12                                                            
 1312 DUM(I,J)=0.0D0                                                            
      WRITE(6,9211)TITLE                                                        
      WRITE(6,1311)                                                             
71311 CONTINUE
      WRITE(6,81311)
      WRITE(6,91311)
 1330 CONTINUE                                                                  
      IF(KTHRSH.GT.0) GO TO 13000                                               
 1342 TNEW=TINC                                                                 
      IF(IZ1.EQ.0)    GO TO 240                                                 
      ITER=0
  241 OLDT=TNEW                                                                 
      CALL BFROMT(AMU,RM1,OLDT,BTA,BSQ)                                         
      CALL UU(AMU,EPS1,RM1,BSQ,EXPOT,U1)                                        
      CALL EFCHRG(AK1,RM1, OLDT,Z1,AKE,ZETASQ)                                  
      ZE  =DSQRT(ZETASQ)*Z1                                                     
      CALL CKZ(BSQ,U1,Z2,Z42,Z52,XZ2,UZ2,RZ2,CKBYZ)                             
      CALL DTBYDX(AKE,AMU,CKBYZ,RM1,BSQ,U1,DTDX)                                
      TT1=DTDX*TDX/DCOS(TANGL*FF)                                               
      TT1=DABS(TT1)                                                             
      TNEW=TINC-TT1/2.0D0                                                       
      ITER=ITER+1
      IF(ITER.GT.20)GO TO 240
      IF(DABS(OLDT/TNEW-1.0D0) -  1.0D-4)240,240,241                            
 1343 AK3=0.0D0                                                                 
      GO TO 1341                                                                
 1349 TT1=-1.0D0                                                                
      TNEW=TINC                                                                 
  240 IF(KDTDX.GT.0) GO TO 1350                                                 
      WRITE(6,20222)AK3                                                         
20222 FORMAT('0  AK3 = ',D12.6)                                                 
      IF(KSPRED.GT.0) TINC=TNEW                                                 
      IF(KSPRED.GT.0) WRITE(6,26210)TINC                                        
26210 FORMAT('0'/' QUANTITIES CORRECTED FOR TARGET ENERGY LOSS'/' ______        
     1_____________________________________'/                                   
     2/'   INCIDENT ENERGY    = ',F12.6,' MEV. ')                               
      WRITE(6,222)        U1,ZE    ,DTDX,TT1                                    
  222 FORMAT(                                      '0  REDUCED ENERGY           
     1= ',F10.4/'   EFFECTIVE CHARGE = ',F10.4/'   STOPPING POWER   = ',        
     2F10.4,'  MEV/CM.'/'   TGT. ENERGY LOSS = ',F10.4,'  MEV.')                
      IF(KRCOIL)260,260,17                                                      
C                                                                               
C                                                                               
C     THIS SECTION CALCULATES THE VALUES OF STOPPING POWER,RANGE,ENERGY         
C     SPREAD,REDUCED THIRD MOMENT, AND EFFECTIVE IONIC CHARGE                   
C                                                                               
C                                                                               
 1350 KLOKA=KLOKA+1                                                             
      IF(IZ1.EQ.0)    GO TO 1399                                                
      FWHM=0.0D0                                                                
      GAM3=0.0D0                                                                
      IF(KLOKA.GT.1) GO TO 1353                                                 
      TMIN=TNEW/U1*2.0D0                                                        
 1353 CONTINUE                                                                  
      IF(TT1.EQ.-1.0D0) GO TO 1358                                              
      IF(TT1.EQ.0.0D0) GO TO 1351                                               
      CALL SIG2(BSQ,C1,TT1,EPS1,EPSP,TNEW,U1,SIGI)                              
      CALL GAMMA3(C2,TT1,EPS1,EPSP,TNEW,U1,GAM3,TRUFAC)                         
      FWHM=2.355D0*DSQRT(SIGI)*TRUFAC                                           
 1351 CONTINUE                                                                  
      DTB=DTDX/1.0D1                                                            
      DTA=(TNEW-TMIN )/2.99D2                                                   
      I=1                                                                       
      J=7                                                                       
      IF(TNEW.LE.TMIN) GO TO 1358                                               
      T=TMIN                                                                    
 1361 CALL BFROMT(AMU,RM1,T,BTA,BSQ)                                            
      CALL UU(AMU,EPS1,RM1,BSQ,EXPOT,U1)                                        
      IF(U1.LT.2.0D0) GO TO 1349                                                
      CALL EFCHRG(AK1,RM1,T,Z1,AKE,ZETASQ)                                      
      CALL CKZ(BSQ,U1,Z2,Z42,Z52,XZ2,UZ2,RZ2,CKBYZ)                             
      CALL DTBYDX(AKE,AMU,CKBYZ,RM1,BSQ,U1,DTDX)                                
      DUM(I,J)=1.0D0/DTDX                                                       
      FAC=(1.0D0-0.5D0*BSQ)/( 1.0D0-BSQ)                                        
      DUM(I,J-6)=AKE*2.0D0*RME/RM1*DUM(I,J)**3*FAC                              
       T=T+DTA                                                                  
      I=I+1                                                                     
      IF(T.GT.TNEW) GO TO 1370                                                  
      IF(I.GT.50) GO TO 1362                                                    
       GO TO 1361                                                               
 1358 CONTINUE                                                                  
      WRITE(6,1359)TINC                                                         
 1359 FORMAT(' ',F9.4,2X,'INVALID CALCULATION-PARTICLE ENERGY TOO LOW FOR PRECIS
     1R PRECISE CALCULATION OF DT/DX AND RANGE')                                
      GO TO 1360                                                                
 1399 WRITE(6,9211)TITLE                                                        
      WRITE(6,1398 )                                                            
 1398 FORMAT('0INVALID CALCULATION-INCIDENT PARTICLE HAS ZERO CHARGE')          
      GO TO 999                                                                 
 1362 I=1                                                                       
      J=J+1                                                                     
      GO TO 1361                                                                
 1370 DTA=DTA*0.33333333333333D0                                                
      STRAG=0.0D0                                                               
      RAYNG=0.0D0                                                               
      DO 1371 I=1,297,2                                                         
      J=(I-1)/50+7                                                              
      L=I-50*(J-7)                                                              
      STRAG=STRAG+DTA*(DUM(L,J-6)+4.0D0*DUM(L+1,J-6)+DUM(L+2,J-6))              
      RAYNG=RAYNG+DTA*(DUM(L,J)+4.0D0*DUM(L+1,J)+DUM(L+2,J))                    
 1371 CONTINUE                                                                  
        STRAG=DSQRT(STRAG)/RAYNG*1.0D2                                          
      RAYNGB=RAYNG*1.0D+3*DENS                                                  
      RAYNG=RAYNG*1.0D1                                                         
      DTA=DTB*1.0D-2 /DENS                                                      
      WRITE(6,1357)TINC,TNEW,TT1,U1,DTB,DTA,RAYNG,RAYNGB,STRAG,FWHM,       
     1GAM3,ZE                                                                   
 1357 FORMAT(' ',3(F9.4,1X),F9.3,3X,D10.4,3X,3(D10.4,2X),2X,F7.3,3X,F9.4        
     1,2X,F8.3,2X,F8.4)                                                         
 1360 CONTINUE                                                                  
      TFAC=DELT                                                                 
      IF(DELT.LT.0.0D0) TFAC=-DELT*TINC                                         
      TINC=TINC+TFAC                                                            
      IF(TINC.GT.TSTOP) GO TO 1389                                              
      IF(KLOKA.EQ.50) GO TO 1355                                                
      GO TO 1342                                                                
 1355 CONTINUE                                                                  
      IF(TINC.GT.TSTOP) GO TO 1389                                              
      WRITE(6,9211)TITLE                                                        
      WRITE(6,1311)                                                             
      WRITE(6,81311)
      WRITE(6,91311)
      KLOKA=1                                                                   
 1311 FORMAT('0  INCIDENT   MEAN', 5X,6HENERGY,5X,7HREDUCED,4X,                
     18HSTOPPING,5X,8HSTOPPING,5X,5HRANGE,7X,5HRANGE,4X,10HSTRAGGLING,          
     23X,6HENERGY,4X,6HGAMMA3,2X,9HEFFECTIVE )                                  
81311 FORMAT('   ENERGY',4X,'ENERGY',4X,
     3'LOSS',7X,'ENERGY',5X,'POWER',8X,'POWER',5X,'DIFFERENCE',2X,'DIFFE        
     4RENCE',2X,'STD. DEVN',4X,'SPREAD',12X,'CHARGE' )                          
91311 FORMAT('    MEV.',6X,'MEV.',
     5  6X,'MEV.',  17X,'MEV/MM.   MEV-CMSQ/MG.     MM.       MG/CMSQ',         
     74X,'PER CENT',5X,'MEV.')                                                  
      GO TO 1342                                                                
 1389 WRITE(6,1390)TMIN                                                         
 1390 FORMAT('0RANGE DIFFERENCES ARE FROM  ',F8.2,' MEV.')                      
      GO TO 999                                                                 
C                                                                               
C     THIS SECTION CALCULATES FURTHER BASIC QUANTITIES                          
C                                                                               
C                                                                               
   17 IF(IZ4.EQ.0)    GO TO 1317                                                
      CALL EP (RME,RM4,EPS4)                                                    
      CALL AKAY(PERM,AVO,DENS,E4,PI,ANMOL,RME,RM4,RM2,Z4,Z2, AK4)                    
      GO TO 1318                                                                
 1317 AK4=0.0D0                                                                 
 1318 WRITE(6,223)AK4                                                           
  223 FORMAT('0  AK4 = ',D12.6)                                                 
  260 CONTINUE                                                                  
      IF(IMP.EQ.0) GO TO 3260                                                   
      RM2=ORM2                                                                  
      Z2=OZ2                                                                    
 3260 CONTINUE                                                                  
      CALL PFROMT(AMU,RM1,TINC,P)                                               
      CALL PFROMT(AMU,RM2,T2L,P2)                                               
      CALL BE(AMU,RM1,RM2,TINC,T2L,P,P2,B)                                      
      CALL GAMMA(B,GAM)                                                         
      CALL WT(AMU,RM1,RM2,TINC,T2L,GAM,WTLAB,WTCM)                              
      CALL WNCOM(AMU,RM1,RM2,WTCM,W1CM)                                         
      CALL PFROMW(AMU,RM1,W1CM,P1CM)                                            
      AK=P1CM/(C*H)                                                             
      DU=RM1+RM2                                                                
      CALL TFROMW(AMU,DU ,WTCM,TCM)                                             
      CALL QVAL(AMU,RM1,RM2,RM3,RM4,QR)                                         
      WTLAB=WTLAB+QE3+QE4                                                       
      WTCM=WTLAB/GAM                                                            
       WRITE(6,210)TCM,P1CM,AK,QR,B,GAM                                         
  210 FORMAT('0  CM. KINETIC ENERGY = ',F12.6,' MEV. '/'   CM. MOMENTUM         
     1      = ',F12.6,' MEV/C '/'   CM. WAVE NUMBER    = ',D12.6/'   Q-V        
     3ALUE OF REACTION= ',F12.6,' MEV. '/'   CM. BETA           = ',F12.6       
     16/'   CM. GAMMA          = ',F14.8)                                       
       IF(TCM.LE.0.0D0) WRITE(6,11211)                                          
11211 FORMAT('0PROBABLE FAULT IN INPUT DATA'/' __________________________       
     1___'//)                                                                   
      CALL WNCOM(AMU,RM3,RM4,WTCM,W3CM)                                         
      CALL BETAW(AMU,RM3,W3CM,BETA3C)                                           
      WRITE(6,233)W3CM,BETA3C                                                   
  233 FORMAT(' CM. TOTAL ENERGY OF PARTICLE THREE  = ', D14.8,'  MEV'/          
     1 ' CM. BETA OF PARTICLE THREE          = ',F12.6)                         
      CALL ANGMAX(B,GAM,BETA3C,ANGM3)                                           
      WRITE(6,231)ANGM3                                                         
  231 FORMAT( ' MAXIMUM LAB ANGLE OF PARTICLE THREE = ',F12.6)                  
      IF(KRCOIL)280,280,270                                                     
  270 CONTINUE                                                                  
      CALL WNCOM(AMU,RM4,RM3,WTCM,W4CM)                                         
      CALL BETAW(AMU,RM4,W4CM,BETA4C)                                           
      WRITE(6,234)W4CM,BETA4C                                                   
  234 FORMAT(' CM. TOTAL ENERGY OF PARTICLE FOUR  = ', D14.8,'  MEV'/           
     1 ' CM. BETA OF PARTICLE FOUR           = ',F12.6)                         
      CALL ANGMAX(B,GAM,BETA4C,ANGM4)                                           
      WRITE(6,232)ANGM4                                                         
  232 FORMAT( ' MAXIMUM LAB ANGLE OF PARTICLE FOUR  = ',F12.6)                  
  280 CONTINUE                                                                  
C                                                                               
C                                                                               
C     THIS SECTION CALCULATES AND PRINTS THE BASIC QUANTITIES FOR THE           
C     SECOND SET OF PARTICLES INVOLVED IN THE KINEMATIC CROSSOVER CALCUL        
C     -ATIONS.                                                                  
C                                                                               
C                                                                               
      IF(KRSOVR.LT.2) GO TO 11280                                               
      WRITE(6,1223)                                                             
      IF(KSPRED)11311,11260,11311                                               
11311 IF(ITB  )21312,21311,21311                                                
21312 TARGB=TARG                                                                
      DENSB=DENS                                                                
      TANGLB=TANGL                                                              
      ANMOLB=ANMOL                                                              
      EXPOTB=EXPOT                                                              
      RMB=RM2                                                                   
      ZB=Z2                                                                     
      GO TO 21313                                                               
21311 RMB=RM2B                                                                  
      ZB=Z2B                                                                    
21313 IF(IEXB  )11312,11314,11313                                               
11312 IF(IZ1B.EQ.0)    GO TO 21314                                              
      CALL BFROMT(AMU,RM1B,TINCB,BTA,BSQB)                                      
      CALL EP(RME,RM1B,EPS1B)                                                   
      CALL AKAY(PERM,AVO,DENSB,E4,PI,ANMOLB,RME,RM1B,RMB ,Z1B,ZB ,AK1B)              
      IF(IEXB.EQ.-1)       DTDXB=DTDXB*1.0D1                                    
      IF(IEXB.EQ.-2)       DTDXB=DTDXB*DENSB                                    
      CALL EFCHRG(AK1B,RM1B,TINCB,Z1B,AKE,ZTASQ)                                
      U=0.0D0                                                                   
      EXPOTB=0.0D0                                                              
      CALL CKP(ZB,EXPOTB,RZ2B,Z42B,Z52B,XZ2B,UZ2B)                              
      CALL CKZ(BSQB,U,ZB,Z42B,Z52B,XZ2B,UZ2B,RZ2B,CKBYZB)                       
      CALL EXFRDT(AK1B,AMU,BSQB,CKBYZB,DTDXB,EPS1B,RM1B,EXPOTB,U)               
      WRITE(6,220)EXPOTB,U                                                      
      GO TO 11316                                                               
21314 AK1B=0.0D0                                                                
      U1B=0.0D0                                                                 
      ZE=0.0D0                                                                  
      DTDX=0.0D0                                                                
      TT1B=0.0D0                                                                
      GO TO 11316                                                               
11314 CALL EXFRZ(ZB ,EXPOTB)                                                    
      WRITE(6,221)EXPOTB                                                        
11313 IF(IZ1B.EQ.0)    GO TO 21314                                              
      CALL EP(RME,RM1B,EPS1B)                                                   
      CALL AKAY(PERM,AVO,DENSB,E4,PI,ANMOLB,RME,RM1B,RMB ,Z1B,ZB ,AK1B)              
11316 IF(IZ3B.EQ.0)    GO TO 21343                                              
      CALL EP(RME,RM3B,EPS3B)                                                   
      CALL AKAY(PERM,AVO,DENSB,E4,PI,ANMOLB,RME,RM3B,RMB ,Z3B,ZB ,AK3B)              
21341 TDXB=TARGB*1.0D-3/DENSB                                                   
      CALL CKP(ZB,EXPOTB,RZ2B,Z42B,Z52B,XZ2B,UZ2B)                              
      IF(IZ1B.EQ.0)    GO TO 11240                                              
      TNEW=TINCB                                                                
      ITER=0
11241 OLDT=TNEW                                                                 
      CALL BFROMT(AMU,RM1B,OLDT,BTAB,BSQB)                                      
      CALL UU(AMU,EPS1B,RM1B,BSQB,EXPOTB,U1B)                                   
      CALL EFCHRG(AK1B,RM1B,OLDT,Z1B,AKE,ZETASQ)                                
      ZE=DSQRT(ZETASQ)*Z1B                                                      
      IF(U1B.LT.2.0D0) GO TO 21349                                              
      CALL CKZ(BSQB,U1B,ZB,Z42B,Z52B,XZ2B,UZ2B,RZ2B,CKBYZB)                     
      CALL DTBYDX(AKE,AMU,CKBYZB,RM1B,BSQB,U1B,DTDXB)                           
      TT1B=DTDXB*TDXB/DCOS(TANGLB*FF)                                           
      TT1B=DABS(TT1B)                                                           
       TNEW=TINCB-TT1/2.0D0                                                     
       ITER=ITER+1
       IF(ITER.GT.20) GO TO 11240
      IF(DABS(OLDT/TNEW-1.0D0) -1.0D-4) 11240,11240,11241                       
21343 AK3=0.0D0                                                                 
      GO TO 21341                                                               
21349 TT1B=-1.0D0                                                               
      TNEW=TINCB                                                                
11240 CONTINUE                                                                  
      IF(KSPRED.GT.0) TINCB=TNEW                                                
      IF(KSPRED.GT.0) WRITE(6,26210 )TINCB                                      
      WRITE(6,222)AK1B,AK3B,U1B,ZE,DTDXB,TT1B                                   
11260 CONTINUE                                                                  
      CALL PFROMT(AMU,RM1B,TINCB,P)                                             
      CALL BE(AMU,RM1B,RM2B,TINCB,0.0D0,P,0.0D0,BB)                             
      CALL GAMMA(BB,GAMB)                                                       
      CALL WT(AMU,RM1B,RM2B,TINCB,0.0D0,GAMB,WTLABB,WTCMB)                      
      CALL WNCOM(AMU,RM1B,RM2B,WTCMB,W1CMB)                                     
      CALL PFROMW(AMU,RM1B,W1CMB,P1CMB)                                         
      AKB=P1CMB/(C*H)                                                           
      DU=RM1B+RM2B                                                              
      CALL TFROMW(AMU,DU,WTCMB,TCMB)                                            
      CALL QVAL(AMU,RM1B,RM2B,RM3B,RM4B,QRB)                                    
      WTLABB =WTLABB+QE3B+QE4B                                                  
      WTCMB=WTLABB/GAMB                                                         
      WRITE(6,210)TCMB,P1CMB,AKB,QRB,BB,GAMB                                    
      IF(TCMB.LE.0.0D0) WRITE(6,11211 )                                         
11212 CONTINUE                                                                  
      CALL WNCOM(AMU,RM3B,RM4B,WTCMB,W3CMB)                                     
      CALL BETAW(AMU,RM3B,W3CMB,BTA3CB)                                         
      WRITE(6,233)W3CMB,BTA3CB                                                  
      CALL ANGMAX(BB,GAMB,BTA3CB,ANGM3B)                                        
      WRITE(6,231)ANGM3B                                                        
11280 CONTINUE                                                                  
      IF(KINEM.EQ.0) GO TO 2000                                                 
C                                                                               
C                                                                               
C     START OF CALCULATION OF KINEMATICS                                        
C                                                                               
C                                                                               
      ANG=ANGS                                                                  
      MAXK3=0                                                                   
      MAXK4=0                                                                   
      WRITE(6,9211)TITLE                                                        
      WRITE(6,9501 )                                                            
      WRITE(6,89501)
      WRITE(6,99501)
 9501 FORMAT('0 N  PART.  LAB.      CM.       LAB. ENERGY  ENERGY LOSS          
     1  XSEC RATIO   LAB. MOMENTUM  MAGNETIC   RECIPROCAL  REDUCED EFFEC        
     2TIVE')                                                                    
89501 FORMAT(' ',8X,'  ANGLE     ANGLE      CORRECTED   ' ,14X,'  CM/LAB        
     4       CORRECTED      RIGIDITY   VELOCITY    ENERGY  CHARGE')             
99501 FORMAT(' ',10X,'DEG.',6X,'DEG.',           10X,'MEV.',10X,'MEV.',2        
     64X,'MEV/C',4X,'  GAUSS-CM.  NSEC/METRE')                                  
      KLOKA=0                                                                   
      KLOKB=0                                                                   
  500 CONTINUE                                                                  
      IF(CMORLB) 51,51,62                                                       
C                                                                               
C                                                                               
C     THIS SECTION FOR LAB ANGLES                                               
C                                                                               
   51 IF(MAXK3.GT.0) GO TO 56                                                   
      ANGL=ANG                                                                  
      ANGM=ANGM3                                                                
      RM=RM3                                                                    
      BETANC=BETA3C                                                             
      WNCM=W3CM                                                                 
      EPS=EPS3                                                                  
      AK=AK3                                                                    
      Z=Z3                                                                      
      IZ=IZ3                                                                    
      N=3                                                                       
      KSCAT=1                                                                   
      KLOKA=KLOKA+1                                                             
      KLOKB=KLOKB+1                                                             
      GO TO 100                                                                 
   52 IF(MAXANG.EQ.1) GO TO 521                                                 
C                                                                               
C     PRINT  ,OUT CALCULATED KINEMATICS                                         
 8000 FORMAT(' ',I3,2X,I1,3X,2(F8.3,2X),4(F10.4,4X),D12.6,3X,F6.2,4X,F7.2,4X    
     12,4X,F6.3)                                                                
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      GO TO 1998                                                                
  521 MAXK3=MAXK3+1                                                             
      IF(MAXK3.GT.1) GO TO 1998                                                 
      WRITE(6,522)KLOKB                                                         
  522 FORMAT(' ',I3,3X,'MAXIMUM LAB. ANGLE EXCEEDED FOR PARTICLE THREE')        
 1998 IF(KRCOIL-2)56,1997,56                                                    
 1997 IF(ANGM3-9.0D1) 53,58,58                                                  
C                                                                               
C                                                                               
C     THIS SUBSECTION APPLIES WHEN THE KINEMATIC EQUATIONS FOR PARTICLE         
C     THREE HAVE TWO SOLUTIONS FOR THE SAME LAB ANGLE.                          
C                                                                               
C                                                                               
   53 IF(MAXANG.EQ.1) GO TO 58                                                  
      ANGL=ANG                                                                  
      ANGM=ANGM3                                                                
      RM=RM3                                                                    
      BETANC=BETA3C                                                             
      WNCM=W3CM                                                                 
      EPS=EPS3                                                                  
      AK=AK3                                                                    
      Z=Z3                                                                      
      IZ=IZ3                                                                    
      N=3                                                                       
      KSCAT=2                                                                   
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
      GO TO 100                                                                 
   54 IF(MAXANG.EQ.1) GO TO 541                                                 
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      GO TO 56                                                                  
  541 WRITE(6,522)KLOKB                                                         
      GO TO 56                                                                  
C                                                                               
C     NEXT SECTION FOR RECOIL OF ABOVE                                          
C                                                                               
   56 IF(KRCOIL.EQ.0) GO TO 1999                                                
      IF(KRCOIL.EQ.2) GO TO 58                                                  
      IF(MAXK3.GT.0 ) GO TO 58                                                  
      ANGC=ANGC-1.8D2                                                           
      RM=RM4                                                                    
      BETANC=BETA4C                                                             
       WNCM=W4CM                                                                
      EPS=EPS4                                                                  
      AK=AK4                                                                    
      Z=Z4                                                                      
      IZ=IZ4                                                                    
      N=4                                                                       
      KRCL=1                                                                    
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
      GO TO 200                                                                 
   57 CONTINUE                                                                  
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      IF(KSCAT.EQ.1) GO TO 1997                                                 
C                                                                               
C                                                                               
C     THIS SUBSECTION FOR RECOIL AT SAME LAB ANGLE AS PARTICLE THREE            
C     NOTE THAT THERE ARE USUALLY TWO SOLUTIONS                                 
C                                                                               
   58 IF(KRCOIL.EQ.1) GO TO 1999                                                
      IF(MAXK4.GT.0) GO TO 1999                                                 
      ANGL=ANG                                                                  
      ANGM=ANGM4                                                                
      RM=RM4                                                                    
      BETANC=BETA4C                                                             
      WNCM=W4CM                                                                 
      EPS=EPS4                                                                  
      AK=AK4                                                                    
      N=4                                                                       
      Z=Z4                                                                      
      IZ=IZ4                                                                    
      KSCAT=3                                                                   
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
       GO TO 100                                                                
   59 IF(MAXANG.EQ.1) GO TO 581                                                 
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      GO TO 60                                                                  
  581 MAXK4=MAXK4+1                                                             
      IF(MAXK4.GT.1) GO TO 1999                                                 
      WRITE(6,582)KLOKB                                                         
  582 FORMAT(' ',I3,3X,'MAXIMUM LAB. ANGLE EXCEEDED FOR PARTICLE FOUR')         
      GO TO 1999                                                                
   60 IF(ANGM4.GE.9.0D1) GO TO 1999                                             
      KSCAT=4                                                                   
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
      GO TO 100                                                                 
   61 IF(MAXANG.EQ.1) GO TO 611                                                 
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      GO TO 1999                                                                
  611 WRITE(6,582)KLOKB                                                         
      GO TO 1999                                                                
C                                                                               
C     THIS SECTION FOR CM. ANGLES                                               
C                                                                               
   62 ANGC=ANG                                                                  
      RM=RM3                                                                    
      BETANC=BETA3C                                                             
      WNCM=W3CM                                                                 
      EPS=EPS3                                                                  
      N=3                                                                       
      Z=Z3                                                                      
      IZ=IZ3                                                                    
      AK=AK3                                                                    
      KLOKA=KLOKA+1                                                             
      KLOKB=KLOKB+1                                                             
      KRCL=2                                                                    
      GO TO 200                                                                 
   63 CONTINUE                                                                  
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
C                                                                               
C                                                                               
C     THIS SECTION FOR RECOIL OF ABOVE                                          
C                                                                               
      IF(KRCOIL.EQ.0) GO TO 1999                                                
      IF(KRCOIL.EQ.2) GO TO 65                                                  
      ANGC=ANGC-1.8D2                                                           
      RM=RM4                                                                    
      BETANC=BETA4C                                                             
      WNCM=W4CM                                                                 
      EPS=EPS4                                                                  
      AK=AK4                                                                    
      N=4                                                                       
      Z=Z4                                                                      
      IZ=IZ4                                                                    
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
      KRCL=3                                                                    
      GO TO 200                                                                 
   64 CONTINUE                                                                  
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
   65 IF(KRCOIL.EQ.1) GO TO 1999                                                
      ANGC=ANGC                                                                 
      KRCL=4                                                                    
      KLOKB=KLOKB+1                                                             
      KLOKA=KLOKA+1                                                             
       GO TO 200                                                                
   66 CONTINUE                                                                  
      WRITE(6,8000)KLOKB,N,ANGL,ANGC,  TBAR,TMINUS,RATIO,P,BRHO,RV,US,ZE        
      IF(MAXK3.GT.0.AND.MAXK4.GT.0) GO TO 10000                                 
 1999 CONTINUE                                                                  
      IF(ANG.GE.ANGF) GO TO 10000                                               
      IF(KLOKA.GE.50) GO TO 10000                                               
  501 CONTINUE                                                                  
      ANG=ANG+ANGI                                                              
      GO TO 500                                                                 
10000 IF(KSPRED.EQ.0.AND.KSMEAR.EQ.0) GO TO 9015                                
      IF(KLOKA.LT.25) GO TO 9021                                                
      WRITE(6,9211)TITLE                                                        
 9021 WRITE(6,9001 )                                                            
      WRITE(6,89001)
      WRITE(6,99001)
 9001 FORMAT('0',8X,'|',37X,'|',42X,'|',19X,'|',1X/                             
     1' N  PART.|',8X,'ENERGY SPREAD IN TARGET',6X,'|',10X,'KINEMATIC EN      
     2ERGY SPREAD',9X,'|', 'TOTAL ENERGY SPREAD|  MULTIPLE   TOTAL')            
89001 FORMAT(
     4' ',8X,'|FLAT-TOP  GAUSSIAN ',8X,'TOTAL',5X,'|  ANGULAR   PLUS BEA        
     5M MULTIPLE   BEAM     |INCL. BEAM+DETEC   |  SCATT.   SMEARING')          
99001 FORMAT(
     6' ',8X,'|',9X,'   FWHM',3X,'|',2X,'FWHM',6X,'FWTM | ACCEPTANCE   S        
     7POT   SCATTERING  DIVERG.  |FWHM      FWTM     |  ANGLE',6X,'ANGLE        
     8'/' ',8X,'|  MEV.',6X,'MEV.',3X,'|  MEV.',6X,'MEV. |    MEV.',3(6X        
     9,'MEV.'),4X,'| MEV.',6X,'MEV.',4X,'|  DEG.',7X,'DEG.')                    
      DO 9010 I=1,KLOKA                                                         
      K=KLOKB-KLOKA+I                                                           
      IF(IDUM(I)+1) 9006,9004,9003                                              
 9006 IF(IDUM(I)+3) 9011,9009,9008                                              
 9008 WRITE(6,9007)K                                                            
      GO TO 9003                                                                
 9007 FORMAT(' ',I3,'   INVALID CALCULATION OF ENERGY SPREAD-PARTICLE EN        
     1ERGY TOO LOW')                                                            
 9009 WRITE(6,9012)K                                                            
      GO TO 9003                                                                
 9012 FORMAT(' ',I3,'   INVALID CALCULATION OF MULTIPLE SCATTERING-ENERGY       
     1Y OF INCIDENT PARTICLE OUTSIDE RANGE OF THEORY')                          
 9011 WRITE(6,9013)K                                                            
 9013 FORMAT(' ',I3,'   INVALID CALCULATION OF MULTIPLE SCATTERING-ENERGY       
     1Y OF SCATTERED PARTICLE OUTSIDE RANGE OF THEORY')                         
      GO TO 9003                                                                
 9004 WRITE(6,9005)K                                                            
 9005 FORMAT(' ',I3,'   INVALID CALCULATION-PARTICLES SCATTERED THROUGH         
     1PLANE OF TARGET')                                                         
      GO TO 9010                                                                
 9003 CONTINUE                                                                  
      WRITE(6,9002)K,IDUM(I),(DUM(I,J),J=1,12)                                  
 9002 FORMAT(' ',I3,2X,I1,3X,10(F7.4,3X),1X,F7.4,4X,F7.4)                       
 9010 CONTINUE                                                                  
 9015 CONTINUE                                                                  
      IF(MAXK3.GT.0.AND.MAXK4.GT.0) GO TO 999                                   
      IF(ANG.GE.ANGF) GO TO 999                                                 
      WRITE(6,9211)TITLE                                                        
      WRITE(6,9501 )                                                            
      WRITE(6,89501)
      WRITE(6,99501)
      KLOKA=0                                                                   
      GO TO 501                                                                 
C                                                                               
C     NEXT SECTION FOR CALCULATION FOR CM. ANGLES                               
C                                                                               
  200 CALL WCTOL1(B,GAM,BETANC,WNCM,ANGC,WNL)                                   
      CALL ANCTOL(B,GAM,BETANC,ANGC,ANGL)                                       
      CALL SIGMA(B,GAM,BETANC,ANGC,RATIO)                                       
      CALL TFROMW(AMU,RM,WNL,T)                                                 
      IDUM(KLOKA)=0                                                             
      DO 21401 JJ=1,12                                                          
21401 DUM(KLOKA,JJ)=0.0D0                                                       
      BRHO=1.0D0                                                                
      ZE=-1.0D0                                                                 
      US=-1.0D0                                                                 
      IF(KSPRED.NE.1) GO  TO 21400                                              
      KSPR=2                                                                    
      GO TO 3000                                                                
C                                                                               
C     ENERGY SPREAD   ROUTINE IN HERE                                           
C                                                                               
21400 TBAR=T                                                                    
      TMINUS=0.0D0                                                              
  214 CONTINUE                                                                  
       CALL PFROMT(AMU,RM,TBAR,P)                                               
      RV=RM*AMU/(P*C*1.0D-11)                                                   
      IF(IZ.EQ.0)    GO TO 21214                                                
      CALL BERHO(BRO,Z,P,BRHO)                                                  
21214 IF(KSMEAR.NE.1) GO TO 215                                                 
      IF(IDUM(KLOKA).LE.-2) GO TO 215                                           
C     ANGULAR SMEARING ROUTINES IN HERE                                         
       KSMR=2                                                                   
       GO TO 4000                                                               
  215 GO TO (57,63,64,66),KRCL                                                  
C                                                                               
C     NEXT SECTION FOR CALCULATION FO RLAB. ANGLES                              
C                                                                               
  100 IF(ANGL.GT.ANGM) GO TO 120                                                
      MAXANG=0                                                                  
      CALL WCTOL2 (AMU,RM,B,GAM,WNCM,ANGL,WNLA,WNLB)                            
      IL=(-1)**KSCAT                                                            
      IF(IL.GE.0) WL=WNLB                                                       
      IF(IL.LT.0) WL=WNLA                                                       
      CALL BETAW(AMU,RM,WL,BETA)                                                
      CALL ALTOC2(B,GAM,BETA,ANGL,ANGC)                                         
      CALL SIGMA(B,GAM,BETANC,ANGC,RATIO)                                       
      CALL TFROMW(AMU,RM,WL,T)                                                  
      IDUM(KLOKA)=0                                                             
      DO 11401 JJ=1,12                                                          
11401 DUM(KLOKA,JJ)=0.0D0                                                       
      BRHO=-1.0D0                                                               
      ZE=-1.0D0                                                                 
      US=-1.0D0                                                                 
      IF(KSPRED.NE.1) GO TO 11400                                               
      KSPR=1                                                                    
       GO TO 3000                                                               
C                                                                               
C     ENERGY SPREAD   ROUTINE IN HERE                                           
C                                                                               
11400 TBAR=T                                                                    
      TMINUS=0.0D0                                                              
  114 CONTINUE                                                                  
       CALL PFROMT(AMU,RM,TBAR,P)                                               
      RV=RM*AMU/(P*C*1.0D-11)                                                   
      IF(IZ.EQ.0)    GO TO 21114                                                
      CALL BERHO(BRO,Z,P,BRHO)                                                  
21114 IF(KSMEAR.NE.1) GO TO 115                                                 
      IF(IDUM(KLOKA).LE.-2) GO TO 115                                           
C       ANGULAR SMEARING ROUTINES IN HERE                                       
      KSMR=1                                                                    
      GO TO 4000                                                                
  120 MAXANG=1                                                                  
      IDUM(KLOKA)=N                                                             
      DO 21120 JJ=1,12                                                          
21120 DUM(KLOKA,JJ)=0.0D0                                                       
  115 GO TO (52,54,59,61),KSCAT                                                 
C                                                                               
C     ENRGY SPREAD   THROUGH TARGET                                             
C                                                                               
 3000 ALPHA=ANGL-TANGL                                                          
      ALP=DABS(ALPHA)                                                           
      IF(IMP.EQ.0) GO TO 3905                                                   
      ORM2=RM2                                                                  
      OZ2=Z2                                                                    
      RM2=RRM2                                                                  
      Z2=ZZ2                                                                    
 3905 CONTINUE                                                                  
      TR=T/TINC                                                                 
      TX=(TINC+TT1/2.0D0)*TR                                                    
       IF(ALP  .GE.9.0D1) TX=(TINC-TT1/2.0D0)*TR                                
      IF(IZ.EQ.0)    GO TO 3005                                                 
      IF(DABS(ALP-90.0D0).LT.1.0D-4) GO TO 3002                                 
      GO TO 3001                                                                
 3005 TTS=0.0D0                                                                 
      GO TO 3010                                                                
 3002 IDUM(KLOKA)=-1                                                            
 3003 TBAR=T                                                                    
      TMINUS=-1.0D0                                                             
      GO TO 701                                                                 
 3001 CONTINUE                                                                  
      TN=TX                                                                     
      ITER=0
  700 TO=TN                                                                     
      CALL BFROMT(AMU,RM,TO,BETA,BESQ)                                          
      CALL UU(AMU,EPS,RM,BESQ,EXPOT,US)                                         
      IF(US.LT.2.0D0) GO TO 709                                                 
      CALL EFCHRG(AK,RM,TO,Z,AKE,ZSQ)                                           
      ZE=DSQRT(ZSQ)*Z                                                           
      CALL CKZ(BESQ,US,Z2,Z42,Z52,XZ2,UZ2,RZ2,CKBYZ)                            
      CALL DTBYDX(AKE,AMU,CKBYZ,RM,BESQ,US,DT)                                  
      TTS=DABS(DT*TDX/DCOS(ALP  *FF))                                           
      TN=TX-TTS/2.0D0                                                           
      ITER=ITER+1
      IF(ITER.GT.20) GO TO 701
      IF(DABS(TO/TN-1.0D0)-1.0D-4) 701,701,700                                  
  709 IDUM(KLOKA)=-2                                                            
      GO TO 3003                                                                
  701 CONTINUE                                                                  
      TBAR=TINC*TR-TTS/2.0D0                                                    
      TMINUS=(TINC+TT1/2.0D0)*TR-TBAR                                           
      IF(KRSOVR.GT.0.OR.KTHRSH.GT.0) GO TO 3010                                 
      TFTB=DABS((TINC-TT1/2.0D0)*TR-TX+TTS)                                     
C                                                                               
C     TFTB IS ENERGY OF PARTICLES THAT SCATTER FROM THE FRONT OF THE            
C     TARGET MINUS ENERGY OF PARTICLES THAT SCATTER FROM THE BACK               
C                                                                               
      IF(ALP  .GE.9.0D1) TFTB=DABS((TINC+TT1/2.0D0)*TR-TX+TTS)                  
      DTA=TT1/2.0D0                                                             
      DTB=TTS/2.0D0                                                             
      TA=TINC+TT1/4.0D0                                                         
      TB=TINC*TR-TTS/4.0D0                                                      
      IF(DTA.GT.5.0D-1*TA) GO TO 709                                            
      IF(DTB.GT.5.0D-1*TB) GO TO 709                                            
      IF(TA.LT.0.0D0.OR.TB.LT.0.0D0) GO TO 709                                  
      IF(IZ1.EQ.0)    GO TO 730                                                 
      CALL GAMMA3(C2,DTA,EPS1,EPSP,TA,U1,GAM3,TRUFAC)                           
      CALL SIG2(BSQ,C1,DTA,EPS1,EPSP,TA,U1,SIGI)                                
      SIGI=SIGI*TRUFAC**2                                                       
  732 IF(IZ.EQ.0)    GO TO 731                                                  
      CALL GAMMA3(C2,DTB,EPS,EPSP,TB,US,GAM3,TRUFAC)                            
      CALL SIG2(BESQ,C1,DTB,EPS,EPSP,TB,US,SIGS)                                
      SIGS=SIGS*TRUFAC**2                                                       
      GO TO 733                                                                 
  730 SIGI=0.0D0                                                                
      GO TO 732                                                                 
  731 SIGS=0.0D0                                                                
  733 SIGT=SIGI+SIGS                                                            
      IF(SIGT.LE.0.0D0) GO TO 734                                               
      SIG=DSQRT(SIGT)                                                           
      GO TO 735                                                                 
  734 SIG=0.0D0                                                                 
  735 FLTOP=TFTB                                                                
      KREF=0                                                                    
      IF(ALP  .GE.90.D0) KREF=1                                                 
      CALL FOLD(FLTOP,SIG ,KREF,FWHM,FWTM)                                      
      IDUM(KLOKA)=N                                                             
      DUM(KLOKA,1 )=FLTOP                                                       
      DUM(KLOKA,2 )=2.355D0*SIG                                                 
      DUM(KLOKA,3 )=FWHM                                                        
      DUM(KLOKA,4 )=FWTM                                                        
 3010 IF(IMP.EQ.0) GO TO 3910                                                   
      RM2=ORM2                                                                  
      Z2=OZ2                                                                    
 3910 GO TO (114,214,6003,6006,6008,6011),KSPR                                  
C                                                                               
C     ANGULAR SMEARING                                                          
C                                                                               
C     CALCULATES EFFECTIVE BEAM SPOT SIZE LOOKING FROM APERTURE                 
C                                                                               
 4000 CONTINUE                                                                  
      IF(IMP.EQ.0) GO TO 4905                                                   
      ORM2=RM2                                                                  
      OZ2=Z2                                                                    
      RM2=RRM2                                                                  
      Z2=ZZ2                                                                    
 4905 CONTINUE                                                                  
      ALPHA=ANGL-TANGL                                                          
      ALP=DABS(ALPHA)                                                           
      IF(IZ.EQ.0) GO TO 4021                                                    
      IF(DABS(ALP-90.0D0).LT.1.0D-4) GO TO 4022                                 
      GO TO 4021                                                                
 4022 IDUM(KLOKA)=-1                                                            
      TBAR=T                                                                    
      TMINUS=-1.0D0                                                             
      GO TO 4030                                                                
 4021 TSPRD=0.0D0                                                               
      DTHETA=0.0D0                                                              
      OSPOT=SPOT                                                                
      IF(HORIZ.EQ.0.0D0) GO TO 4007                                             
      KSPOT=1                                                                   
      SPOT=0.0D0                                                                
 4005 CONTINUE                                                                  
      EFSPOT=SPOT/DCOS(TANGL*FF)                                                
      EFSPOT=DABS(EFSPOT*DCOS(ALPHA*FF))                                        
      DTHETA= DATAN((HORIZ/2.0D0+EFSPOT/2.0D0)/RADIUS)/FF                       
      IKTRL=0                                                                   
 4001 ANGL=ANGL+DTHETA                                                          
      CALL WCTOL2(AMU,RM,B,GAM,WNCM,ANGL,WNLA,WNLB)                             
      IF(IL.GE.0) WL=WNLB                                                       
      IF(IL.LT.0) WL=WNLA                                                       
      CALL TFROMW(AMU,RM,WL,T)                                                  
      IF(IKTRL.GE.1) GO TO 4003                                                 
      TPLUS=T                                                                   
      ANGL=ANGL-3.0D0*DTHETA                                                    
      IKTRL=1                                                                   
      GO TO 4001                                                                
 4003 TSPRD=DABS(TPLUS-T)                                                       
C                                                                               
C     TSPRD IS THE KINEMATIC SPREAD   DUE TO HORIZONTAL APERTURE (WHEN          
C     KSPOT=1) AND DUE TO HORIZONTAL APERTURE AND BEAMSPOT SIZE (WHEN KSPOT     
C     KSPOT=2), ASSUMING FLAT TOP DISTRIBUTIONS                                 
C                                                                               
      ANGL=ANGL+DTHETA                                                          
      GO TO (4007,4004,4006,4008),KSPOT                                         
 4007 SPOT=OSPOT                                                                
      TSPRDA=TSPRD                                                              
      IF(SPOT.EQ.0.0D0) GO TO 4004                                              
      KSPOT=2                                                                   
      GO TO 4005                                                                
 4004 TSPRDB=TSPRD                                                              
      DTHF=DTHETA                                                               
      KSPOT=3                                                                   
      SMB=0.0D0                                                                 
      SMA=0.0D0                                                                 
      IF(IZ1.EQ.0)    GO TO 24004                                               
      TA=TINC+TT1/4.0D0                                                         
      TARGA=DABS(TARG/DCOS(TANGL*FF))                                           
      IF(TARGA.LE.0.0D0) GO TO 24004                                            
      CALL PVFRMT(AMU,RM1,TA,BETSQA,PVA,PVSQA)                                  
      CALL MULSCT(AMU,BETSQA,PVSQA,TARGA,RM2,Z1,Z2,THMA,SMA,IMULT)              
      IF(IMULT.EQ.-1) GO TO 4031                                                
24004 IF(IZ.EQ.0)    GO TO 24005                                                
      TARGB=DABS(TARG/DCOS(ALPHA*FF))                                           
      IF(TARGB.LE.0.0D0) GO TO 24005                                            
      TB=TINC*TR-TTS/4.0D0                                                      
      CALL PVFRMT(AMU,RM ,TB,BETSQB,PVB,PVSQB)                                  
      CALL MULSCT(AMU,BETSQB,PVSQB,TARGB,RM2,Z,Z2,THMB,SMB,IMULT)               
      IF(IMULT.EQ.-1) GO TO 4032                                                
24005 IF(SMA.EQ.0.0D0.AND.SMB.EQ.0.0D0) GO TO 24006                             
      SMT=DSQRT((SMA**2+SMB**2)/2.0D0)                                          
      THMULT=1.1762D0*SMT                                                       
      DTHETA=THMULT                                                             
      IKTRL=0                                                                   
      GO TO 4001                                                                
24006 TSPRD=0.0D0                                                               
      DTHETA=0.0D0                                                              
 4006 TSPRDC=TSPRD                                                              
      DTHM=DTHETA                                                               
      KSPOT=4                                                                   
      IF(DIVERG.EQ.0.0D0) GO TO 4009                                            
      DTHETA=DIVERG/2.0D0                                                       
      IKTRL=0                                                                   
      GO TO 4001                                                                
 4009 TSPRD=0.0D0                                                               
 4008 TSPRDD=TSPRD                                                              
      DTHD=DTHETA                                                               
      FLTSM=2.0D0*DTHF                                                          
      SIGSM=DSQRT(SMT**2+(DIVERG/2.3524D0)**2)                                  
      IF(FLTSM.EQ.0.0D0) GO TO 4010                                             
      KREF=0                                                                    
      CALL FOLD(FLTSM,SIGSM,KREF,FWHMSM,FWTMSM)                                 
      GO TO 4011                                                                
 4010 FWHMSM=SIGSM                                                              
 4011 THET SM=FWHMSM                                                            
      IDUM(KLOKA)=N                                                             
      DUM(KLOKA,5 )=TSPRDA                                                      
      DUM(KLOKA,6 )=TSPRDB                                                      
      DUM(KLOKA,7 )=TSPRDC                                                      
      DUM(KLOKA,8 )=TSPRDD                                                      
      DUM(KLOKA,11)=THMULT                                                      
      DUM(KLOKA,12)=THETSM                                                      
      TOTFLT=DUM(KLOKA,1)+DUM(KLOKA,6)                                          
      TOTG=DSQRT(DUM(KLOKA,2)**2+DUM(KLOKA,7)**2+DUM(KLOKA,8)**2+BEAM**2        
     1+DETEC**2)/2.355D0                                                        
      KREF=0                                                                    
      CALL FOLD(TOTFLT,TOTG,KREF,TOTHM,TOTTM)                                   
      DUM(KLOKA,9 )=TOTHM                                                       
      DUM(KLOKA,10)=TOTTM                                                       
 4030 IF(IMP.EQ.0) GO TO 4930                                                   
      RM2=ORM2                                                                  
      Z2=OZ2                                                                    
 4930 GO TO(115,215),KSMR                                                       
 4031 IDUM(KLOKA)=-3                                                            
      SMA=0.0D0
      GO TO 24004                                                               
 4032 IDUM(KLOKA)=-4                                                            
      SMB=0.0D0
      GO TO 24005                                                               
 2000 CONTINUE                                                                  
      IF(KRSOVR.EQ.0) GO TO 13000                                               
      GO TO (7000,6000,6000,7500,7500),KRSOVR                                   
C                                                                               
C                                                                               
C     THIS SECTION CALCULATES SIMPLE KINEMATIC CROSSOVER FOR ELASTIC            
C     OR INELASTIC SCATTERING OF A PARTICLE FROM TWO DIFFERENT TARGETS          
C     WITHOUT CORRECTIONS FOR ENERGY LOSSES IN THE TARGET.                      
C                                                                               
C                                                                               
 7000 T=TINC                                                                    
      WRITE(6,9211)TITLE                                                        
      WRITE(6,7010 )                                                            
 7010 FORMAT('0    ENERGY    CROSSOVER ANGLE'/'      MEV.',11X,'DEG.'/)         
      I=1                                                                       
 7001 CALL AXOVER(AMU,RM1,RM4,RM4B,QE4,QE4B,T,ANGL)                             
      WRITE(6,7002)T,ANGL                                                       
 7002 FORMAT(' ',F12.4,3X,F9.3)                                                 
      TFAC=DELT                                                                 
      IF(DELT.LT.0.0D0) TFAC=-DELT*T                                            
      T=T+TFAC                                                                  
      I=I+1                                                                     
      IF(T.GT.TSTOP) GO TO 7007                                                 
      IF(I.GT.50) GO TO 7005                                                    
      GO TO 7001                                                                
 7005 I=1                                                                       
      WRITE(6,9211)TITLE                                                        
      WRITE(6,7010 )                                                            
      GO TO 7001                                                                
 7007 GO TO 999                                                                 
 6000 CONTINUE                                                                  
C                                                                               
C     THIS CALCULATION CAN CALCULATE THE CROSSOVER ANGLE FOR ENERGY OR          
C     MAGNETIC RIGIDTY FOR ANY TWO REACTIONS OR ELASTIC OR INELASTIC SCATTERING 
C     OF TWO DIFFERENT PARTICLES FROM TWO DIFFERENT TARGETS                     
C                                                                               
C                                                                               
      IKTRL=1                                                                   
      MSOL=0                                                                    
      NSOL=1                                                                    
      KTRL=1                                                                    
      IF(ANGM3.LT.9.0D1.OR.ANGM3B.LT.9.0D1) NSOL=2                              
      IF(ANGM3.LT.9.0D1) KTRL=2                                                 
      IF(ANGM3B.LT.9.0D1) KTRL=3                                                
      IF(ANGM3.LT.9.0D1.AND.ANGM3B.LT.9.0D1) KTRL=4                             
      NSTEP=18                                                                  
 6001 CONTINUE                                                                  
      I=1                                                                       
      ANGL=0.0D0                                                                
      DELA=DMIN1(ANGM3,ANGM3B)/DFLOAT(NSTEP-1)                                  
 6002 CONTINUE                                                                  
      J=(I-1)/50+1                                                              
      L=I-50*(J-1)                                                              
      CALL WCTOL2(AMU,RM3,B,GAM,W3CM,ANGL,WNLA,WNLB )                           
       CALL TFROMW(AMU,RM3,WNLA,T1)                                             
      IF(KSPRED.EQ.0) GO TO 6004                                                
      T=T1                                                                      
      RM=RM3                                                                    
      EPS=EPS3                                                                  
      Z=Z3                                                                      
      IZ=IZ3                                                                    
      AK=AK3                                                                    
      KSPR=3                                                                    
      GO TO 3000                                                                
 6003 T1=TBAR                                                                   
 6004 GO TO (6007,6005,6007,6005),KTRL                                          
 6005 CALL TFROMW(AMU,RM3,WNLB,T2)                                              
      IF(KSPRED.EQ.0) GO TO 6007                                                
      T=T2                                                                      
      KSPR=4                                                                    
      GO TO 3000                                                                
 6006 T2=TBAR                                                                   
 6007 CALL WCTOL2(AMU,RM3B,BB,GAMB,W3CMB,ANGL,WNLA,WNLB)                        
      CALL TFROMW(AMU,RM3B,WNLA,T3)                                             
      IF(KSPRED.EQ.0) GO TO 6009                                                
      TDXO=TDX                                                                  
      TANGLO=TANGL                                                              
      TINCO=TINC                                                                
      EXPOTO=EXPOT                                                              
      TT1O=TT1                                                                  
      Z2O=Z2                                                                    
      ZZ4=Z42                                                                   
      ZZ5=Z52                                                                   
      XZZ=XZ2                                                                   
      RZZ=RZ2                                                                   
      Z2=ZB                                                                     
      Z42=Z42B                                                                  
      Z52=Z52B                                                                  
      XZ2=XZ2B                                                                  
      RZ2=RZ2B                                                                  
      UZZ=UZ2                                                                   
      UZ2=UZ2B                                                                  
       TDX=TDXB                                                                 
      TANGL=TANGLB                                                              
      TINC=TINCB                                                                
      EXPOT=EXPOTB                                                              
      TT1=TT1B                                                                  
      Z=Z3B                                                                     
      IZ=IZ3B                                                                   
      RM=RM3B                                                                   
      AK=AK3B                                                                   
      EPS=EPS3B                                                                 
      T=T3                                                                      
      KSPR=5                                                                    
      GO TO 3000                                                                
 6008 T3=TBAR                                                                   
 6009 GO TO (6012,6012,6010,6010),KTRL                                          
 6010 CALL TFROMW(AMU,RM3B,WNLB,T4)                                             
      IF(KSPRED.EQ.0) GO TO 6012                                                
      T=T4                                                                      
      KSPR=6                                                                    
      GO TO 3000                                                                
 6011 T4=TBAR                                                                   
      TDX=TDXO                                                                  
      TANGL=TANGLO                                                              
      TINC=TINCO                                                                
      EXPOT=EXPOTO                                                              
      TT1=TT1O                                                                  
      Z2=Z2O                                                                    
      Z42=ZZ4                                                                   
      Z52=ZZ5                                                                   
      XZ2=XZZ                                                                   
      RZ2=RZZ                                                                   
      UZ2=UZZ                                                                   
 6012 CONTINUE                                                                  
 6013 IF(KRSOVR.EQ.2) GO TO 6018                                                
      GO TO (6017,6016,6015,6015),KTRL                                          
 6015 CALL PFROMT(AMU,RM3,T4,P)                                                 
      CALL BERHO(BRO,Z3B,P,T4)                                                  
      IF(KTRL.EQ.3) GO TO 6017                                                  
 6016 CALL PFROMT(AMU,RM3,T2,P)                                                 
      CALL BERHO(BRO,Z3,P,T2)                                                   
 6017 CALL PFROMT(AMU,RM3B,T3,P)                                                
      CALL BERHO(BRO,Z3B,P,T3)                                                  
      CALL PFROMT(AMU,RM3,T1,P)                                                 
      CALL BERHO(BRO,Z3,P,T1)                                                   
 6018 CONTINUE                                                                  
      GO TO (6022,6021,6020,6019),KTRL                                          
 6019 DUM(L,J+9)=T2-T4                                                          
 6020 DUM(L,J+6)=T1-T4                                                          
      IF(KTRL.EQ.3) GO TO 6022                                                  
 6021 DUM(L,J+3)=T2-T3                                                          
 6022 DUM(L,J)=T1-T3                                                            
      ANGL=ANGL+DELA                                                            
      IF(I.EQ.NSTEP) GO TO 6023                                                 
      I=I+1                                                                     
      GO TO (6002,6035,6036     ),IKTRL                                         
 6023 CONTINUE                                                                  
C                                                                               
C                                                                               
C     PREVIOUS SECTION HAS CALCULATED THE DIFFERENCES BETWEEN ALL THE POSSIBLE S
C     SOLUTIONS OF THE OF THE PARTICLE ENERGIES AS A FUNCTION OF ANGLE          
C                                                                               
C     THE NEXT SECTION LOOKS FOR DIFFERENCES OF OPPOSITE SIGN IE. THE           
C     KINEMATIC CROSSOVER POINT                                                 
C                                                                               
C                                                                               
      N=0                                                                       
      M=0                                                                       
 6026 DO 6025 I=2,NSTEP                                                         
      J=(I-1)/50+1+M                                                            
      L=I-50*(J-1-M)                                                            
      JJ=(I-2)/50+1+M                                                           
      LL= I-1-50*(JJ-1-M)                                                       
      IF(DUM(LL,JJ).EQ.0.0D0) GO TO 6024                                        
      IF(DUM(L,J)/DUM(LL,JJ)) 6024,6025,6025                                    
 6024 N=N+3                                                                     
      IF(N.GT.12) GO TO 6025                                                    
      DUM(52,N)=DELA*DFLOAT(I-2)                                                
      DUM(52,N-1)=(DUM(L,J)-DUM(LL,JJ))/DELA                                    
      DUM(52,N-2)=DUM(LL,JJ)                                                    
      IDUM(N)=M                                                                 
 6025 CONTINUE                                                                  
      IF(M.GE.3*(KTRL-1)) GO TO 6030                                            
      GO TO (6030,6027,6028,6027),KTRL                                          
 6027 M=M+3                                                                     
      GO TO 6026                                                                
 6028 M=M+6                                                                     
      GO TO 6026                                                                
 6030 CONTINUE                                                                  
C                                                                               
C     NEXT SECTION COMPARES THE NUMBER OF SOLUTIONS OBTAINED WITH THE           
C     THEORETICAL MAXIMUM NUMBER . IF THESE ARE NOT EQUAL, THE STEP SIZE        
C     IN ANGLE IS HALVED AND THE DIFFERENCES AND CROSSOVER POINTS ARE           
C     RECALCULATED. IF THE SAME NUMBER OF SOLUTIONS ARE OBTAINED  THE           
C     CALCULATION CONTINUES                                                     
C                                                                               
      IF(N/3.GE.NSOL) GO TO 6031                                                
      IF(NSTEP.EQ.144) GO TO 6031                                               
      IF(N.EQ.0) GO TO 36030                                                    
      IF(MSOL.EQ.N/3) GO TO 6031                                                
36030 CONTINUE                                                                  
      MSOL=N/3                                                                  
      NSTEP=NSTEP+NSTEP                                                         
      GO TO 6001                                                                
 6031 NSOL=N/3                                                                  
      IF(N.EQ.0) GO TO 6501                                                     
 6032 M=0                                                                       
 6033 CONTINUE                                                                  
C                                                                               
C     THE NEXT SECTION USES NEWTONS METHOD TO OBTAIN A BETTER VALUE FOR         
C     THE CROSSOVER ANGLES                                                      
C                                                                               
      KLOKA=0                                                                   
 6034 CONTINUE                                                                  
      ANG=DUM(52,M+3)                                                           
      ANGL=ANG-DUM(52,M+1)/DUM(52,M+2)                                          
      IF(ANGL.EQ.0.0D0) TERM=ANGL/ANG                                           
      IF(ANG .EQ.0.0D0) TERM=ANG /ANGL                                          
      IF(DABS(1.0D0-TERM    ).LT.1.0D-6) GO TO 6040                             
      KLOKA=KLOKA+1                                                             
      IF(KLOKA.GT.10) GO TO 6040                                                
      I=1                                                                       
      IKTRL=2                                                                   
      GO TO 6002                                                                
 6035 IKTRL=3                                                                   
      ANGL=ANGL-0.9D0*DELA                                                      
      GO TO 6002                                                                
 6036 KK=IDUM(M+3)                                                              
      DUM(52,M+1)=DUM(I-2,KK+1)                                                 
      DUM(52,M+2)=(DUM(I-1,KK+1)-DUM(I-2,KK+1))/(DELA*1.0D-1)                   
      ANGL=ANGL-1.1D0*DELA                                                      
      DUM(52,M+3)=ANGL                                                          
      GO TO 6034                                                                
 6040 IF(M.GE.3*(NSOL-1)) GO TO 6050                                            
      M=M+3                                                                     
      GO TO 6033                                                                
 6050 CONTINUE                                                                  
      WRITE(6,9211)TITLE                                                        
      WRITE(6,6051  )                                                           
 6051 FORMAT('0    CROSSOVER ANGLES'/'     ________________')                   
      DO 6052 I=1,NSOL                                                          
 6052 WRITE(6,6053)DUM(52,3*(I-1)+3)                                            
 6053 FORMAT('0',6X,F7.3)                                                       
      GO TO 6500                                                                
 6501 WRITE(6,9211)TITLE                                                        
      WRITE(6,6051 )                                                            
      WRITE(6,6502   )                                                          
 6502 FORMAT('0    INVALID CALCULATION-NO CROSSOVER ANGLE(S) FOUND'/'           
     1  PROBABLE ERROR IN INPUT DATA')                                          
 6500 CONTINUE                                                                  
      GO TO 999                                                                 
 7500 CONTINUE                                                                  
13000 CONTINUE                                                                  
      IF(KTHRSH.EQ.0) GO TO 14000                                               
      CALL THRESH(AMU,RM1,RM2,TINC,T2L,QR,TS,IER)                               
      WRITE(6,9211)TITLE                                                        
      WRITE(6,13001 )                                                           
13001 FORMAT('0REACTION THRESHOLD'/' __________________')                       
      WRITE(6,13002)TS                                                          
13002 FORMAT('0',F10.4,' MEV.')                                                 
      CALL TNS(AMU,RM1,RM2,RM3,TS,T2L,QR,T3)                                    
      CALL TNS(AMU,RM1,RM2,RM4,TS,T2L,QR,T4)                                    
      WRITE(6,13004 )                                                           
13004 FORMAT('0'//' LAB. KINETIC ENERGIES AT THRESHOLD'/' PARTICLE 3            
     1 PARTICLE4')                                                              
      WRITE(6,13005)T3,T4                                                       
13005 FORMAT('0',2(F10.4,3X))                                                   
14000 GO TO 999                                                                 
 1000 STOP                                                                      
       END                                                                      
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE MULSCT(AMU,BETSQ,PVSQ,TARG,RM2,Z1,Z2,THETA,SIGM,IMULT)         
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE MULTIPLE SCATTERING USING AN APPROXIMATE THEORY            
C                                                                               
      PI=3.14159265D0                                                           
      IMULT=0                                                                   
       CHICSQ=1.569D-4*Z2*(Z2+1.0D0)*Z1**2*TARG/(RM2*PVSQ)                      
      B=DLOG(2.73D0*(Z2+1.0D0)*Z2**(1.0D0/3.0D0)*Z1**2*TARG/(BETSQ*RM2))        
     1-1.544D-1                                                                 
C                                                                               
C        SOLVE FOR BB   BB-DLOG(BB)=B                                           
C                                                                               
C        TRY  BB=1.5+1.1*B                                                      
C        THEN USE NEWTONS APPROXM.                                              
      BB1= 1.5D0+1.1D0*B                                                        
      IF(BB1.LE.0.0D0)  GO TO 22                                                
    1 BB2=BB1-(BB1-DLOG(BB1)-B)/(1.0D0-1.0D0/BB1)                               
      IF(DABS(1.0D0-BB2/BB1).LT.2.0D-4) GO TO 2                                 
      BB1=BB2                                                                   
      GO TO 1                                                                   
   22 IMULT=-1                                                                  
    2 IF(BB2.LT.3.9D0) IMULT=-1                                                 
      IF(BB2.GT.1.5D1) IMULT=-1                                                 
      IF(IMULT.EQ.-1) GO TO 30                                                  
      X=DSQRT((5.037D-4*BB2**2+0.9956246*BB2-1.273716)/BB2)                     
      THETA=X*DSQRT(CHICSQ*BB2) *1.8D2/PI                                       
      SIGM=THETA/DSQRT(2.0D0)                                                   
   30 RETURN                                                                    
      END                                                                       
      SUBROUTINE AKAY(PERM,AVO,DENS,E4,PI,ANMOL,
     >RME,RM1,RM2,Z1,Z2,AK)                
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE CONSTANT IN THE STOPPING POWER EQUATION                    
C                                                                               
      AA=2.0D0*PI*AVO*DENS*Z2*E4/(RME*RM2*ANMOL)                                
      AK  =AA*Z1**2*RM1/PERM                                                         
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE  ALTOC2(B,GAM,BETANL,ANGL,ANGC)                                
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE CM. ANGLE FROM THE LAB. ANGLE                              
C                                                                               
      ANG=ANGL*1.74532925D-02                                                   
      IF(BETANL.LT.1.0D-08) BETANL=1.0D-08                                      
      X=B/BETANL                                                                
      ANGC=DATAN2(DSIN(ANG),(GAM*(DCOS(ANG)-X)))/1.74532925D-02                 
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE ANCTOL(B,GAM,BETANC,ANGC,ANGL)                                 
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE LAB ANGLE FROM THE CM. ANGLE                               
C                                                                               
      RAD=1.74532925D-02                                                        
      ANG=ANGC*RAD                                                              
      X=B/BETANC                                                                
      ANGL=DATAN2(DSIN(ANG),(GAM*(DCOS(ANG)+X)))/1.74532925D-02                 
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE ANGMAX(B,GAM,BETANC,ANGM)                                      
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE MAXIMUM LAB. ANGLE OF THE PARTICLE                         
C                                                                               
      X=B/BETANC                                                                
      XF=X**2-1.0D0                                                             
      IF(DABS(XF).LE.1.0D-08) XF=0.0D0                                          
      IF(XF.LT.0.0D0) GO TO 1                                                   
      IF(XF.EQ.0.0D0) GO TO 2                                                   
      TANANG=1.0D0/(GAM*DSQRT(XF))                                              
      ANGM=DATAN(TANANG)/1.74532925D-02                                         
      GO TO 3                                                                   
    1 ANGM=1.8D2                                                                
      GO TO 3                                                                   
    2 ANGM=9.0D1                                                                
    3 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE AXOVER(AMU,RM,RM2A,RM2B,QA,QB,T,ANGL)                          
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE KINEMATIC CROSSOVER ANGLE IN ENERGY FOR THE ELASTIC        
C     OR INELASTIC SCATTERING OF ONE PARTICLE FROM TWO DIFFERENT TARGETS        
C                                                                               
      RA=RM*AMU                                                                 
      RB=RM2A*AMU                                                               
      RC=RM2B*AMU                                                               
      P4SQ=(QB-QA)*(2.0D0*RB-2.0D0*RC+QB-QA)*(4.0D0*RB*RC  -(QA+QB)*(2.0        
     1D0*(RB +RC )-QA-QB))/(4.0D0*(RB -RC )**2)                                 
      T4=DSQRT(P4SQ+(RB  -QA)**2) -RB                                           
      DTF=DSQRT((T-T4)*(2.0D0*RA+T-T4)*T*(T+2.0D0*RA))                          
      COSANG=((T-T4)*(2.0D0*RA+T-T4)-P4SQ+T*(T+2.0D0*RA))/(2.0D0*DTF)           
      ANGL= DACOS(COSANG)/1.74532925D-02                                        
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE BE(AMU,RM1,RM2,TI1,TI2,PI1,PI2,B)                              
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE CM. BETA FOR A SYSTEM OF PARTICLES                         
C                                                                               
      SIGN=1.0D0                                                                
      IF(TI2.LT.0.0D0) SIGN=-1.0D0                                              
      B=(PI1+SIGN*PI2)/(TI1+RM1*AMU+RM2*AMU+DABS(TI2))                          
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE BETAW(AMU,RM,W,BETA)                                           
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE BETA OF A PARTICLE GIVEN THE TOTAL ENERGY                  
C                                                                               
      DD=1.0D0 -(AMU*RM/W)**2                                                   
      IF(DD.LT.0.0D0) DD=0.0D0                                                  
      BETA=DSQRT(DD)                                                            
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE BERHO(BRO,Z,P,BRHO)                                            
      IMPLICIT REAL *8 (A-H,O-Z)                                                
C                                                                               
C     CALCULATES THE MAGNETIC RIGIDTY                                           
C                                                                               
      BRHO=BRO*P/Z                                                              
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE BFROMT(AMU,RM,T,B,DD)                                          
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE BETA OF A PARTICLE GIVEN THE KINETIC ENERGY                
C                                                                               
      W=AMU*RM+T                                                                
      DD=1.0D0-(AMU*RM/W)**2                                                    
      IF(DD.LT.0.0D0)  DD=0.0D0                                                 
      B=DSQRT(DD)                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE CKP(Z,EXPOT,RZ,Z4,Z5,XZ,UZ)                                    
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE CONSTANTS FOR SUBROUTINE CKZ                               
C                                                                               
      RZ=1.37D2**2/Z                                                            
      Z4=Z**0.4D0                                                               
      Z5=DSQRT(Z)                                                               
      IF(EXPOT.EQ.0.0D0) EXPOT=9.1D0*Z*(1.0D0+1.9D0/Z**(2.0D0/3.0D0))           
      XZ=0.0D0                                                                  
      UZ=(1.0D0-XZ)*1.022D6/(EXPOT*RZ)                                          
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE CKZ(BETSQ,U,Z,Z4,Z5,XZ,UZ,RZ,CKBYZ)                            
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE INNER SHELL CORRECTION FACTORS                             
C     NOTE THAT THESE HAVE THE CORRECT ASYMPTOTIC BEHAVIOUR AS BETSQ TENDS TO 1 
C     SEE U. FANO'S ARTICLE IN  REPORT ' STUDIES OF PENETRATION OF CHARGED PARTI
C     CLES IM MATTER'                                                           
C                                                                               
      X=RZ*BETSQ                                                                
      X=X+XZ                                                                    
      X6=X**0.6D0                                                               
      IF(X.LT.1.0D0) GO TO 1                                                    
      IF(X.LE.20.0D0) GO TO 2                                                   
      IF(X.LE.200.0D0) GO TO 3                                                  
C     ASYMPTOTIC VALUES                                                         
      CKBYZ=Z4/X+8.0D0*Z/X**2                                                   
      GO TO 5                                                                   
    1 IF(X.EQ.0.0D0) GO TO 10                                                   
      IF(U.EQ.0.0D0) U=UZ*RZ*BETSQ/(1.0D0-XZ)                                   
      X=X-XZ                                                                    
      CKBYZ=-0.1D0*DLOG(1.0D0/X)                                                
      GO TO 5                                                                   
    2 ZFAC=Z5                                                                   
      IF(Z.LT.13.0D0) ZFAC= Z5+4.0D0/Z                                          
      CKBYZ=ZFAC/5.7D0*(1.0D0/X6-1.0D0/X)                                       
      GO TO 5                                                                   
    3 CKBYZ=Z5/(8.15D0*X6)                                                      
      GO TO 5                                                                   
   10 CKBYZ=0.0D0                                                               
    5 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE DTBYDX(AKAY,AMU,CKBYZ,RM,BETSQ,U,DTDX)                         
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE STOPPING POWER                                             
C                                                                               
        TCOR=0.5D0*AMU*RM*BETSQ                                                 
      UCOR=U/DSQRT(1.0D0-BETSQ)                                                 
      DTDX=AKAY*(DLOG(UCOR)-BETSQ-CKBYZ)/TCOR                                   
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE EFCHRG(AK,RM,T,Z,AKEF,ZETASQ)                                  
      IMPLICIT REAL *8 (A-H,O-Z)                                                
C                                                                               
C     CALCULATES THE EFFECTIVE CHARGE OF AN ION AT LOW VELOCITIES               
C                                                                               
      DIMENSION F(31),G(18)                                                     
      DATA F/0.0D0,3.73D-1,5.07D-1,5.76D-1,6.25D-1,6.67D-1,7.05D-1,             
     17.41D-1,7.72D-1,7.99D-1,8.23D-1,8.43D-1,8.63D-1,8.78D-1,8.94D-1,9.        
     206D-1,9.19D-1,9.30D-1,9.41D-1,9.48D-1,9.55D-1,9.61D-1,9.67D-1,            
     39.71D-1,9.74D-1,9.77D-1,9.8D-1,9.83D-1,9.85D-1,9.88D-1,9.9D-1/            
      DATA G/0.0D0,3.82D-1,5.07D-1,5.88D-1,6.5D-1,7.03D-1,7.45D-1,7.85D-1       
     11,8.18D-1,8.5D-1,8.8D-1,9.07D-1,9.3D-1,9.53D-1,9.7D-1,9.83D-1,9.93D-1     
     2D-1,1.0D0/                                                                
      X=T/(  RM*Z**(4.0D0/3.0D0))                                               
      ZETASQ=1.0D0                                                              
      IF(Z.LT.3.0D0) GO TO 2                                                    
      IF(X.GT.2.1D-1) GO TO 1                                                   
      ZETASQ=1.0D0-DEXP(-2.473D1*X+2.476D2*X**2-1.131D3*X**3)                   
      GO TO 1                                                                   
    2 IF(X.GT.0.3D0) GO TO 1                                                    
      L=IDINT(X/1.0D-2)                                                         
      ZETASQ=F(L+1)+(X/1.0D-2-DFLOAT(L))*(F(L+2)-F(L+1))                        
      IF(Z.LT.2.0D0) GO TO 1                                                    
      FAC=1.0D0                                                                 
      IF(L.GT.16) GO TO 4                                                       
      FAC=G(L+1)+(X/1.0D-2-DFLOAT(L))*(G(L+2)-G(L+1))                           
    4 ZETASQ=ZETASQ*FAC                                                         
    1 AKEF=ZETASQ*AK                                                            
    5 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE EP (RME,RM,EPS)                                                
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE EPSILON FACTORS USED IN CALCULATING THE REDUCED ENERGY     
C                                                                               
      EPS=4.0D0*RME/(RM*(1.0D0+RME/RM)**2)                                      
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE EXFRDT(AKAY,AMU,BETSQ,CKBYZ,DTDX,EPS,RM,EXPOT,U)               
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE MEAN EXCITATION POTENTIAL FROM THE VALUE OF DT/DX          
C                                                                               
      DD=DSQRT(1.0D0-BETSQ)                                                     
      TCOR=0.5D0*AMU*RM*BETSQ                                                   
      TREL=TCOR/DD                                                              
      UCOR=DEXP(DTDX*TCOR/AKAY+BETSQ+CKBYZ)                                     
      U=UCOR*DD                                                                 
      EXPOT=EPS*TREL*1.0D6/U                                                    
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE EXFRZ(Z,EXPOT)                                                 
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE MEAN EXCITATION POTENTIAL FROM THE CHARGE OF THE TARGET    
C                                                                               
      EXPOT=9.1D0*Z*(1.0D0+1.9D0/Z**(2.0D0/3.0D0))                              
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE FOLD(FLTOP,SIGT,KREF,FWHM,FWTM)                                
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     ******************************************************************        
C     PERFORMS A CONVOLUTION INTEGRAL FOR A FLAT-TOP DISTRIBUTION AND A         
C     GAUSSIAN AND THEN FINDS THE FULL WIDTH AT HALF MAXIMUM AND THE            
C     FULL WIDTH AT ONE TENTH MAXIMUM BY AN ITERATIVE PROCEEDURE USING          
C     NEWTON'S APPROXIMATION                                                    
C     ******************************************************************        
C                                                                               
C       THIS CALCULATION DOES NOT REPRODUCE THE CORRECT SHAPE OF THE            
C     ENERGY DISTRIBUTIONS DUE TO SCATTERING OF A PARTICLE FROM A TARGET        
C     BUT IT IS A GOOD APPROXIMATION UNLESS THE TARGET IS VERY THIN WHEN        
C     THE DISTRIBUTION TENDS TOWARDS THE LANDAU-SYMON SHAPE OR UNLESS           
C     THE TARGET IS VERY THICK AT LOW INCIDENT ENERGIES WHEN DT/DX IS           
C     A VERY RAPIDLY VARYING FUNCTION OF ENERGY                                 
C                                                                               
      IF(SIGT.LE.0.0D0.AND.FLTOP.LE.0.0D0) GO TO 101                            
      IF(SIGT.LT.1.0D-2*FLTOP) GO TO 102                                        
      R=FLTOP/SIGT                                                              
      IF(R.LT.1.0D-2) GO TO 100                                                 
      PI=3.14159265D0                                                           
      GAUSS=1.1762D0                                                            
      GAUSS2=GAUSS**2                                                           
      CON=DSQRT(2.0D0/PI)                                                       
      G=DSQRT(2.0D0)                                                            
      S=R/2.0D0                                                                 
      IF(KREF.EQ.1) S=R                                                         
      Y=S/G                                                                     
  200  ERFB=DERF(Y)                                                             
C                                                                               
C     ERFB IS THE HEIGHT AT HALF MAXIMUM                                        
C                                                                               
C     CRUDE ESTIMATES OF THE POSITION OF THE FWHM                               
C                                                                               
   10 IF(R-0.2D0)11,12,12                                                       
   11 Z=1.18D0                                                                  
      GO TO 20                                                                  
   12 IF(R-0.4D0)13,14,14                                                       
   13 F=0.8D0                                                                   
      GO TO 18                                                                  
   14 IF(R-1.5D0) 15,16,16                                                      
   15 F=0.7D0                                                                   
      GO TO 18                                                                  
   16 IF(R-5.0D0) 17,19,19                                                      
   17 F=0.8D0                                                                   
   18 Z=DSQRT(F**2*R**2/4.0D0+GAUSS2)                                           
      GO TO 20                                                                  
   19 Z=R/2.0D0                                                                 
C                                                                               
C     CALCULATION OF CONVOLUTION INTEGRAL                                       
C                                                                               
   20 CONTINUE                                                                  
      IF(Z.LT.S) Z=S                                                            
      KLOK=0                                                                    
  120 CONTINUE                                                                  
      Y=(Z+S)/G                                                                 
   21 ERFU=DERF(Y)                                                              
      DERFU=(DERF(Y+1.0D-6)-ERFU)                                               
      Y=0.0D0                                                                   
      IF(Z.EQ.S) GO TO 222                                                      
      Y=(Z-S)/G                                                                 
  222 CONTINUE                                                                  
   22 ERFD=DERF(Y)                                                              
      DERFD=DERF(Y+1.0D-6)-ERFD                                                 
C                                                                               
C     USE NEWTON'S APPROX.                                                      
C                                                                               
      DUM=(DERFU-DERFD)*1.0D6                                                   
      IF(DABS(DUM).LT.1.0D-8) DUM= DSIGN(1.0D-8,DUM)                            
      Z1=Z-(ERFU-ERFD-ERFB)/DUM                                                 
      IF(DABS(1.0D0-Z1/Z).LT.1.0D-6) GO TO 30                                   
      KLOK=KLOK+1                                                               
      Z=Z1                                                                      
      IF(KLOK.GT.10) GO TO 30                                                   
      GO TO 120                                                                 
   30 FWHM=2.0D0*Z1*SIGT                                                        
      IF(KREF.EQ.1) FWHM=FWHM/2.0D0                                             
C                                                                               
C     NOW FIND FWTM                                                             
C     TRY Z = Z1+1.0.  ABOUT RIGHT FOR GAUSSIAN                                 
C                                                                               
      Z=Z1+1.0D0                                                                
      KLOK=0                                                                    
   33 CONTINUE                                                                  
      Y=(Z+S)/G                                                                 
   31 ERFU=DERF(Y)                                                              
      DERFU=DERF(Y+1.0D-6)-ERFU                                                 
      Y=0.0D0                                                                   
      IF(Z.EQ.S) GO TO 232                                                      
      Y=(Z-S)/G                                                                 
  232 CONTINUE                                                                  
   32 ERFD=DERF(Y)                                                              
      DERFD=DERF(Y+1.0D-6)-ERFD                                                 
C                                                                               
C     USE NEWTON'S APPROX.                                                      
C                                                                               
      DUM=(DERFU-DERFD)*1.0D6                                                   
      IF(DABS(DUM).LT.1.0D-8) DUM= DSIGN(1.0D-8,DUM)                            
      Z1=Z-(ERFU-ERFD-ERFB/5.0D0)/DUM                                           
      IF(DABS(1.0D0-Z1/Z).LT.1.0D-6) GO TO 40                                   
      KLOK=KLOK+1                                                               
      Z=Z1                                                                      
      IF(KLOK.GT.10) GO TO 40                                                   
      GO TO 33                                                                  
   40 FWTM=2.0D0*Z1*SIGT                                                        
      IF(KREF.EQ.1) FWTM=FWTM/2.0D0                                             
      GO TO 101                                                                 
  100 FWHM=2.355D0*SIGT                                                         
      FWTM=1.824D0*FWHM                                                         
      GO TO 101                                                                 
  102 FWHM=FLTOP                                                                
      FWTM=FWHM                                                                 
  101 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE GAMMA(B,GAM)                                                   
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE RELATIVISTIC FACTOR GAMMA                                  
C                                                                               
      DD=1.0D0 -B**2                                                            
      IF(DD.LT.0.0D0) DD=0.0D0                                                  
      GAM=1.0D0/DSQRT(DD)                                                       
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE GAMMA3(C2,DT ,EPS,EPSP,T   ,U,GAM3,TRUFAC)                     
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE REDUCED THIRD MOMENT GAMMA 3 FOR THE ENERGY LOSS           
C     DISTRIBUTIONS AND THE RATIO OF TRUE TOCALCULATED FWHM                     
C                                                                               
      DIMENSION FACTOR(20)                                                      
      DATA FACTOR      /0.99D0,0.985D0,0.98D0,0.97D0,0.96D0,0.94D0,0.925        
     1D0,0.91D0,0.88D0,0.83D0,0.78D0,0.74D0,0.68D0,0.62D0,0.575D0,0.53D0        
     2,0.51D0,0.48D0,0.46D0,0.43D0/                                             
      GAM3=DSQRT(C2*EPS*DLOG10(U)*5.0D-3*T/(DT*EPSP))                           
      IF(GAM3.GE.2.0D0) GO TO 2                                                 
      IF(GAM3.LT.1.0D-1) GO TO 1                                                
      L=IDINT(10.D0*GAM3)                                                       
      TRUFAC=FACTOR(L)+(10.0D0*GAM3-DFLOAT(L))*(FACTOR(L+1)-FACTOR(L))          
      GO TO 5                                                                   
    1 TRUFAC=1.0D0                                                              
      GO TO 5                                                                   
    2 TRUFAC=1.0D0/(1.1776D0*GAM3)                                              
    5 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE PFROMT(AMU,RM,T,P)                                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C     CALCULATES THE MOMENTUM FROM THE KINETIC ENERGY                           
C                                                                               
      P=DSQRT(2.0D0*AMU*RM*DABS(T)+T**2)                                        
      IF(T.EQ.0.0D0) P=0.0D0                                                    
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE PFROMW(AMU,RM,W,P)                                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE MOMENTUM FROM THE TOTAL ENERGY                             
C                                                                               
      DD=W**2-(AMU*RM)**2                                                       
      IF(DD.LT.0.0D0) DD=0.0D0                                                  
      P=DSQRT(DD)                                                               
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE PVFRMT(AMU,RM,T,BETSQ,PV,PVSQ)                                 
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
CC     CALCULATES THE MOMENTUM-VELOCITY PRODUCT FROM THE KINETIC ENERGY         
C                                                                               
      BETSQ=1.0D0-1.0D0/(1.0D0+T/(AMU*RM))**2                                   
       PVSQ=BETSQ*(T**2+2.0D0*T*AMU*RM)                                         
      PV=DSQRT(PVSQ)                                                            
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE QVAL(AMU,RM1,RM2,RM3,RM4,QR)                                   
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C     CALCULATES THE Q-VALUE OF A REACTION                                      
C                                                                               
      QR=AMU*(RM1+RM2-RM3-RM4)                                                  
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE SIGMA(B,GAM,BETANC,ANGC,RATIO)                                 
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE RATIO OF THE CM. TO LAB. DIFFERENTIAL CROSS SECTIONS       
C                                                                               
      ANG=ANGC*1.74532925D-02                                                   
      X=B/BETANC                                                                
      DF=GAM*(1.0D0+X*DCOS(ANG))                                                
      D=(DSIN(ANG)**2+(GAM*(X+DCOS(ANG)))**2)**3                                
      IF(D.LE.0.0D0) GO TO 1                                                    
      RATIO=DF/DSQRT(D)                                                         
      GO TO 2                                                                   
    1 RATIO=0.0D0                                                               
    2 RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE TFROMP(AMU,RM,P,T)                                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE KINETIC ENERGY FROM THE MOMENTUM                           
C                                                                               
      T=-RM*AMU+DSQRT((AMU*RM)**2+P**2)                                         
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE SIG2(BETSQ,C1,DTX,EPS,EPSP,T,U,SIG)                            
      IMPLICIT REAL *8 (A-H,O-Z)                                                
C                                                                               
C     CALCULATES THE VARIANCE OF THE ENERGY LOSS DISTRIBUTIONS ACCORDING        
C     TO THE WORK OF TSCHALAR                                                   
C                                                                               
      X=DTX/T                                                                   
      DL=DLOG10(U)                                                              
C                                                                               
C                                                                               
C     THE FOLLOWING CORRECTION FACTORS ARE TAKEN FROM TSCHALAR( REF.  ).        
C     ALFAC IS THE CORRECTION TO  THE SECOND MOMENT FROM RELATIVISTIC           
C     EFFECTS SEE FIG20. NOTE THAT ALPHA'=SIGMA'/E**2 AND ALPHA=SIGMA/T**2      
C     THEREFORE SIGMA'/SIGMA=ALPHA'*E**2/(ALPHA*T**2)                           
C     SINCE  E/T=0.5*BETASQ/(1.0-SQRT(1.0-BETASQ)) WE HAVE A FACTOR EBYTSQ      
C     DTFAC IS A CORRECTION FACTOR FOR C1 WHICH DEPENDS UPON THE FRACTION OF    
C     ENERGY LOST. SEE TSCHALAR FIG  .                                          
C     EBYTSQ IS EXACT. ALFAC AND DTFAC ARE SIMPLE POLYNOMIAL FITS TO THE        
C     APPROPRIATE CURVES.                                                       
C                                                                               
      ALFAC=1.0D0+  (0.7D0-X)*(0.95D0*BETSQ+BETSQ**3)                           
      EBYTSQ=0.5D0*BETSQ/(1.0D0-DSQRT(1.0D0-BETSQ))                             
       EBYTSQ=EBYTSQ*EBYTSQ                                                     
      Y=1.0D-4*X                                                                
      DTFAC=C1+(3.3D0-3.8D0/(DL+1.05D0))*Y+(2.512D1-3.019D1/(DL+0.5D0))*        
     1Y*X*X                                                                     
      SIG=DTFAC*T*DTX*EPS/(EPSP*DL)*ALFAC*EBYTSQ                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE TFROMW(AMU,RM,W,T)                                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE KINETIC ENERGY FROM THE TOTAL ENERGY                       
C                                                                               
      T=W-RM*AMU                                                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE THRESH(AMU,RM1,RM2,T1L,T2L,QR,TS,IF)                           
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES REACTION THRESHOLDS                                            
C                                                                               
      IF=0                                                                      
       SIGN=1.0D0                                                               
      IF(T2L.LE.0.0D0) SIGN=-1.0D0                                              
      QT=-QR                                                                    
      AY=2.0D0*QT*AMU*(RM1+RM2)+QT**2                                           
      TTI=DABS(T2L)                                                             
      W2=RM2*AMU+TTI                                                            
      A=4.0D0*RM2**2*AMU**2                                                     
      B=-4.0D0*(2.0*RM1*RM2*TTI*AMU**2+AY*W2)                                   
      C=(2.0D0*RM1*TTI*AMU-AY)**2                                               
      DD=B**2-4.0D0*A*C                                                         
      IF(DD.LT.0.0D0) DD=0.0D0                                                  
       TS=(-B+SIGN*DSQRT(DD  ))/(2.0D0*A)                                       
      IF(T1L.LT.TS) IF=1                                                        
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE TNS(AMU,RM1,RM2,RM ,TS,T2L,QR,TN)                              
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE LAB. ENERGIES OF REACTION PRODUCTS AT THE                  
C     REACTION THRESHOLD                                                        
C                                                                               
      TN =AMU*RM *(TS+T2L-QR)/((RM1+RM2)*AMU+QR)                                
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
       SUBROUTINE UU(AMU,EPS,RM,BETSQ,EXPOT,U)                                  
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
C                                                                               
C     CALCULATES THE REDUCED ENERGY OF THE PARTICLE                             
C                                                                               
      TREL=0.5D0*RM*AMU*BETSQ/DSQRT(1.0D0-BETSQ)                                
      U=EPS*TREL/EXPOT*1.0D6                                                    
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE WNCOM(AMU,RMA,RMB,WTCM,WNCM)                                   
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE CM. TOTAL ENERGY OF A PARTICLE GIVEN THE CM. TOTAL         
C     ENERGY OF THE SYSTEM                                                      
      WNCM=(WTCM**2+AMU**2*(RMA**2-RMB**2))/(2.0D0*WTCM)                        
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE WCTOL1(B,GAM,BETANC,WNCM,ANGC,WNL)                             
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE LAB. TOTAL ENERGY FROM THE CM. ANGLE                       
C                                                                               
      ANG=ANGC*1.74532925D-02                                                   
      WNL=WNCM*GAM*(1.0D0+B*BETANC*DCOS(ANG))                                   
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE WCTOL2(AMU,RM,B,GAM,WNCM,ANGL,WNLA,WNLB)                       
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE LAB. TOTAL ENERGIES FROM THE LAB. ANGLE                    
C                                                                               
      SIGN=1.0D0                                                                
      IF(DABS(ANGL).GT.9.0D1) SIGN=-1.0D0                                       
      ANG=ANGL*1.74532925D-02                                                   
      CS=DCOS(ANG)                                                              
      AA=B**2*CS**2-1.0D0                                                       
      BB=-2.0D0*WNCM                                                            
      CC=WNCM**2+(RM*AMU*GAM*B*CS)**2                                           
      DD=BB**2+4.0D0*AA*CC                                                      
      IF(DD.LT.0.0D0) DD=0.0D0                                                  
      WNLA=(BB-SIGN*DSQRT(DD))/(2.0D0*GAM*AA)                                   
      WNLB=(BB+SIGN*DSQRT(DD))/(2.0D0*GAM*AA)                                   
      RETURN                                                                    
      END                                                                       
C                                                                               
C     ******************************************************************        
C                                                                               
      SUBROUTINE WT(AMU,RM1,RM2,T1L,T2L,GAM,WTLAB,WTCM)                         
      IMPLICIT REAL*8(A-H,O-Z)                                                  
C                                                                               
C     CALCULATES THE LAB.AND CM. TOTAL ENERGIES OF A SYSTEM                     
C                                                                               
      WTLAB=AMU*RM1+DABS(T1L)+AMU*RM2+DABS(T2L)                                 
      WTCM=WTLAB/GAM                                                            
      RETURN                                                                    
      END                                                                       
      DOUBLE PRECISION FUNCTION DERF(X)
      IMPLICIT REAL*8 (A-H,O-Z)
      P=0.327590011D0
      A=0.254829592D0
      B=-0.284496736D0
      C=+1.421413741D0
      D=-1.453152027D0
      E=1.061405429D0
       X2=X*X
      T=1.0D0/(1.0D0+P*X)
      DERF=1.0D0-(T*(A+T*(B+T*(C+T*(D+T*E)))))*DEXP(-X2)
      RETURN
      END
