* spiral.f - estimate parameters of a log spiral
      PROGRAM SPIRAL
       IMPLICIT NONE
       INTEGER N
       DOUBLE PRECISION ZERO, ONE, TWOPI
       PARAMETER (N = 992,
     $            ZERO = 0.0D0,
     $            ONE = 1.0D0,
     $            TWOPI = 6.2831853071795865D0)

       INTEGER J
       DOUBLE PRECISION AMP, PHASE, FREQ, RATE, Y0R, Y0I,
     $                   E, YAVGR, YAVGI, SOS, RESID, RESULT
       DOUBLE PRECISION X(N), YRI(2*N), W(2*N,4), ZRI(2*N)
       LOGICAL FEXIST
       INTEGER IU, IERR
       CHARACTER*32 IFIL
       CHARACTER*32 HEADER

       INTEGER IOUNIT

       DATA IFIL / '../data/spiral.csv' /

*             Read input.
       INQUIRE (FILE=IFIL, EXIST=FEXIST, IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble inquiring of file'
       IF (.NOT. FEXIST) PRINT *, 'input file not found'

       IU = IOUNIT()
       IF (IU .EQ. -1) PRINT *, 'trouble getting i/o unit!'
       OPEN (UNIT=IU, FILE=IFIL, STATUS='OLD', ACCESS='SEQUENTIAL',
     $        IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble opening ', IFIL

       READ (UNIT=IU, FMT=*, IOSTAT=IERR) HEADER
       IF (IERR .NE. 0) PRINT *, 'trouble reading header'

       DO J = 1, N
         READ (UNIT=IU, FMT=*, IOSTAT=IERR) X(J), YRI(J), YRI(N+J)
         IF (IERR .NE. 0) PRINT *, 'trouble reading line ', J
       END DO

       CLOSE (UNIT=IU, IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble closing file'

*             Call subroutine.
       PRINT *, 'calling spiest...'
       CALL SPIEST (X,YRI,N,RATE,FREQ,AMP,PHASE,Y0R,Y0I,W,IERR)
       PRINT *, 'spiest returns #', IERR
       PRINT *, 'decay rate: ', RATE, ' (-0.3)'
       PRINT *, 'frequency: ', FREQ, ' (0.7)'
       PRINT *, 'amplitude: ', AMP, ' (5.0)'
       PRINT *, 'phase shift: ', PHASE, ' (2.3)'
       PRINT *, 'center:  (', Y0R, ', ', Y0I, ')  (-3.8, 4.2)'

       DO J = 1, N
         E = AMP * EXP(RATE * X(J))
         ZRI(J) = E * COS(TWOPI * FREQ * X(J) - PHASE) + Y0R
         ZRI(J+N) = E * SIN(TWOPI * FREQ * X(J) - PHASE) + Y0I
!         PRINT *, X(J), YRI(J), YRI(N+J), ZRI(J), ZRI(N+J)
       END DO

       YAVGR = ZERO
       YAVGI = ZERO
       DO J = 1, N
         YAVGR = YAVGR + YRI(J)
         YAVGI = YAVGI + YRI(N+J)
       END DO
       YAVGR = YAVGR / DBLE(N)
       YAVGI = YAVGI / DBLE(N)

       SOS = ZERO
       RESID = ZERO
       DO J = 1, N
         SOS = SOS + (YRI(J) - YAVGR)**2 + (YRI(N+J) - YAVGI)**2
         RESID = RESID + (ZRI(J) - YRI(J))**2 + (ZRI(N+J) - YRI(N+J))**2
       END DO
       RESULT = ONE - RESID / SOS
       PRINT *, 'coefficient of determination: ', RESULT
      END