
      SUBROUTINE XCORR(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 TEMPS,POS,IDEPL,IC,PL
      DOUBLE PRECISION TQ(NQ),HELJD(NORUN),HJD,TZERO,TWOPI,EPOCH,EPER 
      DOUBLE PRECISION O(1000),C(1000),ttt,parms(7)
      REAL YMIN,YMAX,CR,NORM,N
      REAL VAL(10000),CORR(10000),BI(10000) 
      DIMENSION Q(NQ),LIM(2,NORUN),DUM(4)
      DIMENSION FREQ(1),AMP(1),PHASE(1)   
      CHARACTER TITLE*80   
      CHARACTER TITHJD*14,YTITLE*20,XOPT*20,YOPT*20,INTITL*(*)
      LOGICAL SPECIAL
C
C     READ IN THE PARAMETERs: IFILE=1 WRITE THE RESULT TO A FILE, 
C                             IFILE=0 DOES NOT. 
C                             IMODE=0 COMPUTES A CIRCULAR CORRELATION 
C                                     FUNCTION (IN THIS CASE ISHIFT MUST 
C                                     BE ZERO). 
C                             IMODE=1 COMPUTES A SHIFTED CORRELATION 
C                                     FUNCTION, AND ISHIFT GIVES THE NUMBER 
C                                     OF LAGS THE 2 TIME SERIES MUST 
C                                     BE SHIFTED. 
C                             PL=0 DOES NOT PLOT THE RESULTS. 
C                             PL=1 DOES, WITH YMIN AND YMAX AS THE Y-AXIS 
C                                     COORDINATES. IF YMIN=YMAX=0, THEN
C                                     AUTOSCALING IS BEING USED.
C
      READ *,IFILE,IMODE,ISHIFT,PL,YMIN,YMAX      
C     
C     COMPUTE THE SYNTHETIC LIGHT CURVE WITH THE FITTED 
C     FREQUENCIES IN FREQ, AMP AND PHASE  
C     
      IF (IFILE.EQ.1) THEN 
        OPEN (30,FILE='xcorr.dat',STATUS='UNKNOWN')
      END IF
      IF ((IFILE.NE.1).AND.(IFILE.NE.0)) THEN
        PRINT *,'NON-VALID INPUT VALUE FOR PARAMETER IFILE'
        PRINT *,'PROGRAM THUS CONSIDERS IFILE=0'
        IFILE=0
      END IF
      IF (IMODE.EQ.0) PRINT *,'MODE=0'
      IF ((IMODE.NE.0).AND.(IMODE.NE.1)) THEN
        PRINT *,'IMODE HAS UNPOSSIBLE VALUE. PROGRAM CONSIDERS IMODE=0'
        IMODE=0
      END IF
C
      FMAX=0.0  
      DO J=1,NFREQ
        FMAX=MAX(FREQ(J),FMAX)  
      ENDDO
      TWOPI=6.2831853
      DT=TQ(2)-TQ(1)
      DO TEMPS=1,NQ
        IF ((TEMPS.GT.1).AND.(TQ(TEMPS)-TQ(TEMPS-1).NE.DT)) THEN
          PRINT *,'TIME SPACING IS NOT REGULAR'
          PRINT *,'CORRELATION FUNCTION CANNOT BE COMPUTED'
          STOP
        END IF
        VAL(TEMPS)=AFIT+(BFIT*TQ(TEMPS))+(CFIT*TQ(TEMPS)*TQ(TEMPS))
        DO J=1,NFREQ
          ARG=DMOD(TWOPI*FREQ(J)*TQ(TEMPS)+PHASE(J),TWOPI) 
          VAL(TEMPS)=VAL(TEMPS)+AMP(J)*COS(ARG)   
        ENDDO
      ENDDO
      
      IF (IMODE.EQ.1) GOTO 100
C
C     COMPUTE THE CORRELATION FUNCTION CORR FOR IMODE=0
C
      DO IC=1,NQ+1
        CR=0
        DO L=1,NQ
          POS=L+IC-1
          IF (POS.GT.NQ) POS=POS-NQ
          IF (POS.GT.NQ) POS=POS-NQ
          CR=CR+(Q(L)*VAL(POS))
        ENDDO
        NORM=NORM+(CR*CR)
        CORR(IC)=CR/NQ
        BI(IC)=IC  
      ENDDO
      N=SQRT(NORM)
C
C     WRITE OUT THE RESULT INTO FILE 'xcorr.dat'
C
      IF (IFILE.EQ.1) THEN
        DO IWR=1,NQ+1
          CORR(IWR)=CORR(IWR)/N
          WRITE (30,'(F7.0,3X,E12.5)') BI(IWR),CORR(IWR)
        ENDDO
      END IF
      GOTO 200
C
C     COMPUTE THE CORRELATION FUNCTION CORR FOR IMODE=1
C
100   IF (ISHIFT.GT.NQ) ISHIFT=NQ-1
      DO IC=-ISHIFT,ISHIFT
        CR=0
        DO L=1,NQ
          IF ((L+IC.GT.0).AND.(L+IC.LT.NQ-ISHIFT+1)) THEN
            CR=CR+(Q(L)*VAL(L+IC))
          END IF
        ENDDO
        NORM=NORM+(CR*CR)
        CORR(IC+ISHIFT+1)=CR/(NQ-ABS(IC))
        BI(IC+ISHIFT+1)=IC  
      ENDDO
      N=SQRT(NORM)
C
C     WRITE OUT THE RESULT INTO FILE 'xcorr.dat'
C
      IF (IFILE.EQ.1) THEN
        DO IWR=-ISHIFT,ISHIFT
          CORR(IWR+ISHIFT+1)=CORR(IWR+ISHIFT+1)/N
      	  WRITE(30,'(F7.0,3X,E12.5)')BI(IWR+ISHIFT+1),CORR(IWR+ISHIFT+1)
        ENDDO
      END IF
C
C     PLOT THE RESULT IF REQUIRED (IF PL=1)
C
200   IF (PL.EQ.0) GOTO 500 
C
C     INITIALISE THE PLOTTING IF NECESSARY
C
      IF(IPGPL.EQ.0)THEN   
        IPGPL=1
        CALL PGBEG(0,PGPLDEV,1,1)
      END IF
C     SET THE CHARACTER HEIGHT
      CALL PGSCH(1.2)
C     SET THE LINE WIDTH (2=DARKER THAN DEFAULT)
      CALL PGSLW(2) 
C
C     ARRANGE FOR VIEWPORT TO BE CENTRED IN THE PLOTTING DEVICE AREA
C     AND TO OCCUPY:
C        IF /VPS DEVICE, 50 % OF WIDTH IN X, 30 % OF WIDTH IN Y.
C        IF /PS DEVICE , 50 % OF WIDTH IN X, 60 % OF WIDTH IN Y.
C
      IF(INDEX(PGPLDEV,'/VPS') .NE.0.OR.
     *   INDEX(PGPLDEV,'/XWIN').NE.0)THEN
        CALL PGSVP(0.20,0.85,0.2,0.6)
      END IF
      IF(INDEX(PGPLDEV,'/PS').NE.0)THEN
        CALL PGSVP(0.3,0.8,0.25,0.85)
      END IF
C
C     X SCALING
C
      IF (IMODE.EQ.0) THEN
        XLO=1
        XHI=NQ+1   
      ELSE
        XLO=1        
        XHI=2*ISHIFT+1
      END IF          
      XMIN=1.0E30
      XMAX=-1.0E30
      DO J=XLO,XHI
        IF(BI(J).LT.XMIN)XMIN=BI(J)
        IF(BI(J).GT.XMAX)XMAX=BI(J)
      END DO
C
C     Y SCALING
C
      IF(YMIN.EQ.0.0.AND.YMAX.EQ.0.0)THEN   
        YMIN=1.0E30
        YMAX=-1.0E30
        DO J=XLO,XHI
          IF(CORR(J).GT.YMAX)YMAX=CORR(J)
          IF(CORR(J).LT.YMIN)YMIN=CORR(J)
        END DO
        IF (YMAX.GT.0) THEN  
          YMAX=YMAX*1.02
        ELSE
          YMAX=YMAX*0.95
        END IF
        IF (YMIN.LT.0) THEN
          YMIN=YMIN*1.02
        ELSE
          YMIN=YMIN*0.95
        END IF  
      END IF
C
C     ESTABLISH MAPPING BETWEEN VIEWPORT AND WORLD COORDINATES
C
      CALL PGSWIN(XMIN,XMAX,YMIN,YMAX)
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 PGMTXT( 'T', 3.0, -0.2, 0.0, '.')     
      CALL PGMTXT( 'B', 4.0, 1.1, 0.0, '.')    
      CALL PGSCH(1.2)
C
C     ANNOTATE VIEWPORT WITH FRAME, TICK MARKS AND NUMERIC LABELS
C
      CALL PGBOX('BCNTS1',0.0,0,'BCNTS1',0.0,0)
C
C     PLOT SPECTRUM
C
      KD=XHI-XLO+1      
      CALL PGLINE(KD,BI,CORR)
      IF (IMODE.EQ.1) THEN
        CALL PGMOVE(0,YMIN)
	CALL PGDRAW(0,YMAX)
      END IF 
C     USE PGLAB TO DRAW X- AND Y-AXIS LABELS
      IF(LEN(INTITL).GT.0)
     *CALL PGLAB('Shift','Cross correlation function',INTITL)
      IF(LEN(INTITL).EQ.0)
     *CALL PGLAB('Shift','Cross correlation function',' ')

500   IF (IFILE.EQ.1) CLOSE(30)
      CALL PGPAGE
      RETURN
      END   

