      program redi_wave
c
c     Program takes input time series and performs various operations
c	Command line driven
c
c     redi_wave -w -g -r -v [-l lat lon -t time] file1 .... filen
c
c     -w computes Wood Anderson amplitude
c     -g computes ground acceleration, velocity, and displacement
c     -e computes energy magnitude - requires 3-C data
c     -p computes characteristic function and picks - requires 3-C data
c     -r computes ramp function
c
c     -l inputs event location (lat, lon in decimal degrees)
c     -t inputs event origin time
c     -d inputs event depth (in km)
c
c     -f inputs filter information (npole, fl, fh in hz)
c
c     -coord inputs file for instrument coordinates
c     -resp inputs file for instrument response
c
c     -v sets the verbose output option
c     -s sets the save flag
c
      include '/usr/local/include/qlib2.inc'
c
c     maxp is the dimension of the data array
c     maxp2 is the dimension of the array for the fft
c     maxps is the parameter for the maximum request
c     -- this is the data array to be doubled to the next power of 2
c     maxfl is the parameter for the maximum number of files
c
      parameter (maxp  = 65536)
      parameter (maxp2 = 32768)
      parameter (maxps = 32767)
      parameter (maxfl = 75)
      parameter (mphase = 60)
c
c structure for sdr header
c
      include 'station_hdr.inc'
      include 'pick_hdr.inc'
c
      integer*4 idata(maxp)
      real*4 d1(maxp)
      complex*8 dw(maxp2)
c
      real*4 dt
      real*4 slat, slon, sele
      real*4 elat, elon, edep
      real*4 fl0, fh0
      real*4 dtr, dtc
      real*4 tip(mphase),pip(mphase),dtdh(mphase),dddp(mphase)
	real*4 pi, deg2km, deg2rad
      real*8 td1, td2
      integer*4 it1s, it1t, it2s, it2t
      integer*4 npts, iflag, nfile, ierr, icomp
      integer*4 ifile(maxfl), np0, idc
      integer*4 ic(3)
      character*1 cid
      character*3 chan, chan_iem
      character*4 stat, stat_iem
      character*8 plist, slist, phip(mphase)
      character*256 fileo(maxfl), filen, argv
      character*256 file_coord, file_resp, file_sens
      character*256 file_save, file_pick
      logical lwood, lground, lresp, lexist, lloc, ltime, lverbose
      logical lenergy, ldep, lsave, lpick, lramp, l3, lsignal
      logical parse_date, lclip
c
      external myhandler, parse_date
c
      record /INT_TIME/ org_time, time_iem, t1_time, t2_time,
     .  p_time
      record /EXT_TIME/ ext_time
      record /PICK/ pick(100),iasp(100)
c
c     common blocks
c
      common /filter/ np0, fl0, fh0
c
c     data arrays
c
      data io /6/
      data io_out /7/
      data io_save /8/
      data io_in /2/
      data nfile /0/
      data lwood,lground,lresp,lexist,lloc,ltime,lverbose /7*.false./
      data lenergy, lpick, lramp, l3, ldep, lclip /6*.false./
      data icomp /0/
      data file_coord /'/home/redi/stat/bdsn.stat.all'/
      data file_resp /'/home/redi/stat/bdsn.resp'/
c     data file_resp /'/usr/contrib/data/bdsn/instr.resp'/
      data file_sens /'/usr/contrib/data/bdsn/instr.sens'/
      data file_pick /'pickfile.pf'/
      data np0,fl0,fh0 /0, 0.0, 0.0/
      data plist,slist /'P','S'/
      data ipick,itime /0,0/
      data s2n /0.0/
c
c     turn on ieee error handling
c
c     ie = ieee_handler ('set', 'all', myhandler)
c
c     conversion from degrees to km
c
      rnk = 6371.0
      pi = 4.0*atan(1.0)
      deg2rad = pi/180.
      deg2km = rnk * deg2rad
c
c     get command line arguments
c
      inum = iargc()
      if (inum .eq. 0 .or. inum .eq. 1) then
        write(io,*)' redi_wave {-w -g -r} [-l lat lon -t time] '
        write(io,*)'           {-v -s} file1 .... filen'
        write(io,*)' -w          computes WA amplitude'
        write(io,*)' -g          computes ground motion'
        write(io,*)' -r          computes ramp function'
        write(io,*)' -e          computes energy magnitude'
        write(io,*)'               {requires -l and -d and 3-comp data}'
        write(io,*)' -p          computes characteristic function'
        write(io,*)'               {requires 3-comp data}'
        write(io,*)' '
        write(io,*)' -n s/n      input signal/noise ratio for testing'
        write(io,*)' -t time     input event origin time'
        write(io,*)' -l lat lon  input event location'
        write(io,*)' -d depth    input event depth'
        write(io,*)' '
        write(io,*)' -f np fl fh input filter parameters'
        write(io,*)' -c dtr dtc  input chr/detector parameters'
        write(io,*)' '
        write(io,*)' -coord      input station coordinate file'
        write(io,*)' -resp       input station response file'
        write(io,*)' '
        write(io,*)' -pick       input pickfile'
        write(io,*)' '
        write(io,*)' -v          verbose output option'
        write(io,*)' -s          save output data'
        write(io,*)' file1...    input SDR format file(s)'
        write(io,*)' '
        stop
      else
        i = 1
c
c       This large do while loops over the command line arguments
c
        do while (i .le. inum)
c
cccc      Check for options
c
          call getarg(i,argv)
          if (argv .eq. '-w') then
            lwood = .true.
          elseif (argv .eq. '-g') then
            lground = .true.
          elseif (argv .eq. '-e') then
            lenergy = .true.
            l3 = .true.
          elseif (argv .eq. '-p') then
            lpick = .true.
            l3 = .true.
          elseif (argv .eq. '-r') then
            lramp = .true.
            l3 = .true.
          elseif (argv .eq. '-v') then
            lverbose = .true.
          elseif (argv .eq. '-s') then
            lsave = .true.
          elseif (argv .eq. '-n') then
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, s2n, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode S/N ratio'
              call exit(1)
            end if
            lsignal = .true.
          elseif (argv .eq. '-l') then
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, elat, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode latitude'
              call exit(1)
            end if
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, elon, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode longitude'
              call exit(1)
            end if
		if (lverbose) then
		  write(io,*)' Coordinates: ', elat, elon
		end if
            lloc = .true.
          elseif (argv .eq. '-t') then
            i = i + 1
            call getarg(i,argv)
            ltime = f_parse_date(org_time, argv)
            if (.not. ltime) then
              write(io,'(a)') '  Redi_wave - Unable to decode time'
              call exit(1)
            end if
		if (lverbose) then
		  write(io,*) ' Time: ', org_time.year, org_time.second, 
     &          org_time.usec
		end if
          elseif (argv .eq. '-d') then
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, edep, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode depth'
              call exit(1)
            end if
		if (lverbose) then
		  write(io,*)' Depth: ', edep
		end if
            ldep = .true.
          elseif (argv .eq. '-f') then
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, rnum, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode np0'
              call exit(1)
            end if
            np0 = rnum
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, fl0, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode fl'
              call exit(1)
            end if
            i = i + 1 
            call getarg(i,argv)  
            call cipnum(argv, fh0, iflag) 
            if (iflag .ne. 0) then 
              write(io,'(a)') '  Redi_wave - Unable to decode fh'
              call exit(1) 
            end if 
		if (lverbose) then
		  write(io,*)' Filter: ', np0, fl0, fh0
		end if
          elseif (argv .eq. '-c') then
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, rnum, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode chr'
              call exit(1)
            end if
            idc = rnum
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, dtr, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode dtr'
              call exit(1)
            end if
            i = i + 1
            call getarg(i,argv)
            call cipnum(argv, dtc, iflag)
            if (iflag .ne. 0) then
              write(io,'(a)') '  Redi_wave - Unable to decode dtc'
              call exit(1)
            end if
            if (lverbose) then
              write(io,*)' Characteristic function: ', idc
              write(io,*)' Detection constants: ', dtr, dtc
            end if
          elseif (argv .eq. '-coord') then
            i = i + 1
            call getarg(i,argv)
            file_coord = argv
            if (lverbose) then
              call klen(file_coord,kc)
              write(io,*)' Coordinate file: ', file_coord(1:kc)
            end if
          elseif (argv .eq. '-pick') then
            i = i + 1
            call getarg(i,argv)
            file_pick = argv
            if (lverbose) then
              call klen(file_pick,kc)
              write(io,*)' Pick file: ', file_pick(1:kc)
            end if
          elseif (argv .eq. '-resp') then
            i = i + 1
            call getarg(i,argv)
            file_resp = argv
            if (lverbose) then
              call klen(file_resp,kc)
              write(io,*)' Response file: ', file_resp(1:kc)
            end if
          else
c
c           get a list of all the file names
c
            nfile = nfile + 1
            if (nfile .gt. maxfl) then
              nfile = nfile - 1
              write(io,*)'  Redi_wave - too many data files'
            else
              fileo(nfile) = argv
            end if
          end if
          i = i + 1
        end do
      end if
c
      if (nfile .eq. 0) then
        write(io,*) '  Redi_wave - no input data files '
        call exit(1)
      end if
c
c     sort the files names
c
      call sort_c(nfile,fileo,ifile)
c
c     tell us what to expect
c
      if (lverbose) then
        write(io,*) ' '
        if (lwood) then
          write(io,'(a)') '  Computing Wood-Anderson synthetics'
        end if
        if (lground) then
          write(io,'(a)') '  Computing ground motion values'
        end if
        if (l3) then
          if (lenergy) then
            write(io,'(a)') '  Computing energy magnitude'
          end if
          if (lpick) then
            write(io,'(a)') '  Computing characteristic function'
          end if
          if (lramp) then
            write(io,'(a)') '  Computing ramp growth rate'
          end if
        end if
      end if
c
c     Now begin processing data
c
      if (lverbose) then
          write(io,*) ' '
          write(io,*) ' Looping over files'
      end if
      do jj = 1, nfile
        lclip = .false.
        filen = fileo(ifile(jj))
        call klen(filen,kf)
        if (lverbose) then
          write(io,*) ' ',filen(1:kf)
        end if
c
c       read data in  - the data is originally integer
c
        call sdr_read(filen,maxps,idata,head,ierr)
        if (ierr .ge. 999) then
          write(io,*) '  Redi_wave - error opening file ',ierr
          go to 999
        end if
        stat  = head.station_id
        chan  = head.channel_id
        npts  = head.num_pts
        dt  = 1.0/head.sample_rate
c
c       check for clipping - stolen from Bob
c
        call check_clip(chan,npts,idata,ierr)
        if (ierr .ne. 0) then
          write(io,*) '  Redi_wave - data may be clipped ',ierr
          if (.not. lramp) then
            go to 999
          else
            write(io,*) '  Using data anyway'
            lclip = .true.
          end if
        end if
c
c       get instrument information
c
        inquire(file=file_coord,exist=lexist)
	  if (.not. lexist) then
          write(io,*) '  Redi_wave - coordinate file does not exist'
          go to 999
        end if
c
        call fetch_coord(io_in,file_coord,stat,slat,slon,sele,ierr)
	  head.slat = slat
	  head.slon = slon
	  head.sele = sele
        if (ierr .ne. 0) then
          write(io,*) '  Redi_wave - error with coordinates ',ierr
          go to 999
        end if
c
        if (lramp) then
          inquire(file=file_sens,exist=lexist)
  	    if (.not. lexist) then
            write(io,*) '  Redi_wave - sensitivity file does not exist'
            go to 999
          end if
          call fetch_sens(io_in,file_sens,head,ierr) 
          if (ierr .ne. 0) then
            write(io,*) '  Redi_wave - error with sensitivity ',ierr
            go to 999
          end if
        end if
c
        inquire(file=file_resp,exist=lexist)
	  if (.not. lexist) then
          write(io,*) '  Redi_wave - response file does not exist'
          go to 999
        end if
        call fetch_response(io_in,file_resp,head,ierr)
        if (ierr .ne. 0) then
          write(io,*) '   Redi_wave - error with response ',ierr
          go to 999
        end if
c
c       read the pickfile, if it exists
c
        inquire(file=file_pick,exist=lexist)
        if (lexist) then 
          call fetch_pick(io_in,file_pick,head,pick,ipick,ierr) 
          if (ierr .ne. 0) then
            write(io,*) '   Redi_wave - error reading pickfile', ierr
          end if
c
          if (lverbose) then
            do ip = 1, ipick
              call f_int_to_ext(pick(ip).ptime,ext_time)
              write(io,*) pick(ip).pid,' ',
     .        ext_time.year, ext_time.doy, ext_time.hour,
     .        ext_time.minute, ext_time.second, ext_time.usec
            end do
          end if
c
        end if
c
c       calculate the epicentral distance, azimuth, and focal distance
c
        if (lloc) then
          call azimth(elat,elon,slat,slon,delta,azim,bazim)
          dist = delta * deg2km
          if (ldep) then
            rdist = sqrt(dist*dist + edep*edep)
          end if
        else
          delta = 0.0
          dist = 0.0
          rdist = 0.0
        end if
c 
c       calculate some travel times 
c 
        if (lenergy .or. lramp) then
        if (lloc .and. ltime .and. ldep) then 
          call ttimes(edep,delta,plist,nip,tip,pip,dtdh,dddp,phip) 
          t1 = tip(1) 
          it1s = int(t1)  
          it1t = (t1 - it1s)*usecs_per_sec  
          call f_add_time(org_time, it1s, it1t, t1_time)  
          iasp(1).pid = phip(1)
          iasp(1).ptime = t1_time
c
          call ttimes(edep,delta,slist,nip,tip,pip,dtdh,dddp,phip)
          t1 = tip(1)
          it1s = int(t1)
          it1t = (t1 - it1s)*usecs_per_sec
          call f_add_time(org_time, it1s, it1t, t1_time)
          iasp(2).pid = phip(1)
          iasp(2).ptime = t1_time
c
          itime = 2
        end if 
        end if 
c
c       check for signal to noise ratio
c
c       if (lsignal) then
c         call check_s2n(head,idata,iasp,ierr)
c         if (ierr .ne. 0) then
c           write(io,*) '  Redi_wave - data may be clipped ',ierr
c           go to 999
c         end if
c       end if
c
c       extract desired window length in some cases
c
        if (lenergy .or. lramp) then
          t1 = -30.
          t2 = 60.
          it1s = int(t1)
          it1t = (t1 - it1s)*usecs_per_sec
          it2s = int(t2)
          it2t = (t2 - it2s)*usecs_per_sec
c
          if (lverbose) then
            call f_int_to_ext(head.begtime,ext_time)
            write(io,*) ' Before window extraction: ',
     .        stat,' ',chan,' ',head.num_pts,
     .        ext_time.year, ext_time.doy, ext_time.hour,
     .        ext_time.minute, ext_time.second, ext_time.usec
          end if
c
          if (ipick .ne. 0) then
            call f_add_time(pick(1).ptime, it1s, it1t, t1_time)
            call f_add_time(pick(1).ptime, it2s, it2t, t2_time)
            call extra_window(head,idata,t1_time,t2_time,ierr)
          elseif (itime .ne. 0) then
            call f_add_time(iasp(1).ptime, it1s, it1t, t1_time)
            call f_add_time(iasp(1).ptime, it2s, it2t, t2_time)
            call extra_window(head,idata,t1_time,t2_time,ierr)
          else
            write(io,*) '  Redi_wave - unable to extract window'
          end if
50        continue
          if (ierr .ne. 0) then
            write(io,*) '  Redi_wave - insufficient data ',ierr
            go to 999
          else
            npts = head.num_pts
          end if
          if (lverbose) then 
            call f_int_to_ext(head.begtime,ext_time) 
            write(io,*) ' After window extraction:  ',
     .        stat,' ',chan,' ',head.num_pts,
     .      ext_time.year, ext_time.doy, ext_time.hour, 
     .      ext_time.minute, ext_time.second, ext_time.usec 
          end if 
        end if
c
c       check length of timeseries for arrivals of interest
c
        if (lwood .or. lground) then
          if (lloc .and. ltime) then
            call check_window(head.begtime,head.endtime,org_time,
     .        dist,ierr)
            if (ierr .ne. 0) then
              write(io,*) '  Redi_wave - insufficient window ',ierr
              go to 999
            end if
          end if
        end if
c
c       check length of timeseries for FFT
c
        n2 = ipow(npts*2)
        if (n2 .gt. maxp) then
          write(io,*) '  Redi_wave - FFT exceeded '
          go to 999
        end if  
c
c       summarize what we know about this data
c
c       if (lverbose) then
c         write(io,*) stat, ' ', chan, ' ', npts, ' ',dt
c         write(io,*) head.slat,' ',head.slon,' ',head.sele
c         write(io,*) head.azi,' ',head.dip,' ',head.iunit
c         write(io,*) head.ds,' ',head.gain
c         write(io,*) head.npoles,' ',head.nzeros
c       endif
c
c ********************************************************************
c We are through with preparations - now let's do something
c ********************************************************************
c
c       Calculate Wood Anderson response
c
        if (lwood) then
          do ii = 1, npts
            d1(ii) = idata(ii)
          end do
          call est_wood(head,d1,dw,io_out,ierr)
          if (ierr .ne. 0) then
            write(io,*) '  Redi_wave - error in WA est ',ierr
            go to 999
          end if
          if (lsave) then
            file_save = filen(1:kf)//'.wood'
            call save_to_ascii(head,d1,io_save,file_save,ierr)
          end if
        end if
c
c       Calculate ground displacement/velocity/acceleration
c
        if (lground) then
          do ii = 1, npts
            d1(ii) = idata(ii)
          end do
c
c         changed from dist to rdist 9/23/98
c
          call est_ground(head,d1,dw,rdist,io_out,ierr)
          if (ierr .ne. 0) then
            write(io,*) '  Redi_wave - error in GM est ',ierr
            go to 999
          end if
          if (lsave) then
            file_save = filen(1:kf)//'.ground'
            call save_to_ascii(head,d1,io_save,file_save,ierr)
          end if
        end if
c
c       Calculate response spectra
c
        if (lresp) then
          write(io,*) '  Redi_wave - option not implemented'
        end if
ccccc
c
c       These options require 3 component data
c       -- do the overhead now in common
c
       if (l3) then
c
          if (icomp .ne. 0) then
            if (stat .eq. stat_iem) then
              if (chan(1:2) .eq. chan_iem(1:2)) then
                td1 = f_tdiff(head.begtime, time_iem)
                if (td1 .eq. 0.d0) then
                  icomp = icomp + 1
                  go to 1000
                end if
              end if
            end if
            write(io,*) '  Redi_wave - unexpected channel ',td1
          end if
c
c         try again
c
          icomp = 1
          stat_iem = stat
          chan_iem = chan
          time_iem = head.begtime
          do kk = 1, 3
            ic(kk) = 0
          end do
c
1000      continue
          cid = chan(3:3)
          if (cid .eq. 'Z' .or. cid .eq. '1') then
            ic(1) = ic(1) + 1
          elseif (cid .eq. 'N' .or. cid .eq. '2') then
            ic(2) = ic(2) + 1
          elseif (cid .eq. 'E' .or. cid .eq. '3') then
            ic(3) = ic(3) + 1
          end if
c
c         Calculate energy magnitude
c
          if (lenergy) then
            do ii = 1, npts
              d1(ii) = idata(ii)
            end do
            call est_energy(head,d1,dw,dist,icomp,ic,io_out,ierr)
            if (ierr .ne. 0) then
              write(io,*) '  Redi_wave - error in EN est ',ierr
              go to 399
            end if
            if (lsave) then
               file_save = filen(1:kf)//'.energy'
              call save_to_ascii(head,d1,io_save,file_save,ierr)
            end if
399         continue
          end if
c
c         Calculate characteristic function and pick arrivals
c
          if (lpick) then
c            do ii = 1, npts
c              d1(ii) = idata(ii)
c            end do
c            if (idc .eq. 0) then
c              idc = 1
c              dtr = 9.
c              dtc = 15.
c            end if
c            call est_cfr(head,idc,dtr,dtc,d1,dw,icomp,ic,
c     .          file_pick,io_out,ierr)
c            if (ierr .ne. 0) then
c              write(io,*) '  Redi_wave - error in CF est ',ierr
c              go to 499
c            end if
c            if (lsave) then
c              file_save = filen(1:kf)//'.pick'
c              call save_to_ascii(head,d1,io_save,file_save,ierr)
c            end if
c499         continue
          end if
c
c         Calculate ramp function and parameters
c
          if (lramp) then
            if (ipick .eq. 0) then
              write(io,*) '  Redi_wave - tt required'
              go to 599
            end if
            do ii = 1, npts
              d1(ii) = idata(ii)
            end do
            td1 = 0.d0
            td2 = 0.d0
c
c           use the first pick only - so assume it is good
c
            ip = 1
            p_time = pick(ip).ptime
c           call est_ramp(head,d1,p_time,dist,icomp,ic,lclip,
c     .          io_out,ierr)
            if (ierr .ne. 0) then
              write(io,*) '  Redi_wave - error in RF est ',ierr
              go to 599
            end if
c           if (lsave) then
c             file_save = filen(1:kf)//'.ramp'
c             call save_to_ascii(head,d1,io_save,file_save,ierr)
c           end if
599         continue
          end if
c
c         reinitialize the variables if we are through
c
          if (icomp .ge. 3) then
            icomp = 0
            do kk = 1, 3
              ic(kk) = 0
            end do
          end if
c
c       end of if statement over 3 component data
c
        end if
c
999     continue
c
c       end of loop over stuff with filen
c
       end do
c
       stop
       end

