      PROGRAM RPS
C Modified 3/10/95 W Rae
C THIS PROGRAM ALSO NEEDS THE COMMON BLOCK IN .rpscom WHICH IS
C INCLUDED AT COMPILE STAGE BY 'INCLUDE' STATEMENTS
C----------------------------------------------------------------------
C      Written Jan 1985   SCA
C
C      Last modified :
C       10 July 1985 to use same sign conventions as version on PSFGENDAT
C       28 Oct  1985 to allow absorber on either detector
C                    and choice of angles in log file
C        9 Nov  1985 to allow absorbers for protons,deuterons & tritons
C       19 Feb  1987 moved to PSFGENDAT
C          July 1987 converted to GEC Fortran-77 (see CGEC comments
C                    for non-standard bits) and RPS2 incorporated
C                    (allows restrictions on thetastar & psi)
C
C      Function: three-body kinematics for sequential breakup reactions
C      Given reaction & detector settings, converts from thetastar, psi
C      to lab.angles
C
C      Output  : log file and/or line printer plots
C
CGEC   Output is formatted for 132 columns x 80 lines/page
C
C-----------------------------------------------------------------------
C
C   Common data for program RPS
C
CGEC   NB: This is 'rpscom' rearranged to suit GEC
CGEC        - must change anything here as well as in 'rpscom'
C
      PARAMETER (NXSET=64, NYSET=45, NTH=40, MAXANG=20, MAXEN=20)
C
      CHARACTER*1 CHARS(15)
      CHARACTER*6 PTCLE(3)
      CHARACTER*20 CONTF, LOGF, PLOTF
      CHARACTER*50 NEWF
      CHARACTER*1 SYMBOL, ICR, OPT, IZER, IELIM, OVWRT
C
C    Output arrays:
      CHARACTER*1 HIT(NXSET*NYSET), NOHIT(NXSET*NYSET)
      CHARACTER*1 TTHIT(0:NTH*NTH-1), TTNOHT(0:NTH*NTH-1)
C
      REAL ABSRBS(3)
C
      INTEGER IDET(2), MASSES(3)
      INTEGER PLOTU, LOGU, CONTU, NCONTU, SCRTCU
      INTEGER ICHO
C
      DIMENSION READNO(2*MAXANG)
      DIMENSION AMIN(MAXANG), AMAX(MAXANG)
      DIMENSION BMIN(MAXANG), BMAX(MAXANG)
C
      DIMENSION ENERGS(MAXEN)
      DIMENSION ENERGY(3), THETA(3), ABSOLD(3)
C
      LOGICAL LOVL, LZER, LLOPEN
      LOGICAL TEST(4), FLAG, LCORR(3), LFIRST
      LOGICAL LPAGE, OK, LEXIST, LPOPEN, LLIST, LOPN, LELIM, THPOK
C
C
      EQUIVALENCE (MASSES(1),MHI)
C
      EQUIVALENCE (ENERGY(1),EHICOR), (ENERGY(2),ELICOR),               
     +(ENERGY(3),ERECOR)
      EQUIVALENCE (THETA(1),THHI), (THETA(2),THLI),                     
     +(THETA(3),THREC)
C
      COMMON/PARMS/EBEAM,MB,MT,MHI,MLI,MREC,Q2,EX,EXI,EXINC,EXF,        
     +THRESH,DISTA,DIAMA,DISTB,DIAMB,                                   
     +EAMIN,EBMIN,EAMAX,EBMAX,ABSRBS,THSMIN,PSIMIN,                     
     +IVAR,ELIMAX,IDET,THMINA,THMINB,                                   
     +DTH,DPSI,DTHA,DTHB,LOVL,                                          
     +THAMIN,THAMAX,THBMIN,THBMAX
C
      DATA PTCLE/'    HI','    LI','RECOIL'/
      DATA CHARS/'A','B','C','D','0','1','2','3','4','5','6',           
     +'7','8','9','I'/
C
      DATA ABSOLD/0.0,0.0,0.0/
      DATA LPOPEN/.FALSE./
      DATA LLOPEN/.FALSE./
C
C--------------------------------------------------------------
C
      RADIAN = 90.0/ASIN(1.0)
C
C    Defaults:
      DTH = 1.0
      DPSI = 2.0
      THMIN = 0.0
      PSIMIN = 0.0
C
CGEC Unit numbers: (plot file, log file, control file, new control
CGEC                file and scratch work file repectively)
      PLOTU = 12
      LOGU  = 11
      CONTU = 7
      NCONTU= 9
      SCRTCU= 10
C
C    Set up parameters
C
CGEC
    4 CONTINUE
      WRITE(6,5)
    5 FORMAT('Create or Read a control file ? ',$)
      READ(5,10) ICR
   10 FORMAT(A1)
CGEC
   11 CONTINUE
      WRITE(6,12)
   12 FORMAT('File name : ',$)
      READ(5,'(A)') CONTF
C
      INQUIRE(FILE=CONTF,EXIST=LEXIST,OPENED=LOPN) 
C
      IF(ICR.EQ.'C' .OR. ICR.EQ.'c') THEN
         IF (LEXIST) THEN
CGEC  Have to check whether file already exists - no version numbers
            CONTINUE
            WRITE(6,13) CONTF
   13       FORMAT('File ',A20,' already exists; overwrite? (Y/N) : ',$)
            READ(5,10) OVWRT
            WRITE(6,*) ' '
            IF (OVWRT.EQ.'Y' .OR. OVWRT.EQ.'y') THEN
               IF (.NOT.LOPN)
     +           OPEN (UNIT=CONTU, FILE=CONTF, STATUS = 'OLD')
               CLOSE (UNIT=CONTU, STATUS = 'DELETE')
            ELSE
               CONTF = '                    '
               GOTO 11
            ENDIF
         ENDIF
C
         CALL CREATF 
* create new control file
         CALL WRITEF(CONTU,CONTF,LOGF) 
* save it
         ELSE IF(ICR.EQ.'R' .OR. ICR.EQ.'r') THEN 
         IF (LEXIST) THEN 
            CALL READF(CONTF) 
* read existing file
            WRITE(6,46) 
            CALL WRITEF(6,CONTF,LOGF) 
* write parms on screen
            ELSE 
            WRITE(6,14) CONTF 
   14       FORMAT(' File ',A20,' does not exist;',                           
     +       ' please try again')
            CONTF = '                    '
            GOTO 11 
         ENDIF 
         ELSE 
         GOTO 4 
      ENDIF
C
C
C   Options to change parameters or perform calculation
C
   20 WRITE(6,*) ' '
CGEC
      CONTINUE
      WRITE(6,21)
   21 FORMAT('OPTION (type Z for list) : ',$)
      READ(5,10) OPT
      WRITE(6,*) ' '
C
      IF(OPT.EQ.'Z' .OR. OPT.EQ.'z') THEN
         WRITE(6,*)'The responses to OPTION are :' 
         WRITE(6,*)'       beam (E)nergy and initial masses'
         WRITE(6,*)'      (M)asses of decay products'
         WRITE(6,*)'      (Q)-values'
         WRITE(6,*)'      (D)etected particles'
         WRITE(6,*)'      (A) distance & diameter'
         WRITE(6,*)'      (B) distance & diameter'
         WRITE(6,*)'      (R)ange-energy corrections'
         WRITE(6,*)'      (T)heta star and psi min. & bin sizes'
         WRITE(6,*)'      (L)ab. detector angles displayed on plot'
         WRITE(6,*)'      (P)arameter listing'
         WRITE(6,*)'      (W)rite out new control file'
         WRITE(6,*)'      (C)alculation'
         WRITE(6,*)'      e(X)it'
C
         ELSE IF (OPT.EQ.'X' .OR. OPT.EQ.'x') THEN 
         GOTO 350 
C
         ELSE IF (OPT.EQ.'E' .OR. OPT.EQ.'e') THEN 
CGEC
         CONTINUE
         WRITE(6,24) 
   24    FORMAT('Ebeam, Mbeam, Mtarget : ',$)
         READ(5,*) EBEAM,MB,MT 
C
         ELSE IF (OPT.EQ.'M' .OR. OPT.EQ.'m') THEN 
CGEC
   25    CONTINUE
         WRITE(6,26) 
   26    FORMAT('MHI, MLI : ',$)
         READ(5,*) MHI,MLI 
C        Check recoil mass and set up MASSES for absorber determination
         MREC = MB + MT - MHI - MLI 
         IF (MREC.LE.0) THEN 
            WRITE(6,*) ' '
            WRITE(6,*) 'Zero or negative recoil mass ; please try again'
            GOTO 25 
         ENDIF 
C
         ELSE IF (OPT.EQ.'Q' .OR. OPT.EQ.'q') THEN 
CGEC
         CONTINUE
         WRITE(6,28) 
   28    FORMAT('Two-body Q-value, breakup threshold : ',$)
         READ(5,*) Q2,THRESH 
C          Q2 is Q-value for beam+target->R.P.(g.s.)+recoil
C
         ELSE IF (OPT.EQ.'D' .OR. OPT.EQ.'d') THEN 
CGEC
  450    CONTINUE
         WRITE(6,29) 
   29    FORMAT('Particles detected in telescopes A & B (1=HI, 2=LI,',    
     +    ' 3=recoil) : ',$)
         READ(5,*) IDET 
C          IDET(1) is ptcle in telescope A, IDET(2) ptcle in B
         DO 1004 I=1,2 
            IF (IDET(I).LT.1 .OR. IDET(I).GT.3) THEN
               WRITE(6,*) 'Invalid value'
               GOTO 450
            ENDIF
            IF (MASSES(IDET(I)).GT.4) THEN 
               ABSRBS(IDET(I))=0.0 
               WRITE(6,501) IDET(I)
  501          FORMAT(' Absorber ',I1,' set to zero')
            ENDIF 
 1004    CONTINUE 
C
         ELSE IF (OPT.EQ.'A' .OR. OPT.EQ.'a') THEN 
CGEC
         CONTINUE
         WRITE(6,30) 
   30    FORMAT('Telescope A distance from target & diameter : ',$)
         READ(5,*) DISTA,DIAMA 
C
         ELSE IF (OPT.EQ.'B' .OR. OPT.EQ.'b') THEN 
CGEC
         CONTINUE
         WRITE(6,32) 
   32    FORMAT('Telescope B distance from target & diameter : ',$)
         READ(5,*) DISTB,DIAMB 
C
         ELSE IF (OPT.EQ.'R' .OR. OPT.EQ.'r') THEN 
CGEC
         CONTINUE
         WRITE(6,34) 
   34    FORMAT('Minimum & maximum A energies : ',$)
         READ(5,*) EAMIN,EAMAX 
CGEC
         CONTINUE
         WRITE(6,35) 
   35    FORMAT('Minimum & maximum B energies : ',$)
         READ(5,*) EBMIN,EBMAX 
C
         IF (MASSES(IDET(1)).LE.4) THEN 
C        i.e. detected particle A is an alpha,triton,deuteron or proton
CGEC
            CONTINUE
            WRITE(6,36) 
   36       FORMAT('A absorber thickness in microns : ',$)
            READ(5,*) ABSRBS(IDET(1)) 
         ENDIF 
         IF (MASSES(IDET(2)).LE.4) THEN 
C        i.e. detected particle B is an alpha,triton,deuteron or proton
CGEC
            CONTINUE
            WRITE(6,37) 
   37       FORMAT('B absorber thickness in microns : ',$)
            READ(5,*) ABSRBS(IDET(2)) 
         ENDIF 
C
         ELSE IF (OPT.EQ.'T' .OR. OPT.EQ.'t') THEN 
CGEC
         CONTINUE
         WRITE(6,38) 
   38    FORMAT('Thetastar minimum and increment : ',$)
         READ(5,*) THSMIN,DTH 
CGEC
         CONTINUE
         WRITE(6,39) 
   39    FORMAT('Psi minimum and increment : ',$)
         READ(5,*) PSIMIN,DPSI 
C
         ELSE IF (OPT.EQ.'L' .OR. OPT.EQ.'l') THEN 
CGEC
         CONTINUE
         WRITE(6,40) 
   40    FORMAT('A lab angles as minimum and increment : ',$)
         READ(5,*) THMINA,DTHA 
CGEC
         CONTINUE
         WRITE(6,41) 
   41    FORMAT('B lab angles as minimum and increment : ',$)
         READ(5,*) THMINB,DTHB 
C          Angles for TH(A) vs TH(B) plot
C
         ELSE IF (OPT.EQ.'P' .OR. OPT.EQ.'p') THEN 
         WRITE(6,46) 
   46    FORMAT(' Present parameters are :'/)
         CALL WRITEF(6,CONTF,LOGF) 
* write parms on screen
         WRITE(6,*)' ' 
C
         ELSE IF (OPT.EQ.'W' .OR. OPT.EQ.'w') THEN 
C          Write out new control file
CGEC
   47    CONTINUE
         WRITE(6,12) 
         READ(5,'(A)') CONTF
         INQUIRE(FILE=CONTF,EXIST=LEXIST)
         IF (LEXIST) THEN
CGEC
            CONTINUE
            WRITE(6,48) CONTF
   48       FORMAT('File ',A20,' already exists; overwrite? (Y/N) : ',$)
            READ(5,10) OVWRT
            WRITE(6,*) ' '
            IF (OVWRT.EQ.'Y' .OR. OVWRT.EQ.'y') THEN
               OPEN (UNIT=NCONTU, FILE=CONTF, STATUS='OLD')
               CLOSE (UNIT=NCONTU, STATUS='DELETE')
            ELSE
               CONTF = '                    '
               GOTO 47
            ENDIF
         ENDIF
         CALL WRITEF(NCONTU,CONTF,LOGF)
C
         ELSE IF (OPT.EQ.'C' .OR. OPT.EQ.'c') THEN 
         GOTO 50 
C
      ENDIF
C
      GOTO 20
C
   50 CONTINUE
C
C
C   Type of output
C
      WRITE(6,51)
   51 FORMAT(' Calculation option :'/' 0=original, 1=all angles on',    
     + ' same plot, 2=auto-detector, 3=log file,')
CGEC
      CONTINUE
      WRITE(6,52)
   52 FORMAT(' 4=back to main option level : ',$)
      READ(5,*) IVAR
      IF (IVAR.EQ.4) GOTO 20
      IF (IVAR.LT.0 .OR. IVAR.GT.3) GOTO 50
C
C    Angular radii of detectors (for options 0 & 1)
C
      ANGA2 = RADIAN * ATAN2(0.5*DIAMA,DISTA)
      ANGB2 = RADIAN * ATAN2(0.5*DIAMB,DISTB)
C
C    Ex and detector angles
C
      LLIST = .FALSE.
* for angles
C
C      Ex goes into ENERGS
C
      IF (IVAR.NE.3) THEN
* not log file
         WRITE(6,*) ' '
         IF (.NOT.LPOPEN) THEN
C          LPOPEN is .TRUE. if plot file is open
CGEC
  300       CONTINUE
            WRITE(6,301)
  301       FORMAT('Plot file name : ',$)
            READ(5,'(A)') PLOTF
            INQUIRE(FILE=PLOTF,EXIST=LEXIST)
            IF (LEXIST) THEN
CGEC
               CONTINUE
               WRITE(6,13) PLOTF
               READ(5,10) OVWRT
               WRITE(6,*) ' '
               IF (OVWRT.EQ.'Y' .OR. OVWRT.EQ.'y') THEN
                  OPEN(UNIT=PLOTU,FILE=PLOTF,STATUS='OLD')
                  CLOSE(UNIT=PLOTU,STATUS='DELETE')
                  ELSE
                  PLOTF = '                    '
                  GOTO 300
               ENDIF
            ENDIF
         ENDIF
C
CGEC
         CONTINUE
         WRITE(6,53) 
   53    FORMAT('Ex as initial,increment,final (MeV) ',               
     +    '(-1 for list) : ',$)
         CALL FREEFM(EXI,EXINC,EXF) 
         IF (EXI.EQ.-1.0) THEN 
* enter list of Ex
            DO 1005 I=1,MAXEN 
* (0 marks end)
               READNO(I)=0.0 
 1005       CONTINUE 
CGEC
            CONTINUE
            WRITE(6,54) 
   54       FORMAT('List of excitation energies : ',$)
            CALL FREEF2(READNO) 
            NEN=0 
* no. of Ex
            DO 1006 I=1,MAXEN 
               IF (READNO(I).EQ.0.0) GOTO 55 
* end
               ENERGS(I)=READNO(I) 
               NEN=NEN+1 
 1006       CONTINUE 
   55       CONTINUE 
            ELSE 
* i,inc,f
            IF (EXINC.LE.0.0) EXINC = 1.0 
            IF (EXF.LT.EXI) EXF = EXI 
            NEN = (EXF - EXI)/EXINC + 1 
* no. of Ex
C           Put Ex values to be used into ENERGS
            DO 1007 I=1,NEN 
               ENERGS(I) = EXI + (I-1)*EXINC 
 1007       CONTINUE 
         ENDIF 
C
C        Detector limits go into AMIN & AMAX, BMIN & BMAX
C        Use absolute values of input
C
         IF (IVAR.NE.2) THEN 
* not auto-detector
CGEC
            CONTINUE
            WRITE(6,56) 
   56       FORMAT('Angle A as initial,increment,final  ',
     +        '(-1 for list) : ',$)
            CALL FREEFM(THAI,THAINC,THAF) 
            IF (THAI.EQ.-1.0) THEN 
* input list of angle pairs
               LLIST = .TRUE. 
               DO 1008 I=1,2*MAXANG 
                  READNO(I)=0.0 
 1008          CONTINUE 
CGEC
               CONTINUE
               WRITE(6,57) 
   57          FORMAT('List of A & B angle pairs : ',$)
C              in order (a,b), (a,b,), etc.
C              Assumes you won't enter zero - this is used to mark end
               CALL FREEF2(READNO) 
               NANGA=0 
* no. of A angles
               DO 1009 I=1,MAXANG 
                  IF (READNO(2*I-1).EQ.0.0) GOTO 60 
                  AMIN(I) = ABS(READNO(2*I-1)) - ANGA2 
                  AMAX(I) = ABS(READNO(2*I-1)) + ANGA2 
                  BMIN(I) = ABS(READNO(2*I)) - ANGB2 
                  BMAX(I) = ABS(READNO(2*I)) + ANGB2 
                  NANGA=NANGA+1 
 1009          CONTINUE 
   60          CONTINUE 
               NANGB=1 
* to give MAXTH=NANGA
               ELSE 
               THAI = ABS(THAI) 
               THAF = ABS(THAF) 
               IF (THAINC.LE.0.0) THAINC = 1.0 
               IF (THAF.LT.THAI) THAF = THAI 
               NANGA = (THAF - THAI)/THAINC + 1 
               DO 1010 I=1,NANGA 
                  AMIN(I) = THAI + (I-1)*THAINC - ANGA2 
                  AMAX(I) = AMIN(I) + 2*ANGA2 
 1010          CONTINUE 
CGEC
               CONTINUE
               WRITE(6,64) 
   64          FORMAT('Angle B as initial,increment,final     ',
     +             '           : ',$)
               CALL FREEFM(THBI,THBINC,THBF) 
               THBI = ABS(THBI) 
               THBF = ABS(THBF) 
               IF (THBINC.LE.0.0) THBINC = 1.0 
               IF (THBF.LT.THBI) THBF = THBI 
               NANGB = (THBF - THBI)/THBINC + 1 
               DO 1011 I=1,NANGB 
                  BMIN(I) = THBI + (I-1)*THBINC - ANGB2 
                  BMAX(I) = BMIN(I) + 2*ANGB2 
 1011          CONTINUE 
            ENDIF 
C
            ELSE 
* auto-detector
            DO 1012 I=1,2*MAXANG 
               READNO(I)=0.0 
 1012       CONTINUE 
CGEC
            CONTINUE
            WRITE(6,76) 
   76       FORMAT('A detector boundaries : ',$)
C           Zero used to mark end of input again
            CALL FREEF2(READNO) 
            NANGA=0 
            DO 1013 I=1,MAXANG 
               IF (READNO(I+1).EQ.0.0) GOTO 80 
* end of input
               AMIN(I) = ABS(READNO(I)) 
               AMAX(I) = ABS(READNO(I+1)) 
               NANGA=NANGA+1 
 1013       CONTINUE 
   80       CONTINUE 
C
            DO 1014 I=1,2*MAXANG 
               READNO(I)=0.0 
 1014       CONTINUE 
CGEC
            CONTINUE
            WRITE(6,86) 
   86       FORMAT('B detector boundaries : ',$)
            CALL FREEF2(READNO) 
            NANGB=0 
            DO 1015 I=1,MAXANG-1 
               IF (READNO(I+1).EQ.0.0) GOTO 90 
               BMIN(I) = ABS(READNO(I)) 
               BMAX(I) = ABS(READNO(I+1)) 
               NANGB=NANGB+1 
 1015       CONTINUE 
   90       CONTINUE 
         ENDIF 
C
         WRITE(6,*) ' '
CGEC
         CONTINUE
         WRITE(6,92) 
   92    FORMAT('Reasons for elimination or not (Y/N) : ',$)
         READ(5,10) IELIM 
         IF (IELIM.EQ.'Y' .OR. IELIM.EQ.'y') THEN
           LELIM = .TRUE.
           ELSE
           LELIM = .FALSE.
         ENDIF
* second plot?
         ELSE 
* log file
         WRITE(6,*) ' '
         IF (.NOT.LLOPEN) THEN
CGEC
  303       CONTINUE
            WRITE(6,305)
  305       FORMAT('Log file name : ',$)
            READ(5,'(A)') LOGF
            INQUIRE(FILE=LOGF,EXIST=LEXIST)
            IF (LEXIST) THEN
CGEC
               CONTINUE
               WRITE(6,13) LOGF
               READ(5,10) OVWRT
               WRITE(6,*) ' '
               IF (OVWRT.EQ.'Y' .OR. OVWRT.EQ.'y') THEN
                  OPEN(UNIT=LOGU,FILE=LOGF,STATUS='OLD')
                  CLOSE(UNIT=LOGU,STATUS='DELETE')
                  ELSE
                  LOGF = '                    '
                  GOTO 303
               ENDIF
            ENDIF
            ELSE
            WRITE(LOGU,230)
C           new page
         ENDIF
C
CGEC
         CONTINUE
         WRITE(6,94) 
   94    FORMAT('Ex (MeV) : ',$)
         READ(5,*) ENERGS(1) 
         NEN=1 
         NANGA=1 
         NANGB=1 
C          Choose angles for log file
CGEC
         CONTINUE
         WRITE(6,96) 
   96    FORMAT('Thetastar initial, increment, final : ',$)
         CALL FREEFM(THSMIN,DTH,THSMAX) 
         IF (DTH.LE.0.0) DTH = 1.0 
CGEC
         CONTINUE
         WRITE(6,98) 
   98    FORMAT('Psi       initial, increment, final : ',$)
         CALL FREEFM(PSIMIN,DPSI,PSIMAX) 
         IF (DPSI.LE.0.0) DPSI = 1.0 
      ENDIF
C
C    Choose restriction on thetastar, psi
   99 WRITE(6,*)' '
      WRITE(6,101) 
  101 FORMAT('No restrictions on angles, or restrict calculation to')
CGEC
      CONTINUE
      WRITE(6,102)
  102 FORMAT('psi = -thstar, psi = 180 deg, thstar = 0 (0/1/2/3) : ',$)
      READ(5,*) ICHO
      IF (ICHO.LT.0 .OR. ICHO.GT.3) THEN
         WRITE(6,*) 'Invalid value'
         GOTO 99
      ENDIF
C
      IF (ICHO.EQ.1 .OR. ICHO.EQ.2) THEN
C     Restrict values of psi
         IPSI1 = -1
         IPSI2 = 1
      ELSE IF (IVAR.EQ.3) THEN
         IPSI1 = 1
         IPSI2 = (PSIMAX - PSIMIN)/DPSI + 1
         PSIMIN = PSIMIN - 0.5*DPSI
      ELSE
         IPSI1 = 1
         IPSI2 = NXSET
      ENDIF
      IF (ICHO.EQ.3) THEN
C     Restrict values of thetastar
         ITH1 = -1
         ITH2 = 1
      ELSE IF (IVAR.EQ.3) THEN
         ITH1 = 1
         ITH2 = (THSMAX - THSMIN)/DTH + 1
         THSMIN = THSMIN - 0.5*DTH
      ELSE
         ITH1 = 1
         ITH2 = NYSET
      ENDIF
C
C    Should we print empty plots? - should only arise for option 0
      LZER=.FALSE.
      IF (IVAR.EQ.0) THEN
         WRITE(6,*)' '
CGEC
         CONTINUE
         WRITE(6,59) 
   59    FORMAT('Print zero hits or not (Y/N) : ',$)
         READ(5,10) IZER 
         IF (IZER.EQ.'Y' .OR. IZER.EQ.'y') LZER=.TRUE.
      ENDIF
C
      IF (IVAR.EQ.3) THEN
* open log file 
         IF (.NOT.LLOPEN) THEN
            CALL RENAME(LOGF,NEWF)
CGEC        Listing filetype to get carriage control
            OPEN(UNIT=LOGU,FILE=LOGF,ERR=401,FORM='PRINT')
            LLOPEN = .TRUE.
         ENDIF
         CALL WRITEF(LOGU,LOGF,LOGF) 
* write out parms
         WRITE(6,*) ' ' 
         WRITE(6,*) 'Now just wait a couple of minutes...' 
         WRITE(6,*) ' ' 
      ENDIF
C
C
      M = MB + MT
* total mass
C
      MRP = MHI + MLI
* mass of resonant particle
C
      MREC = M - MRP
* mass of recoil nucleus
C
C----------------------Loop over Ex--------------------------------
C
      DO 1016 IEX=1,NEN
         EX = ENERGS(IEX) 
C
         HITSUM=0 
* Total no. of hits
         LOVL=.FALSE. 
* Overlaps?
C
C    Erel
         EREL = EX - THRESH 
C
         IF (EREL.LE.0.0) THEN
            WRITE(6,*) ' ' 
            IF (EREL.LT.0.0) WRITE(6,*) 'Erel negative' 
            IF (EREL.EQ.0.0) WRITE(6,*) 'Erel = 0'
            IF (IVAR.EQ.3) WRITE(LOGU,*)' Ex below breakup threshold;',         
     +       ' calculation aborted'
            GOTO 20 
         ENDIF 
C
         VREL = SQRT ( 2.0 * EREL * MRP/(MLI * MHI) ) 
C
         VRELHI = VREL * MLI/MRP 
C
         VRELLI = VREL * MHI/MRP 
C
C    Velocity of c.m. frame w.r.t. lab.
C
         IF (EBEAM.LE.0.0) THEN 
            WRITE(6,*) ' '
            WRITE(6,*) ' Negative beam energy' 
            IF (IVAR.EQ.3) WRITE(LOGU,*) ' Beam energy negative: ',             
     +       'calculation aborted'
            GOTO 20 
         ENDIF 
C
         VCM = (SQRT(2.*MB*EBEAM)) / M 
C
C    Velocity of resonant particle in c.m. frame
C
          EICM = MT * EBEAM / M 
* initial energy in c.m.
         EFCM = EICM + Q2 - EX 
* final energy in c.m.
         IF(EFCM.LE.0.0) THEN 
            WRITE(6,*)' '
            WRITE(6,*)'Not enough c.m. energy',
     +        ' to excite this state in r.p.'
            IF (IVAR.EQ.3) WRITE(LOGU,*)' Not enough c.m. energy to ',          
     +        'excite this state in r.p.'
            GOTO 20 
         ENDIF 
C
         VRL2CM = SQRT(2.0 * EFCM * M/(MRP * MREC)) 
C      Relative velocity of R.P. & recoil in c.m. frame
C
         VRPCM = VRL2CM * MREC/M 
* vel. of R.P. in c.m.
C
C    If absorber, calculate lookup tables
C
         DO 1017 I=1,2 
            II = IDET(I) 
            IF (MASSES(II).LE.4 .AND. ABSRBS(II).GT.0.0) THEN 
               LCORR(II)=.TRUE. 
               IF (ABSRBS(II).NE.ABSOLD(II)) THEN 
                  CALL ABSORB(ABSRBS(II),EBEAM,II,MASSES(II)) 
                  ABSOLD(II) = ABSRBS(II) 
               ENDIF 
            ENDIF 
 1017    CONTINUE 
C
C    Write out results to log file
C
         IF (IVAR.EQ.3) THEN 
            WRITE(LOGU,105) 
  105       FORMAT(1X,//)
            WRITE(LOGU,110) M,MRP,MREC 
  110       FORMAT(' Total mass =',I4,
     +        '; mass of resonant particle =',I4,'; mass of recoil =',
     +        I4/)
            WRITE(LOGU,120) EREL,VREL 
  120       FORMAT(' Erel =',F10.3,'; Vrel =',F10.3/)
            WRITE(LOGU,130) VRELHI, VRELLI 
  130       FORMAT(' Vel of HI & LI in rest frame of r.p. =',
     +        F10.3,' ,',F10.3)
            WRITE(LOGU,140) VCM 
  140       FORMAT(' Velocity of centre of mass Vcm =',F10.3/)
            WRITE(LOGU,150) VRPCM 
  150       FORMAT(' Velocity of  r.p.  in 2-body c.m. frame =',F10.3)
            WRITE(LOGU,155) VRPCM*MRP/MREC 
  155       FORMAT(' Velocity of recoil in 2-body c.m. frame =',F10.3/)
            WRITE(LOGU,160) DTH,DPSI 
  160       FORMAT(' Dth =',F6.2,'; dpsi =',F6.2/)
            DO 1018 I=1,3 
               IF (ABSRBS(I).GT.0.0) WRITE(LOGU,165) PTCLE(I),ABSRBS(I) 
  165          FORMAT(' ',A6,' absorber thickness =',F10.1,' um'/)
 1018       CONTINUE 
            WRITE(LOGU,170) ENERGS(1)
  170       FORMAT(' Ex (resonant particle) =',F10.1,' MeV')
C
            WRITE(LOGU,175) 
  175       FORMAT('1',5X,'THSTAR',8X,'PSI',5X,'TH(HI)',5X,'TH(LI)',          
     +       4X,'TH(REC)',6X,'E(HI)',6X,'E(LI)',5X,'E(REC)',                   
     +       6X,'E(HI)',6X,'E(LI)',5X,'E(REC)',2X,'E(A)+E(B)'                  
     +       /94X,'uncorr     uncorr     uncorr'/)
            IOUT=3 
* no. of lines printed in log file
            ELSE 
            MAXTH = NANGA*NANGB 
C          No. of loops over angle for one Ex
         ENDIF 
C
C
         ITIMTH = 0 
* no. of iterations over angle
         ICOUNT=0 
* no. of plots with some hits
         ELIMAX = 0.0 
* max. energy deposited in LI detector
C
C----------------------Loops over A & B detector angles----------------
C
         DO 1019 ITHA=1,NANGA 
            DO 1020 ITHB=1,NANGB 
C
               ITIMTH = ITIMTH + 1 
C
               IF (IVAR.NE.0 .AND. ITIMTH.NE.1) GOTO 179 
C      Only zero HIT for option 0, or first time around
               DO 1021 I=1,NXSET*NYSET 
                  HIT(I) = ' ' 
                  NOHIT(I) = '.'
                  IF (I.LE.NTH*NTH) THEN 
                     TTHIT(I-1) = ' ' 
                     TTNOHT(I-1) = '.' 
                  ENDIF 
 1021          CONTINUE 
  179          CONTINUE 
C
C    Angular range covered by each detector
C
               IF (LLIST) THEN 
                  IITHB = ITHA 
                  ELSE 
                  IITHB = ITHB 
               ENDIF 
               THAMAX = AMAX(ITHA) 
               THAMIN = AMIN(ITHA) 
               THBMAX = BMAX(IITHB) 
               THBMIN = BMIN(IITHB) 
C
               IHITS = 0 
* no. of hits
C
C------------------------Loop over thetastar---------------------------
C
               DO 1022 ITH=ITH1,ITH2
C
                  IF (ICHO.EQ.3) THEN
C                 Thetastar = 0 +/- 2 degrees
                     THSTAR = ITH * 2.0
                     IITH2  = (THSTAR - THSMIN)/DTH + 1.0
C                    altered from '+2'
                     IITH   = (IITH2 - 1) * NXSET
                  ELSE
                     THSTAR = THSMIN + (ITH-0.5)*DTH 
                     IITH = (ITH-1)*NXSET 
                  ENDIF
C          Thstar defined as +ve on HI side of beam
C
                  IF (IVAR.EQ.3) THEN 
* log file
                     WRITE(LOGU,*)' ' 
* space between each thstar
                     IOUT = IOUT+1 
                     IF (IOUT.EQ.80) THEN 
* new page - need titles
                        WRITE(LOGU,175) 
                        IOUT=3 
                     ENDIF 
                  ENDIF 
C
                  THSRAD = THSTAR/RADIAN 
C
C    Components of R.P. & recoil velocity in c.m. frame
C
                  VRPCMX = VRPCM * COS(THSRAD) 
                  VRPCMY = VRPCM * SIN(THSRAD) 
                  VRECMX = -VRPCMX * MRP/MREC 
                  VRECMY = -VRPCMY * MRP/MREC 
C
C    Convert to lab. frame
C
                  VRPLX= VRPCMX + VCM 
                  VRPLY= VRPCMY 
                  VRECLX = VRECMX + VCM 
                  VRECLY = VRECMY 
C
C    Recoil angle in the lab. (on opposite side of beam to H.I.)
C
                  THREC = RADIAN * ATAN2(VRECLY,VRECLX)
C
C    Recoil energy
                  EREC = 0.5 * MREC * (VRECLX**2 + VRECLY**2) 
                  IF (LCORR(3)) THEN 
* absorber on recoil det
                     CALL LOOKUP(EREC,ERECOR,3) 
                     ELSE 
                     ERECOR = EREC 
                  ENDIF 
C
C
C-----------------------Loop over psi--------------------------------
Coct
                  DO 1023 IPSI=IPSI1,IPSI2
                     IF (ICHO.EQ.1) THEN
C                    psi = - thetastar +/- 2 degrees
                        PSI = -THSTAR + IPSI*2.0
                        IIPSI2 = (PSI - PSIMIN)/DPSI + 1.5
                        II = IITH + IIPSI2
                     ELSE IF (ICHO.EQ.2) THEN
C                    psi = 180 +/- 2 degrees
                        PSI = 180.0 + IPSI*2.0
                        IIPSI2 = (PSI - PSIMIN)/DPSI + 1.5
                        II = IITH + IIPSI2
                     ELSE
                        PSI = PSIMIN + (IPSI-0.5)*DPSI 
                        II = IITH + IPSI 
                     ENDIF
C
                     PSIRAD = PSI/RADIAN 
C
C           NB: psi=0 when LI decays forwards
C
C          Components of Vrel(HI) & Vrel(LI)
C
                     VRHIX = -VRELHI * COS(PSIRAD) 
                     VRHIY = VRELHI * SIN(PSIRAD) 
                     VRLIX = VRELLI * COS(PSIRAD) 
                     VRLIY = -VRELLI * SIN(PSIRAD) 
C
C          Lab. velocities of HI & LI
C
                     VHIX = VRPLX + VRHIX 
                     VHIY = VRPLY + VRHIY 
                     VLIX = VRPLX + VRLIX 
                     VLIY = VRPLY + VRLIY 
C
                     DO 1024 I=1,4 
                        TEST(I) = .TRUE. 
* one element for each test
* energy & angle, A & B
 1024                CONTINUE 
C
C          Calculate theta(HI), theta(LI), E(HI), E(LI)
C
                     THHI = RADIAN * ATAN2(VHIY,VHIX)
                     THLI = RADIAN * ATAN2(VLIY,VLIX)
                     EHI = 0.5 * MHI * (VHIX**2 + VHIY**2) 
                     ELI = 0.5 * MLI * (VLIX**2 + VLIY**2) 
C
C          Correct energies for absorber if necessary
C
                     IF(LCORR(1)) THEN 
* absorber on "HI" detector
                        CALL LOOKUP(EHI,EHICOR,1) 
                        ELSE 
                        EHICOR = EHI 
                     ENDIF 
                     IF (LCORR(2)) THEN 
* absorber on LI detector
                        CALL LOOKUP(ELI,ELICOR,2) 
                        ELSE 
                        ELICOR = ELI 
                     ENDIF 
C
C          Assign required energies & angles
                     EA = ENERGY(IDET(1)) 
                     EB = ENERGY(IDET(2)) 
                     THA = THETA(IDET(1)) 
                     THB = THETA(IDET(2)) 
C	type *, ea,eb,tha,thb
C	type *, eamin,eamax,ebmin,ebmax
C	type *, thamin,thamax,thbmin,thbmax
C
C          Log file
                     IF (IVAR.EQ.3) THEN 
                        WRITE(LOGU,180) THSTAR,PSI,THHI,THLI,THREC,       
     +                  EHICOR,ELICOR,ERECOR,EHI,ELI, EREC,EA+EB 
  180                   FORMAT(1X,12F11.3)
                        IOUT = IOUT+1 
                        IF (IOUT.EQ.80) THEN 
* new page
                           WRITE(LOGU,175) 
                           IOUT=3 
                        ENDIF 
                        GOTO 195 
* no tests for log file
                     ENDIF 
C
C
C          Check each parameter against constraints
C
C
                     IF (ABS(THA).GT.ABS(THAMAX) .OR. ABS(THA).LT.
     +                   ABS(THAMIN)) TEST(1) = .FALSE. 
C
                     IF (ABS(THB).GT.ABS(THBMAX) .OR. ABS(THB).LT.
     +                   ABS(THBMIN)) TEST(2) = .FALSE. 
C
                     IF (EA.GT.EAMAX .OR. EA.LT.EAMIN) 
     +                   TEST(3) = .FALSE.
C
C          Set up for plot of TH(A) vs TH(B)
C
                     IF (EB.GT.EBMAX .OR. EB.LT.EBMIN) 
     +                   TEST(4) = .FALSE.
C	type *, (test(iwr),iwr=1,4)
                     OK=.FALSE. 
                     THPOK=.FALSE.
                     IA = (ABS(THA)-ABS(THMINA))/DTHA 
                     IB = (ABS(THB)-ABS(THMINB))/DTHB
                     ITT = IA*NTH + IB 
* position in array
C
                     IF (IA.GE.0 .AND. IA.LE.(NTH-1) .AND. IB.GE.0 
     +                  .AND. IB.LE.(NTH-1)) OK=.TRUE. 
* within bounds of theta(A)-theta(B) plot
                     IF(THSTAR.GE.THSMIN .AND. THSTAR.LE.
     +                   (THSMIN+NYSET*DTH).AND. PSI.GE.PSIMIN 
     +                    .AND. PSI.LE.(PSIMIN+NXSET*DPSI)) THPOK=.TRUE.
* within bounds of thetastar-psi plot
C
                     IF (TEST(1).AND.TEST(2).AND.TEST(3).AND.TEST(4))
     +               THEN 
* a hit
                      IF (THPOK .AND. OK) THEN
                        IHITS = IHITS+1 
                        IF (MLI.LE.4 .AND. (IDET(1).EQ.2 .OR. 
     +                     IDET(2).EQ.2)) THEN 
                           IF (ELICOR.GT.ELIMAX) ELIMAX = ELICOR 
                        ENDIF 
C
                        IF (LELIM) THEN 
                           NOHIT(II)=' ' 
                           TTNOHT(ITT)=' ' 
                        ENDIF 
C
                        IF (IVAR.EQ.0) THEN 
                           HIT(II) = '*' 
                           TTHIT(ITT) = '*' 
                           WRITE(100,*) THSTAR,-1*PSI
                           ELSE 
                           IF (MAXTH.LE.36) THEN 
C                          26 letters + 10 digits
                              ISYMB2=ITIMTH 
                              ELSE 
                              ISYMB2=ICOUNT 
                           ENDIF 
                           ICH=54
* A-Z
                           IF (ISYMB2.LT.11) ICH=47 
* 0-9
                           ISYMB = ICH+ISYMB2 
                           SYMBOL = CHAR(ISYMB) 
                           IF (HIT(II).EQ.' ') THEN 
                              HIT(II) = SYMBOL
*XVGR output 
                              WRITE(100,*) THSTAR,-1*PSI,THREC
                              WRITE(200,*) THB,THA,THREC
                              ELSE 
                              WRITE(100,*) THSTAR,-1*PSI
                             HIT(II) = '+' 
* overlap
                              LOVL = .TRUE. 
                           ENDIF 
                           IF (TTHIT(ITT).EQ.' ') THEN 
                              TTHIT(ITT) = SYMBOL 
                              ELSE 
                              IF (TTHIT(ITT).NE.SYMBOL) 
     +                              TTHIT(ITT)='+' 
                              LOVL = .TRUE. 
                           ENDIF 
                        ENDIF 
                       ENDIF
                       ELSE 
                        IF (LELIM) THEN 
                           IF (IVAR.EQ.0 .OR. ITIMTH.EQ.1) THEN 
                              SYMBOL = CHARS(ICIND(TEST)) 
                              IF (THPOK) NOHIT(II) = SYMBOL 
                              IF (OK) THEN
                                IF (TTNOHT(ITT).NE.' ')
     +                              TTNOHT(ITT) = SYMBOL
                              ENDIF
                           ENDIF 
C                  Only first time for given Ex
                        ENDIF 
                     ENDIF 
C
C
  195                CONTINUE 
 1023             CONTINUE 
C-----------------------End of loop over psi------------------------
C
 1022          CONTINUE 
C----------------------End of loop over thetastar------------------
C
               IF (IVAR.EQ.3) GOTO 20 
* finished log file
C
               LFIRST=.TRUE. 
C
C    Open printer file
C
               IF (.NOT.LPOPEN) THEN
                  CALL RENAME(PLOTF,NEWF)
CGEC
                  OPEN(UNIT=PLOTU,FILE=PLOTF,STATUS='NEW',FORM='PRINT')
                  LPAGE=.FALSE. 
* no new page 1st time
                  LPOPEN=.TRUE. 
               ENDIF 
C
               HITSUM = HITSUM + IHITS 
               IF (IHITS.EQ.0) THEN 
* no hits
                  WRITE(6,*)' ' 
                  WRITE(6,205) EX,(THAMAX+THAMIN)/2.0,
     +             (THBMAX+THBMIN)/2.0 
  205             FORMAT(' No hits for Ex=',F6.2,' MeV ; A angle =',
     +              F6.2,', B angle =',F6.2/)
                  IF (IVAR.EQ.0 .AND. .NOT.LZER) GOTO 280 
                  IF (ITIMTH.EQ.MAXTH .AND. HITSUM.EQ.0 .AND. 
     +               IEX.EQ.NEN .AND. NEN.GT.1) THEN 
                     WRITE(6,210) 
  210                FORMAT(' No hits at all......nothing printed')
                     GOTO 20 
                  ENDIF 
C          Not printing zero hits
                  ELSE 
* some hits
                  ICOUNT=ICOUNT+1 
                  IF (ICOUNT.GT.36) THEN 
* only 36 characters used
                     WRITE(PLOTU,227) (THAMAX+THAMIN)/2.0, 
     +                (THBMAX+THBMIN)/2.0 
  227                FORMAT(5X,'Sorry, ran out of characters at',/
     +                'TH(A)=',F7.3,'; TH(B)=',F7.3)
                     FLAG = .TRUE.
                  ENDIF 
               ENDIF 
C
               LPAGE=.TRUE. 
* flag for new page
C
               IF (IVAR.NE.0) THEN 
                  IF (ITIMTH.LT.MAXTH) GOTO 290 
C          Keep going for options 1 & 2
               ENDIF 
C
  228          CONTINUE 
C
C    Print first page
C
               OPEN(UNIT=SCRTCU,STATUS='SCRATCH',ERR=421) 
C
               CALL PRNTF2(SCRTCU,TTHIT,LFIRST) 
C      Put TTHIT into temporary file
C
               IF (LPAGE) WRITE(PLOTU,230) 
  230          FORMAT('1')
* new page
               CALL PRINTF(PLOTU,HIT,SCRTCU,LFIRST) 
C      Write out parameters, HIT & TTHIT
               LFIRST=.FALSE. 
               LPAGE=.TRUE. 
C
               CLOSE(UNIT=SCRTCU) 
C
               IF (.NOT.LELIM) GOTO 280 
* not second page
C
C    Print second page
C
               OPEN(UNIT=SCRTCU,STATUS='SCRATCH',ERR=431) 
               CALL PRNTF2(SCRTCU,TTNOHT,LFIRST) 
C      Put TTNOHT into temporary file
               IF (LPAGE) WRITE(PLOTU,230) 
               CALL PRINTF(PLOTU,NOHIT,SCRTCU,LFIRST) 
C      Write out key, NOHIT & TTNOHT
               IF (LFIRST) THEN 
                  WRITE(PLOTU,230) 
                  CALL PRINT2(PLOTU,1) 
* key on separate page
               ENDIF 
               CLOSE(UNIT=SCRTCU) 
C
C
  280          CONTINUE 
               IF (FLAG) GOTO 350 
* ran out of characters - stop
C
  290          CONTINUE 
		WRITE(100,*) ' '
		WRITE(200,*) ' '
 1020       CONTINUE 
C---------------------End of loop over B angle--------------------------
 1019    CONTINUE 
C---------------------End of loop over A angle--------------------------
 1016 CONTINUE
C-----------------------End of loop over Ex-----------------------------
C
      GOTO 20
* back to OPTION
C
C
  350 CONTINUE
      IF (IVAR.EQ.3) THEN
         WRITE(6,360) LOGF
  360    FORMAT('Log file is ',A20)
         CLOSE(UNIT=LOGU)
         WRITE(6,*) ' '
      ENDIF
* log file
      IF (LPOPEN) THEN
         WRITE(6,370) PLOTF
  370    FORMAT('Plot file is ',A20)
         CLOSE(UNIT=PLOTU)
         WRITE(6,*) ' '
      ENDIF
* printer file
C
      STOP
C[
C   I/O errors:
C
  401 STOP 'Error in opening new log file in RPS'
  411 STOP 'Error in opening new plot file in RPS'
  421 STOP 'Error in opening first scratch work file in RPS'
  431 STOP 'Error in opening second scratch work file in RPS'
C
  399 END
C
C----------------------------------------------------------------------
C----------------------------------------------------------------------
C
      INTEGER FUNCTION ICIND(TEST)
C
C      Finds out which tests were failed for given thstar,psi
C      and gives index of corresponding code in CHARS
C
      LOGICAL TEST(4)
C
      J=0
* no. of failures
      DO 1025 I=1,4
         IF (.NOT.TEST(I)) J=J+1 
 1025 CONTINUE
C
      IF (J.EQ.1) THEN
         DO 1026 I=1,4 
            IF(.NOT.TEST(I)) ICIND=I 
 1026    CONTINUE 
         ELSE IF (J.EQ.2) THEN 
         DO 1027 I=1,4 
            DO 1028 K=I+1,4 
               IF((.NOT.TEST(I)).AND.(.NOT.TEST(K))) THEN 
                  ICIND = I+K+2 
                  IF(I.GE.2) ICIND=ICIND+1 
                  GOTO 400 
               ENDIF 
 1028       CONTINUE 
 1027    CONTINUE 
  400    CONTINUE 
         ELSE IF (J.EQ.3) THEN 
         DO 1029 I=1,4 
            IF (TEST(I)) ICIND=I+10 
 1029    CONTINUE 
         ELSE 
         ICIND=15 
      ENDIF
C
      END
C-----------------------------------------------------------------------
C      KEY:
C
C      Failed test            Eliminated by                  ICIND
C      -----------            -------------                  -----
C           1                   TH(A)                          1
C           2                   TH(B)                          2
C           3                   E(A)                           3
C           4                   E(B)                           4
C          1+2                  TH(A)+TH(B)                    5
C          1+3                  TH(A)+ E(A)                    6
C          1+4                  TH(A)+ E(B)                    7
C          2+3                  TH(B)+ E(A)                    8
C          2+4                  TH(B)+ E(B)                    9
C          3+4                  E(A) + E(B)                   10
C         2+3+4                 TH(B)+ E(A)+ E(B)             11
C         1+3+4                 TH(A)+ E(A)+ E(B)             12
C         1+2+4                 TH(A)+TH(B)+ E(B)             13
C         1+2+3                 TH(A)+TH(B)+ E(A)             14
C        1+2+3+4                TH(A)+TH(B)+ E(A)+ E(B)       15
C
C--------------------------------------------------------------------------
C
C
      SUBROUTINE CREATF
C
C     Set up parameters for control file
C     Format statements assume NXSET=64
C
CGEC
      INCLUDE 'rpscom'
C
CGEC
      CONTINUE
      WRITE(6,414)
  414 FORMAT('Ebeam, Mbeam, Mtarget : ',$)
      READ(5,*) EBEAM,MB,MT
C
CGEC
      CONTINUE
      WRITE(6,416)
  416 FORMAT('MHI, MLI : ',$)
      READ(5,*) MHI,MLI
      MREC = MB + MT - MHI - MLI
C
CGEC
      CONTINUE
      WRITE(6,418)
  418 FORMAT('Two-body Q-value, breakup threshold (MeV) : ',$)
      READ(5,*) Q2,THRESH
C
CGEC
  450 CONTINUE
      WRITE(6,419)
  419 FORMAT('Particles detected in telescopes A & B (1=HI, 2=LI,',    
     +' 3=recoil) : ',$)
      READ(5,*) IDET
      DO 460 I=1,2
         IF (IDET(I).LT.1 .OR. IDET(I).GT.3) THEN
             WRITE(6,*) 'Invalid value'
             GOTO 450
         ENDIF
  460 CONTINUE
C
CGEC
      CONTINUE
      WRITE(6,420)
  420 FORMAT('Telescope A distance from target & diameter : ',$)
      READ(5,*) DISTA,DIAMA
C
CGEC
      CONTINUE
      WRITE(6,422)
  422 FORMAT('Telescope B distance from target & diameter : ',$)
      READ(5,*) DISTB,DIAMB
C
CGEC
      CONTINUE
      WRITE(6,424)
  424 FORMAT('Minimum & maximum A energies : ',$)
      READ(5,*) EAMIN,EAMAX
CGEC
      CONTINUE
      WRITE(6,425)
  425 FORMAT('Minimum & maximum B energies : ',$)
      READ(5,*) EBMIN,EBMAX
C
      IF(MASSES(IDET(1)).LE.4) THEN
CGEC
         CONTINUE
         WRITE(6,426) 
  426    FORMAT('A absorber thickness in microns : ',$)
         READ(5,*) ABSRBS(IDET(1)) 
      ENDIF
      IF(MASSES(IDET(2)).LE.4) THEN
CGEC
         CONTINUE
         WRITE(6,427) 
  427    FORMAT('B absorber thickness in microns : ',$)
         READ(5,*) ABSRBS(IDET(2)) 
      ENDIF
C
CGEC
      CONTINUE
      WRITE(6,430)
  430 FORMAT('Thetastar & psi min. : ',$)
      CALL FREEFM(THSMIN,PSIMIN,TEMP)
C     Default is 0,0
C
CGEC
      CONTINUE
      WRITE(6,432)
  432 FORMAT('Thetastar & psi increments : ',$)
      CALL FREEFM(DTH,DPSI,TEMP)
C     Defaults:
      IF (DTH.EQ.0.0) DTH=1.0
      IF (DPSI.EQ.0.0) DPSI=2.0
C
CGEC
      CONTINUE
      WRITE(6,434)
  434 FORMAT('Min. A & B lab angles : ',$)
      CALL FREEFM(THMINA,THMINB,TEMP)
C     Default 0,0
C
CGEC
      CONTINUE
      WRITE(6,436)
  436 FORMAT('A & B lab angle increments : ',$)
      CALL FREEFM(DTHA,DTHB,TEMP)
C     Defaults:
      IF (DTHA.EQ.0.0) DTHA=1.0
      IF (DTHB.EQ.0.0) DTHB=1.0
C
      RETURN
C
      END
C---------------------------------------------------------------
C
      SUBROUTINE WRITEF(IUNIT,FNAME,LOGF)
C
C     Writes out control file FNAME on unit IUNIT
C     Leaves output file on IUNIT open/closed as it finds it
C
      CHARACTER*20 FNAME, LOGF
      LOGICAL LFOPEN
C
CGEC
      INCLUDE 'rpscom'
C
      IF (IUNIT.NE.6) THEN
* not for tty
         INQUIRE(FILE=FNAME,OPENED=LFOPEN)
         IF (.NOT. LFOPEN) OPEN(UNIT=IUNIT,FILE=FNAME,STATUS='NEW',
     +        ERR=599)
      ENDIF
C
      WRITE(IUNIT,519) EBEAM
  519 FORMAT(1X,F10.3,16X,'! Ebeam')
C
      WRITE(IUNIT,521) MB,MT
  521 FORMAT(1X,2I10,6X,'! Mbeam, Mtarget')
C
      WRITE(IUNIT,523) MHI,MLI
  523 FORMAT(1X,2I10,6X,'! MHI, MLI')
C
      WRITE(IUNIT,525) Q2
  525 FORMAT(1X,F10.3,16X,'! Two-body Q-value')
C
      WRITE(IUNIT,526) THRESH
  526 FORMAT(1X,F10.3,16X,'! Breakup threshold')
C
      IF (IUNIT.EQ.6.OR. FNAME.EQ.LOGF) THEN
         WRITE(IUNIT,527) (PTCLE(IDET(I)),I=1,2)
  527    FORMAT(5X,A6,4X,A6,6X,'! Detected particles A & B')
         ELSE
         WRITE(IUNIT,528) IDET
  528    FORMAT(1X,2I10,6X,'! Detected particles A & B')
      ENDIF
C
      WRITE(IUNIT,529) DISTA, DIAMA
  529 FORMAT(1X,2F10.3,6X,'! A dist & diam')
C
      WRITE(IUNIT,530) DISTB, DIAMB
  530 FORMAT(1X,2F10.3,6X,'! B dist & diam')
C
      WRITE(IUNIT,531) EAMIN, EAMAX
  531 FORMAT(1X,2F10.3,6X,'! EA min, EA max')
      WRITE(IUNIT,532) EBMIN, EBMAX
  532 FORMAT(1X,2F10.3,6X,'! EB min, EB max')
C
      IF (MASSES(IDET(1)).LE.4) THEN
         WRITE(IUNIT,533) ABSRBS(IDET(1))
  533    FORMAT(1X,F10.3,16X,'! A absorber thickness')
      ENDIF
      IF (MASSES(IDET(2)).LE.4) THEN
         WRITE(IUNIT,534) ABSRBS(IDET(2))
  534    FORMAT(1X,F10.3,16X,'! B absorber thickness')
      ENDIF
C
      WRITE(IUNIT,535) THSMIN,PSIMIN
  535 FORMAT(1X,2F10.3,6X,'! Thstar min, psi min')
C
      WRITE(IUNIT,536) DTH,DPSI
  536 FORMAT(1X,2F10.3,6X,'! Thstar & psi increments')
C
      WRITE(IUNIT,537) THMINA,THMINB
  537 FORMAT(1X,2F10.3,6X,'! Min. A & B lab. angles')
C
      WRITE(IUNIT,538) DTHA,DTHB
  538 FORMAT(1X,2F10.3,6X,'! A & B lab. angle increments')
C
C
      IF (.NOT.(IUNIT.EQ.6 .OR. LFOPEN)) CLOSE(UNIT=IUNIT)
C
      RETURN
C
  599 WRITE(6,*) 'Error in opening control file in WRITEF'
      RETURN
C
      END
C
C---------------------------------------------------------------
C
      SUBROUTINE READF(FNAME)
C
C     Read control file
C
      CHARACTER*20 FNAME
      INTEGER RCONTU
C
CGEC
      INCLUDE 'rpscom'
C
      RCONTU = 8
      OPEN(UNIT=RCONTU,FILE=FNAME,STATUS='OLD',ERR=699)
C
      READ(RCONTU,*) EBEAM
      READ(RCONTU,*) MB,MT
      READ(RCONTU,*) MHI,MLI
      MREC = MB + MT - MHI - MLI
      READ(RCONTU,*) Q2
      READ(RCONTU,*) THRESH
      READ(RCONTU,*) IDET
      READ(RCONTU,*) DISTA,DIAMA
      READ(RCONTU,*) DISTB,DIAMB
      READ(RCONTU,*) EAMIN,EAMAX
      READ(RCONTU,*) EBMIN,EBMAX
      IF (MASSES(IDET(1)).LE.4) READ(RCONTU,*) ABSRBS(IDET(1))
      IF (MASSES(IDET(2)).LE.4) READ(RCONTU,*) ABSRBS(IDET(2))
      READ(RCONTU,*) THSMIN,PSIMIN
      READ(RCONTU,*) DTH,DPSI
      READ(RCONTU,*) THMINA,THMINB
      READ(RCONTU,*) DTHA,DTHB
C
      CLOSE(UNIT=RCONTU)
      RETURN
C
  699 WRITE(6,*) 'Error in opening control file in READF'
      RETURN
      END
C
C---------------------------------------------------------------
C
      SUBROUTINE PRINTF(IUNIT,AIN,IN,LFIRST)
C
C     To print out array AIN and parameters
C     Writes out on IUNIT, reads in on IN
C     Assumes IUNIT & IN open
C
      CHARACTER*51 LINE
C     CHARACTER*1 AIN(NXSET*NYSET)
      CHARACTER*1 AIN(64*45)
      LOGICAL LFOPEN, LFIRST
C
CGEC
      INCLUDE 'rpscom'
C
      INQUIRE(UNIT=IUNIT,OPENED=LFOPEN)
      IF (.NOT.LFOPEN) THEN
         WRITE(6,*) 'Unit ', IUNIT, ' not open in PRINTF'
         RETURN
      ENDIF
C
      INQUIRE(UNIT=IN,OPENED=LFOPEN)
      IF (.NOT.LFOPEN) THEN
         WRITE(6,*) 'Unit ', IN, ' not open in PRINTF'
         RETURN
      ENDIF
C
      REWIND(UNIT=IN)
C
      WRITE(IUNIT,*) ' '
      IF (LFIRST) THEN
        IF (IVAR.EQ.2) THEN
           WRITE(IUNIT,700) EX
           ELSE
           WRITE(IUNIT,702) EX,DISTA,DIAMA
           IF (IVAR.EQ.0) THEN
              WRITE(IUNIT,703) (THAMAX+THAMIN)/2.0,
     +           (THBMAX+THBMIN)/2.0,DISTB,DIAMB
              ELSE
              WRITE(IUNIT,704) DISTB,DIAMB
           ENDIF
         ENDIF
         WRITE(IUNIT,*) ' ' 
         WRITE(IUNIT,705) EBEAM,EAMIN,EAMAX 
         WRITE(IUNIT,706) MB,MT,EBMIN,EBMAX 
         WRITE(IUNIT,707) MHI,MLI,(ABSRBS(IDET(I)),I=1,2) 
         WRITE(IUNIT,708) Q2,THRESH,(PTCLE(IDET(I)),I=1,2) 
         ELSE 
         CALL PRINT2(IUNIT,2) 
      ENDIF
C
  700 FORMAT(1X,F10.3,16X,'! Ex')
  701 FORMAT(1X,2F10.3,6X,'! A size, B size')
  702 FORMAT(1X,F10.3,16X,'! Ex',33X,2F10.3,5X,'!A dist & diam')
  703 FORMAT(1X,2F10.3,6X,'!Angle A, Angle B',19X,2F10.3,6X,            
     +'!B dist & diam')
  704 FORMAT(63X,2F10.3,6X,'!B dist & diam')
  705 FORMAT(1X,F10.3,16X,'!Ebeam',30X,2F10.3,6X,'!EA min, EA max')
  706 FORMAT(1X,2I10,6X,'!Mbeam, Mtarget',21X,2F10.3,6X,                
     +'!EB min, EB max')
  707 FORMAT(1X,2I10,6X,'!MHI, MLI',27X,2F10.3,6X,                      
     +'!Absorber thicknesses A & B')
  708 FORMAT(1X,2F10.3,6X,'!Two-body Q-value, breakup threshold',       
     +4X,A6,4X,A6,6X,'!Detected particles A & B')
  711 FORMAT(1X,2I10,6X,'!MHI, MLI',31X,A6,4X,A6,6X,'!Detected',        
     +' particles A & B')
  713 FORMAT(1X,2F10.3,6X,'!Two-body Q-value, breakup threshold')
C
      WRITE(IUNIT,*)' '
      WRITE(IUNIT,*)' '
      WRITE(IUNIT,*)' '
      IF (.NOT.LFIRST) WRITE(IUNIT,*)' '
C
      IF (MLI.LE.4 .AND. (IDET(1).EQ.2 .OR. IDET(2).EQ.2) .AND.         
     +  LFIRST) WRITE(IUNIT,714) ELIMAX
  714 FORMAT(86X,'Max. LI energy in detector =',F9.3,' MeV')
      WRITE(IUNIT,*)' '
      WRITE(IUNIT,715) ((PSIMIN + (I-1)*DPSI),I=1,NXSET+1,10) 
  715 FORMAT(1X,7(F9.1,1X),9X,A51)
C
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,718) LINE
  718 FORMAT(8X,7('!',9X),2X,A51)
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,717) (THSMIN+NYSET*DTH,I=1,2), LINE
  717 FORMAT(1X,F5.1,68('-'),F5.1,1X,A51)
C
      DO 1004 IY=2,NYSET
        READ(IN,722,ERR=799) LINE
  722   FORMAT(A51)
        IF(MOD(IY,5).EQ.0) THEN
          WRITE(IUNIT,725) (THSMIN+(NYSET-IY)*DTH), (AIN((NYSET-IY)*     
     +     NXSET+IX),IX=1,NXSET), (THSMIN+(NYSET-IY)*DTH),LINE 
  725     FORMAT(1X,F5.1,'-!',64(A),'!-',F5.1,1X,A51)
          ELSE 
          IF(IY.EQ.18) THEN 
            WRITE(IUNIT,726) (AIN((NYSET-IY)*NXSET+IX), IX=1,NXSET),    
     +        LINE 
  726       FORMAT(1X,'THSTAR!',64(A),'!',7X,A51)
            ELSE 
            WRITE(IUNIT,727) (AIN((NYSET-IY)*NXSET+IX), IX=1,NXSET),    
     +        LINE 
  727       FORMAT(7X,'!',64(A),'!',7X,A51)
          ENDIF 
        ENDIF
 1004 CONTINUE
C
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,740) LINE
  740 FORMAT(7X,'!',64('-'),'!',7X,A51)
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,741) LINE
  741 FORMAT(8X,6('!',9X),'!',11X,A51)
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,715) ((PSIMIN+(I-1)*DPSI),I=1,NXSET+1,10),LINE
      READ(IN,722,ERR=799) LINE
      WRITE(IUNIT,745) LINE
  745 FORMAT(42X,'PSI',35X,A51)
C
  770 RETURN
C
  799 WRITE(6,*) 'Error in reading in PRINTF'
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PRINT2(IUNIT,INDEX)
C
C     Subroutine to print out key for plot showing what causes limits
C
CGEC
      INCLUDE 'rpscom'
C
      IF (INDEX.EQ.1) THEN
* all in one column down page
         WRITE(IUNIT,805) 
  805    FORMAT(1X,//)
         WRITE(IUNIT,810) 
  810    FORMAT(11X,'KEY:',/)
         DO 1005 I=1,15  
            GOTO (812,814,816,818,820,822,824,826,828,830,832,
     +        834,836,838,840),I 
  812       WRITE(IUNIT,853) CHARS(I) 
            GOTO 845 
  814       WRITE(IUNIT,854) CHARS(I) 
            GOTO 845 
  816       WRITE(IUNIT,855) CHARS(I) 
            GOTO 845 
  818       WRITE(IUNIT,856) CHARS(I) 
            GOTO 845 
  820       WRITE(IUNIT,857) CHARS(I) 
            GOTO 845 
  822       WRITE(IUNIT,858) CHARS(I) 
            GOTO 845 
  824       WRITE(IUNIT,859) CHARS(I) 
            GOTO 845 
  826       WRITE(IUNIT,860) CHARS(I) 
            GOTO 845 
  828       WRITE(IUNIT,861) CHARS(I) 
            GOTO 845 
  830       WRITE(IUNIT,862) CHARS(I) 
            GOTO 845 
  832       WRITE(IUNIT,863) CHARS(I) 
            GOTO 845 
  834       WRITE(IUNIT,864) CHARS(I) 
            GOTO 845 
  836       WRITE(IUNIT,865) CHARS(I) 
            GOTO 845 
  838       WRITE(IUNIT,866) CHARS(I) 
            GOTO 845 
  840       WRITE(IUNIT,867) CHARS(I) 
            GOTO 845 
  845       CONTINUE                          
 1005    CONTINUE 
         WRITE(IUNIT,846) 
C
  846    FORMAT(5X,'''.'' denotes point not accessible for this',          
     +     ' range of thstar,psi')
         ELSE 
* three columns
         DO 1006 I=1,5  
            GOTO (870,872,874,876,878),I 
  870       WRITE(IUNIT,882) CHARS(I),CHARS(I+5),CHARS(I+10) 
            GOTO 880 
  872       WRITE(IUNIT,883) CHARS(I),CHARS(I+5),CHARS(I+10) 
            GOTO 880 
  874       WRITE(IUNIT,884) CHARS(I),CHARS(I+5),CHARS(I+10) 
            GOTO 880 
  876       WRITE(IUNIT,885) CHARS(I),CHARS(I+5),CHARS(I+10) 
            GOTO 880 
  878       WRITE(IUNIT,886) CHARS(I),CHARS(I+5),CHARS(I+10) 
            GOTO 880 
  880       CONTINUE                          
1006     CONTINUE 
      ENDIF
      RETURN
C
  853 FORMAT(5X,'''',A,'''',3X,'TH(A)')
  854 FORMAT(5X,'''',A,'''',3X,'TH(B)')
  855 FORMAT(5X,'''',A,'''',3X,'E(A)')
  856 FORMAT(5X,'''',A,'''',3X,'E(B)')
  857 FORMAT(5X,'''',A,'''',3X,'TH(A)+TH(B)')
  858 FORMAT(5X,'''',A,'''',3X,'TH(A)+E(A)')
  859 FORMAT(5X,'''',A,'''',3X,'TH(A)+E(B)')
  860 FORMAT(5X,'''',A,'''',3X,'TH(B)+E(A)')
  861 FORMAT(5X,'''',A,'''',3X,'TH(B)+E(B)')
  862 FORMAT(5X,'''',A,'''',3X,'E(A) +E(B)')
  863 FORMAT(5X,'''',A,'''',3X,'TH(B)+E(A) +E(B)')
  864 FORMAT(5X,'''',A,'''',3X,'TH(A)+E(A) +E(B)')
  865 FORMAT(5X,'''',A,'''',3X,'TH(A)+TH(B)+E(B)')
  866 FORMAT(5X,'''',A,'''',3X,'TH(A)+TH(B)+E(A)')
  867 FORMAT(5X,'''',A,'''',3X,'ALL FOUR')
C
  882 FORMAT(1X,'KEY:',5X,'''',A,'''',3X,'TH(A)',11X,'''',A,'''',3X,    
     +'TH(A)+E(A)',5X,'''',A,'''',5X,'TH(B)+ E(A)+E(B)',                
     +10X,'''.'' denotes point not accessible')
  883 FORMAT(10X,'''',A,'''',3X,'TH(B)',11X,'''',A,'''',3X,             
     +'TH(A)+E(B)',5X,'''',A,'''',5X,'TH(A)+ E(A)+E(B)',                
     +14X,'for this range of thstar,psi')
  884 FORMAT(10X,'''',A,'''',4X,'E(A)',11X,'''',A,'''',3X,              
     +'TH(B)+E(A)',5X,'''',A,'''',5X,'TH(A)+TH(B)+E(B)')
  885 FORMAT(10X,'''',A,'''',4X,'E(B)',11X,'''',A,'''',3X,              
     +'TH(B)+E(B)',5X,'''',A,'''',5X,'TH(A)+TH(B)+E(A)')
  886 FORMAT(10X,'''',A,'''',3X,'TH(A)+TH(B)',5X,'''',A,'''',4X,        
     +'E(A)+E(B)',5X,'''',A,'''',5X,'ALL FOUR')
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE PRNTF2(IUNIT,AIN,LFIRST)
C
C     Plot TH(A) vs TH(B)
C
CGEC     Format statements assume NTH=40
C
C
C     CHARACTER*1 AIN(0:NTH*NTH-1)
      CHARACTER*1 AIN(0:40*40-1)
      CHARACTER*5 TH(3)
      LOGICAL LFOPEN, LFIRST
C
CGEC
      INCLUDE 'rpscom'
      DATA TH/'TH HI','TH LI','TH RE'/
C
      INQUIRE(UNIT=IUNIT,OPENED=LFOPEN)
      IF (.NOT.LFOPEN) THEN
         WRITE(6,*) 'Unit ', IUNIT, ' not open in PRNTF2'
         RETURN
      ENDIF
C
      IF (LFIRST .AND. LOVL) THEN
        WRITE(IUNIT,902)
  902   FORMAT(6X,'''+'' denotes overlap'//)
        ELSE
        WRITE(IUNIT,903)
  903   FORMAT(1X,//)
      ENDIF
C
      WRITE(IUNIT,905) (THMINB+I*DTHB,I=0,NTH,10)
  905 FORMAT(5X,10(F6.1,4X))
C
      WRITE(IUNIT,906)
  906 FORMAT(9X,5('!',9X))
C
      WRITE(IUNIT,908) THMINA+NTH*DTHA
  908 FORMAT(1X,F6.1,43('-'))
C
      DO 1007 I=1,NTH 
      IF (MOD(I,5).EQ.0) THEN
         WRITE(IUNIT,910) (THMINA+(NTH-I)*DTHA), 
     +     (AIN((NTH-I)*NTH+J),J=0,NTH-1) 
  910    FORMAT(1X,F6.1,'-!',40(A),'!')
         ELSE 
         IF (I.EQ.22) THEN 
            WRITE(IUNIT,912) TH(IDET(1)), (AIN((NTH-I)*NTH+J),J=0,NTH-1)
  912       FORMAT(1X,A6,' !',40(A),'!')
            ELSE IF (I.EQ.23) THEN 
            WRITE(IUNIT,913) (AIN((NTH-I)*NTH+J),J=0,NTH-1) 
  913       FORMAT(2X,'(A)   !',40(A),'!')
            ELSE 
            WRITE(IUNIT,914) (AIN((NTH-I)*NTH+J),J=0,NTH-1) 
  914       FORMAT(7X,' !',40(A),'!')
         ENDIF 
      ENDIF                         
1007  CONTINUE
C
      WRITE(IUNIT,915)
  915 FORMAT(8X,42('-'))
      WRITE(IUNIT,916)
  916 FORMAT(9X,4('!',9X),'!')
      WRITE(IUNIT,905) (THMINB+I*DTHB,I=0,NTH,10)
      WRITE(IUNIT,917) TH(IDET(2))
  917 FORMAT(21X,A5,' (B)')
C
      RETURN
      END
C
C---------------------------------------------------------------
C
      SUBROUTINE RENAME(OLD,NEW)
C
CGEC     Adds necessary GEC file attributes to file name
CSUN     Not necessary for SUN I think therefore return
C
      CHARACTER*20 OLD
      CHARACTER*30 ATTR
      CHARACTER*50 NEW
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE ABSORB(ABSBR,EF,ITAB,IA)
C
C     Calculates remaining energy of alpha-particles after passing
C         through mylar absorber
C
C     Procedure:
C         (a) produce arrays of energy before (EIN) & after (EOUT)
C      passing through absorber
C         (b) produce lookup table (TAB) by quadratic interpolation
C      from these arrays
C         (c) do the lookup by linear interpolation between points
C      in table
C
C     Taken from program ABSORBER & PRJ.FOR
C     All data from Northcliffe & Schilling
C
C     Table produced in TAB with NPTS points, from input arrays
C     of up to NPTSIN points
C     Modified 28 Oct 1985 to produce two tables         SCA
C        9 Nov 1985 to do protons, deuterons & tritons   SCA
C
C*************************************************************************
C
      PARAMETER (NPTSIN=200, NPTS=500)
      PARAMETER (EI=1.0,DE=1.0)
      PARAMETER (RHO=1.395)
* density of mylar in g/cm3
C
      DIMENSION TAB(2*NPTS,3)
      DIMENSION EIN(NPTSIN+2), EOUT(NPTSIN+2)
C
C-------------------------------------------------------------------------
C
      IZ = 1
      IF (IA.EQ.4) IZ = 2
C
      ABSMG=ABSBR*RHO/10.
* thickness in mg/cm2
C
      IF(EF.GT.200.) EF=200.
      IF(((EF-EI)/DE).GE.NPTSIN) EF = EI + DE*(NPTSIN-1)
C
      CALL ENERGY(ABSMG,EMIN,IZ,IA)
C
      NP = 0
      DO 1004 E=EI,EF,DE
         NP=NP+1 
         IF(E.LE.EMIN) THEN 
            EX = 0.0 
            GOTO 700 
         ENDIF 
C     Get range R(mg/cm2) for input energy E
         CALL RANGE(E,R,IZ,IA) 
C     Find excess of this range over ABSBRer thickness
         RX = R - ABSMG 
C     Get energy for reduced range
         CALL ENERGY(RX,EX,IZ,IA) 
C
  700    EIN(NP+1) = E 
         EOUT(NP+1) = EX 
 1004 CONTINUE
C
C**********************************************************************
C
C     Extend
      EIN(1) = EIN(2) - DE
      EOUT(1) = EOUT(2)
      EIN(NP+2) = EIN(NP+1) + DE
      EOUT(NP+2) = EOUT(NP+1)
C
      ESTEP = (EIN(NP+2) - EIN(1))/(NPTS-1)
C
C     Fill in table
      DO 1005 I=1,NPTS
         E0 = EIN(1) + (I-1)*ESTEP 
C         Check if outside range
         IF(E0.LE.EIN(1)) THEN 
            TAB(2*I-1,ITAB) = E0 
            TAB(2*I,ITAB) = EOUT(1) 
            GOTO 710 
            ELSE IF(E0.GE.EIN(NP+2)) THEN 
            TAB(2*I-1,ITAB) = E0 
            TAB(2*I,ITAB) = EOUT(NP+2) 
            GOTO 710 
         ENDIF 
C
C         Quadratic interpolation
         DO 1006 J=1,NP 
            IF(EIN(J).GE.E0) GOTO 705 
 1006    CONTINUE 
C
  705    E1 = EIN(J-1) 
         E2 = EIN(J) 
         E3 = EIN(J+1) 
C
C
         Y0 = EOUT(J-1)*(E0-E2)*(E0-E3)/((E1-E2)*(E1-E3)) + EOUT(J)*(E0-
     +   E1)*(E0-E3)/((E2-E1)*(E2-E3)) + EOUT(J+1)*(E0-E1)*(E0-E2)/((E3-
     +   E1)*(E3-E2)) 
         IF(Y0.LT.0.0) Y0=0.0 
         TAB(2*I-1,ITAB) = E0 
         TAB(2*I,ITAB) = Y0 
  710    CONTINUE 
 1005 CONTINUE
C
      RETURN
C
C*******************************************************************
C
      ENTRY LOOKUP(XIN,XOUT,ITAB)
C
C     Use the table to look up EIN
C
      IX = (XIN - TAB(1,ITAB))/ESTEP + 1
C
C     Linear interpolation between points in TAB; IX is x-value
C     just below XIN
C
      IF(IX.LT.0) THEN
         XOUT = TAB(2,ITAB) 
         RETURN 
         ELSE IF(IX.GE.NPTS) THEN 
         XOUT = TAB(2*NPTS,ITAB) 
         RETURN 
         ELSE 
         X1 = TAB(2*IX-1,ITAB) 
         X2 = TAB(2*IX+1,ITAB) 
         Y1 = TAB(2*IX,ITAB) 
         Y2 = TAB(2*IX+2,ITAB) 
         XOUT = ((Y2-Y1)*XIN + Y1*X2 - X1*Y2)/ESTEP 
      ENDIF
C
      RETURN
C
      END
C*************************************************************************
      SUBROUTINE RANGE(E,R,IZ,IA)
C
C     Outputs range R(mg/cm2 mylar) for input energy E/MeV
C
      DIMENSION C(0:7,2)
C     Coefficients of 7th order polynomial used to fit N&S data
C        (log(R/A) vs log(E/A))
C     Row 1 for protons, row 2 for alphas
      DIMENSION S(0:1,2)
C     Coefficients of straight line used for extrapolation
      DIMENSION AMASS(4)
C
      DATA C/.4152, 1.6040, .1586,-.1506, .1321,.00420,-.06778,.02403,  
     +-.1771, 1.3710, .6403,-.1305,-.3354,.00350, .1097, .0306/
      DATA S/.3275,1.7933,-.2438,1.7573/
C
      DATA AMASS/1.007825,2.014102,3.016049,4.002603/
C     Atomic masses of 1H, 2H, 3H, 4He
C
C
      IF (E.LE.0.0) THEN
         R=0.0 
         RETURN 
      ENDIF
C
      IF (IA.EQ.4) THEN
         EA = E/AMASS(4) 
         ELSE IF (IA.EQ.1) THEN 
         EA = E 
         ELSE 
         EA = E * AMASS(1)/AMASS(IA) 
      ENDIF
      EN = ALOG10(EA)
C
C     Evaluate polynomial
      IF (EN.EQ.0.0) THEN
         POLY = C(0,IZ) 
         GOTO 745 
      ENDIF
      POLY = 0.0
      DO 1007 I=0,7
         POLY = POLY + C(I,IZ)*EN**I 
 1007 CONTINUE
C
C     Evaluate straight line
  745 STLIN = S(0,IZ) + S(1,IZ)*EN
C
      IF (EA.LT.2.0) THEN
C     Pure polynomial
         R = POLY 
         GOTO 749 
         ELSE IF (EA.LT.12.0) THEN 
C     Mix polynomial & st.line
         R = (POLY*(12.-EA) + STLIN*(EA-2.))/10. 
         GOTO 749 
         ELSE 
C     Pure straight line
         R = STLIN 
      ENDIF
C
  749 R = 10.**R
      IF (IA.EQ.2 .OR. IA.EQ.3) THEN
         R = R * AMASS(IA)/AMASS(1) 
C         Correct proton data for deuterons or tritons
         ELSE 
         R = R * IA 
      ENDIF
C
      RETURN
C
      END
C
C*************************************************************************
      SUBROUTINE ENERGY(R,E,IZ,IA)
C
C     Outputs energy E/MeV corresponding to range R(mg/cm2 mylar)
C
      DIMENSION D(0:7,2)
      DIMENSION T(0:1,2)
      DIMENSION AMASS(4)
C
      DATA D/-.2676,.6758,-.1109,.0995,-.0490,.00435,.00418,-.00099,    
     +.1265,.6252,-.2094,.2805,-.0436,-.0922,.0167,  .0104/
      DATA T/-.1826,.5576,.1388,.5690/
C
      DATA AMASS/1.007825,2.014102,3.016049,4.002603/
C
      IF (R.LE.0.0) THEN
         E=0.0 
         RETURN 
      ENDIF
      IF (IA.EQ.2 .OR. IA.EQ.3) THEN
         RA = R * AMASS(1)/AMASS(IA) 
         ELSE 
         RA = R/IA 
      ENDIF
      RN = ALOG10(RA)
C
      IF (RN.EQ.0.0) THEN
         POLY = D(0,IZ) 
         GOTO 755 
      ENDIF
      POLY = 0.0
      DO 1008 I=0,7
         POLY = POLY + D(I,IZ)*RN**I 
 1008 CONTINUE
C
  755 STLIN = T(0,IZ) + T(1,IZ)*RN
C
      IF (IZ.EQ.1) THEN
         IF (RA.LT.8.102) THEN 
            E = POLY 
            ELSE IF (RA.LT.183.0) THEN 
            E = (POLY*(183.0-RA) + STLIN*(RA-8.102))/174.9 
            ELSE 
            E = STLIN 
         ENDIF 
         ELSE 
         IF (RA.LT.1.957) THEN 
            E = POLY 
            ELSE IF (RA.LT.45.26) THEN 
            E = (POLY*(45.26-RA) + STLIN*(RA-1.957))/43.3 
            ELSE 
            E = STLIN 
         ENDIF 
      ENDIF
C
      IF (IA.EQ.1) THEN
         E = 10.**E 
         ELSE IF (IA.EQ.4) THEN 
         E = AMASS(4)*(10.**E) 
         ELSE 
         E = AMASS(IA)*(10.**E)/AMASS(1) 
      ENDIF
      RETURN
C
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE FREEFM(F1,F2,F3)
      COMMON/FMT/ISIZE
      DIMENSION READ(3)
      ISIZE=3
      READ(1)=F1
      READ(2)=F2
      READ(3)=F3
      CALL GOBBLE(READ)
      F1=READ(1)
      F2=READ(2)
      F3=READ(3)
      RETURN
      END
C
C---------------------------------------------------
C
      SUBROUTINE FREEF2(READ)
CGEC  DIMENSION READ(2*MAXANG)
      DIMENSION READ(2*20)
      COMMON/FMT/ISIZE
      INCLUDE 'rpscom'
      ISIZE=2*MAXANG
      CALL GOBBLE(READ)
      RETURN
      END
C
C---------------------------------------------------
C
      SUBROUTINE GOBBLE(READ)
C
C ** TO READ <ISIZE> (FLOATING POINT) NUMBERS IN "FREE FORMAT".
C ** COPIED FROM [NAJWINFIE] & MODIFIED 25/1/85       SCA
      LOGICAL BLANK
      CHARACTER *80 LINE
      COMMON/DIGT2/LINE
      COMMON/FMT/ISIZE
      DIMENSION READ(ISIZE)
    1 READ(5,10)LINE
      IF(LINE(1:1).EQ.';') GO TO 1
   10 FORMAT(A80)
      DO 30 I=1,ISIZE
   30 READ(I)=0.0
C
C ** KHAR points to position in line, KNUM keeps track of which number
C ** we're dealing with.
C
      KHAR=1
      DO 100 KNUM=1,ISIZE
         FNUM=0.0 
         TEN=10. 
         TENTH=1. 
         SIGNF=1. 
   50    IF(LINE(KHAR:KHAR).EQ.'.') THEN 
            TEN=1. 
            TENTH=0.1 
            KHAR=KHAR+1 
         END IF 
         CALL DIGIT(BLANK,KHAR,DGT) 
         KHAR=KHAR+1 
         IF(.NOT.BLANK) THEN 
            IF(DGT.GT.9.5) SIGNF=-1. 
            IF(DGT.GT.9.5) GO TO 50 
            FNUM= FNUM*TEN + DGT*TENTH 
            IF(TENTH.LT.1.) TENTH = TENTH*0.1 
            IF(KHAR.LE.70) GO TO 50 
         END IF 
         READ(KNUM)=FNUM*SIGNF 
  100 CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C
      SUBROUTINE DIGIT(BLANK,K,D)
C
C ** Returns numerical value of character of line when called by FREEFM.
C ** If a digit (0...9) or a minus sign BLANK = .FALSE. , any other
C ** character (i.e. a delimiter) BLANK = .TRUE.
C
      LOGICAL BLANK
      CHARACTER *80 LINE
      CHARACTER *11 NUMBER
      COMMON/DIGT2/LINE
      DATA NUMBER/'0123456789-'/
      BLANK=.TRUE.
      DO 10 J=1,11
         IF(LINE(K:K).EQ.NUMBER(J:J)) THEN 
            BLANK=.FALSE. 
            D=J-1 
C  n.b. a minus sign flagged by D = 10.
            GO TO 20 
         END IF 
   10 CONTINUE
   20 RETURN
      END
