*-----------------------------------------------------------------------
* star.f - frequency of a variable star
*-----------------------------------------------------------------------
      PROGRAM STAR
       IMPLICIT NONE
       INTEGER N, P
       PARAMETER (N = 116,        ! length of data
     $            P = 3)          ! sinusoid parameters to determine

       DOUBLE PRECISION ZERO, ONE, TWOPI
       PARAMETER (ZERO = 0.D0,
     $            ONE = 1.D0,
     $            TWOPI = 6.2831853071795865D0)

       INTEGER IERR, IOU, J, K
       DOUBLE PRECISION F, X(N), Y(N), S(N), X1, Z(N), RESULT
       INTEGER IND(P)
       DOUBLE PRECISION ARR(N,P), SOLN(P), WORK(P), COLS(P), RELTOL
       CHARACTER*32 FNAME, LINDAT

       INTEGER IOUNIT

       DATA RELTOL / 1.D-6 /
       DATA FNAME / '../data/var_star.csv' /
*-----------------------------------------------------------------------
* Variable star 1SWASPJ141941.31055341.9
* https://www.superwasp.org/vespa/source/1SWASPJ141941.31+055341.9/
*
* This data is an extract of the file v0.9_1SWASPJ141941.31055341.9.fits
* The 116 observations are the longest subsequence without gaps longer
* than one hour.
*
* Time is in hours from the beginning of the .fits file
* Flux is the TAMFLUX2 field, the value given after correction with
* the Tamuz algorithm.
*-----------------------------------------------------------------------

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       IOU = IOUNIT()
       IF (IOU < 1) THEN
         PRINT *, 'trouble getting i/o unit!'
         STOP
       END IF
       OPEN (UNIT=IOU, FILE=FNAME, ACCESS='SEQUENTIAL', STATUS='OLD',
     $           IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble opening ', FNAME

       READ (UNIT=IOU, FMT=*, IOSTAT=IERR) LINDAT  ! burn header
       IF (IERR .NE. 0) PRINT *, 'trouble reading header'

       DO J = 1, N
         READ (UNIT=IOU, FMT=*, IOSTAT=IERR) X(J), Y(J)
         IF (IERR .NE. 0) PRINT *, 'trouble reading data record ', J
       END DO
       CLOSE (UNIT=IOU, IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble closing file!'

       X1 = X(1)
       DO J = 1, N               ! start time at zero
         X(J) = X(J) - X1
       END DO

*-----------------------------------------------------------------------
*             Call subroutine.
*-----------------------------------------------------------------------
       CALL FREEST (X, Y, N, F, S, IERR)
       PRINT *, 'freest returns error# ', IERR
       PRINT *, 'period in hours ', 1.D0 / F

*-----------------------------------------------------------------------
*             Evaluate fit.
*-----------------------------------------------------------------------
       DO J = 1, N
         ARR(J,1) = COS(TWOPI * F * X(J))
         ARR(J,2) = SIN(TWOPI * F * X(J))
         ARR(J,3) = ONE
       END DO
       CALL OSOLVD (N,P,ARR,N,Y,RELTOL,K,SOLN,WORK,COLS,IND,IERR)
       PRINT *, 'osolvd returns #', IERR, ' with rank ', K

       DO J = 1, N
         Z(J) =   SOLN(1) * COS(TWOPI * F * X(J))
     $          + SOLN(2) * SIN(TWOPI * F * X(J))
     $          + SOLN(3)
       END DO
       CALL CODD (N, Y, Z, RESULT)
       PRINT *, 'a: ', SOLN(1), ' b: ', SOLN(2), ' c: ', SOLN(3)
       PRINT *, 'coefficient of determination: ', RESULT
       
      END