      SUBROUTINE PLOTVEL(A,TOBS,TCURV,V1OBS,V1CURV,V2OBS,V2CURV,V1MOD,
     1                   V2MOD,DELV1,DELV2,DELVSQ,
     2                   NVPTS,STARNAME,FK5,DO_1,DO_2)

C PLOTS RADIAL VELOCITY DATA, MODEL, AND RESIDUALS ON THE SCREEN
C
C 3 Sept 1991     J. T. Armstrong

      IMPLICIT UNDEFINED (A-Z)

      INCLUDE 'VELPARM.INC'

      REAL*8 A1_SINI,A2_SINI,V0,DV0DT
      REAL*8 TOBS(MAXVDATA),V1OBS(MAXVDATA),V2OBS(MAXVDATA)
      REAL*8 V1MOD(MAXVDATA),V2MOD(MAXVDATA),DELV1(MAXVDATA)
      REAL*8 DELV2(MAXVDATA),DELVSQ
      REAL*8 TCURV(MAXVMOD),V1CURV(MAXVMOD),V2CURV(MAXVMOD)
      REAL*8 T,A(11),V1,V2
      REAL*8 E_TRUE,PERIOD,EPOCH,ARG_PERI

      INTEGER*4 I,FK5,NVPTS,IER

      CHARACTER*1 ANSWER
      CHARACTER*128 PLOTFILE
      CHARACTER*80 STARNAME, COMMENT
      CHARACTER*100 TITLE(7)

      LOGICAL DO_1,DO_2

C VARIABLES FOR PLOTTING: NOTE THAT PGPLOT WANTS REAL*4

      REAL*4 TMAX,TMIN,VMAX,VMIN,PHMIN,PHMAX,VLO,VHI,TLO,THI
      REAL*4 TPLOT(MAXVMOD),V1CUR4(MAXVMOD),V2CUR4(MAXVMOD)
      REAL*4 TPL,VPLOT
      LOGICAL HCOPY
      INTEGER*4 IPLOT,XR,YR

C
C    +++++++++++++++++++PLOT.INC+++++++++++++++++++
C  To make hardcopies with the cal-tek graphics package, the user
C  must specify a file to store the graphics commands and then
C  copy /b those files to the printer.  For each hardcopy plot a 
C  new file name must be specified.  This include statement sets
C  up a file name that will be incremented each time a plot is
C  made in any subroutine.
C
      CHARACTER*4 PLOT_DEV
      INTEGER*4   NPLOTS
      COMMON /HCPBLK/ NPLOTS
      COMMON /PFILE  / PLOT_DEV

C    ++++++++++++++++++++++++++++++++++++++++++++++

      EXTERNAL PGBEGIN
      

      E_TRUE   = A(2) 
      PERIOD   = A(3)
      EPOCH    = A(4)
      ARG_PERI = A(7)
      A1_SINI  = A(8)
      A2_SINI  = A(9)
      V0       = A(10)
      DV0DT    = A(11)

C SET UP GRAPHICS ENVIRONMENT

C MAXIMA AND MINIMA
      TMIN = 0.
      TMAX = PERIOD
      PHMIN = 0.
      PHMAX = 1.
      VMAX = -1000.
      VMIN = 1000.
D     WRITE(OUTC,*) ' VMIN, VMAX = ',VMIN,VMAX
D     WRITE(OUTC,*) ' TMIN, TMAX = ',TMIN,TMAX
      DO I=1,NVPTS
        IF(V1OBS(I).LT.999.) THEN
          IF(V1OBS(I).GT.VMAX) VMAX=V1OBS(I)
          IF(V1MOD(I).GT.VMAX) VMAX=V1MOD(I)
          IF(V1OBS(I).LT.VMIN) VMIN=V1OBS(I)
          IF(V1MOD(I).LT.VMIN) VMIN=V1MOD(I)
        END IF
        IF(V2OBS(I).LT.999.) THEN
          IF(V2OBS(I).GT.VMAX) VMAX=V2OBS(I)
          IF(V2MOD(I).GT.VMAX) VMAX=V2MOD(I)
          IF(V2OBS(I).LT.VMIN) VMIN=V2OBS(I)
          IF(V2MOD(I).LT.VMIN) VMIN=V2MOD(I)
        END IF
      END DO
      DO I=1,MAXVMOD
          IF(V1CURV(I).GT.VMAX) VMAX=V1CURV(I)
          IF(V1CURV(I).LT.VMIN) VMIN=V1CURV(I)
          IF(V2CURV(I).GT.VMAX) VMAX=V2CURV(I)
          IF(V2CURV(I).LT.VMIN) VMIN=V2CURV(I)
      END DO
D     WRITE(OUTC,*) ' VMIN, VMAX = ',VMIN,VMAX
C TITLES
      IF(DO_1.AND..NOT.DO_2) THEN
          TITLE(1) = 'RADIAL VELOCITY FOR '//STARNAME//', STAR 1'
      ELSE IF(DO_2.AND..NOT.DO_1) THEN
          TITLE(1) = 'RADIAL VELOCITY FOR '//STARNAME//', STAR 2'
      ELSE
          TITLE(1) = 'RADIAL VELOCITIES FOR '//STARNAME
      END IF
      TITLE(2) = 'PHASE'
      TITLE(3) = 'PHASE (days)'
      TITLE(4) = 'RADIAL VELOCITY (km/s)'
      WRITE(TITLE(5),1001) A1_SINI,A2_SINI,E_TRUE,PERIOD
      WRITE(TITLE(6),1002) EPOCH,ARG_PERI/RPDEG,V0,DV0DT,DELVSQ

C HARDCOPY?

      WRITE(OUTC,*) ' HARD COPY?'
      READ(INC,*) ANSWER
      CALL CAPS(ANSWER)
      IF(ANSWER.EQ.'Y') THEN
          HCOPY=.TRUE.
      ELSE
          HCOPY=.FALSE.
      END IF

   40 IF(HCOPY) THEN
          PLOT_DEV = '/LPM'
          IER = PGBEGIN(0,PLOTFILE//PLOT_DEV,1,1)
      ELSE
          IER = PGBEGIN(0,'/EGA',1,1)
      END IF

      CALL PGSLW(1)
      CALL PGSCF(1)
      IPLOT = 1
      IF(IER.NE.1) GO TO 900

C MAKE THE GRAPHICS BOX

      CALL PGSCI(5)           ! CYAN
      CALL PGRNGE(VMIN,VMAX,VLO,VHI)
D     WRITE(OUTC,*) ' VMIN,VMAX,VLO,VHI = ',VMIN,VMAX,VLO,VHI
      CALL PGRNGE(TMIN,TMAX,TLO,THI)
D     WRITE(OUTC,*) ' TMIN,TMAX,TLO,THI = ',TMIN,TMAX,TLO,THI
D     WRITE(OUTC,*) ' ANY CHARACTER TO CONTINUE'
D     READ(INC,*) ANSWER
      CALL PGSCH(1.5)
      CALL PGENV(TLO,THI,VLO,VHI,0,1)
      CALL EGA_RESTORE_DEFAULT
      CALL PGSCI(1)           ! WHITE
      CALL PGSCH(1.5)
      CALL PGMTEXT('T',1.5,0.5,0.5,TITLE(1))
      CALL PGSCH(1.0)
      CALL PGMTEXT('B',3.0,0.5,0.5,TITLE(2))
      CALL PGMTEXT('L',3.0,0.5,0.5,TITLE(4))
      CALL PGMTEXT('B',4.0,0.5,0.5,TITLE(5))
      CALL PGMTEXT('B',5.2,0.5,0.5,TITLE(6))

C PLOT THE OBSERVED POINTS, MODEL POINTS, AND RESIDUAL POINTS

      CALL EGA_RESTORE_DEFAULT
      DO I=1,NVPTS
          IF(TOBS(I).LT.EPOCH) THEN
              TOBS(I)=TOBS(I)+PERIOD*(1+INT((EPOCH-TOBS(I))/PERIOD))
D             WRITE(OUTC,*) ' ADDING TO TOBS ',I,
D    1                       1+INT((EPOCH-TOBS(I))/PERIOD)
          END IF
      END DO
      DO I=1,NVPTS
          TPLOT(I) = MOD(TOBS(I)-EPOCH,PERIOD)
D         WRITE(OUTC,*) ' TPLOT = ',TPLOT(I)
      END DO
      CALL PGSCI(1)
      DO I=1,NVPTS
          TPL = TPLOT(I)
          IF(DO_1.AND.V1OBS(I).LT.999.) THEN
              VPLOT = V1OBS(I)
              CALL PGPOINT(1,TPL,VPLOT,17)
          END IF
          IF(DO_2.AND.V2OBS(I).LT.999.) THEN
              VPLOT = V2OBS(I)
              CALL PGPOINT(1,TPL,VPLOT,17)
          END IF
      END DO
C      CALL PGSCI(2)                   ! RED
C      DO I=1,NVPTS
C          TPL = TPLOT(I)
C          VPLOT = V1MOD(I)
C          CALL PGPOINT(1,TPL,VPLOT,22)
C          VPLOT = V2MOD(I)
C          CALL PGPOINT(1,TPL,VPLOT,22)
C      END DO
      DO I=1,NVPTS
          TPL = TPLOT(I)
          IF(DO_1) THEN
              VPLOT = DELV1(I)
              CALL PGSCI(3)               ! GREEN
              CALL PGPOINT(1,TPL,VPLOT,22)
          END IF
          IF(DO_2) THEN
              VPLOT = DELV2(I)
              CALL PGSCI(2)               ! RED
              CALL PGPOINT(1,TPL,VPLOT,22)
          END IF
      END DO

C PLOT THE MODEL

      CALL PGSCI(1)
C      CALL PGSLS(4)                   ! DOTTED
      DO I=1,MAXVMOD
C          TPLOT(I)=(I-1)*PERIOD/(MAXVMOD-1)
          TPLOT(I) = TCURV(I)-EPOCH
          V1CUR4(I) = V1CURV(I)
          V2CUR4(I) = V2CURV(I)
      END DO
      IF(DO_1) THEN
          CALL PGPOINT(MAXVMOD,TPLOT,V1CUR4,1)
      END IF
      IF(DO_2) THEN
          CALL PGPOINT(MAXVMOD,TPLOT,V2CUR4,1)
      END IF

      IF (HCOPY) THEN
          CALL PGEND
          HCOPY = .FALSE.
          GO TO 40
      END IF

      CALL EGA_RESTORE_DEFAULT
C      CALL LOCATE (0,1)
C      CALL WRITE_STRING('Any character to continue')
C      CALL PGCURSE(XR,YR,ANSWER)

      CALL PGEND
      RETURN

  900 WRITE(OUTC,*) ' Some kind of problem with PGBEGIN'
      RETURN

 1001 FORMAT(' a1,2 sin i = ',G12.3,', ',G12.3,' km;  e = ',F6.4,
     1       ';  P = ',F11.6)                             ! 70 CHARS
 1002 FORMAT(' T=JD',F11.3,'; |gw=',F7.2,'; V|d0|u=',F7.2,
     1       ' km/s; dV/dt=',F7.4,' km/s/yr; dV|u2|d=',F12.3) ! 86 CHARS

      END

