program cmplxdemod
c
c     Copyright (C) 1995
c
c     Michael Mann
c
c     perform envelope analysis in a moving loop through a given time
c     series
c
c     v1.0
c
c  uses  stripped version of original  multitaper spectral code of
c  Jeffrey Park (c)
c
c   mmax     - maximum number of spatial dofs desired (ie, max # of time series)
c   ntapmax  - maximum number of spectral dofs desired
c   nscanmax - maximum length of time series (may have to be increased!)
c
      parameter (nscanmax=4000,ntapmax=3)
      real*8 ell,env0,sum,cs,sn
      real aaa(nscanmax),atot1(nscanmax,180),atot2(nscanmax,180)
      integer icount(nscanmax)
c
c
      character*4 chead(158)
      complex zed
c
      common/npiprol/anpi
      common/nnaamm/iwflag
      common/work/b1(8194)
      common/taperz/ta(4100,16),tai(4100,16),dcf(4100,16),amu(4100,4)
      common/staple/tad(8200)
      common/stap2/tas(8192,16)
      common/stap4/cs(8200),sn(8200)
      common/envel1/ell(20)
      common/data/a(50000)
c
      common/pltopt/npad,nst,mpts,nmin,nmax,t(3),fmin,fmax,nf,
     x                                ddt,nwin,inc,tt(3)
c
      equivalence (zed,reim),(iy,dum),(iah,ahead),(iah,chead),
     x          (tad,env0)
c
c
      dimension dum(35),ahead(158),env0(2,4100)
      iwflag=0
      zero = 0.0
      one = 1.0
      pi=3.14159265358979
c
c
c
c     dt    = time interval spacing in units of years
c
c     here we're assuming monthly data
c
      dt = 0.0833333
c
c     set frequency range for DFT, frequency spacing, zero
c     padding, etc.
c     ** NOTE: maximum padded DFT length is npad=8192 **
c
c
      npad = 8192
c
c     f1 = minimum frequency for analysis = 0
c     f2 = maximum frequency  (in cyc/year -- this must always be
c          less than/equal to the nyquist sampling frequency = 0.5/dt)
c     ( presently set for the range 0 to 0.5 cycle year, or, periods greater
c
      f1 = zero
      f2 = 0.5/dt
c
c
c
      ddf=1./(npad*dt)
      eps=.1*ddf
      nf1=(f1+eps)/ddf
      nf2=(f2+eps)/ddf
      nf=nf2-nf1+1
      nf = nf-1
      nfmax = nf
c
      t(1)=f1
      t(2)=ddf
      t(3)=ddf
      nn1=(f1+eps)/ddf
      fmin=n1*ddf
      nn2=(f2+eps)/ddf
      fmax=n2*ddf
      nmin=nn1*2+1
      nmax=nn2*2+2
c
c
c     read in data
c
c     *** this part of the program may have to modified depending
c         on data format ***
c
c     here we assume input file format with dummy x column,
c     and y value taken from 2nd column
c
      open (unit=15,file='evol.in',status='old')
      do i=1,nscanmax
         read (15,*,end=99) adum,aaa(i)
      end do
99    close (unit=15)
c
c     determine length of series
c
      nscan = i-1
c
222   format (4x,f6.2)
      write (6,*) 'read in: '
      write (6,*) ' time series of length: ',nscan
      write (6,*) '(',nscan*dt, ' years long)'
c
c
c     get frequency, determine associated DFT index
c
      write (6,*) 'frequency for complex demod (cycles/year):'
      read (5,*) freq
      iif = 1+(freq+eps)/ddf
c
c     default settings
c
      npi = 2
      anpi=float(npi)
      nwin=3
      imode = 0
      nmode = 1
      iseg = 0
      n1 = 1
      n2 = nscan
      npts = n2-n1+1
      iformat = 0
      iref = 1
      jsmoo = 0
      iaverage = 0
      irecon = 0
      imoderecon = 0
      ioutput = 0
      itype = 0
      nbeg=1
      nend = nscan
      neval = nscan
      nmove = nscan
      ninterv = nmove

444   write (6,*) 'present settings'
      write (6,*) '1) use',nwin, ' x ',npi,' pi tapers'
      write (6,*) '2) data interval (',nbeg,' to ',nend,' )'
      if (itype.eq.0) then
        write (6,*) 
     $    '3) fixed interval (change this for time-dependent analysis)'
      else
         write (6,*) '3) evolutive analysis'
         write (6,*) '    --',nmove,' year moving window'
         write (6,*) '    --',ninterv,' year time step'
      endif
      if (jsmoo.eq.0) then
        write (6,*) '4) inversion type :minimum norm (gen advised)'
      else
        if (jsmoo.eq.1) then
          write (6,*) '4) inversion type: minimum slope'
        else
           write (6,*) '4) inversion type: minimum roughness'
        endif
      endif
      write (6,*) 'item to change (0=continue)'
      read (5,*) icont
      if (icont.eq.1) then
       goto 1111
      else
       if (icont.eq.2) then
         goto 2222
       else
         if (icont.eq.3) then
          goto 3333
         else
            if (icont.eq.4) then
              goto 4444
            else
              goto 9999
            endif
         endif
       endif
      endif
c
c     anpi = desired time/frequency bandwith "p" (ie, "p pi tapers")
c     nwin = number of tapers used (maximum is 2p-1)
c
1111  write (6,*) 'time-frequency bandwidth product "p"'
      read (5,*) npinew
      if (ntapmax.lt.2*npinew-1) then
         write (6,*) 'sorry--exceeds maximum bandwidth'
      else
         write (6,*) 'number of tapers'
         read (5,*) nwinnew
         if (nwinnew.gt.2*npinew-1) then
           write (6,*) 'sorry--cannot exceed 2*p-1'
         else
           npi = npinew
           anpi = float(npi)
           nwin=nwinnew
         endif
      endif
      goto 444

2222  write (6,*) 'beginning time index'
      read (5,*) nbeg
      write (6,*) 'ending time index'
      read (5,*) nend
      nmove = nend-nbeg+1
      neval = nend-nbeg+1
      ninterv = nmove
      goto 444

3333  write (6,*) 'choose'
      write (6,*) '(0) fixed window'
      write (6,*) '(1) evolutive analysis'
      read (5,*) itype
      if (itype.eq.1) then
         write (6,*) '(WARNING--for some reason certain choices'
         write (6,*) 'of window and interval of parameters lead'
         write (6,*) 'to an error (zero output). If this happens,'
         write (6,*) 'make another choice. sorry!)'
         write (6,*) 'window width (in months e.g., 120 for 10 year)'
         read (5,*) nmove
         write (6,*) 
     $ 'interval between evaluations (in months, e.g., 12 for 1 year)'
         read (5,*) ninterv
      endif
      goto 444

c
c     "minimum norm" : minimize the overall size of the envelope of the
c              oscillation. Generally the most faithful constraint, but
c              may artificially minimize amplitude ends of time interval.
c     "minimum slope" :minimize the average slope of the envelope of
c              the oscillation. Particularly useful for reconstructing
c              secular trends in the series which do not have the periodic
c              property of non-secular oscillations.
c
c     "minimum roughness" : minimize the average 2nd derivative of the
c              envelope of the oscillation. Useful if the amplitude of the
c              envelope is believed to be changing abruptly near the ends of
c              the time interval
c
4444  write (6,*) 'inversion for time domain'
      write (6,*) '(0) minimum norm'
      write (6,*) '(1) minimum slope'
      write (6,*) '(2) minimum roughness'
      read (5,*) jsmoo
      goto 444
c
c     we're happy with the settings. Continue...
c
9999  open (unit=4,file='cmplxdemod-amp.out',status='unknown')
      open (unit=7,file='cmplxdemod-phase.out',status='unknown')
c
c     construct tapers appropriate for window width
c
      npts = nmove
      call taper(npts,nwin,ell)
c
      write (6,*) 'nscan ',nscan
      write (6,*) 'neval ',neval
      write (6,*) 'nbeg ',nbeg
      write (6,*) 'nend ',nend
      write (6,*) 'dt    ',dt
      write (6,*) 'npts  ',npts
      write (6,*) 'nmove ',nmove
      write (6,*) 'ninterv ',ninterv
      niter = (neval+0.5-nmove)/ninterv + 1
      write (6,*) 'beginning 1st of ',niter, ' iteration/s...'
      do i=1,nscan
         icount(i)=0
      end do
c
c     begin moving time window loop
c
      do it = 1,niter
        n1 = 1+(it-1)*ninterv
        n2 = n1+npts-1
        ncent = (n1+n2-1)/2
        write (6,*) 'n1  n2 ',n1,n2
c
c       calculate DFTs
c
c       demean gridpoint series locally
c
        sum = 0.
        do i=n1,n2
           sum = sum+aaa(i)
        end do
        sum = sum/float(npts)
c
c
        do iwin=1,nwin
          do i=1,npad
             b1(i)=0.
          end do
          do i=n1,n2
            b1(i-n1+1)=(aaa(i)-sum)*tas(i+1-n1,iwin)
          end do
          call realft(b1,npad/2,1)
          b1(npad-1)=b1(2)
          b1(npad)=0.
          b1(2)=0.
          j=0
          do i=1,npad,2
            j=j+1
            ta(j,iwin)=b1(i)
            tai(j,iwin)=b1(i+1)
          end do
        end do
c
c       do time reconstructions
c
        do i=1,2*nscanmax
           a(i)=0.0
        end do
c
c       decimate tapers
c
1660    inc=(npts-1)/500+1
        j=0
        do i=1,npts,inc
          j=j+1
          do k=1,nwin
           tas(j,k)=tas(i,k)
          end do
        end do
        mpts=j
        ddt=dt*inc
c
        ifreq = iif
        ff0 = freq
        if (ifreq.eq.1) ta(ifreq,iwin)=ta(ifreq,iwin)/2.0
c
c       invert spectral envelope to get time envelope
c
        call envel(ff0,jsmoo)
        call cossin(npts,cs,sn,1.d0,0.d0,dble(ff0),dble(dt))
        finc=float(inc)
c
        do i=1,mpts-1
           ii=(i-1)*inc
           do j=1,inc
            im = ii+j
            a1=(env0(1,i)*(inc+1-j)+env0(1,i+1)*(j-1))/finc
            a2=(env0(2,i)*(inc+1-j)+env0(2,i+1)*(j-1))/finc
            envmag = sqrt(a1**2+a2**2)
            envphase = -180.0*atan2(a2,a1)/pi
            index = n1+im-1
            imonth = mod(index-1,12)+1
c            phase0 = (-float(imonth-1)+0.5)*dt*360.0
c            phase0 = dt*360.0
c            phase0 = 0.0
            envphase = envphase + 360.0
            if (envphase.gt.360.0) envphase=envphase-360.0
c            phase0 = 15.0
            phase0 = 0.0
            envphase = envphase-phase0
c            if (envphase.gt.360.0) envphase=envphase-360.0
c            if (envphase.gt.180.0) envphase=envphase-360.0
c            if (envphase.lt.-360.0) envphase=envphase+360.0
c            if (envphase.lt.-180.0) envphase=envphase+360.0
            if (index.eq.ncent) then
             icount(index)=icount(index)+1
             atot1(index,icount(index))=envmag
             atot2(index,icount(index))=envphase
            endif
           end do
        end do
c
c       linear extrapolate reconstructions one point
c
        ii=(mpts-1)*inc
        e1=2.*env0(1,mpts)-env0(1,mpts-1)
        e2=2.*env0(2,mpts)-env0(2,mpts-1)
        do j=1,inc
            im = ii+j
            a1=(env0(1,mpts)*(inc+1-j)+e1*(j-1))/finc
            a2=(env0(2,mpts)*(inc+1-j)+e2*(j-1))/finc
            envmag = sqrt(a1**2+a2**2)
            envphase = -180.0*atan2(a2,a1)/pi
            index = n1+im-1
            imonth = mod(index-1,12)+1
            phase0 = (float(imonth)+0.5)*dt*360.0
            envphase = envphase+phase0
            if (envphase.gt.360) envphase=envphase-360.0
            if (index.eq.ncent) then
             icount(index)=icount(index)+1
             atot1(index,icount(index))=envmag
             atot2(index,icount(index))=envphase
            endif
        end do
c
c       provide progress update
c
        open (unit=12,file='svdevals-standard.status',
     $         status='unknown')
        if ((niter.ge.10).and.(mod(it-1,niter/10).eq.0)) then
         write (12,*) it,' iterations of ',niter,' done...'
         write (6,*) it,' iterations of ',niter,' done...'
        endif
        close (unit=12)
c
c       terminate time loop
c
      end do
c
c     now determine ensemble averaged amplitude and phase
c     information
c
      do i=1,nscan
       time = float(i-1)*dt
       envmag = 0.0
       envphase = 0.0
       do j=1,icount(i)
          envmag = envmag+atot1(i,j)
          envphase = envphase+atot2(i,j)
       end do
       if (icount(i).gt.0) then
       envmag = envmag/float(icount(i))
       envphase = envphase/float(icount(i))
       if (envphase.gt.180.0) envphase=envphase-360.0
       write (4,888) time,envmag
       write (7,888) time,envphase
       endif
      end do
      write (6,*) 'done.'
c
c     done -- close files
c
      close (unit=4)
      close (unit=7)
888   format (2f12.2,i4)
c888   format (f6.2,f12.4)
      stop
      end