*-----------------------------------------------------------------------
* damp.f - analyze the damped vibration data of Mohamed Gharib
*-----------------------------------------------------------------------
      PROGRAM DAMP
       IMPLICIT NONE
       INTEGER N, P
       PARAMETER (N = 6000, P=3)
       DOUBLE PRECISION X(N), Y(N), W1(N), W2(N), F, R,
     $    TOL, CARR(N,P), WK(P), COL(P), SOLN(P), AMP, PHI, Z(N), RESULT
       DOUBLE PRECISION X1
       INTEGER IND(P)
       INTEGER IERR, IOU, I, J, K
       INTEGER DIV, NSUB, OFFS
       CHARACTER*32 FNAME, HEADER
       INTEGER IOUNIT
       DOUBLE PRECISION ONE, TWOPI
       PARAMETER (ONE = 1.D0, TWOPI = 6.283185307179586477)

       DATA FNAME / '../data/damp.csv' /
       DATA TOL / 1.D-6 /

*-----------------------------------------------------------------------
*             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) 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 y line ', J
       END DO
       CLOSE (UNIT=IOU, IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble closing file!'

*-----------------------------------------------------------------------
*             Call subroutine.
*-----------------------------------------------------------------------
       CALL DAMEST (X, Y, N, R, F, W1, W2, IERR)
       PRINT *, 'damest returns error# ', IERR
       PRINT*, 'growth rate: ', R, ' frequency: ', F

       DO J = 1, N
         CARR(J,1) = EXP(R * X(J)) * COS(TWOPI * F * X(J))
         CARR(J,2) = EXP(R * X(J)) * SIN(TWOPI * F * X(J))
         CARR(J,3) = ONE
       END DO

       CALL OSOLVD (N, P, CARR, N, Y, TOL, I, SOLN, WK, COL, IND, IERR)
       IF (IERR .NE. 0) PRINT*, 'osolvd returns #', IERR
       IF (I .NE. P) PRINT *, 'rank: ', I

       AMP = SQRT(SOLN(1)**2 + SOLN(2)**2)
       PHI = ATAN2(SOLN(2), SOLN(1))
       PRINT *, 'amplitude: ', AMP, ' phase shift: ', PHI

       DO J = 1, N
         Z(J) = AMP * EXP(R*X(J)) * COS(TWOPI*F*X(J) - PHI) + SOLN(3)
       END DO
       CALL CODD (N, Y, Z, RESULT)
       PRINT *, 'coefficient of determination: ', RESULT

*-----------------------------------------------------------------------
*             Subdivide domain
*-----------------------------------------------------------------------
  20   FORMAT ()
       PRINT 20
       PRINT *, 'in ten pieces:'
       DIV = 10
       NSUB = 600
       DO 30 K = 1, DIV
        PRINT *, 'piece ', K
        OFFS = NSUB * (K-1) + 1
        X1 = X(OFFS)
        DO J = 1, NSUB
          X(OFFS+J-1) = X(OFFS+J-1) - X1
        END DO

        CALL DAMEST (X(OFFS), Y(OFFS), NSUB, R, F, W1, W2, IERR)
         IF (IERR .NE. 0) THEN
           PRINT *, 'damest returns error# ', IERR
           GO TO 30
         END IF
         PRINT*, 'growth rate: ', R, ' frequency: ', F

         DO J = 1, NSUB
           X1 = X(OFFS+J-1)
           CARR(J,1) = EXP(R * X1) * COS(TWOPI * F * X1)
           CARR(J,2) = EXP(R * X1) * SIN(TWOPI * F * X1)
           CARR(J,3) = ONE
         END DO

         CALL OSOLVD (NSUB, P, CARR, N, Y(OFFS), TOL, I, SOLN,
     $               WK, COL, IND, IERR)
         IF (IERR .NE. 0) PRINT *, 'osolvd returns #', IERR
         IF (I .NE. P) PRINT *, 'rank: ', I

         AMP = SQRT(SOLN(1)**2 + SOLN(2)**2)
         PHI = ATAN2(SOLN(2), SOLN(1))
         PRINT *, 'amplitude: ', AMP, ' phase shift: ', PHI

         DO J = 1, NSUB
           X1 = X(OFFS+J-1)
           Z(J) = AMP * EXP(R * X1) * COS(TWOPI * F * X1 - PHI)+ SOLN(3)
         END DO

         CALL CODD (NSUB, Y(OFFS), Z, RESULT)
         PRINT *, 'coefficient of determination: ', RESULT
         PRINT 20
  30   CONTINUE

      END