
      SUBROUTINE PLOTQ(Q,TQ,NQ,NORUN,LIM,HELJD,NFREQ,FREQ,AMP,PHASE,
     *                 AFIT,BFIT,CFIT,INTITL,ITPHAS,O,C,NSHIFT)     
C     
C     ROUTINE TO PLOT ARRAY Q AT TIMES TQ ON MULTIPLE FRAMES
C     
      COMMON /PGPL/PGPLDEV,IPGPL
      CHARACTER PGPLDEV*80
      COMMON /PLOTFL/PLNAME
      common /ecl/modee,param(30),nparam,ai
      CHARACTER PLNAME*12
      INTEGER TRMLEN
      EXTERNAL TRMLEN
      DOUBLE PRECISION TQ(NQ),HELJD(NORUN),HJD,TZERO,TWOPI,EPOCH,EPER 
      DOUBLE PRECISION O(1000),C(1000),ttt,parms(30),TBASE,TSHIFT
      DIMENSION Q(NQ),LIM(2,NORUN),DUM(4)
      DIMENSION FREQ(1),AMP(1),PHASE(1)   
      DIMENSION YYMIN(36),YYMAX(36),YYSCAL(36),YYAXIS(36)
      CHARACTER TITLE*80/' '/,XTITLE*6/' '/,INTITL*80   
      CHARACTER TITHJD*14,YTITLE*20,XOPT*20,YOPT*20
      LOGICAL SPECIAL
C
C     PARAMETERS:
C     -----------
C     XL:               LENGTH OF X-AXIS AS FRACTION OF WIDTH OF PAGE
C     YL:               LENGTH OF Y-AXIS OF EACH PANEL AS A FRACTION OF
C                       HEIGHT OF PAGE
C   
C     TZERO:            TIME OF START OF FIRST FRAME. THIS IS USUALLY ZERO
C                       IF THE TIME SCALE IS IN SEC. IF 999.0, THE PROGRAM
C                       SEARCHES FOR THE MINIMUM AND MAXIMUM FRACTIONAL
C                       JULIAN DATE IN EACH RUN AND USES THESE AS THE
C                       PLOTTING LIMITS FOR EACH PANEL
C     TFRAME:           TIME SPANNED BY EACH FRAME (IF 0.0, TFRAME=TIME
C                       SPANNED BY DATA TO BE PLOTTED). IF NEGATIVE, TFRAME
C                       IS SET POSITIVE AND TSHIFT IS SET TO 86400.0 TSHIFT
C                       IS USUALLY = TFRAME BUT IN THIS CASE IT IS SET TO
C                       1 DAY SO THAT EACH PANEL COVERS THE SAME FRACTIONAL
C                       JULIAN DAY
C
C     YMIN:             MINIMUM OF Y-AXIS OF EACH PANEL  (IF BOTH ARE ZERO,
C     YMAX:             MAXIMUM OF Y-AXIS OF EACH PANEL   AUTOSCALING IS DONE)
C
C     IFIT:          1  DRAW IN FIT (LSTSQ, NONLIN ETC.) AS A SOLID LINE
C
C
C                   99  REREAD IFIT AND THE FOLLOWING ADDITIONAL PARAMETERS:
C
C     IFIT:          1  DRAW IN FIT (LSTSQ, NONLIN ETC.) AS A SOLID LINE
C
C     I1:               INDEX IN Q OF FIRST POINT TO BE PLOTTED (IF 0, I1=1)
C     I2:               INDEX IN Q OF LAST  POINT TO BE PLOTTED (IF 0, I2=NQ)
C
C     LINTYP:        0  PLOT DATA AS FILLED CIRCLES
C                    1            AS CROSSES
C                    2            USING DIFFERENT SYMBOLS FOR EACH RUN
C                                               AND COLOURS 
C                   -1            AS CONTINUOUS LINES ("JOIN THE DOTS")
C
C     OVERLP:       0.0 NO VERTICAL OVERLAP BETWEEN EACH PANEL
C               NON-0.0 VERTICAL OVERLAP OF OVERLP*YL
C
C     IFIXY:         1  USE DIFFERENT VERTICAL SCALE FOR EACH PANEL
C                    0      SAME VERTICAL SCALE FOR EACH PANEL
C 
C     IFILE:         1  WRITE DATA OUT TO FILE data.dat AND FIT OUT TO FILE
C                       model.dat. TIMES ARE WRITTEN OUT IN SEC.
C                   -1  WRITE DATA OUT TO FILE data.dat AND FIT OUT TO FILE
C                       model.dat. TIMES ARE WRITTEN OUT AS HJD.
C 

      I1=0
      I2=0
      LINTYP=0    
      OVERLP=0.0
      IFIXY=0
      IFILE=0
      LABRIGHT=0

C
C     READ IN THE PARAMETERS
C
      READ *,XL,YL,TZERO,TFRAME,YMIN,YMAX,IFIT        
      IF(IFIT.EQ.99)THEN
C
C       READ THE EXTRA LINE OF PARAMETERS
C
        READ *,IFIT,I1,I2,LINTYP,OVERLP,IFIXY,IFILE
      END IF
      IF(IFIT.EQ.2)THEN
C
C       READ IN THE EPHEMERIS
C
        READ *,NFREQ,(FREQ(J),AMP(J),PHASE(J),J=1,NFREQ)
      END IF
  
      IF(IFIT.EQ.-3.OR.IFIT.EQ.-4)THEN
        READ *,NMOMEN,(O(J),J=1,NMOMEN)
      END IF

C     
C     CHECK THE INPUT
C     
      IF(LINTYP.LT.-1.OR.LINTYP.GT.2)THEN
        PRINT *,'### LINTYP ONLY ALLOWED VALUES IN RANGE'   
        PRINT *,'### -1 TO 2. RESETTING TO ZERO.'
        LINTYP=0
      END IF

      IF(IFILE.NE.0)THEN
C
C       WRITE THE DATA AND THE MODEL OUT TO A FILE
C  
        READ *,YADD
        OPEN(30,FILE='data.dat',STATUS='UNKNOWN')
        OPEN(31,FILE='model.dat',STATUS='UNKNOWN')
      END IF

      IF(IFIXY.EQ.1)THEN
C
C       THIS OPTION REQUIRES THE NUMBER OF FRAMES TO BE SPECIFIED
C       AS WELL AS THE YAXIS SIZE, THE Y AXIS MINIMUM AND THE Y
C       AXIS MAXIMUM FOR EACH FRAME
C
        READ *,NYYFRM
        DO J=1,NYYFRM
          READ *,YYAXIS(J),YYMIN(J),YYMAX(J)
        END DO
      END IF

      IF(XL.LE.0.0.OR.XL.GT.0.8)XL=0.8
      IF(YL.EQ.0.0)YL=0.4     

C
C     THIS ALLOWS A SUBSECTION OF Q TO BE PLOTTED EITHER
C     BY SPECIFYING WHICH POINTS ARE TO BE PLOTTED OR
C     BY SPECIFYING WHICH RUNS ARE TO BE PLOTTED (I1,I2<0).
C
      IF(I1.LT.0)THEN   
        I1=-I1    
        I2=-I2    
        IF(I1.LE.0)I1=1 
        IF(I2.GT.NORUN)I2=NORUN     
        NI=LIM(1,I1)    
        NF=LIM(2,I2)    
      ELSE  
        NI=I1     
        NF=I2     
      END IF
      IF(NI.LE.0)NI=1   
      IF(NF.GT.NQ)NF=NQ 
      IF(NF.EQ.0)NF=NQ  

C     
C     THIS SECTION DEALS WITH X-AXIS (TIME) SCALING
C     
      DELT=TQ(NI+1)-TQ(NI)    
      DO J=NI,NF-1
        IF(DELT.GT.ABS(TQ(J+1)-TQ(J)))DELT=ABS(TQ(J+1)-TQ(J))
      END DO
      TMAX=-1.0E30
      TMIN=1.0E30
      DO J=NI,NF
        IF(TMAX.LT.TQ(J))TMAX=TQ(J)       ! Minimum and maximum times
        IF(TMIN.GT.TQ(J))TMIN=TQ(J)       ! (in sec) of first and last pt
      END DO
      EPSILON=MIN(3.0*DELT,1000.0)        ! Slight offset to improve look
      EPSIL  =EPSILON                     ! of x-axes
      TBASE=TMIN                          ! Time at left edge of each frame
      TRANGE=TMAX-TMIN                    ! Time span of data to be plotted
      HJD=HELJD(1)+TMIN/86400.0D0         ! HJD at start of plot
      JD=HJD                              ! JD is integer Julian Date
                                          ! Additional variables include:
                                          ! TFRAME=time span of each frame
                                          ! TMINLAB, TMAXLAB are the times
                                          ! with which to label the x-axis

      IF(ITPHAS.EQ.0)THEN     
C
C       NORMAL TIME PLOTTING (AS OPPOSED TO PHASE PLOTTING)
C
        XTITLE='(Sec)'    
        LXTITL=5  
        TSHIFT=0.0
        IF(TZERO.NE.0.0)THEN
C
C         IT IS ASSUMED THAT TZERO IS THE HJD FROM WHICH PLOTTING
C         WILL BEGIN. THIS MEANS A RECOMPUTATION OF TMIN, TBASE ETC.
C
          IF(ABS(TZERO-999.0).LT.0.01)THEN
C
C           SPECIAL CASE WHERE PLOTTING IS DONE ON PANELS COVERING
C           THE SAME FRACTIONAL JULIAN DATE RANGE. THIS RANGE IS
C           DEFINED BY THE SMALLEST AND LARGEST FRACTIONAL JULIAN
C           DATE IN EACH RUN. THE VARIABLE TSHIFT IS THE TIME TO SHIFT
C           ALONG BY TO GO TO THE NEXT FRAME. THIS IS USUALLY = TFRAME
C           (I.E. FRAMES ARE "CONTIGUOUS"). HOWEVER, IN THIS CASE, TSHIFT
C           IS 86400.0 SEC = 1 DAY.
C
            TFRMAX=-1.0E30
            TFRMIN=1.0E30
            DO J=1,NORUN
              IF(TFRMIN.GT.HELJD(J)-INT(HELJD(J)))
     *          TFRMIN=HELJD(J)-INT(HELJD(J))
              IF(TFRMAX.LT.HELJD(J)-INT(HELJD(J)) +
     *                     (TQ(LIM(2,J))-TQ(LIM(1,J)))/86400.0)
     *          TFRMAX=HELJD(J)-INT(HELJD(J)) +
     *                 (TQ(LIM(2,J))-TQ(LIM(1,J)))/86400.0
            END DO
            TZERO=INT(HELJD(1))+INT(TFRMIN*100.0)/100.0
            TFRAME=INT(TFRMAX*100.0+1.0)/100.0-INT(TFRMIN*100.0)/100.0
            TSHIFT=86400.0
          END IF
          IF(ABS(TFRAME).GE.15.0)THEN
C           TFRAME IS TAKEN TO BE IN SEC
            IF(TFRAME.LT.0.0)THEN
C             TIME BETWEEN PANELS IS 86400.0 S (1 DAY)
              TFRAME=-TFRAME
              TSHIFT=86400.0   ! See TZERO=999 for explanation of TSHIFT
            END IF
            XTITLE='(Sec)'
            LXTITL=5
            TMINLAB=0.0
            TMAXLAB=TMINLAB+TFRAME
            HJD=TZERO
            JD=HJD
          ELSE
C           TFRAME IS TAKEN TO BE IN DAYS
            IF(TFRAME.LT.0.0)THEN
C             TIME BETWEEN PANELS IS 86400.0 S (1 DAY)
              TFRAME=-TFRAME
              TSHIFT=86400.0   ! See TZERO=999 for explanation of TSHIFT
            END IF
            XTITLE='(Days)'  
            LXTITL=6
            TMINLAB=TZERO-INT(TZERO)
            TMAXLAB=TMINLAB+TFRAME
            EPSIL=EPSIL/86400.0D0
            TFRAME=TFRAME*86400.0D0   
            JD=TZERO    
            HJD=JD
          END IF
          TMIN=(TZERO-HELJD(1))*86400.0D0  
          TRANGE=TMAX-TMIN    
          TBASE=TMIN    
          IF(TRANGE.GT.86400.0)LABRIGHT=1
        ELSE
C
C         IT IS ASSUMED THAT PLOTTING BEGINS FROM TIME OF FIRST POINT.
C         3 CASES NEED TO BE DISTINGUISHED: TFRAME=0, 0<TFRAME<=15  AND
C         15 <TFRAME
C
          IF(TFRAME.EQ.0.0)THEN
            TFRAME=TRANGE
            XTITLE='(Sec)'
            LXTITL=5
            TMINLAB=0.0
            TMAXLAB=TMINLAB+TFRAME
          ELSE IF(TFRAME.GE.15.0)THEN
C           TFRAME IS TAKEN TO BE IN SEC
            XTITLE='(Sec)'
            LXTITL=5
            TMINLAB=0.0
            TMAXLAB=TMINLAB+TFRAME
          ELSE IF(TFRAME.GT.0.0.AND.TFRAME.LE.15.0)THEN
C           TFRAME IS TAKEN TO BE IN DAYS
            XTITLE='(Days)'  
            LXTITL=6
            TMINLAB=0.0              
            TMAXLAB=TMINLAB+TFRAME
            EPSIL=EPSIL/86400.0D0
            TFRAME=TFRAME*86400.0D0   
          ELSE
            PRINT *,'### TFRAME UNDEFINED OR NEGATIVE.'
            PRINT *,'### FATAL ERROR - HALTING'
            STOP
          END IF
        END IF
      ELSE  
C
C       IT IS ASSUMED THAT SOME COMMAND LIKE TIMBAS TO PHASE FOLD
C       THE DATA HAS BEEN USED SO THAT ITPHAS=1 AND ALL THE NUMBERS IN
C       TQ ARE PHASES, NOT TIMES.
C
        XTITLE='PHASE'  
        LXTITL=5  
        IF(TFRAME.EQ.0.0)THEN 
          TFRAME=TRANGE
          TMINLAB=0.0
          TMAXLAB=TMINLAB+TFRAME
        ELSE IF(TFRAME.GT.0.0)THEN  
          TMIN=TZERO
c         TMAX=TZERO+TFRAME
c         TRANGE=TMAX-TMIN
          TBASE=TMIN
          TMINLAB=TMIN
          TMAXLAB=TMINLAB+TFRAME
        ELSE
          PRINT *,'### TFRAME UNDEFINED OR NEGATIVE.'
          PRINT *,'### FATAL ERROR - HALTING'
          STOP
        END IF    
        TSHIFT=TFRAME
      END IF
      IF(TSHIFT.EQ.0.0)TSHIFT=TFRAME   ! See TZERO=999 for definition of TSHIFT
      NFRAME=(TRANGE/TSHIFT)+0.99999
C     
C     CHECK NOW FOR FRAMES TO BE OMITTED DUE TO LONG TIME GAPS IN DATA  
C     
      IFRAME=0    
      MISS=0                  ! This is set to 1 if there are any missing
      TEND=TMIN+TFRAME-TSHIFT ! frames
      IPT=NI-1    
    5 IPT=IPT+1   
      IF(IPT.GT.NF)GO TO 15   
      NSUB=0
   10 IF(TQ(IPT).GT.TEND)THEN 
        IFRAME=IFRAME+1 
        TEND=TEND+TSHIFT   ! See TZERO=999 for definition of TSHIFT  
        NSUB=NSUB+1     
        GO TO 10  
      END IF
      IF(NSUB.GT.1)MISS=1     
      IF(NSUB.GT.0)NFRAME=NFRAME-NSUB+1   
      GO TO 5     
C     
   15 CONTINUE    
      IF(IFIXY.EQ.0)THEN
        YFRAME=NFRAME*YL*(1.0-OVERLP) + YL*OVERLP
      ELSE
        YFRAME=0.0
        DO J=1,NFRAME
          YFRAME=YFRAME+YYAXIS(J)
        END DO
      END IF
      IF(YFRAME.GT.0.8)THEN      
        PRINT *,'TOO MANY FRAMES (>PAGE HEIGHT)'  
        PRINT *,'NUMBER OF FRAMES IS ',NFRAME   
        PRINT *,'TFRAME IS ',TFRAME 
        PRINT *,'TSHIFT IS ',TSHIFT 
        PRINT *,'TRANGE IS ',TRANGE 
        PRINT *,'X-AXIS LENGTH IS ',XL    
        PRINT *,'Y-AXIS LENGTH IS ',YL    
        PRINT *,'TOTAL Y AXIS LENGTH IS ',YFRAME
        STOP
      END IF
      IF(NFRAME.LE.0)THEN     
        PRINT *,'CALCULATED NUMBER OF FRAMES <=0'     
        PRINT *,'RETHINK YOUR HORIZONTAL SCALING'     
        PRINT *,'NSUB= ',NSUB
        PRINT *,'TSHIFT= ',TSHIFT
        PRINT *,'TRANGE= ',TRANGE
        PRINT *,'TFRAME= ',TFRAME
        PRINT *,'NFRAME= ',NFRAME
        STOP
      END IF
      IF(IFIXY.EQ.1.AND.NYYFRM.NE.NFRAME)THEN
        PRINT *,'### WARNING WARNING WARNING ###'
        PRINT *,'CALCULATED NUMBER OF FRAMES FROM TIME RANGE: ',NFRAME
        PRINT *,'NOT EQUAL TO NO. SPECIFIED BY 99 OPTION: ',NYYFRM
        PRINT *,'WILL USE THE CALCULATED NUMBER FOR PLOTTING'
      END IF

C
C     THIS SECTION DEALS WITH Y AXIS SCALING
C
      QMAX=-1.0E30
      QMIN=1.0E30
      DO 20 J=NI,NF   
        IF(Q(J).GT.QMAX)QMAX=Q(J) 
        IF(Q(J).LT.QMIN)QMIN=Q(J) 
   20 CONTINUE  
C     REPLACE QMAX, QMIN IF YMAX, YMIN SPECIFIED 
      IF(YMIN.NE.0.0.OR.YMAX.NE.0.0)THEN
        QMIN=YMIN
        QMAX=YMAX
      END IF

      PRINT 25,NI,NF,NFRAME,OVERLP,XL,TFRAME,TMIN,HJD,LINTYP
   25 FORMAT(12X,'PLOTTING OF Q FROM POINTS  ',I6,' TO ',I6/
     *       12X,'ON ',I5,' FRAMES WITH OVERLAP FRACTION ',F5.2/
     *       12X,'X-AXIS ',0PF5.2,'  TIME SCALE ',1PE12.4/
     *       12X,'TIME AT START OF AXIS IS ',1PE12.4/
     *       12X,'CORRESPONDING TO HJD  ',0PF15.5/
     *       12X,'LINTYP                ',I15)
      IF(IFIXY.EQ.0)THEN
        PRINT 30,YL,QMAX,QMIN
   30   FORMAT(12X,'Y-AXIS ',0PF5.2,'     YMAX ',1PE12.4,
     *         '  YMIN ',E12.4)
      ELSE
        PRINT 35
   35   FORMAT(12X,'USER SPECIFIED SCALING FOR EACH FRAME:')
        PRINT 30,(YYAXIS(J),YYMAX(J),YYMIN(J),J=1,NFRAME)
      END IF
      PRINT 36,INTITL
   36 FORMAT(12X,'TITLE : ',A)
C     
C     START PLOTTING - INITIALISE IF NECESSARY
C
      IF(IPGPL.EQ.0)THEN
        IPGPL=1
        CALL PGLDEV
        CALL PGBEG(0,PGPLDEV,1,1)
      END IF
      CALL PGSCH(1.2)                     ! Set the character height
      CALL PGSLW(3)                       ! Set the line width
      VPXLO=(1.0-XL)/2.0
      VPXHI=VPXLO+XL  
      VPYLO=0.1 
      VPYHI=VPYLO+NFRAME*YL
      CALL PGSVP(VPXLO,VPXHI,VPYLO,VPYHI) ! Define viewport on view surface
C
C     WRITE A TINY DOT AT THE UPPER LEFT CORNER AND LOWER RIGHT CORNER
C     OF THE PLOT SO THE WINDOW SHOWS NICELY IN GHOSTVIEW OF THE PGPLOT.PS FILE
C
      CALL PGSCH(0.01)
c     CALL PGMTXT( 'T', 600.0, -0.95*VPXLO/(VPXHI-VPXLO), 0.0, '.' )
c     CALL PGMTXT( 'B', 200.0, 1.0+0.95*VPXLO/(VPXHI-VPXLO), 0.0, '.' )
c     CALL PGMTXT( 'T', 600.0, -0.5*VPXLO/(VPXHI-VPXLO), 0.0, '.' )
c     CALL PGMTXT( 'B', 200.0, 1.0+0.5*VPXLO/(VPXHI-VPXLO), 0.0, '.' )
      CALL PGSCH(1.2)
     
C     WRITE OUT THE TITLE. INCLUDE THE JULIAN DATE IF TZERO /= 0.0 OR
C     TITLE HAS THE CHARACTER @ IN IT AND PROVIDED LABRIGHT=0
      TITLE=INTITL
      JDLAB=0
      IF(TZERO.NE.0.0.AND.LABRIGHT.EQ.0)JDLAB=1
      DO J=1,TRMLEN(TITLE)
        IF(TITLE(J:J).EQ.'@')THEN
          TITLE(J:J)=' '
          IF(LABRIGHT.EQ.0)JDLAB=1
        END IF
      END DO
      CALL PGMTXT( 'T', 3.0, 0.0, 0.0, TITLE ) 

      IF(JDLAB.EQ.1.AND.ITPHAS.EQ.0)THEN     
        INTHJD=HJD
        IF(ABS(HJD-INTHJD).GT.1.0E-5)THEN
          IF(HJD.GT.100000.0)THEN
            WRITE(TITHJD,'(F13.5,''+'')')HJD   
          ELSE IF(HJD.GE.10000.0)THEN
            WRITE(TITHJD,'(F11.5,''+'')')HJD
          ELSE
            WRITE(TITHJD,50)HJD
   50       FORMAT(F10.5,'+')  
          END IF
c         CALL PGMTXT( 'T', 4.0, 1.0, 1.0, TITHJD )
        ELSE
          WRITE(TITHJD,51)INTHJD   
   51     FORMAT('JD ',I5,'+')  
          CALL PGMTXT( 'T', 1.5, 0.0, 0.0, TITHJD )
        END IF
      END IF
      CALL PGSCH(1.0) 
C     
C     NOW LOOP OVER THE FRAMES
C
C     FIRST REDEFINE THE VIEWPORT TO OCCUPY A VERTICAL HEIGHT EQUAL
C     TO THE FULL VIEWPORT DIVIDED BY THE NUMBER OF FRAMES
C       
      IFRAME=1    
      IPT=NI-1    
      TWOPI=3.141592653D0*2.0D0     
      DO 500 K=1,NFRAME 
        NEWFR=1   
        IF(IFIXY.EQ.1)THEN
C
C         SET THE VERTICAL SCALE FACTORS ON A FRAME BY FRAME BASIS
C         FROM THE ONES THAT WERE READ IN
C
          YL=YYAXIS(K)
          QMIN=YYMIN(K)
          QMAX=YYMAX(K)
        END IF

        IF(K.EQ.1)THEN
          YHI=VPYHI
        ELSE
          IF(IFIXY.EQ.0)THEN
            YHI=YHI-(1.0-OVERLP)*YL
          ELSE
            YHI=YHI-(1.0-OVERLP)*YYAXIS(K-1)
          END IF
        END IF
        YLO=YHI-YL
C
C       REDEFINE THE VIEWPORT (REGION OCCUPIED BY THE FRAME OF THE GRAPH)
C       TO BE THE NEXT FRAME TO BE PLOTTED   
C
        CALL PGSVP(VPXLO,VPXHI,YLO,YHI) 
C
C       DEFINE THE MAPPING OF USER UNITS TO THE VIEWPORT IN TERMS
C       OF TMINLAB AND TMAXLAB (FOR LABELLING PURPOSES ONLY)
C
        CALL PGSWIN(TMINLAB-EPSIL,TMAXLAB+EPSIL,QMIN,QMAX)
C
C       NOW DRAW PART OF THE BOX: Y-AXES ON BOTH SIDES AND FOR THE
C       FIRST FRAME THE TOP X-AXIS OR FOR THE LAST FRAME THE BOTTOM
C       X-AXIS.
C
        XOPT=' '
        IF(K.EQ.1)THEN
          XOPT='CMTS1'
        ELSE IF(K.EQ.NFRAME)THEN
          XOPT='BNTS1'
        END IF
        CALL PGBOX( XOPT, 0.0, 0, ' ', 0.0, 0 )
 
        YOPT=' '
        IF(OVERLP.EQ.0.0.OR.K.EQ.1)THEN
          YOPT='BCT'
          IF(IFIXY.EQ.0.AND.K.EQ.1)YOPT=YOPT(1:4)//'NV' 
          IF(IFIXY.EQ.1.AND.MOD(K,2).EQ.0)YOPT=YOPT(1:4)//'NV' 
          IF(IFIXY.EQ.1.AND.MOD(K,2).EQ.1)YOPT=YOPT(1:4)//'MV' 
        END IF
c       CALL MYBOX( ' ', XTICK, 0, YOPT, YTICK, 0 )
        CALL PGBOX( ' ', 0.0, 0, YOPT, 0.0, 0 )
C
C       IF OVERLP /= 0.0, MAKE SPECIAL PROVISION FOR DRAWING Y AXES AROUND
C       THE WHOLE FRAME (DO IT ONCE I.E. WHEN K=1)
C       
        IF(OVERLP.GT.0.0.AND.K.EQ.1)THEN
          YHI1=VPYHI
          YLO1=YHI1-YFRAME
          CALL PGSVP(VPXLO,VPXHI,YLO1,YHI1)
C         QMIN1 IS THE VALUE OF THE USER UNITS AT THE BOTTOM OF THE
C         WHOLE PANEL ASSUMING QMIN AND QMAX SCALE TO THE TOP FRAME
          QMIN1=QMAX-YFRAME/YL*(QMAX-QMIN)
          CALL PGSWIN(TMINLAB-EPSIL,TMAXLAB+EPSIL,QMIN1,QMAX)
C         NOW DRAW PART OF THE BOX: ONLY THE Y-AXES ON BOTH SIDES
          XOPT=' '
          YOPT=' '
C         YOPT='BCTSNV'
          YOPT='BC'
c         CALL MYBOX( XOPT, XTICK, 0, YOPT, YTICK, 0 )
          CALL PGBOX( XOPT, 0.0, 0, YOPT, 0.0, 0 )    
        END IF

C
C       DRAW THE X AXIS TITLE IF IT'S THE FIRST FRAME
C
        IF(K.EQ.1)THEN
          CALL PGSCH(0.8)
          CALL PGMTXT( 'T', 3.0, 0.5, 0.5, XTITLE )       
          CALL PGSCH(1.0)
        END IF
C
C       IF OVERLP=0.0, DRAW INTERIOR X-AXIS
C
        IF(OVERLP.EQ.0.0)THEN
          XOPT=' '
          YOPT=' '
          IF(K.EQ.1)THEN
            XOPT='BS'
          ELSE IF(K.EQ.NFRAME)THEN
            XOPT='CS'
          ELSE
            XOPT='BCS'
          END IF
          CALL PGBOX( XOPT, 0.0, 0, YOPT, 0.0, 0 )
        END IF
C
C       REDEFINE THE VIEWPORT ON THE VIEWING SURFACE AS WELL AS THE
C       MAPPING OF USER UNITS TO THE VIEWPORT IN TERMS OF THE ACTUAL 
C       UNITS USED (TBASE AND TBASE+TFRAME)
C
        CALL PGSVP(VPXLO,VPXHI,YLO,YHI) 
        TB=TBASE
        CALL PGSWIN(TB-EPSILON,TB+TFRAME+EPSILON,QMIN,QMAX)

        JD=HELJD(1)+TBASE/86400.0D0   
        WRITE(YTITLE,'(A3,I5)')'JD ',JD           
        IF(LABRIGHT.EQ.1)CALL PGMTXT( 'RV', 4.0, 0.5, 0.5, YTITLE )       
C     
  250   IF(IFIT.EQ.0)GO TO 295    
        IF(IFIT.EQ.-1)GO TO 280     
        IF(IFIT.EQ.-2)GO TO 290
        IF(IFIT.EQ.-4)GO TO 292
C     
C       PLOT THE FITTED FREQUENCIES IN FREQ, AMP AND PHASE  
C     
        FMAX=0.0  
        DO 255 J=1,NFREQ
          FMAX=MAX(FREQ(J),FMAX)  
  255   CONTINUE  
        IF(NFREQ.GT.0)THEN
          DELFIT=1.0/FMAX/20.0  
        ELSE
          DELFIT=TFRAME/20.0
          IF(IFIT.EQ.-5)THEN
            DELFIT=1.0
            PARMS(1)=AFIT
            PARMS(2)=BFIT
            PARMS(3)=CFIT
          END IF
          IF(IFIT.EQ.-6.OR.IFIT.EQ.-7.OR.  
     *       IFIT.EQ.-8)THEN
            DELFIT=(TQ(2)-TQ(1))
            DO J=1,NPARAM
              PARMS(J)=PARAM(J)
            END DO
          END IF
        END IF
        TIME=TBASE-DELFIT     
        IC=1
  260   TIME=TIME+DELFIT
        IF(TIME.GT.TBASE+TFRAME)THEN
          IF(IFIT.EQ.-3)THEN
            GO TO 292
          ELSE
            GO TO 295 
          END IF
        END IF
        VAL=AFIT+(BFIT*TIME)+(CFIT*TIME*TIME)
        IF(IFIT.EQ.-5)THEN
          TTT=TIME
          VAL=FUNC(TTT,3,PARMS,MODEE)
        END IF
        IF(IFIT.EQ.-6)THEN
          TTT=TIME
          VAL=FUNC1(TTT,NPARAM,PARMS,AI)
        END IF
        IF(IFIT.EQ.-7)THEN
          TTT=TIME
          VAL=FUNC2(TTT,NPARAM,PARMS,AI)
        END IF
        IF(IFIT.EQ.-8)THEN
          TTT=TIME
          VAL=FUNC3(TTT,NPARAM,PARMS,AI)
        END IF
        DO 275 J=1,NFREQ
          ARG=DMOD(TWOPI*FREQ(J)*TIME+PHASE(J),TWOPI) 
          VAL=VAL+AMP(J)*COS(ARG)   
  275   CONTINUE  
        IF(IFILE.EQ.1)WRITE(31,276)TIME,VAL+YADD
        IF(IFILE.EQ.-1)
     *    WRITE(31,276)(HELJD(1)+TIME/86400.0D0),VAL+YADD
  276   FORMAT(F20.7,1PE20.7)
        IF(IC.EQ.1)THEN     
          IC=0  
          CALL PGMOVE(TIME, VAL)    
        END IF  
        CALL PGDRAW(TIME, VAL)
        GO TO 260 
C     
  280   CONTINUE  
C     
C       THIS SECTION IS TO DRAW ARROWS OVER THE PLOT AT INTERVALS 
C       SPECIFIED BY THE FREQUENCY FREQ(1) STARTING AT TIME 
C       T=-PHASE(1)/2*PI*FREQ(1)    
C     
        TIME=-PHASE(1)/(TWOPI*FREQ(1))    
        TB=TBASE
        TIME=TIME+IFIX(TB*FREQ(1)-2)/FREQ(1) 
  285   TIME=TIME+1.0/FREQ(1) 
        IF(TIME.GT.TBASE+TFRAME)GO TO 295 
        IF(TIME.LT.TBASE)GO TO 285  
        T=(TIME-TBASE)/TFRAME 
        CALL PGSCH(0.5)
        CALL PGSAH(2,30.0,1.0)
        CALL PGARRO(TIME,QMAX,TIME,QMAX-0.035*NFRAME*(QMAX-QMIN))
        CALL PGSCH(1.0)
        GO TO 285 
C     
C       THIS SECTION IS TO MARK IN OBSERVED AND CALCULATED TIMES  
C       COMPUTED BY ROUTINE SOMINC. 
C
  290   CONTINUE
        DO 291 J=1,NSHIFT
          CC=(C(J)-HELJD(1))*86400.0D0
          OO=(O(J)-HELJD(1))*86400.0D0
          IF(CC.GE.TBASE.AND.CC.LE.TBASE+TFRAME)THEN
            T=CC
            CALL PGSCH(0.5)
            CALL PGSAH(2,30.0,1.0)
            CALL PGARRO(T,QMAX,T,QMAX-0.035*NFRAME*(QMAX-QMIN))
            CALL PGSCH(1.0)
          END IF
          IF(OO.GE.TBASE.AND.OO.LE.TBASE+TFRAME)THEN
            T=OO
            CALL PGSCH(0.5)
            CALL PGSAH(2,30.0,1.0)
            CALL PGARRO(T,QMAX,T,QMAX-0.035*NFRAME*(QMAX-QMIN))
            CALL PGSCH(1.0)
          END IF
  291   CONTINUE
C     
C       THIS SECTION IS TO MARK IN ARROWS OVER SOME OBSERVED
C       TIMES STORED IN O(J). THERE ARE NMOMEN OF THEM.
C
  292   CONTINUE
        DO 293 J=1,NMOMEN
          OO=(O(J)-HELJD(1))*86400.0D0
          IF(OO.GE.TBASE.AND.OO.LE.TBASE+TFRAME)THEN
            T=OO
            CALL PGSCH(0.5)
            CALL PGSAH(2,30.0,1.0)
            CALL PGARRO(T,QMAX,T,QMAX-0.035*NFRAME*(QMAX-QMIN))
            CALL PGSCH(1.0)
          END IF
  293   CONTINUE
C
C       NOW PLOT THE DATA     
C
  295   CONTINUE
        IF(LINTYP.EQ.0)THEN
          CALL PGSCH(0.4)
          INDEX=17
        ELSE IF(LINTYP.EQ.1.OR.LINTYP.EQ.2)THEN
          CALL PGSCH(0.8)
          INDEX=2
        END IF
        IPT=IPT+1 
        IF(IPT.GT.NF)GO TO 600
        IF(IFILE.EQ.1)
     *    WRITE(30,*)TQ(IPT),Q(IPT)+YADD
        IF(IFILE.EQ.-1)
     *    WRITE(30,*)(HELJD(1)+TQ(IPT)/86400.0D0),Q(IPT)+YADD
        T=TQ(IPT)
        IF(T.GT.TBASE+TFRAME)THEN 
  300     IFRAME=IFRAME+1     
          TBASE=TBASE+TSHIFT  
          IF(T.GT.TBASE+TFRAME)GO TO 300
          IPT=IPT-1     
          CALL PGSCH(1.0)
          GO TO 500     
        END IF    
        QPOS=Q(IPT)
        IF(LINTYP.EQ.0.OR.LINTYP.EQ.1)CALL PGPT(1,T,QPOS,INDEX)  
        IF(LINTYP.EQ.2)THEN
C         DECIDE WHICH RUN IT BELONGS TO
          DO KK=1,NORUN
            IF(IPT.GE.LIM(KK,1).AND.IPT.LE.LIM(KK,2))THEN
              IF(KK.EQ.1)INDEX=17
              IF(KK.EQ.2)INDEX=4
              IF(KK.GT.2)INDEX=KK-1
              CALL PGSCI(MOD(KK,5))
              GO TO 310
            END IF
          END DO
          INDEX=2
  310     CALL PGPT(1,T,QPOS,INDEX)
        END IF
        IF(LINTYP.EQ.-1)THEN
          IF(NEWFR.EQ.1)CALL PGMOVE(T,QPOS)
          CALL PGDRAW(T,QPOS)
        END IF
        NEWFR=0   
        CALL PGSCI(1)
        GO TO 295 
C     
  500 CONTINUE    
  600 CONTINUE
      CALL PGSCH(1.0)
      CALL PGPAGE      
      IF(IFILE.NE.0)CLOSE(30)
      IF(IFILE.NE.0)CLOSE(31)
      RETURN
      END   
