 
C ----------------------------------------------------------------------
C		     Subroutine VINTERP
C
C ----------------------------------------------------------------------
 
      SUBROUTINE VINTERP(NDEV,EL,RNGES,GTIM,NPTS,NHOP,NPR)
 
	implicit double precision (a-h,o-z)
 
      REAL*8 NPOS,EL(NPTS),RNGES(NPTS),GTIM(NPTS)
C     COMMON DECK "CRAYCOM" INSERTED HERE
      COMMON/RAYCOM/NRAYS,MXRYS,MXHOP,MXXHOP,NRSLIM,ICODE
      COMMON/RAYCOM/NHOPS(333),MNHOPS(333)
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
	character id*80
      common /ident/ id
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      REAL*8 NRSLMI
      EQUIVALENCE (W(80),APOG) ,(W(401),PTYP) ,(W(403),ELMIN)
     1   ,(W(404),ELMAX) ,(W(405),ELVINC) ,(W(406),RGINC)
     2   ,(W(407),HOPLME) ,(W(408),NRSLME)
     3   ,(W(409),GAPLO), (W(410),GAPHI), (W(411),SYMRED)
     4   ,(W(428),GRPMIN) ,(W(429),GRPMAX) ,(W(430),GRPVINC)
     5   ,(W(432),HOPLMG) ,(W(433),NRSLMG)
     6   ,(W(434),RFSPEED) ,(W(435),TMFT) ,(W(436),RMAX)
     7   ,(W(453),RCRNG) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
      REAL*8 MAXSTP,MAXERR,INTYP,LLAT,LLON
      EQUIVALENCE (EARTHR,W(1)),(RAY,W(2)),(XMTRH,W(3)),(TLAT,W(4)),
     1 (TLON,W(5)),(OW,W(6)),(FBEG,W(7)),(FEND,W(8)),(FSTEP,W(9)),
     2 (AZ1,W(10)),(AZBEG,W(11)),(AZEND,W(12)),(AZSTEP,W(13)),
     3 (BETA,W(14)),(ELBEG,W(15)),(ELEND,W(16)),(ELSTEP,W(17)),
     8   (RUNSUP,W(18)),(RCVRH,W(20)),
     4 (ONLY,W(21)),(HOP,W(22)),(MAXSTP,W(23)),(PLAT,W(24)),(PLON,W(25))
     5,(HMAX,W(26)),(RAYFNC,W(29)),(EXTINC,W(33)),
     6 (HMIN,W(27)),(RGMAX,W(28)),
     8 (INTYP,W(41)),(MAXERR,W(42)),(ERATIO,W(43)),
     6 (STEP1,W(44)),(STPMAX,W(45)),(STPMIN,W(46)),(FACTR,W(47)),
     7 (SKIP,W(71)),(RAYSET,W(72)),(PRTSRP,W(74)),(HITLET,W(75))
     9 ,(BINRAY,W(76)),(PAGLN,W(77)),(PLT,W(81)),(PFACTR,W(82)),
     1 (LLAT,W(83)),(LLON,W(84)),(RLAT,W(85)),(RLON,W(86))
     2,(TIC,W(87)),(HB,W(88)),(HT,W(89)),(TICV,W(96))
C
C     COMMON DECK "CEFILES" INSERTED HERE
      COMMON/EFILES/INFILE,IOTFIL,LOGFIL,NCNTRL,NRYSFL
C
      external binary
C
      DATA KOUNT/0/
C
      KOUNT=KOUNT+1
      IF(KOUNT.EQ.1) WRITE(LOGFIL,35)
35    FORMAT('  ENTRY HOP-NO RAY-NO  EST EL      ELEV(LOW/HIGH)   '
     1   ,'   RANGES(LOW/HIGH)     SLOPE  GROUP TIMES')
C
      OPOS=BINARY(RNGES(1)-RCRNG)
      DO 100 I=2,NPTS
         NPOS=BINARY(RNGES(I)-RCRNG)
         IF(OPOS*NPOS.EQ.1.0) GO TO 100
         OPOS=NPOS
         IM1=I-1
         IF(RNGES(I).GT.RMAX .OR. RNGES(IM1).GT.RMAX) GO TO 100
         IF(RNGES(I).EQ.RNGES(IM1)) GO TO 100
         SLOPE=(EL(I)-EL(IM1))/(RNGES(I)-RNGES(IM1))
C
C     USE RANGE SENSITIVITY AS CUTOFF ON RAYS
         ELF=EL(IM1)+SLOPE*(RCRNG-RNGES(IM1))
30       IF(DABS(SLOPE).LE.SLPCUT) GO TO 60
C
crb   There are several lines added here by Ron Brent to write
crb   information to the file HARPO.IN.
ccccc NOTE:  The lines added by Ron Brent have now been commented out,
ccccc        since NDEV is now passed into this routine as 43.  Thus,
ccccc        RB's correct modifications are no longer needed. jah.

         IF(NPR.NE.0) WRITE(NDEV,40) '  0'
         IF(NPR.NE.0) WRITE(NDEV,40) ID
         NPR=NPR+1
         WRITE(NDEV,45) sngl(ELF)
40    FORMAT(10A8)
45    FORMAT(' 15',T4,F8.3,T18,'1',T25,'ELEVATION ANGLE')
      FCTR=1.0
      IF(RNGES(I).NE.RNGES(IM1))
     1   FCTR=(RCRNG-RNGES(IM1))/(RNGES(I)-RNGES(IM1))
60    WRITE(LOGFIL,50)
     1  KOUNT,NHOP,I,sngl(ELF),sngl(EL(IM1)),sngl(EL(I))
     2  ,sngl(RNGES(IM1)),sngl(RNGES(I)),sngl(SLOPE)
     3  ,sngl(GTIM(IM1)),sngl(GTIM(I))
     4  ,sngl(GTIM(IM1)+FCTR*(GTIM(I)-GTIM(IM1)))
50    FORMAT(3I6,3F11.5,2F11.4,F8.4,3F11.5)
100   CONTINUE
C
      END
 
C ----------------------------------------------------------------------
C			Function BINARY
C
C ----------------------------------------------------------------------
 
       double precision FUNCTION BINARY(X)
	
 	real*8 x
 
       BINARY=1.0
       IF(X.GE.0.0) RETURN
       BINARY=-1.0
       END
 
C ----------------------------------------------------------------------
C			Subroutine KURVE
C
C ----------------------------------------------------------------------
 
      SUBROUTINE KURVE(X,Y,NPTS,MARK,XMIN,XMAX)
 
	implicit double precision (a-h,o-z)
 
C+P   (KURVE) MODIFIED VRS OF DISSPLA 'CURVE'. DOES CLIPPING AND
C+    MARKER SELECTION.
      REAL*8 XW(333),YW(333),X(1),Y(1)
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      equivalence (w(436),rmax)
c
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      XL=XMIN
      XR=XMAX
      YB=0.0
      YT=RMAX
      NW=333
      CALL CPLOT(X,Y,NPTS,XW,YW,NW)
      IF(NW.LE.0) RETURN
C
      IF(MARK.NE.0) CALL MARKER(MARK)
      N=0
      IF(MARK.NE.0) N=-1
      CALL CURVE (XW,YW,NW,N)
      END
 
C ----------------------------------------------------------------------
C			Subroutine PLOTINI
C
C ----------------------------------------------------------------------
 
      SUBROUTINE PLOTINI(DHITE,SYMRED)
	
	implicit double precision (a-h,o-z)
 
C+P   (PLOTTING) INITIALIZE PLOTTING FOR EIGEN PROGRAM
C
C     COMMON DECK "HDR" INSERTED HERE
      CHARACTER*10 INITID*80,DAT,TOD
      COMMON/HDR/SEC
      COMMON/HDRC/INITID,DAT,TOD
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
C
C.....Write a code to TXTOUT indicating that eigen plots are being made.
C
      write (ndevtxt,111) 40
111   format (i2)
 
      CALL COMPRS
      CALL GRACE(0.,0.)
      CALL PHYSOR(2.,1.3)
      CALL PAGE(12.,10.)
      CALL HEIGHT(DHITE)
      CALL SETANN
      S=SYMRED
C  DEFAULT TO 15% REDUCTION ON POINT SIZE OVER DISSPLAY'S DEFAULT
      IF(S.EQ.0.0) S=.15
      CALL SCLPIC(S)
      END
 
C ----------------------------------------------------------------------
C			Subroutine PLTTLT
C
C ----------------------------------------------------------------------
 
      SUBROUTINE PLTTLT(NCODE)
	
	implicit double precision (a-h,o-z)
 
C+P   (PLTTLT) PUT TITLES ON PLOTS FOR EIGEN PROGRAM
C
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
	character id*80
      common /ident/id
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      REAL*8 NRSLMI
      EQUIVALENCE (W(80),APOG) ,(W(401),PTYP) ,(W(403),ELMIN)
     1   ,(W(404),ELMAX) ,(W(405),ELVINC) ,(W(406),RGINC)
     2   ,(W(407),HOPLME) ,(W(408),NRSLME)
     3   ,(W(409),GAPLO), (W(410),GAPHI), (W(411),SYMRED)
     4   ,(W(428),GRPMIN) ,(W(429),GRPMAX) ,(W(430),GRPVINC)
     5   ,(W(432),HOPLMG) ,(W(433),NRSLMG)
     6   ,(W(434),RFSPEED) ,(W(435),TMFT) ,(W(436),RMAX)
     &   ,(W(437),rmin)
     7   ,(W(453),RCRNG) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
      COMMON/GRPC/GRPN,GRPX
C
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "CODES" INSERTED HERE
      COMMON/CODES/ICODRE,ICODRT,ICODEI
C
      COMMON/RAYDEV/NRYIND,NDEVTMP,NFRMAT,NDEVBIN,NDEVTXT
 
      CHARACTER TDID*30
      COMMON/CDATIM/TDID
	character label*80	
	equivalence (label,ilabel)
C
      DATA NTCHRS/21/
C
C.....Changed 7.8 by 6.9 to 5.65 by 5.0 for plot dimensions.
 
C.....Write a code to TXTOUT indicating that eigen plots are being made.
C
      write (ndevtxt,111) 40
111   format (i2)
 
      IF(NCODE.EQ.ICODRE) CALL TITLE (TDID,NTCHRS
     1     ,'ELEV. ANGLE OF TRANSMISSION (#DEG!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
      IF(NCODE.EQ.ICODRT.AND.TMFT.EQ.1.) CALL TITLE (TDID,NTCHRS
     1     ,'RELATIVE TRAVEL TIME (#SEC!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
      IF(NCODE.EQ.ICODRT.AND.TMFT.NE.1.) CALL TITLE (TDID,NTCHRS
     1     ,'RELATIVE TRAVEL TIME (#MS!)$',100
     2     ,'RANGE (#KM!)$',100
     3     ,5.65d0,5.0d0)
 
 
CCCCC      CALL FRAME  !!!!!
      ELINC=DEGS*ELVINC
      IF(ELVINC.LE.0.0) ELINC=DEGS*(ELMAX-ELMIN)/10.0
      GRPINC=GRPVINC
      IF(GRPVINC.LE.0.0) GRPINC=(GRPMAX-GRPMIN)/10.0
 
C.....Add the following text to TXTOUT.  This is a subheading for the
C     plot.  Eigen plots are essentially identical to profile plots,
C     therefore there must be a subheading for PSGRAPH.
 
      write (ndevtxt,100) 3, 0, 0
100   format (i2,2i4)
      write (ndevtxt,200) 0,80
200   format (2i3)
      write (ndevtxt,300) id
300   format (a)
 
      IF(NCODE.EQ.ICODRE)
     1   CALL GRAFB(DEGS*ELMIN,ELINC,DEGS*ELMAX,rmin,RGINC,RMAX)
      IF(NCODE.EQ.ICODRT)
     1  CALL GRAFB(GRPMIN*TMFT,GRPINC*TMFT,GRPMAX*TMFT,rmin,RGINC, RMAX)
C
 
      RETURN
      END
 
C ----------------------------------------------------------------------
C			Function F
C
C ----------------------------------------------------------------------
 
 
      character FUNCTION F(C)
 
	implicit double precision (a-h,o-z)
 
C+P   MODIFY ANNOTATIONS FOR EIGEN PROGRAM DEPENDING ON ANNOTATION
C+    MODEL.
      CHARACTER C*(*)
C
      CALL ANNFIL(C,F)
      RETURN
      END
 
C ----------------------------------------------------------------------
C			Function GETRAYS
C
C ----------------------------------------------------------------------
 
      FUNCTION GETRAYS(NDEV,NRUN)
 
	implicit double precision (a-h,o-z)
 
C+
C+P   (GETRAYS) READ NRUN'TH RAYSET FROM FILE UNIT NDEV.
C+    IF LAST READ WENT PAST NRUN, DO RE-WIND
C+
C
      DATA LRUN/1000/
      IF(NRUN.LE.LRUN) THEN
          REWIND NDEV
          LRUN=1
      ENDIF
      LR=LRUN
      LRUN=NRUN+1
      DO 100 I=LR,NRUN
        GETRAYS=GETRAY1(NDEV)
        IF(GETRAYS.NE.0.0) RETURN
100     CONTINUE
        RETURN
        END
C8888888888888
C ----------------------------------------------------------------------
C			Function GETRAY1
C
C ----------------------------------------------------------------------
      FUNCTION GETRAY1(NDEV)
C
	implicit double precision (a-h,o-z)
C
C+
C+P   (GETRAY1) READ SINGLE RAYSET FROM FILE UNIT NDEV.
C+    ASSUME FILE IS POSITIONED AT OR BEFORE TRANSMIT RAY.
C+    READ WILL TERMINATE AT END-OF-RAYSET OR FILE EOF
C+
C
      CHARACTER FID*3
      LOGICAL NEWREAD
      CHARACTER LINE*80,C,CEOF,TYPE,rcvr,TRANS,STATE,MINDIS
      REAL*8 AZIMUTH,ELEVA,RANGE
C     COMMON DECK "CONST" INSERTED HERE
      COMMON/PCONST/CREF,RGAS,GAMMA
      COMMON/MCONST/PI,PIT2,PID2,DEGS,RAD,ALN10
C     COMMON DECK "WWE" INSERTED HERE
      PARAMETER (NWARSZ=1000)
	character id*80
      COMMON/WW/dum(10),MAXW,W(NWARSZ)
      REAL*8 NRSLMI
      EQUIVALENCE (W(80),APOG) ,(W(401),PTYP) ,(W(403),ELMIN)
     1   ,(W(404),ELMAX) ,(W(405),ELVINC) ,(W(406),RGINC)
     2   ,(W(407),HOPLME) ,(W(408),NRSLME)
     3   ,(W(409),GAPLO), (W(410),GAPHI), (W(411),SYMRED)
     4   ,(W(428),GRPMIN) ,(W(429),GRPMAX) ,(W(430),GRPVINC)
     5   ,(W(432),HOPLMG) ,(W(433),NRSLMG)
     6   ,(W(434),RFSPEED) ,(W(435),TMFT) ,(W(436),RMAX)
     7   ,(W(453),RCRNG) ,(W(454),SLPCUT)
     8   ,(W(456),HOPLMI) ,(W(457),NRSLMI) ,(W(458),ESET)
C     COMMON DECK "CRAYCOM" INSERTED HERE
      COMMON/RAYCOM/NRAYS,MXRYS,MXHOP,MXXHOP,NRSLIM,ICODE
      COMMON/RAYCOM/NHOPS(333),MNHOPS(333)
C     COMMON DECK "CRAYDAT" INSERTED HERE
      COMMON/RAYDAT/AZIMUTH(333),ELEVA(333),RANGE(10000)
      COMMON/RAYDAT/GRPTIM(10000)
      COMMON/RAYDAT/NELMTS
C
C     COMMON DECK "CEFILES" INSERTED HERE
c      COMMON/EFILES/INFILE,IOTFIL,LOGFIL,NCNTRL,NRYSFL
C
C
      DATA TYPE,rcvr,MINDIS,TRANS,CEOF/2*'R','M','T',' '/
      DATA NEWREAD/.TRUE./
      DATA UNDEF,BELEVA/2*1.E25/
C
      GEOF=0.0
      NRSETS=0
      NRAYS=0
      MX=0
      BAZIMUT=UNDEF
C+
C+    RESET 'STATE' FLAG USED TO INITIATE NEW RAY SERIES
C+
      STATE=CEOF
C  CONVERT ANGLES TO DEGREES
      ELMIND=DEGS*ELMIN
      ELMAXD=DEGS*ELMAX
      temp1 = fud(elmind)
      temp2 = fud(elmaxd)
C
C+
C+    LOOP THRU ALL INPUT LINES
C+
      DO 500 NLINE=1,100000
C+
C+    ATTEMPT TO FIND NEXT TRANSMIT RAY
C+
C
C.....READ A LINE FROM FILE, EXIT ON EOF
      DO 10 I=1,300
           if(newread) READ(NDEV,'(A)',IOSTAT=ios) LINE
           if (ios.lt.0) then
                geof = 1.0
                go to 900
           endif
C+
           newread = .true.
C+
C+    USING CHARACTER CODE IN 80TH COLUMN LOCATE NEXT TRANSMIT SERIES
C+    C can be 'T', 'R', 'M'.  The usual situation is:  a single line
C     with 'T' in 80th column followed by numerous (NHOP, i believe) lines
C     with 'R' in 80th column.  If 'M' is in the 80th column, this indicates
C     that the line is skipped and there will be a gap in the plot.
C
          C=LINE(80:80)
C+
C+    IF TRANSMIT RAY NOT SEEN, IGNORE OTHER LINES OF INPUT
C+
          IF(STATE.eq.TRANS.and.c.eq.ceof) goto 900
c
          IF(C .NE. TRANS) GO TO 15
          INC=1
          IF(NRAYS.gt.0) then
ccccc              print*,'nrays = ',nrays
              NHOPS(NRAYS)=NHOP+NMIND
              MNHOPS(NRAYS)=MINHOP
              MINHOP=0
              IF(NHOP.LE.0) then
                   nhops(nrays) = 9999
ccccc              INC=0
              endif
              IF(NRAYS.ge.MXRYS) stop' NUMBER OF RAYS EXCEEDS LIMIT'
          endif
          NRAYS=NRAYS+INC
C
C  RESET ALL THE COUNTERS FOR RAYSETS
          NMIND=0
          NHOP=0
          MINHOP=0
C+
C+I   READ AZIMUTH AND ELEVATION OF NEW RAY
C+
          READ(LINE,2222) FID,AZIMUTH(NRAYS),ELEVA(NRAYS)
          IF(BAZIMUT.NE.UNDEF.AND.BAZIMUT.NE.AZIMUTH(NRAYS)) then
               print*,'bazimut undef'
               GO TO 480
          endif
          IF(BELEVA.EQ.ELEVA(NRAYS)) then
              print*,'beleva=eleva(nrays)'
              gO TO 490
          endif
          BELEVA=ELEVA(NRAYS)
          BAZIMUT=AZIMUTH(NRAYS)
C
c
          IF(BELEVA.LT.temp1.OR.BELEVA.GT.temp2) then
              print*,'failed temp tst'
              GO TO 490
          endif
C
          IF(STATE.EQ.TRANS) then
              GO TO 500
          endif
C
15        IF(STATE.NE.TRANS) GO TO 155
          IF(C .EQ. rcvr) then
               GO TO 100
          endif
          IF(C .eq. MINDIS) then
               nmind = nmind + 1
          endif
 
C  SET MINIMUM SETS FLAG
 
          GO TO 10
C
155       IF(C.eq.TRANS) then
             STATE=TRANS
             IF(MX.NE.0) then
                 print*,'mx .ne.0,going to 500'
                 GO TO 500
             endif
             READ(LINE,3333) MX
             PRINT 4444,MXXHOP,MX
             MXHOP=MIN0(MX,MXXHOP)
          endif
 
10    CONTINUE  ! i=1,300
C
      stop'QUESTIONABLE RAYSET FILE FORMAT'
100   IF(NHOP+NMIND.GE.MXXHOP) then
           print*,'nhop,nmind,mxxhop = ',nhop,nmind,mxxhop
           GO TO 490
      endif
 
      NHOP=NHOP+1
      NN=NHOP+NMIND
      IF(NHOP.EQ.1) MINHOP=NN
      LNHOP=(NRAYS-1)*MXXHOP+NN
C+
C+    READ RANGE,DELAY AND HOP NUMBER FOR EACH RAY OF SERIES
C+
      READ (LINE,5555) RANGE(LNHOP), GRPTIM(LNHOP), NHOPIN
      if (inc.eq.0) then
           range(lnhop) = undef
      endif
 
C
C+
C+    COMPLAIN IF HOP NUMBER DOES NOT COMPARE WITH RUNNING VALUE
C+
      IF(NHOPIN.EQ.NN) GO TO 500
c        WRITE(LOGFIL,1200) NRAYS,NHOP,NHOPIN,NMIND
c1200  FORMAT(' FOR RAY NO',I3,',NHOP=',I10
c    1      ,' ,NHOPIN=',I4,' ,NMIND=',I4,'..LOGIC ERROR1')
         STOP 'logic error, see code'
C
C     CHANGE OF AZIMUTH SO HOLD THIS SET
480   NEWREAD=.FALSE.
      NRAYS=NRAYS-1
      RETURN
C
C     ABORT THIS RAY
490   print*,'state = ceof'
      STATE=CEOF
500   CONTINUE
      PRINT *,' LINE LIMIT ON RAYSETS REACHED'
C
900   GETRAY1=GEOF
      IF(NRAYS.LE.0) then
           RETURN
      endif
      IF(NHOP.gt.0) then
           NHOPS(NRAYS)=NHOP+NMIND
           MNHOPS(NRAYS)=NMIND+1
           RETURN
      endif
      NRAYS=NRAYS-1
 
2222      FORMAT(A,T44,2F10.5)
3333         FORMAT(T78,I2)
4444         FORMAT(' HOPLIM=',I3,' MXHOP=',I4)
5555  FORMAT(T10,F9.4,T37,F10.5,T78,I2)
      RETURN
      END
c88888888888
C ----------------------------------------------------------------------
C			Function FUD
C
C ----------------------------------------------------------------------
 
 
      FUNCTION FUD(X)
 
	implicit double precision (a-h,o-z)
 
C     RETURN VALUE SLIGHTLY LOW FOR NEGATIVE VALUES, SLIGHTLY HIGH
C     FOR POSITIVE VALUES. USED WHEN ROUNDOFF PROBLEMS MESS COMPARES.
C
      FUD=0.0
C***      IF(X.LT.0.0) FUD=1.001*X
C***      IF(X.GT.0.0) FUD=.999*X     !   <-- NO!
C*** Changed by A Newhall    Bug?
      FUD=1.001*X
 
      RETURN
      END
 
C ----------------------------------------------------------------------
C			Subroutine CPLOT
C
C ----------------------------------------------------------------------
 
 
      SUBROUTINE CPLOT(X,Y,NP,XW,YW,NW)
 
	implicit double precision (a-h,o-z)
 
C+P   (CPLOT) CLIPPED PLOT LINE PROGRAM. USES STANDARD CLIPPING
C+    ALGORITHM(CLIP) TO DETERMINE BREAKPOINTS AND INSERTS ANY
C+    LINE SEGMENTS TO PRODUCE CONNECTED LINE.
C
      INTEGER CLIP
      LOGICAL OUTSIDE
      REAL*8 X(*),Y(*),XW(*),YW(*)
C
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
      DATA XS,YS/2*0.0/
C
C     HANDLE SINGLES SPECIALLY
      IF(NP.EQ.1) GO TO 180
C
C     FIRST TRIM THE POINTS BASED ON Y VALUES
      NPP=NP+1
      NM=NP
      CALL YTRIM(Y,YT,NM)
      NM=NM-1
      N=1
      IF(NM.LT.1) GO TO 140
C
      DO 100 I=1,NM
      XI=X(I)
      YI=Y(I)
      XX=X(I+1)
      YY=Y(I+1)
      OUTSIDE=CLIP(XI,YI,XX,YY).NE.0
      IF(OUTSIDE) GO TO 100
      IF((XS.EQ.XI.AND.YS.EQ.YI).OR.N.EQ.1) GO TO 50
      IF(N.EQ.NW) then
ccccc      CALL STOPIT('LOGIC ER 1')
           stop 'n = nw' !!??????
       endif
C      WRITE(6,111) I,sngl(XI),sngl(YI),sngl(XS),sngl(YS),
C     1             sngl(XX),sngl(YY)
111   FORMAT(I5,10F8.2)
      XW(N)=XS
      YW(N)=YS
      N=N+1
50    XS=XX
      YS=YY
      IF(N.EQ.NW) then
ccccc      CALL STOPIT('LOGIC ER 1')
           stop 'n = nw' !!??????
      endif      !!!CALL STOPIT('LOGIC ER 2')
      XW(N)=XI
      YW(N)=YI
      N=N+1
100   CONTINUE
C
      IF(.NOT.OUTSIDE) GO TO 150
140   NW=N-1
      CALL YTRIM(YW,YT,NW)
      RETURN
150   XW(N)=XX
      YW(N)=YY
      NW=N
      CALL YTRIM(YW,YT,NW)
C
      RETURN
C
180   NW=0
      IF(X(1).GT.XR) RETURN
      IF(X(1).LT.XL) RETURN
      IF(Y(1).GT.YT) RETURN
      IF(Y(1).LT.YB) RETURN
      XW(1)=X(1)
      YW(1)=Y(1)
      NW=1
      END
 
C ----------------------------------------------------------------------
C			Subroutine YTRIM
C
C ----------------------------------------------------------------------
 
      SUBROUTINE YTRIM(Y,YT,NP)
 
	implicit double precision (a-h,o-z)
 
C+P   (YTRIM) ELIMINATES SEGMENTS OF LINE FROM HIGH END WHOSE
C+    Y VALUES EXCEED GIVEN VALUE.
      REAL*8 Y(*)
      IF(NP.LT.2) RETURN
      NPP=NP+1
      K=NPP
      DO 10 I=1,NP
         IF(Y(NPP-I).LT.YT) GO TO 15
10    K=K-1
15    NP=MIN0(K,NP)
      END
 
C ----------------------------------------------------------------------
C			Function CLIP
C
C ----------------------------------------------------------------------
 
      INTEGER FUNCTION CLIP(X1,Y1,X2,Y2)
 
	implicit double precision (a-h,o-z)
 
C+P   (CLIP) PERFORM STANDARD CLIPPING ON SINGLE LINE SEGMENTS
C+    (X1,Y1;X2,Y2) AGAINST WINDOW (XL,YB;XR,YT). 'CLIP' TAKES
C+    ON NONZERO VALUE IF LINE SEGMENT IS COMPLETELY OUTSIDE WINDOW.
C
      LOGICAL SWAPPED
      INTEGER C1,C2
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      SWAPPED=.FALSE.
      C2=KODE(X2,Y2)
C
C     BEGIN "COMPUTE" LOOP
C
C     WE TREAT PT ONE AS SPECIAL
10    C1=KODE(X1,Y1)
      CLIP=IAND(C1,C2)
ccccc111   FORMAT(4F8.2,3(1X,O4),L5)  !!!!!not used!?
C     IF BOTH END PTS ARE OFF-SCREEN, SO IS LINE
      IF(CLIP.NE.0) GO TO 100
C     IF BOTH END PTS ARE ON-SCREEN SO IS THE LINE
      IF(C1.EQ.0 .AND. C2.EQ.0) GO TO 100
C
      IF(C1.NE.0) GO TO 15
C     AT LEAST ONE END POINT IS OFF-SCREEN, MAKE SURE PT ONE IS ONE
      SWAPPED=.NOT.SWAPPED
      I=C1
      C1=C2
      C2=I
      CALL SWAP2(X1,X2,Y1,Y2)
C
C     SUCCESSIVELY REPLACE LINE BY ITS PROJECTIONS ALONG PRIMARY DIVIDER
15    IF(IAND(C1,1).EQ.0) GO TO 20
        Y1=Y1+(Y2-Y1)*(XL-X1)/(X2-X1)
        X1=XL
        GO TO 10
20    IF(IAND(C1,8).EQ.0) GO TO 30
        Y1=Y1+(Y2-Y1)*(XR-X1)/(X2-X1)
        X1=XR
        GO TO 10
30    IF(IAND(C1,64).EQ.0) GO TO 40
        X1=X1+(X2-X1)*(YB-Y1)/(Y2-Y1)
        Y1=YB
        GO TO 10
40    IF(IAND(C1,512).EQ.0) GO TO 10
        X1=X1+(X2-X1)*(YT-Y1)/(Y2-Y1)
        Y1=YT
C
      GO TO 10
100   IF(SWAPPED) CALL SWAP2(X1,X2,Y1,Y2)
C
      END
 
C ----------------------------------------------------------------------
C			Subroutine SWAP2
C
C ----------------------------------------------------------------------
 
      SUBROUTINE SWAP2(X1,X2,Y1,Y2)
 
	implicit double precision (a-h,o-z)
 
C
        T=X1
        X1=X2
        X2=T
        T=Y1
        Y1=Y2
        Y2=T
      END
 
C ----------------------------------------------------------------------
C			Function KODE
C
C ----------------------------------------------------------------------
 
      FUNCTION KODE(X,Y)
 
	implicit double precision (a-h,o-z)
 
C
      COMMON/CALTC/INI,FX,FY,XL,XR,YB,YT
C
      IF(X.GE.XL) GO TO 10
        K=1
        GO TO 100
10    IF(X.LE.XR) GO TO 50
        K=8
        GO TO 100
C
50    K=0
100   IF(Y.GE.YB) GO TO 110
        K=IOR(64,K)
      GO TO 200
110   IF(Y.LE.YT) GO TO 200
        K=IOR(512,K)
C
200   KODE=K
      RETURN
      END

