c
c
c-- For Fitting Sum of Damped Sines and Cosines to Measured Time Series Data
c
c      i.e.,  for simultaneous fits of the form:
c
c      ( y_2 )   (  Sum[damped sines]  )   (     0    )   (  Sum[sines]  )
c      ( --- ) = ( ------------------- ) + ( -------- ) + ( ------------ )
c      ( y_1 )   ( Sum[damped cosines] )   ( exp(-Bt) )   ( Sum[cosines] )
c
c      with last two terms (pure exponential & pure sines/cosines) optional
c
c**************************************************************************
c
      subroutine calc ( k, nl, n, mmax, iv, c, alf, t, eta )
c
c--------------------------------------------------------------------------
c
c  input :
c ---------
c       k = data point number
c      nl = no. of non-linear parameters
c       n = no. of linear parameters
c    mmax = maximum no. data points allowed
c      iv = no. of independent variables
c       c = n-vector of linear parameters
c     alf = nl-vector of non-linear parameters
c       t = (m x iv) matrix of independent variable values
c
c  output :
c ----------
c      ph = (n or n+1)-vector of function values : 
c           ph(j) = phi( t(k,1), ... , t(k,iv) ; alf(1), ... , alf(nl) ) ,
c              j = 1, 2, ... , n, (n+1)
c           ( returned through  common block /blkph/ )
c     eta = calculated model prediction :  
c           eta = c(1)*ph(1) + ... + c(n)*ph(n) ( + ph(n+1) )
c
c-------------------------------------------------------------------------
c
      PARAMETER ( maxn1 = 101 )
c                            ------  max. allowed no. of phi functions + 1
c
      dimension c(n), alf(nl), t(mmax,iv)
c
      common /blkph/ ph(maxn1)
c
      common /mtotbl/ mtot
c
      common /modelbl/ ncon, ncwav, npexp, npsin
c
c                   ncon = no. constant (baseline correction) terms in model
c                   ncwav = no. coupled sine/cosine waves in model
c                   npexp = no. pure exponential (0 or 1)
c                   npsin = no. pure (undamped) sinusoids
c
c------------------------------------------------------------------------
c
      twopi = 2.0 * acos( -1.0 )
c
      movr2 = mtot / 2
c
      tk = t(k,1)
c
      if ( nl .eq. -2 ) then
c                          ----- Fit two straight lines to the data
             n = 4
c
             if ( k .le. movr2 ) then
                                     ph(1) = 1.0
                                     ph(2) = tk
                                     ph(3) = 0.0
                                     ph(4) = 0.0
                    else
                        ph(1) = 0.0
                        ph(2) = 0.0
                        ph(3) = 1.0
                        ph(4) = tk
             end if
c
             eta = 0.0
             do j = 1, n
                eta = eta + c(j)*ph(j)
             end do
c
c           ++++++++
             return
c           ++++++++
c
      end if            
c
c
c   For nl > 0 ,  Fit sum of damped sines and cosines to the data
c
      eta = 0.0
c
      if ( ncon .gt. 0 ) then
                         iph = 1    
              if ( k .le. movr2 ) then
                                      ph(iph)   = 0.0
                                      ph(iph+1) = 1.0
                  else
                                      ph(iph)   = 1.0
                                      ph(iph+1) = 0.0
              end if
c
              eta = eta + c(iph)*ph(iph) + c(iph+1)*ph(iph+1)
c
      end if
c
      do kwav = 1, ncwav
c
         iph  = ncon + 2 * kwav - 1
         ialf =        2 * kwav - 1
c
         arg = twopi * alf(ialf+1) * tk
c
         if ( k .le. movr2 ) then
c
                 ph(iph)   = exp( -alf(ialf)*tk ) * sin( arg )
                 ph(iph+1) = exp( -alf(ialf)*tk ) * cos( arg )
c
             else
c
                ph(iph)   =  exp( -alf(ialf)*tk ) * cos( arg )
                ph(iph+1) = -exp( -alf(ialf)*tk ) * sin( arg )
c
         end if
c
         eta = eta + c(iph)*ph(iph) + c(iph+1)*ph(iph+1)
c
      end do
c
      if (npexp .gt. 0) then
                            iph  = ncon + 2 * ncwav + 1
                            ialf =        2 * ncwav + 1
c
              if ( k .le. movr2 ) then
                                      ph(iph)   = 0.0
                                      ph(iph+1) = exp( -alf(ialf)*tk )
                  else
                                      ph(iph)   = exp( -alf(ialf)*tk )
                                      ph(iph+1) = 0.0
              end if
c
              eta = eta + c(iph)*ph(iph) + c(iph+1)*ph(iph+1)
c
      end if
c
      if (npsin .gt. 0) then
c
             do isin = 1, npsin
c
                iph  = ncon + 2*ncwav + 2*npexp + (2*isin - 1)
                jalf =        2*ncwav +   npexp +    isin
c
                arg = twopi * alf(jalf) * tk
                if ( k .le. movr2 ) then
                                    ph(iph)   = sin( arg )
                                    ph(iph+1) = cos( arg )
                        else
                            ph(iph) =  cos( arg )
                            ph(iph+1) = -sin( arg )
                end if
c
                eta = eta + c(iph)*ph(iph) + c(iph+1)*ph(iph+1)
c
             end do
c
      end if
c
c    ++++++++
      return 
c    ++++++++
c
      end 
c
c
c***********************************************************************
c
      subroutine deriv ( k, t, mmax, iv, alf, nl, np1 )
c
c----------------------------------------------------------------------
c
c  input :
c ---------
c
c       k = data point number
c       t = (m x iv) matrix of independent variable values
c    mmax = maximum no. data points allowed
c      iv = no. of independent variables
c     alf = nl-vector of non-linear parameters
c      nl = no. of non-linear parameters
c     np1 = ( n + 1 ) , where  n = no. of linear parameters
c
c  output :
c ----------
c     dph = ( (n or n+1) x nl ) matrix of partial derivatives :
c             dph(i,j) = d( ph(i) ) / d( alf(j) )  evaluated at 
c                      ( t(k,1), ... , t(k,iv) ; alf(1), ... , alf(nl) )
c                        i = 1, 2, ... , n, (n+1) ,   j = 1, 2, ... , nl
c           ( returned through  common block /blkpdh/ )
c
c-------------------------------------------------------------------------
c
      PARAMETER ( maxn1 = 101 , maxkg = 50 )
c
c                 maxn1 = max. allowed no. of phi functions
c                 maxkg = max. allowed no. of non-lin. parameters alf
c
      dimension  alf(nl), t(mmax,iv)
c
      common /blkdph/ dph(maxn1,maxkg)
c
      common /mtotbl/ mtot
c
      common /modelbl/ ncon, ncwav, npexp, npsin
c
c                   ncon = no. constant (baseline correction) terms in model
c                   ncwav = no. coupled sine/cosine waves in model
c                   npexp = no. pure exponentials (0 or 1)
c                   npsin = no. pure (undamped) sinusoids
c
c------------------------------------------------------------------------
c
      twopi = 2.0 * acos( -1.0 )
c
      do 10 j=1,nl
        do 10 i=1,np1
          dph(i,j) = 0.0
   10 continue
c
      n = np1 - 1
c
      tk = t(k,1)
c
      movr2 = mtot / 2
c
c           dph(i,j) should be set to the partial derivative
c                    of ph(i) with respect to alf(j)
c
      do i = 1, n, 2
c
         expfac = tk*exp( -alf(i)*tk )
         sinfac = sin( twopi*alf(i+1)*tk )
         cosfac = cos( twopi*alf(i+1)*tk )
c
         if ( k .le. movr2 ) then
c
                       dph(ncon+i  ,i  ) = -expfac * sinfac
                       dph(ncon+i  ,i+1) =  expfac * cosfac * twopi
                       dph(ncon+i+1,i  ) = -expfac * cosfac
                       dph(ncon+i+1,i+1) = -expfac * sinfac * twopi
c
             else
c
                  dph(ncon+i  ,i  ) = -expfac * cosfac
                  dph(ncon+i  ,i+1) = -expfac * sinfac * twopi
                  dph(ncon+i+1,i  ) =  expfac * sinfac
                  dph(ncon+i+1,i+1) = -expfac * cosfac * twopi
c   
         end if
c
      end do
c
      if ( npexp .gt. 0 ) then
c
                             jalf = 2*ncwav + 1
                             iph  = 2*ncwav + 1 + ncon
c
c-----            if ( k .gt. movr2 ) then
            if ( k .le. movr2 ) then
                         dph(iph+1,jalf) = -tk*exp( -alf(jalf)*tk )
                    else
                         dph(iph,  jalf) = -tk*exp( -alf(jalf)*tk )
            end if
c
      end if
c
      if ( npsin .gt. 0 ) then
c
             do isin = 1, npsin
c
                iph  = ncon + 2*ncwav + 2*npexp + (2*isin - 1)
               jalf =        2*ncwav +   npexp +    isin
c
                arg = twopi * alf(jalf) * tk
                twopit = twopi * tk
c
                if ( k .gt. movr2 ) then
c
                            dph(iph,  jalf) =  twopit * cos( arg )
                            dph(iph+1,jalf) = -twopit * sin( arg )
                    else
                        dph(iph,  jalf) = -twopit * sin( arg )
                        dph(iph+1,jalf) = -twopit * cos( arg )
                end if
c
             end do
c
      end if
c
c    ++++++++
      return
c    ++++++++
c
      end 
c
c
c**************************************************************************
c
      subroutine wt ( m, weight, y )
c
c-------------------------------------------------------------------------
c
c  input :
c ----------
c       m = no. of data points
c       y = m-vector of observations of the dependent variable
c
c  output :
c -----------
c  weight = m-vector of weights for the fit
c
c------------------------------------------------------------------------
c
      dimension  y(m), weight(m)
c
       do 10 i=1,m
         weight(i) = 1.0/y(i)**2
   10  continue
c
      return
c
      end 
c
c
c**************************************************************************
c
      subroutine trans ( t, bt, y, by, iv, m, mmax )
c
c--------------------------------------------------------------------------
c
c  input :
c ---------
c
c       t = ( m x iv ) matrix of independent variable values
c       y = m-vector of observed values of dependent variable
c      iv = no. of independent variables
c       m = no. of data points
c    mmax = max. allowed no. of data points
c
c  output :
c ----------
c       t = transformed matrix of independent variable values
c      bt = original, untransformed, independent variable matrix,
c           i.e.,  bt(output) = t(input)
c       y = transformed vector of dependent variable observations
c      by = original, untransformed, dependent variable vector,
c           i.e.,  by(output) = y(input)
c----------------------------------------------------------------------- 
c
      dimension t(mmax,iv), bt(mmax,iv), y(m), by(m)
c
      common /t0blok/ t0
c
c-----------------------------------------------------------------------
c
      t0 = t(1,1)
c
      do i = 1, m
c 
         bt(i,1) = t(i,1)
c
         t(i,1) = t(i,1) - t0
c
         by(i) = y(i)
c
      end do
c
c
c    ++++++++
      return
c    ++++++++
c
      end 
c
c
c***************************************************************************
c
      subroutine baktrn ( eta, beta, m )
c
c---------------------------------------------------------------------------
c
c  input :
c ---------
c     eta = m-vector of predicted values of dependent variable
c       m = no. of data points
c
c  output :
c ----------
c    beta = m-vector of back-transformed model predictions
c
c--------------------------------------------------------------------------
c
      dimension eta(m), beta(m)
c
      common /t0blok/ t0
c
c-----------------------------------------------------------------------
c
      do i = 1, m
c
         beta(i) = eta(i)
c
      end do
c
c    ++++++++
      return
c    ++++++++
c
      end 
