c  anirec_synth_circle - anirec to compute synthetics for a specified model and
c  a specified station location from one seismic record read from in_recfunk
c  the code will take the station location  and generate a ring of 120 events that are 90 degrees
c  away, one for every 3 degrees back-azimuth.  To do this we generate a great circle path
c  with the station at the pole of the great circle
c  modified to accept SACfiles with start time info in the header, in A-marker (ahead(9)) 
c  because the timing information is not used, the value of A-marker is not important.
c  However, in the synthetics, the A-marker and the T1-marker variables in the header are set to 5.0 seconds
c  for easy analysis by recfunk codes that use the header information for timing.
c  the synthetic files are s_NNN.bh[zrt], where NNN is the back azimuth.
c  the code writes in_recpick_circle, a list of the filenames for recfunk_pick and other recfunk codes that
c  expect the start time in the SAC header
c
c  written originally 12/3/03 JJP
c  3/3/04 -- Added Kluge -- I have added code that forces the phase velocity to
c  fluctuate with BZ to mimic the effect of a dipping layer, changing its incidence angle with a 
c  cos(\phi) dependence.  This kluge mimics the P-SV conversion variation caused by a dipping interface
c  but not the P-SH conversion behavior.  This kluge was used for the last check on RF interpretation in
c  Park et al 2004 in JGR, the paper on the Cascadia subduction zone.  Look for the word "kluge" in the code
c  and make certain that the code between begin-kluge and end-kluge comment-lines is commented
c  
c  xf77 -o /Users/jjpark/bin/anirec_synth_circle anirec_synth_circle.f /Users/jjpark/Plotxy/plotlib.a /Users/jjpark/Ritz/eislib.a /Users/jjpark/Ritz/jlib.a
c
c  anirec - program to calculate receiver-function response of 
c  a stack of anisotropic layers to a plane wave incident from below
c  cannibalized from aniprop.f  11/18/95
c  for hexagonally symmetric media
c  reads fast axis orientation, constants A,B,C,D,E from file animodel
c  calculate quadratic eigenvalue problem based on the Christoffel matrix
c  see appendix of P. Shearer's thesis
c
c  read model, phase velocity of incident wave, P, SV, or SH
c
c  calc the eigenvector decomps for the layers
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
c  find distortion of reference wavelet
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*80 name,namo,title,xlabel,ylabel(3),xlabel2
      character *4 chead(158)
      complex*16 pp,u0,ee,z1,z0,xnu,e1,e2,zla,xl
      complex*16  rt,tt,rt0,trc,pfac,u,resp,zz
      real*4 cc4,zz4,frqq,amn,amx,baseline,dat4,a,ahead,dat44
      real*4 d_p,d_pkp,s_p,s_pkp
      common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp
      common/data/a(40000)
      common/header/ahead(158)
      common/stfff/w(3,101),t(3,3),ttl(3,3),s(3,3),stl(3,3),
     x                                           r(3,3),x(3),y(3)
      common/model/z(100),dz(100),rho(101),vp(101),vs(101),vp2(101),
     x               vp4(101),vs2(101),vss(101)
      common/model2/xmu(101),xla(101),xmu2(101),xla2(101),xla4(101)
      common/propag/xnu(6,101),xl(6,100),pfac(6,3),u(3,6)
      common/mstff/qq(6,6),wr(6),wi(6),zr(6,6),zi(6,6),iv(6),fv(6)
      common/pstff/pp(3),u0(3),ee(6,6,101),e1(6,6),e2(6,6),zla(6)
      common/rstff/rt(3,3,101),tt(3,3,101),rt0(3,3),trc(3,3)
      common/nstff/cc4(8200),zz4(8200),dat4(8200,3,3),ccc(8200)
      common/nstfff/dat44(8200,3,3)
      common/disper/resp(3,3,4100),frqq(4100)
      common/disper2/roota(101),rootb(101),jtrval(101),kroots(4100)
      common/evanes/ievan(10000)
      common/rotstuf/pol(3),ev(3),aux(3),evrot(3,3),ev1(3)
      dimension iah(158)
      equivalence (iah,ahead),(chead,ahead)
      data pi/3.14159265358979d0/,eps/1.d-6/,tol/1.d-3/
c  we reduce the condition numbers of matrices by 
c  normalizing physical quantities to make them dimensionless
c  Ill use the normal-mode normalizations
c  which are a little peculiar for the crust, but what the hey!
      rbar=5.515d3
      ren=1.075190645d-3
      radi=6.371d6
      rdi=pi/180.
      vbar=ren*radi
      con=rbar*vbar**2
      pi=3.14159265358979
      rad=180./pi
      z1=dcmplx(1.d0,0.d0)
      z0=dcmplx(0.d0,0.d0)
c  Notes on the angle conventions for w-hat, the axis of symmetry:
c  
c  In the anisotropic reflectivity code, subroutine matget *assumes* a
c  coordinate system in which z is down (anti-vertical), x is the radial
c  direction, and y is anti-transverse.  Therefore, the position angles
c  theta,phi for w-hat are tilt relative to down, and azimuth defined as a
c  rotation from x towards y.  This rotation is CCW if viewed from below,
c  and CW if viewed from above.  Since w-hat and -(w-hat) define the same
c  axis of symmetry, the position angles *also* can be defined as
c  theta=(tilt from vertical) and phi=(rotation from anti-x (anti-radial)
c  towards anti-y (transverse)).  Viewed from above, this phi rotation is
c  CW, and defines the strike of w-hat relative to the arrival azimuth of
c  the wave.  
c  
c  In order to compute seismograms for a variety of back-azimuths, we
c  assume that the default (BAZ=0)  is a wave approaching from the north, so that
c  radial is south and transverse is east.  Due to this orientation, the
c  synthetic code interprets the layered model as having w-hat position
c  angles defined as theta=(tilt from vertical) phi=(strike CW from N).
c  For an event at back-azimuth psi (CW from N), routine matget rotates w-hat
c  from geographic coordinates to ray-based coordinates before computing  
c  reflectivity matrices.  If a wave arrives at back-azimuth psi, the strike
c  of the axis of symmetry w-hat relative to its arrival azimuth is
c  phi'=phi-psi.  The code performs this rotation with this code in
c  subroutine matget, for w-hat azimuth "az":
c  
c        caz=dcos(rdi*az)
c        saz=-dsin(rdi*az)   ! sin(-az)
c        do n=1,nlp
c          ww(3)=w(3,n)
c          ww(1)=w(1,n)*caz-w(2,n)*saz
c          ww(2)=w(1,n)*saz+w(2,n)*caz
c  ...
c  
c  In this manner, the axes of symmetry of the model, saved in array w(.,.),
c  are never modified.
c  in the driver code, "az" is the variable "baz" for back azimuth
c  
c  default baz = 0
      baz=0.
c  set up for 20 sps data
c      dt=0.05
      print *,'enter sample time e.g. 0.05 for 20 sps'
      read(5,*) dt
c  nfrq must be .le.npad/2
      npad=8192
      nfrq=4096
      dur=npad*dt
      df=1./dur
      frqmax=nfrq*df
      print *,'dt,df,duration of record,max frequency:',dt,df,dur,frqmax
c      print *,'NOTE: cosine^2 taper will be applied up to fmax'
      print *,'input model? (space-return: animodel)'
      read(5,102) name
 102  format(a)
      if(name(1:1).eq.' ') then
        open(7,file='animodel',form='formatted')
      else
        open(7,file=name,form='formatted')
      endif
      read(7,102) title
      print *,title
      read(7,*) nl
c  read in theta,phi in degrees - polar coords of fast axis
      nlp=nl+1
      nlm=nl-1
      do i=1,nlp
        read(7,*) theta,phi
        w(1,i)=dsin(rdi*theta)*dcos(rdi*phi)
        w(2,i)=dsin(rdi*theta)*dsin(rdi*phi)
        w(3,i)=dcos(rdi*theta)
        print *,(w(j,i),j=1,3)
c  read depth to ith interface, vp (m/sec), pk-to-pk cos(2th) relative P pert
c  pk-to-pk cos(4th) relative P pert, v_s,  pk-to-pk cos(2th) relative S pert
c  density (kg/m**3)
        read(7,*) z(i),vp(i),vp2(i),vp4(i),vs(i),vs2(i),rho(i)
c  recall that we interpret fractional values of b,c,e 
c  as peak-to-peak relative velocity perts.
c  therefore, e=0.02 is 2% pert to mu from slowest to fastest
        xmu(i)=rho(i)*vs(i)**2/con
        xmu2(i)=vs2(i)*xmu(i)
        xla(i)=rho(i)*vp(i)**2/con
        xla2(i)=vp2(i)*xla(i)
        xla4(i)=vp4(i)*xla(i)
        vs(i)=vs(i)/vbar
        vp(i)=vp(i)/vbar
        rho(i)=rho(i)/rbar
        z(i)=z(i)/radi
      end do
      close(7)
      do i=2,nl
        dz(i)=z(i)-z(i-1)
      end do
      dz(1)=z(1)
c  print the organ-pipe mode count for 1Hz
c  the lowest layer (nl+1) is taken as evanescent region.
      sdelay=0.
      pdelay=0.
      do i=1,nl
        sdelay=sdelay+(dz(i)/vs(i))/ren
        pdelay=pdelay+(dz(i)/vp(i))/ren
      end do
      print *, 'organ-pipe mode count at 1 Hz in stack of layers: S & P'
      print *, 'Also the S & P travel time from the halfspace of model'
      print 104,sdelay,pdelay
  104 format(2f10.1)
c  search for cmin, cmax is halfspace velocity
      cmin=vs(1)
      vss(1)=vs(1)
      do i=2,nlp
        if(cmin.gt.vs(i)) cmin=vs(i)
        vss(i)=vs(i)
      end do
  900 csmin=vs(nlp)*vbar/1000.
      cpmin=vp(nlp)*vbar/1000.
      print *,'minimum phase velocity for S wave (km/sec)',csmin
      print *,'minimum phase velocity for P wave (km/sec)',cpmin
      print *,'either BAZ is determined by the loop count'
      print *,'or, in interactive mode, BAZ is user-input'
      print *,'if BAZ=0, the wave is arriving from N'
      print *,'in anirec.bh? the N component is the -RADIAL'
      print *,'          and the E component is the TRANSVERSE'
      print *,'if BAZ=90, the wave is arriving from E'
      print *,'  in anirec.bh? the N component is the -TRANSVERSE'
      print *,'and the             E component is the -RADIAL'
      print *,'if BAZ=180, the wave is arriving from S'
      print *,'  in anirec.bh? the N component is the RADIAL'
      print *,'and the             E component is the -TRANSVERSE'
      print *,'if BAZ=270, the wave is arriving from W'
      print *,'  in anirec.bh? the N component is the TRANSVERSE'
      print *,'and the             E component is the RADIAL'
c  version for reading data files and using baz and phase velocity of the
c  phases in the observed records 
c  uses a table of phase velocity versus epicentral distance (from rfmigrate)
      call read_ttimes
      print *,'code requires a file of datafile names, end with <stop>'
      print *,'you can <grep bh in_recfunk> to extract the filenames'
      print *,'file of seismic records to read? (space-ret: in_recpick)'
      read(5,102) name
      if(name(1:1).eq.' ') then
        open(10,file='in_recpick',form='formatted')
      else
        open(10,file=name,form='formatted')
      endif
      open(11,file='in_recpick_circle',form='formatted')
   10 print *,'enter filename e.g. 1998.361.bhz'
      read(10,101) name
      do i=80,1,-1
        if(name(i:i).eq.'/') go to 246
      end do
      islash=0
  246 islash=i
c  we assume that the last character is either 'r','t',or 'z' in data filenames
c  we replace with these letters as we read the files
      do i=80,1,-1
        if(name(i:i).ne.' ') then
          namlen=i
          go to 98
        endif
      end do
      print *,name
      print *,'no filename?'
      stop
   98 continue
      name(namlen:namlen)='r'
      print *,name
      call sacread(name,a,ierr)
      dt=ahead(1) 
      nscan=iah(80)
      if(nscan.gt.45000) then
        print *,'Careful! data length greater than 45000:', nscan
        pause
      endif
c  read the station location and generate 120 event locations
c  the event location is placed at DEL=90, baz=0 by adding 90 to polat
      polat=ahead(32)
      polon=ahead(33)
      print *,polat, polon,' = polat, polon'
      pol(3)=sin(polat/rad)
      pol(1)=cos(polat/rad)*cos(polon/rad)
      pol(2)=cos(polat/rad)*sin(polon/rad)
      evlat=polat+90.
      evlon=polon
      if(evlat.gt.90.) then
        evlat=180.-evlat
	evlon=evlon+180.
      endif
      if(evlon.gt.180.) evlon=evlon-360.
      ev(3)=sin(evlat/rad)
      ev(1)=cos(evlat/rad)*cos(evlon/rad)
      ev(2)=cos(evlat/rad)*sin(evlon/rad)
      aux(3)=ev(1)*pol(2)-ev(2)*pol(1)
      aux(1)=ev(2)*pol(3)-ev(3)*pol(2)
      aux(2)=ev(3)*pol(1)-ev(1)*pol(3)
      print *,(pol(i),i=1,3), '= pol'
      print *,(ev(i),i=1,3), '= ev'
      print *,(aux(i),i=1,3), '= aux'
c  for the 3x3 rotation matrix -- we advance the source by 3 degrees around a circle
      sn=sin(3./rad)
      cs=cos(3./rad)
      do i=1,3
        do j=1,3
          evrot(i,j)=pol(i)*pol(j)
	  evrot(i,j)=evrot(i,j)+(cs*ev(i)+sn*aux(i))*ev(j)
	  evrot(i,j)=evrot(i,j)+(-sn*ev(i)+cs*aux(i))*aux(j)
	end do
      end do
      do i=1,3
        print *,(evrot(i,j),j=1,3), '= evrot row', i
      end do
c      pause
      ibaz=0
      ahead(54)=90.
      ahead(51)=pi*6371./2.
 7117 if (ibaz.ge.360) go to 950
      print *,evlat,evlon, ' = evlat,evlon'
      ahead(53)=float(ibaz)
      ahead(36)=evlat
      ahead(37)=evlon
      write(namo,1011),'s_',ibaz,'.bhz'
 1011 format(a,i3,a)
      if(namo(3:3).eq.' ') namo(3:3)='0'
      if(namo(4:4).eq.' ') namo(4:4)='0'
      write(11,101) namo
c      write(11,101) '5 85 85'
      baz=ahead(53)
      epi=ahead(54)
      print *,'BAZ,EPI',baz,epi
      gcarc=epi
c  get p-slowness -- iflag branches on the start time of the analysis window
      igcarc=gcarc
      fac=gcarc-igcarc
      print *,'igcarc,i,gcarc,fac',igcarc,i,gcarc,fac
      slow=(1.-fac)*s_p(igcarc)+fac*s_p(igcarc+1)
c KLUGE  for generating anirec synthetics for a specified epicentral distance
c  that is, a phase velocity that corresponds to a specific angle of incidence
c  then perturbs the angle of incidence to match the effect of a dipping interface
c  choose DELTA=80
c      vpref=8.1
c      igcarc=80
c      slow=s_p(igcarc)
c      sangi=vpref*slow
c      angi=asin(sangi)
c      angid=angi*rad
c      print *,angid,' = angle of incidence in degrees'
cc  perturb angle of incidence by 15 degrees at maz (east-west BAZ)
c      angid=angid+sin(baz/rad)*15.
c      print *,angid,' = perturbed angle of incidence in degrees'
c      sangi2=sin(angid/rad)
c      slow=sangi2/vpref
c  END KLUGE
c  the stack-model velocities are in SI units, s_p is in sec/km
      slw=(slow/1000.)
c  non-dimensionalize cc
      cc=1000./(slow*vbar)
      cvv=1000./slow
      print *,'phase vel (m/sec & dimensionless)',cvv,cc
c  calc the eigenvector decomps for the layers
c  need to identify upgoing P,SV,SH in the halfspace
      call matget(nl,cc,baz)
      print 101,'thru matget'
c  loop over frequency, calc reflection/transmission matrices
c  calc 3-comp transfer fct response at surface
      do jf=1,nfrq
        om=2.d0*pi*jf*df/ren
        frqq(jf)=jf*df
        call respget(nl,om,cc,resp(1,1,jf))
      end do
      print 101,'thru respget'
      cccc=cc*vbar/1000.
      write(xlabel,1002) cccc
 1002 format(' Freq (Hz) -->  phase velocity',f10.2)
      print 101,'past ylabel writing'
      write(xlabel2,1003) cccc
 1003 format(' Time (sec) -->  phase velocity',f10.2)
      ylabel(1)='Incident P wave'
      print 101,ylabel(1)
      ylabel(2)='Incident SV wave'
      ylabel(3)='Incident SH wave'
c version for data-supplied cc and baz
      ipulse=3 
      if(ipulse.le.2) then
        print *,'enter wavelet period in seconds'
        read(5,*) per
      endif
      t1=0.
      dur=npad*dt
      t1=0.
      t2=dur
c version for gridsearch cc and baz
      t1=0.
      t2=100.0
c  end versions
      npts=t2/dt+1
      nst=t1/dt+1
c  we start at dt, with duration 2T -- ONE CYCLE
      npul=per/dt
      do i=1,npad
        zz4(i)=0.d0
      end do
      fac=2.d0*pi/per
      if(ipulse.eq.1) then
        do i=1,npul
          time=i*dt
          zz4(i)=(dsin(fac*time/2.d0))**2
        end do
      elseif(ipulse.eq.2) then
        do i=1,npul
          time=i*dt
          zz4(i)=dsin(fac*time)*(dsin(fac*time/2.d0))**2
        end do
      elseif(ipulse.eq.3) then
        do i=1,160
          xx=0.05*(i-80)
c          zz4(i)=xx*cos(2.*pi*xx/sqrt(1.+xx**2))*exp(-0.4*xx**2)
          zz4(i)=abs(cos(5.*pi*xx/sqrt(1.+xx**2)))*exp(-0.3*xx**2)
c          zz4(i)=abs(cos(2.*pi*xx/sqrt(1.+xx**2)))*exp(-0.2*xx**2)
        end do
      else
        do i=1,160
          xx=0.05*(i-80)
          zz4(i)=cos(5.*pi*xx/sqrt(1.+4*xx**2))*exp(-4.0*xx**2)
        end do
      endif
c      print *,(zz4(i),i=1,npul)
      call refft(zz4,npad,1,1)
c  zero the DC and Nyquist 
c  ick switches the sign of y and z components to transverse & vertical
      zz4(1)=0.
      zz4(2)=0.
      print 101,ylabel(1)
      do itype=1,3
        ick=1
        baseline=0.
        do k=1,3
          cc4(1)=0.
          cc4(2)=0.
          if(k.gt.1) ick=-1
          do jf=1,nfrq
            zz=ick*dcmplx(dble(zz4(2*jf+1)),dble(zz4(2*jf+2)))
            zz=zz*resp(k,itype,jf)
            cc4(2*jf+1)=dreal(zz)
            cc4(2*jf+2)=dimag(zz)
          end do
          do jf=2*nfrq+3,npad
            cc4(jf)=0.
          end do
          call refft(cc4,npad,-1,-1)
          amx=cc4(1)
          amn=cc4(1)
          do i=1,npts
            amx=amax1(amx,cc4(i+nst))  
            amn=amin1(amn,cc4(i+nst))  
          end do
          baseline=baseline-amn
          do i=1,npad
            dat4(i,k,itype)=cc4(i)+baseline
            dat44(i,k,itype)=cc4(i)
          end do
          baseline=baseline+amx
        end do
      end do
      cc4(1)=nst*dt
      cc4(2)=dt
c      do ity=1,3
c        do k=1,3
c          call plotit(cc4,dat4(nst,k,ity),dum,npts,title,xlabel2,
cc     x        ylabel(ity),1,0,0,0,(k/3))
c     x        ylabel(ity),1,0,0,0,(30+ity)*(k/3))
c        end do
c      end do
C        if(ity.eq.3) then  ! INCIDENT SH
C        if(ity.eq.2) then  ! INCIDENT SV
c        if(ity.eq.1) then  ! INCIDENT P -- OUTPUT TO FILE
      print 101,ylabel(1)
      ity=1
c use organ-pipe count to correct for traveltime of Swave thru stack
c          tdelay=sdelay-13.
c use organ-pipe count to correct for traveltime of Pwave thru stack
c          tdelay=pdelay-20.-300000./cvv
	  tdelay=pdelay-10.
c  no correction
c          tdelay=0.
          nst0=tdelay/ahead(1)
c          if(nst0.lt.0) nst0=0
c comp1=radial comp2=transverse comp3=vert
c  we transform radial & transverse into east & north
          ahead(53)=baz
c  here we rotate the horizontals 
c  if we are in the loop, BAZ is determined by the loop count
c  if we are in interactive mode, BAZ=0, and the wave is arriving from N
c  if BAZ=0, the wave is arriving from N
c    in anirec.bh? the N component is the -RADIAL
c              and the E component is the TRANSVERSE
c  if BAZ=90, the wave is arriving from E
c    in anirec.bh? the N component is the -TRANSVERSE
c  and the             E component is the -RADIAL
c  if BAZ=180, the wave is arriving from S
c    in anirec.bh? the N component is the RADIAL
c  and the             E component is the -TRANSVERSE
c  if BAZ=270, the wave is arriving from W
c    in anirec.bh? the N component is the TRANSVERSE
c  and the             E component is the RADIAL
      cs=cos(rdi*baz)
      sn=sin(rdi*baz)
      iah(80)=npad
c  set the start of the time window for analysis (5 seconds after start)
c  is set in the T1 marker of the SAC header, and in the A-marker for older code
      ahead(12)=5.0
      ahead(9)=5.0
c  write BOTH radial/transverse and north/east horizontal components
c  the radial/transverse SAC files are less vulnerable to BAZ confusion in SAC
      write(namo,101) namo(1:8)//'r'
      print *,islash,namo
      ahead(58)=baz+180.
      ahead(59)=90.
      chead(151)='BHR '
c  the synthetic may suffer a wrap-around alias, with late-arriving phases in the 
c  early part of record.  We therefore save the wrap-around
      if(nst0.gt.0) then
        do i=1,npad-nst0
          a(i)=dat44(nst0+i,1,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,nst0
          a(npad-nst0+i)=dat44(i,1,ity)+0.0001*(rand(0)-0.5)
        end do
      else
        do i=1,npad+nst0
          a(-nst0+i)=dat44(i,1,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,-nst0
          a(i)=dat44(npad+nst0+i,1,ity)+0.0001*(rand(0)-0.5)
        end do        
      endif      
c      print 101, (chead(i),i=111,158)
      call sacout(namo,a)
      write(namo,101) namo(1:8)//'t'
      print *,islash,namo
      ahead(58)=baz+90.
      ahead(59)=90.
      chead(151)='BHT '
      if(nst0.gt.0) then
        do i=1,npad-nst0
          a(i)=dat44(nst0+i,2,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,nst0
          a(npad-nst0+i)=dat44(i,2,ity)+0.0001*(rand(0)-0.5)
        end do
      else
        do i=1,npad+nst0
          a(-nst0+i)=dat44(i,2,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,-nst0
          a(i)=dat44(npad+nst0+i,2,ity)+0.0001*(rand(0)-0.5)
        end do        
      endif
c      print 101, (chead(i),i=111,158)
      call sacout(namo,a)
      write(namo,101) namo(1:8)//'z'
      print *,islash,namo
      ahead(58)=0.
      ahead(59)=0.
      chead(151)='BHZ '
      if(nst0.gt.0) then
        do i=1,npad-nst0
          a(i)=dat44(nst0+i,3,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,nst0
          a(npad-nst0+i)=dat44(i,3,ity)+0.0001*(rand(0)-0.5)
        end do
      else
        do i=1,npad+nst0
          a(-nst0+i)=dat44(i,3,ity)+0.0001*(rand(0)-0.5)
        end do
        do i=1,-nst0
          a(i)=dat44(npad+nst0+i,3,ity)+0.0001*(rand(0)-0.5)
        end do        
      endif
c      print 101, (chead(i),i=111,158)
      call sacout(namo,a)
  103 format(a,i2,a,i3,a)  
c  now we rotate the event position around the great circle 90 degrees from station
      do i=1,3
        ev1(i)=0.
	do j=1,3
	  ev1(i)=ev1(i)+evrot(i,j)*ev(j)
	end do
      end do
      ibaz=ibaz+3
      do i=1,3
        ev(i)=ev1(i)
      end do
      print *,ev(1),ev(2),ev(3)
      evlat=rad*asin(ev(3))
      evlon=rad*atan2(ev(2),ev(1))
      print *,evlat,evlon,' = evlat evlon'      
      go to 7117
  950 close(10)
      write(11,101) 'stop'
      close(11)
  101 format(80a)
      stop
      end
      subroutine matget(nl,cc,az)
c  SPECIAL VERSION: az rotates the w-hat vector by -az degrees
c  returns stress-displacement vectors for a stack of anisotropic layers
c  P waves may be evanescent, but S waves are oscillatory in the stack
c  the weirdness seen in the surface wave code should not appear in 
c  a receiver function code
c  however, the iev parameter is retained to avoid leaving timebombs 
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      complex*16 pp,u0,ee,pw,uw,pu,z1,z0,xnu,eye,e1,e2,zla,rtm
      complex*16 pfac,u,xl
      common/stfff/w(3,101),t(3,3),ttl(3,3),s(3,3),stl(3,3),
     x                     			 r(3,3),x(3),y(3)
      common/model/z(100),dz(100),rho(101),vp(101),vs(101),vp2(101),
     x               vp4(101),vs2(101),vss(101)
      common/model2/xmu(101),xla(101),xmu2(101),xla2(101),xla4(101)
      common/propag/xnu(6,101),xl(6,100),pfac(6,3),u(3,6)
      common/rrt/rtm(6,6,100)
      common/mstff/qq(6,6),wr(6),wi(6),zr(6,6),zi(6,6),iv(6),fv(6)
      common/pstff/pp(3),u0(3),ee(6,6,101),e1(6,6),e2(6,6),zla(6)
      common/qstff/qi(6,6),xr(6),xi(6),yr(6),yi(6),ips(3)
      dimension ww(3)
      data pi/3.14159265358979d0/,eps/1.d-6/,tol/1.d-7/
c  set iev=1 ** should be superfluous, save for now
c  toggle to iev=0 if there is a purely propagating wave in the top layer n=1
      iev=1
      z1=dcmplx(1.d0,0.d0)
      z0=dcmplx(0.d0,0.d0)
      eye=dcmplx(0.d0,1.d0)
      rbar=5.515d3
      ren=1.075190645d-3
      radi=6.371d6
      rdi=pi/180.
      vbar=radi*ren
      con=rbar*radi*radi*ren*ren
      nlp=nl+1
      nlm=nl-1
c  first calculate vertical wavenumbers and propagating waves for each layer
c   requires an eigenvector problem be solved
c  in general, the evanescent vertical wavenumbers have nonzero real parts 
c   complex exponential fct is used to avoid endless branching
c  horizontal slowness p_x
      px=1.d0/cc
      caz=dcos(rdi*az)
      saz=-dsin(rdi*az)   ! sin(-az)
      do n=1,nlp
        ww(3)=w(3,n)
        ww(1)=w(1,n)*caz-w(2,n)*saz
        ww(2)=w(1,n)*saz+w(2,n)*caz
        a=xla(n)
        b=xla2(n)
        c=xla4(n)
        d=xmu(n)
        e=xmu2(n)
c        print *,'a,b,c,d,e',a,b,c,d,e
        fact=8.d0*ww(1)*ww(1)*c+2.d0*e
        facs=16.d0*ww(1)*ww(3)*c
        facr=8.d0*ww(3)*ww(3)*c+2.d0*e
c        print *,'a,b,c,d,e',a,b,c,d,e
c        print *,'w(.,n),fact,facs,facr',(ww(l),l=1,3),fact,facs,facr
        do i=1,3
c  first the what-0-what tensor
          do j=1,3
            t(j,i)=fact*ww(j)*ww(i)
            s(j,i)=facs*ww(j)*ww(i)
            r(j,i)=facr*ww(j)*ww(i)
          end do
c  next the identity tensor - correct an error on 7/6/95
          t(i,i)=t(i,i)+d+e*(2.d0*ww(1)*ww(1)-1.d0)
          s(i,i)=s(i,i)+4.d0*e*ww(1)*ww(3)
          r(i,i)=r(i,i)+d+e*(2.d0*ww(3)*ww(3)-1.d0)
        end do 
c        print 101,(ww(i),i=1,3)
c        print 101,fact,facs,facr
c        print *,'t,s,r'
c        print 101,((t(i,j),j=1,3),i=1,3)
c        print 101,((s(i,j),j=1,3),i=1,3)
c        print 101,((r(i,j),j=1,3),i=1,3)
        fac=b-4.d0*c-2.d0*e
c  next the what-0-xhat and what-0-zhat tensors
        do i=1,3
          t(1,i)=t(1,i)+fac*ww(1)*ww(i)
          t(i,1)=t(i,1)+fac*ww(1)*ww(i)
          s(1,i)=s(1,i)+fac*ww(3)*ww(i)
          s(i,1)=s(i,1)+fac*ww(3)*ww(i)
          s(3,i)=s(3,i)+fac*ww(1)*ww(i)
          s(i,3)=s(i,3)+fac*ww(1)*ww(i)
          r(3,i)=r(3,i)+fac*ww(3)*ww(i)
          r(i,3)=r(i,3)+fac*ww(3)*ww(i)
        end do
        fac=a-b+c-d+e
c  finally the xhat-0-xhat, zhat-0-zhat, xhat-0-zhat, zhat-0-xhat tensors
        t(1,1)=t(1,1)+fac
        s(3,1)=s(3,1)+fac
        s(1,3)=s(1,3)+fac
        r(3,3)=r(3,3)+fac
c  mult by horizontal slowness and calc the modified T-matrix 
        do i=1,3
          do j=1,3
            t(j,i)=t(j,i)*px*px
            s(j,i)=s(j,i)*px
          end do
          t(i,i)=t(i,i)-rho(n)
        end do
c  calculate R**(-1).S, R**(-1).T, using routine solve
        nn=3
        do i=1,3
          do j=1,3
            y(j)=s(j,i)
          end do
          call solve(nn,r,x,y)
          do j=1,3
            stl(j,i)=x(j)
          end do
          nn=-3
        end do
        do i=1,3
          do j=1,3
            y(j)=t(j,i)
          end do
          call solve(nn,r,x,y)
          do j=1,3
            ttl(j,i)=x(j)
          end do
        end do
c  fill the 6x6 Q-matrix
        do i=1,3
          do j=1,3
            qq(j,i)=-stl(j,i)
            qq(j,i+3)=-ttl(j,i)
            qq(j+3,i)=0.d0
            qq(j+3,i+3)=0.d0
          end do
          qq(i+3,i)=1.d0
        end do
c  solve eigenvalue problem for polarization vectors and vertical slownesses
c  matrix system is nonsymmetric real valued
c  solution from the eispack guide
        call balanc(6,6,qq,is1,is2,fv)
        call elmhes(6,6,is1,is2,qq,iv)
        call eltran(6,6,is1,is2,qq,iv,zr)
        call hqr2(6,6,is1,is2,qq,wr,wi,zr,ierr)
        if(ierr.ne.0) then
          print *, ierr,'   error!'
          stop
        endif      
        call balbak(6,6,is1,is2,fv,6,zr)
c        print *,'for layer',n
c        print *, 'for phase velocity',cc,'  the vertical slownesses are'
c        print 101,(wr(i),wi(i),i=1,6)
c        pause
  101 format(6g12.4)
c  eigenvector unpacking, see EISPACK guide, page 88
c  bad eigenvector order is flagged by wi(i)>0. for odd i
        iflag=0
        do i=1,6
          if(wi(i).eq.0.d0) then
            if(n.eq.1) iev=0
            do j=1,6
              zi(j,i)=0.d0
            end do
          elseif(wi(i).gt.0.d0) then
c  bad eigenvector order is flagged by wi(i)>0 for even i
            if((i/2)*2.eq.i) then
              iflag=iflag+1
              iv(iflag)=i
            endif
            do j=1,6
              zi(j,i)=zr(j,i+1)
            end do
          else
            do j=1,6
              zi(j,i)=-zi(j,i-1)
              zr(j,i)=zr(j,i-1)
            end do
          endif
c  normalize by the last three indices
          sum=0.d0
          do j=4,6
            sum=sum+zr(j,i)**2+zi(j,i)**2
          end do
          sum=dsqrt(sum)
          do j=1,6
            zr(j,i)=zr(j,i)/sum
            zi(j,i)=zi(j,i)/sum
          end do         
        end do
c  assemble the stress-displacement vectors
c  calculate the traction components, with i removed
        pp(1)=dcmplx(px,0.d0)
        pp(2)=z0
        do k=1,6
          pp(3)=dcmplx(wr(k),wi(k))
          do i=1,3
            u0(i)=dcmplx(zr(i+3,k),zi(i+3,k))
          end do
          pu=z0
          pw=z0
          uw=z0
          abcde=a-b+c-2.d0*d+2.d0*e
          bce=b-4.d0*c-4.d0*e
          de=d-e
          do i=1,3
            pu=pu+pp(i)*u0(i)
            pw=pw+pp(i)*ww(i)
            uw=uw+u0(i)*ww(i)
          end do
          do i=1,3
            e1(i,k)=u0(i)
            e1(i+3,k)=ww(i)*(pu*ww(3)*bce+8.d0*pw*uw*ww(3)*c
     x                       +2.d0*(pw*u0(3)+uw*pp(3))*e)
            e1(i+3,k)=e1(i+3,k)+pp(i)*(u0(3)*de+2.d0*uw*ww(3)*e)
            e1(i+3,k)=e1(i+3,k)+u0(i)*(pp(3)*de+2.d0*pw*ww(3)*e)
          end do
          e1(6,k)=e1(6,k)+pu*abcde+pw*uw*bce
c  almost lastly, mult traction by i
          do i=1,3
            e1(i+3,k)=eye*e1(i+3,k)
          end do
        end do
c  reorder into upgoing and downgoing waves
c  we use the exp(-i*omega*t) convention with z increasing downward
c  so downgoing oscillatory waves have p_z>0, k_z real
c  downgoing evanescent waves have Im(p_z)>0
c  if the axis of symmetry is tilted, there are cases where a pair of 
c  near-horizontal plane waves will be both upgoing or both downgoing
c  since Chen's algorithm depends on a 3,3 split, we must adopt a kluge
c  similarly, there are cases where the EISPACK routines dont return
c  the vertical wavenumbers in ordered pairs, but mix them up a bit
c  this seems to cause problems, so a fix is necessary
c
c  first, test for bad eigenvector order, switch k-1->k+1, k->k-1, k+1->k
c   worst case is iflag=2, real,imag1+,imag1-,imag2+,imag2-,real
        if(iflag.gt.0) then
          do i=1,iflag
            k=iv(i)
            wrr=wr(k-1)
            wii=wi(k-1)
            wr(k-1)=wr(k)
            wi(k-1)=wi(k)
            wr(k)=wr(k+1)
            wi(k)=wi(k+1)
            wr(k+1)=wrr
            wi(k+1)=wii
            do j=1,6
              pu=e1(j,k-1)
              e1(j,k-1)=e1(j,k)
              e1(j,k)=e1(j,k+1)
              e1(j,k+1)=pu
            end do
          end do
        endif
c  second, divide into upgoing and downgoing waves
        isum=0
        do k=1,6
          iv(k)=0
          if(wi(k).eq.0.d0.and.wr(k).gt.0) iv(k)=1
          if(wi(k).gt.0.d0) iv(k)=1
          isum=isum+iv(k)
        end do
c  if up and downgoing cohorts are not equal, switch the sense of the
c  pure-oscillatory wave with smallest wavenumber
  140   continue
        if(isum.ne.3) then
          wr0=0.d0
          do k=1,6
            wr0=dmax1(wr0,dabs(wr(k)))
          end do
          do k=1,6
            if(wi(k).eq.0.d0) then
              if(dabs(wr(k)).lt.wr0) then
                wr0=dabs(wr(k))
                kk=k
              endif
            endif
          end do
          if(iv(kk).eq.0) then
            iv(kk)=1
          else
            iv(kk)=0
          endif
c  check that we have equal up/down cohorts
          isum=0
          do k=1,6
            isum=isum+iv(k)
          end do
          go to 140
        endif
        jdown=1
        jup=4
c        print *,'for layer',n,'  the vert wavenums are (0=up,1=dn)'
        do k=1,6
          if(iv(k).eq.1) then
            ki=jdown
            jdown=jdown+1
          else
            ki=jup
            jup=jup+1
          endif
          do i=1,6
            ee(i,ki,n)=e1(i,k)
          end do
c  incorporate the factor of i into the stored vertical slowness 
          xnu(ki,n)=dcmplx(-wi(k),wr(k))
        end do
 1008 format(a,2g15.6,a,2g15.6)
      end do
c  now, must identify which upgoing waves in the halfspace are P,SV,SH
c  crud, this goes back to array ee
c  3: SH is y-motion
c  2: SV is (-sqrt((1/vs)**2-p_x**2),0,-p_x)  ! recall that z points down
c  1: P is (p_x,0,-sqrt((1/vp)**2-p_x**2)
c  so we branch on size of u_y, and relative sign of u_x and u_z
c      print *,'in the halfspace:'
c      do i=4,6
c        print *,'for i*k_z=',xnu(i,nlp),', the disp-stress vector is'
c        do j=1,6
c          xi(j)=dimag(ee(j,i,nlp))
c          xr(j)=dreal(ee(j,i,nlp))
c        end do
c        print 101,(xr(j),j=1,6),(xi(j),j=1,6)
c      end do
      do i=4,6
        ips(i-3)=3
        if(zabs(ee(2,i,nlp)).lt.dsqrt(tol)) then  ! not SH
          test=dreal(ee(1,i,nlp))/dreal(ee(3,i,nlp))
          if(test.gt.0.d0) then
            ips(i-3)=2
          else
            ips(i-3)=1
          endif
        endif
      end do
      print *,'wave types:',(ips(i),i=1,3)
      if(ips(1).eq.ips(2).or.ips(1).eq.ips(3).or.ips(2).eq.ips(3)) then
        print *,'WHOA! Wave types are not distinct!'
        print *,'Is there anisotropy in the halfspace?  Remove it!!'
        stop  
      endif     
      return
      end
      subroutine matget_old(nl,cc)
c  returns stress-displacement vectors for a stack of anisotropic layers
c  P waves may be evanescent, but S waves are oscillatory in the stack
c  the weirdness seen in the surface wave code should not appear in 
c  a receiver function code
c  however, the iev parameter is retained to avoid leaving timebombs 
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      complex*16 pp,u0,ee,pw,uw,pu,z1,z0,zz,xnu,eye,e1,e2,zla,rtm
      complex*16 pfac,u,xl
      common/stfff/w(3,101),t(3,3),ttl(3,3),s(3,3),stl(3,3),
     x                     			 r(3,3),x(3),y(3)
      common/model/z(100),dz(100),rho(101),vp(101),vs(101),vp2(101),
     x               vp4(101),vs2(101),vss(101)
      common/model2/xmu(101),xla(101),xmu2(101),xla2(101),xla4(101)
      common/propag/xnu(6,101),xl(6,100),pfac(6,3),u(3,6)
      common/defect/idfct(4,101),adf(2,101)
      common/rrt/rtm(6,6,100)
      common/mstff/qq(6,6),wr(6),wi(6),zr(6,6),zi(6,6),iv(6),fv(6)
      common/pstff/pp(3),u0(3),ee(6,6,101),e1(6,6),e2(6,6),zla(6)
      common/qstff/qi(6,6),xr(6),xi(6),yr(6),yi(6),ips(3)
      data pi/3.14159265358979d0/,eps/1.d-6/,tol/1.d-7/
c  set iev=1 ** should be superfluous, but save for now
c  toggle to iev=0 if there is a purely propagating wave in the top layer n=1
      iev=1
      z1=dcmplx(1.d0,0.d0)
      z0=dcmplx(0.d0,0.d0)
      eye=dcmplx(0.d0,1.d0)
      rbar=5.515d3
      ren=1.075190645d-3
      radi=6.371d6
      vbar=radi*ren
      con=rbar*radi*radi*ren*ren
      nlp=nl+1
      nlm=nl-1
c  first calculate vertical wavenumbers and propagating waves for each layer
c   requires an eigenvector problem be solved
c  in general, the evanescent vertical wavenumbers have nonzero real parts 
c   complex exponential fct is used to avoid endless branching
c  horizontal slowness p_x
      px=1.d0/cc
      do n=1,nlp
        a=xla(n)
        b=xla2(n)
        c=xla4(n)
        d=xmu(n)
        e=xmu2(n)
c        print *,'a,b,c,d,e',a,b,c,d,e
        fact=8.d0*w(1,n)*w(1,n)*c+2.d0*e
        facs=16.d0*w(1,n)*w(3,n)*c
        facr=8.d0*w(3,n)*w(3,n)*c+2.d0*e
        do i=1,3
c  first the what-0-what tensor
          do j=1,3
            t(j,i)=fact*w(j,n)*w(i,n)
            s(j,i)=facs*w(j,n)*w(i,n)
            r(j,i)=facr*w(j,n)*w(i,n)
          end do
c  next the identity tensor - correct an error on 7/6/95
          t(i,i)=t(i,i)+d+e*(2.d0*w(1,n)*w(1,n)-1.d0)
          s(i,i)=s(i,i)+4.d0*e*w(1,n)*w(3,n)
          r(i,i)=r(i,i)+d+e*(2.d0*w(3,n)*w(3,n)-1.d0)
        end do 
c        print 101,(w(i,n),i=1,3)
c        print 101,fact,facs,facr
c        print *,'t,s,r'
c        print 101,((t(i,j),j=1,3),i=1,3)
c        print 101,((s(i,j),j=1,3),i=1,3)
c        print 101,((r(i,j),j=1,3),i=1,3)
        fac=b-4.d0*c-2.d0*e
c  next the what-0-xhat and what-0-zhat tensors
        do i=1,3
          t(1,i)=t(1,i)+fac*w(1,n)*w(i,n)
          t(i,1)=t(i,1)+fac*w(1,n)*w(i,n)
          s(1,i)=s(1,i)+fac*w(3,n)*w(i,n)
          s(i,1)=s(i,1)+fac*w(3,n)*w(i,n)
          s(3,i)=s(3,i)+fac*w(1,n)*w(i,n)
          s(i,3)=s(i,3)+fac*w(1,n)*w(i,n)
          r(3,i)=r(3,i)+fac*w(3,n)*w(i,n)
          r(i,3)=r(i,3)+fac*w(3,n)*w(i,n)
        end do
        fac=a-b+c-d+e
c  finally the xhat-0-xhat, zhat-0-zhat, xhat-0-zhat, zhat-0-xhat tensors
        t(1,1)=t(1,1)+fac
        s(3,1)=s(3,1)+fac
        s(1,3)=s(1,3)+fac
        r(3,3)=r(3,3)+fac
c  mult by horizontal slowness and calc the modified T-matrix 
        do i=1,3
          do j=1,3
            t(j,i)=t(j,i)*px*px
            s(j,i)=s(j,i)*px
          end do
          t(i,i)=t(i,i)-rho(n)
        end do
c  calculate R**(-1).S, R**(-1).T, using routine solve
        nn=3
        do i=1,3
          do j=1,3
            y(j)=s(j,i)
          end do
          call solve(nn,r,x,y)
          do j=1,3
            stl(j,i)=x(j)
          end do
          nn=-3
        end do
        do i=1,3
          do j=1,3
            y(j)=t(j,i)
          end do
          call solve(nn,r,x,y)
          do j=1,3
            ttl(j,i)=x(j)
          end do
        end do
c  fill the 6x6 Q-matrix
        do i=1,3
          do j=1,3
            qq(j,i)=-stl(j,i)
            qq(j,i+3)=-ttl(j,i)
            qq(j+3,i)=0.d0
            qq(j+3,i+3)=0.d0
          end do
          qq(i+3,i)=1.d0
        end do
c  solve eigenvalue problem for polarization vectors and vertical slownesses
c  matrix system is nonsymmetric real valued
c  solution from the eispack guide
        call balanc(6,6,qq,is1,is2,fv)
        call elmhes(6,6,is1,is2,qq,iv)
        call eltran(6,6,is1,is2,qq,iv,zr)
        call hqr2(6,6,is1,is2,qq,wr,wi,zr,ierr)
        if(ierr.ne.0) then
          print *, ierr,'   error!'
          stop
        endif      
        call balbak(6,6,is1,is2,fv,6,zr)
c        print *,'for layer',n
c        print *, 'for phase velocity',cc,'  the vertical slownesses are'
c        print 101,(wr(i),wi(i),i=1,6)
c        pause
  101 format(6g12.4)
c  eigenvector unpacking, see EISPACK guide, page 88
c  bad eigenvector order is flagged by wi(i)>0. for odd i
        iflag=0
        do i=1,6
          if(wi(i).eq.0.d0) then
            if(n.eq.1) iev=0
            do j=1,6
              zi(j,i)=0.d0
            end do
          elseif(wi(i).gt.0.d0) then
c  bad eigenvector order is flagged by wi(i)>0 for even i
            if((i/2)*2.eq.i) then
              iflag=iflag+1
              iv(iflag)=i
            endif
            do j=1,6
              zi(j,i)=zr(j,i+1)
            end do
          else
            do j=1,6
              zi(j,i)=-zi(j,i-1)
              zr(j,i)=zr(j,i-1)
            end do
          endif
c  normalize by the last three indices
          sum=0.d0
          do j=4,6
            sum=sum+zr(j,i)**2+zi(j,i)**2
          end do
          sum=dsqrt(sum)
          do j=1,6
            zr(j,i)=zr(j,i)/sum
            zi(j,i)=zi(j,i)/sum
          end do         
        end do
c  assemble the stress-displacement vectors
c  calculate the traction components, with i removed
        pp(1)=dcmplx(px,0.d0)
        pp(2)=z0
        do k=1,6
          pp(3)=dcmplx(wr(k),wi(k))
          do i=1,3
            u0(i)=dcmplx(zr(i+3,k),zi(i+3,k))
          end do
          pu=z0
          pw=z0
          uw=z0
          abcde=a-b+c-2.d0*d+2.d0*e
          bce=b-4.d0*c-4.d0*e
          de=d-e
          do i=1,3
            pu=pu+pp(i)*u0(i)
            pw=pw+pp(i)*w(i,n)
            uw=uw+u0(i)*w(i,n)
          end do
          do i=1,3
            e1(i,k)=u0(i)
            e1(i+3,k)=w(i,n)*(pu*w(3,n)*bce+8.d0*pw*uw*w(3,n)*c
     x                       +2.d0*(pw*u0(3)+uw*pp(3))*e)
            e1(i+3,k)=e1(i+3,k)+pp(i)*(u0(3)*de+2.d0*uw*w(3,n)*e)
            e1(i+3,k)=e1(i+3,k)+u0(i)*(pp(3)*de+2.d0*pw*w(3,n)*e)
          end do
          e1(6,k)=e1(6,k)+pu*abcde+pw*uw*bce
c  almost lastly, mult traction by i
          do i=1,3
            e1(i+3,k)=eye*e1(i+3,k)
          end do
        end do
c  reorder into upgoing and downgoing waves
c  we use the exp(-i*omega*t) convention with z increasing downward
c  so downgoing oscillatory waves have p_z>0, k_z real
c  downgoing evanescent waves have Im(p_z)>0
c  if the axis of symmetry is tilted, there are cases where a pair of 
c  near-horizontal plane waves will be both upgoing or both downgoing
c  since Chen's algorithm depends on a 3,3 split, we must adopt a kluge
c  similarly, there are cases where the EISPACK routines dont return
c  the vertical wavenumbers in ordered pairs, but mix them up a bit
c  this seems to cause problems, so a fix is necessary
c
c  first, test for bad eigenvector order, switch k-1->k+1, k->k-1, k+1->k
c   worst case is iflag=2, real,imag1+,imag1-,imag2+,imag2-,real
        if(iflag.gt.0) then
          do i=1,iflag
            k=iv(i)
            wrr=wr(k-1)
            wii=wi(k-1)
            wr(k-1)=wr(k)
            wi(k-1)=wi(k)
            wr(k)=wr(k+1)
            wi(k)=wi(k+1)
            wr(k+1)=wrr
            wi(k+1)=wii
            do j=1,6
              pu=e1(j,k-1)
              e1(j,k-1)=e1(j,k)
              e1(j,k)=e1(j,k+1)
              e1(j,k+1)=pu
            end do
          end do
        endif
c  second, divide into upgoing and downgoing waves
        isum=0
        do k=1,6
          iv(k)=0
          if(wi(k).eq.0.d0.and.wr(k).gt.0) iv(k)=1
          if(wi(k).gt.0.d0) iv(k)=1
          isum=isum+iv(k)
        end do
c  if up and downgoing cohorts are not equal, switch the sense of the
c  pure-oscillatory wave with smallest wavenumber
  140   continue
        if(isum.ne.3) then
          wr0=0.d0
          do k=1,6
            wr0=dmax1(wr0,dabs(wr(k)))
          end do
          do k=1,6
            if(wi(k).eq.0.d0) then
              if(dabs(wr(k)).lt.wr0) then
                wr0=dabs(wr(k))
                kk=k
              endif
            endif
          end do
          if(iv(kk).eq.0) then
            iv(kk)=1
          else
            iv(kk)=0
          endif
c  check that we have equal up/down cohorts
          isum=0
          do k=1,6
            isum=isum+iv(k)
          end do
          go to 140
        endif
        jdown=1
        jup=4
c        print *,'for layer',n,'  the vert wavenums are (0=up,1=dn)'
 1001 format(i2,2g15.6)
        do k=1,6
          if(iv(k).eq.1) then
            ki=jdown
            jdown=jdown+1
          else
            ki=jup
            jup=jup+1
          endif
          do i=1,6
            ee(i,ki,n)=e1(i,k)
          end do
c  incorporate the factor of i into the stored vertical slowness 
          xnu(ki,n)=dcmplx(-wi(k),wr(k))
        end do
c  OK, here's where we check whether two downgoing stress-disp vectors
c  are nearly parallel - we check the dotproducts of displacement components
        do i=1,4
          idfct(i,n)=0
          adf((i+1)/2,n)=0.d0
        end do
        do i=1,2
          do j=i+1,3
            r1=0.d0
            r2=0.d0
            zz=z0
            do k=1,3
              r1=r1+zabs(ee(k,i,n))**2
              r2=r2+zabs(ee(k,j,n))**2
              zz=zz+ee(k,j,n)*conjg(ee(k,i,n))
            end do
            qqq=1.d0-zabs(zz)/dsqrt(r1*r2)
            if(qqq.lt.tol) then
              ccc=cc*vbar
              idfct(1,n)=i
              idfct(2,n)=j
c              print 1008,'vert slownesses',xnu(i,n),' and',xnu(j,n)
c  we average eigenvalues (vert slownesses)
c  and solve for eigenvector in subroutine defective
              xnu(i,n)=(xnu(i,n)+xnu(j,n))/2.d0
              xnu(j,n)=xnu(i,n)
c  calculate the extravector for defective repeated eigenvalue
              call defective(i,j,n,adf(1,n),a,b,c,d,e,px)
c              print *,i,j,n,ccc,qqq,adf(1,n)
            endif
          end do
        end do
 1008 format(a,2g15.6,a,2g15.6)
c  OK, here's where we check whether two upgoing stress-disp vectors
c  are nearly parallel - we check the dotproducts of displacement components
        do i=4,5
          do j=i+1,6
            r1=0.d0
            r2=0.d0
            zz=z0
            do k=1,3
              r1=r1+zabs(ee(k,i,n))**2
              r2=r2+zabs(ee(k,j,n))**2
              zz=zz+ee(k,j,n)*conjg(ee(k,i,n))
            end do
            qqq=1.d0-zabs(zz)/dsqrt(r1*r2)
            if(qqq.lt.tol) then
              ccc=cc*vbar
              idfct(3,n)=i
              idfct(4,n)=j
c              print 1008,'vert slownesses',xnu(i,n),' and',xnu(j,n)
c  we average the eigenvalues
              xnu(i,n)=(xnu(i,n)+xnu(j,n))/2.d0
              xnu(j,n)=xnu(i,n)
c  calculate the extravector for defective repeated eigenvalue
c  as well as coefficient (adf) of linear (z-z0) term
              call defective(i,j,n,adf(2,n),a,b,c,d,e,px)
c              print *,i,j,n,ccc,qqq,adf(2,n)
            endif
          end do
        end do
      end do
c  now, must identify which upgoing waves in the halfspace are P,SV,SH
c  crud, this goes back to array ee
c  3: SH is y-motion
c  2: SV is (-sqrt((1/vs)**2-p_x**2),0,-p_x)  ! recall that z points down
c  1: P is (p_x,0,-sqrt((1/vp)**2-p_x**2)
c  so we branch on size of u_y, and relative sign of u_x and u_z
      print *,'in the halfspace:'
      do i=4,6
        print *,'for i*k_z=',xnu(i,nlp),', the disp-stress vector is'
        do j=1,6
          xi(j)=dimag(ee(j,i,nlp))
          xr(j)=dreal(ee(j,i,nlp))
        end do
        print 101,(xr(j),j=1,6),(xi(j),j=1,6)
      end do
      do i=4,6
        ips(i-3)=3
        if(zabs(ee(2,i,nlp)).lt.dsqrt(tol)) then  ! not SH
          test=dreal(ee(1,i,nlp))/dreal(ee(3,i,nlp))
          if(test.gt.0.d0) then
            ips(i-3)=2
          else
            ips(i-3)=1
          endif
        endif
      end do
      print *,'wave types:',(ips(i),i=1,3)        
      return
      end
      subroutine respget(nl,om,cc,resp)
c  returns surface response for a stack of anisotropic layers
c  incident p,sv,sh waves with freq om and phase velocity cc
c  iev=1 if the waves are evanescent in the top layer, iev=0 otherwise
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      complex*16 pp,u0,ee,z1,z0,xnu,eye,e1,e2,zla,rtm
      complex*16 rt,tt,rt0,trc,xl,pfac,u,co,resp,ur
      common/stfff/w(3,101),t(3,3),ttl(3,3),s(3,3),stl(3,3),
     x                     			 r(3,3),x(3),y(3)
      common/model/z(100),dz(100),rho(101),vp(101),vs(101),vp2(101),
     x               vp4(101),vs2(101),vss(101)
      common/model2/xmu(101),xla(101),xmu2(101),xla2(101),xla4(101)
      common/propag/xnu(6,101),xl(6,100),pfac(6,3),u(3,6)
      common/defect/idfct(4,101),adf(2,101)
      common/rrt/rtm(6,6,100)
      common/mstff/qq(6,6),wr(6),wi(6),zr(6,6),zi(6,6),iv(6),fv(6)
      common/pstff/pp(3),u0(3),ee(6,6,101),e1(6,6),e2(6,6),zla(6)
      common/pstf2/co(6,101),ur(3)
      common/qstff/qi(6,6),xr(6),xi(6),yr(6),yi(6),ips(3)
      common/rstff/rt(3,3,101),tt(3,3,101),rt0(3,3),trc(3,3)
      dimension resp(3,3)
      data pi/3.14159265358979d0/,eps/1.d-6/,tol/1.d-7/
c  set iev=1
c  toggle to iev=0 if there is a purely propagating wave in the top layer n=1
      iev=1
      z1=dcmplx(1.d0,0.d0)
      z0=dcmplx(0.d0,0.d0)
      eye=dcmplx(0.d0,1.d0)
      rbar=5.515d3
      ren=1.075190645d-3
      radi=6.371d6
      vbar=radi*ren
      con=rbar*radi*radi*ren*ren
      nlp=nl+1
      nlm=nl-1
c  first calculate vertical wavenumbers and propagating waves for each layer
c    an eigenvector problem was solved in prior subroutine
c   results stored in array ee(6,6,101)
c  in general, the evanescent vertical wavenumbers have nonzero real parts 
c   complex exponential fct is used to avoid endless branching
c
c  calculate modified R/T coefficients
c  first calc the propagation factors
c  note that for dipping fast axes the upgoing and downgoing wavenumbers are
c  independent, so we must calc all to be safe
      do n=1,nl
        do k=1,3
          xl(k,n)=zexp(om*xnu(k,n)*dz(n))		! downgoing
          xl(k+3,n)=zexp(-om*xnu(k+3,n)*dz(n))	! upgoing
        end do
      end do
c      do i=1,6
c        print 1002,xnu(i,3),xl(i,3)
c      end do
 1002 format('i*k_z:',2g15.6,',  propfac is',2g15.6)
c  calculate modified R/T coefficients at each interface
      do n=1,nl
c  rearrange to e1: waves approaching and e2: waves leaving an interface
        do k=1,3
          do i=1,6
            e1(i,k)=ee(i,k,n+1)
            e2(i,k)=ee(i,k,n)
            e1(i,k+3)=-ee(i,k+3,n)
            e2(i,k+3)=-ee(i,k+3,n+1)
          end do
          zla(k)=xl(k,n)
          if(n.lt.nl) then
            zla(k+3)=xl(k+3,n+1)   
          else 
c  reference the upcoming wave amplitude to the top of halfspace
c  therefore no propagation factor, not valid for evanescent waves in halfspace
c  in surface wave code this is zero, so that upgoing evanescent waves vanish
            zla(k+3)=1.d0
          endif
        end do
c mult the columns of e2
        do k=1,6
          do i=1,6
            e2(i,k)=e2(i,k)*zla(k)
          end do
        end do
c  the possibility of defective matrices must be contemplated here
c  k=1,2,3 columns are downgoing in nth layer
c  k=4,5,6 columns are upgoing in (n+1)th layer
c  the vector e2(.,k1) has already been multiplied by exponential factor zla
        if(idfct(1,n).ne.0) then
          k1=idfct(1,n)
          k2=idfct(2,n)
          do i=1,6
            e2(i,k2)=e2(i,k2)+adf(1,n)*dz(n)*e2(i,k1)
          end do
        endif
c  the sign change on dz is for upgoing waves
        if(idfct(3,n+1).ne.0) then
          k1=idfct(3,n+1)
          k2=idfct(4,n+1)
          do i=1,6
            e2(i,k2)=e2(i,k2)-adf(2,n+1)*dz(n+1)*e2(i,k1)
          end do
        endif
c  in order to use csolve to invert e1, must separate into real/imag parts
c  its clumsy, but im lazy
c  we calc e1^{-1}\cdot e2\cdot \Gamma one column at a time
        do k=1,6
          do i=1,6
            qq(i,k)=dreal(e1(i,k))
            qi(i,k)=dimag(e1(i,k))
          end do
        end do
        nn=6
        do k=1,6
          do i=1,6
            yr(i)=dreal(e2(i,k))
            yi(i)=dimag(e2(i,k))
          end do
          call csolve(nn,qq,qi,xr,xi,yr,yi)
          nn=-6
          do i=1,6
            rtm(i,k,n)=dcmplx(xr(i),xi(i))
          end do
        end do  
      end do
c  calc R_ud at the free surface
c  note that first two factors in Chen (20) dont collapse
c  mult by inv-matrix one column at a time
      do k=1,3
        do i=1,3
          rt0(i,k)=ee(i+3,k+3,1)*xl(k+3,1)
          s(i,k)=dreal(ee(i+3,k,1))
          t(i,k)=dimag(ee(i+3,k,1))
        end do
      end do
c  the possibility of defective matrices must be contemplated here
c  these waves are upgoing in 1st layer
c  the sign change on dz is for upgoing waves, and xl(k1,1)=xl(k2,1)
      if(idfct(3,1).ne.0) then
        k1=idfct(3,1)
        k2=idfct(4,1)-3
        do i=1,3
          rt0(i,k2)=rt0(i,k2)-adf(2,1)*dz(1)*ee(i+3,k1,1)*xl(k1,1)
        end do
      endif
      nn=3
      do k=1,3
        do i=1,3
          yr(i)=dreal(rt0(i,k))
          yi(i)=dimag(rt0(i,k))
        end do
        call csolve(nn,s,t,xr,xi,yr,yi)
        nn=-3
        do i=1,3
          rt0(i,k)=-dcmplx(xr(i),xi(i))
        end do
      end do  
c  recursive calc of generalized R/T coefs:
c  in contrast to the surface-wave code, we start from the top layer and
c  iterate down to the halfspace
c  we also uses submatrices of generalized R/T matrix in different order
      do n=1,nl
c  first the generalized upward-transmission coef:
        do k=1,3
          do i=1,3
            trc(i,k)=z0
            if(n.gt.1) then
              do j=1,3
                trc(i,k)=trc(i,k)-rtm(i+3,j,n)*rt(j,k,n-1)
              end do
            else
c  use free-surface reflection matrix in top layer (interface "zero")
              do j=1,3
                trc(i,k)=trc(i,k)-rtm(i+3,j,n)*rt0(j,k)
              end do
            endif
          end do
          trc(k,k)=trc(k,k)+z1
        end do  
        do k=1,3
          do i=1,3
            s(i,k)=dreal(trc(i,k))
            t(i,k)=dimag(trc(i,k))
          end do
        end do
        nn=3
        do k=1,3
          do i=1,3
            yr(i)=dreal(rtm(i+3,k+3,n))
            yi(i)=dimag(rtm(i+3,k+3,n))
          end do
          call csolve(nn,s,t,xr,xi,yr,yi)
          nn=-3
          do i=1,3
            tt(i,k,n)=dcmplx(xr(i),xi(i))
          end do
        end do  
c  next the generalized reflection coef:
        do k=1,3
          do i=1,3
            trc(i,k)=z0
            if(n.gt.1) then
              do j=1,3
                trc(i,k)=trc(i,k)+rt(i,j,n-1)*tt(j,k,n)
              end do
            else
c  use free-surface reflection matrix in top layer (interface "zero")
              do j=1,3
                trc(i,k)=trc(i,k)+rt0(i,j)*tt(j,k,n)
              end do
            endif
          end do
        end do  
        do k=1,3
          do i=1,3
            rt(i,k,n)=rtm(i,k+3,n)
            do j=1,3
              rt(i,k,n)=rt(i,k,n)+rtm(i,j,n)*trc(j,k)
            end do
          end do
        end do  
      end do
 1001 format(6f14.6)
c      print *,'free-surface reflection'
c      print 1001,((rt0(i,j),j=1,3),i=1,3)
c      do n=1,nl
c        print *,'interface',n
c        print 1001,((rt(i,j,n),j=1,3),i=1,3)
c        print 1001,((tt(i,j,n),j=1,3),i=1,3)
c      end do
c  using the p,sv,sh identification, we propagate upward to the surface,
c  calculate the wave coefs in the top layer, then the particle displacement
      do iup=1,3
        do i=1,3
          co(i+3,nlp)=z0
        end do
        co(iup+3,nlp)=z1
c  from upgoing coefs in the n+1 layer, calculate
c  upgoing coefs in the nth layer, downgoing coefs in the n+1 layer
        do n=nl,1,-1
          do i=1,3
            co(i+3,n)=z0
            co(i,n+1)=z0
            do j=1,3
              co(i+3,n)=co(i+3,n)+tt(i,j,n)*co(j+3,n+1)
              co(i,n+1)=co(i,n+1)+rt(i,j,n)*co(j+3,n+1)
            end do
          end do
        end do
c  then downgoing coefs in the top layer:
        do i=1,3
          co(i,1)=z0
          do j=1,3
            co(i,1)=co(i,1)+rt0(i,j)*co(j+3,1)
          end do
        end do
c        print *,'upgoing coefs'
c        print 1001,((co(j+3,n),j=1,3),n=1,nlp)
c        print *,'downgoing coefs'
c        print 1001,((co(j,n),j=1,3),n=1,nlp)
c  calc the surface displacement
        h1=0.d0
        h2=z(1)
        do i=1,3
          ur(i)=z0
          do k=1,3
            ur(i)=ur(i)+co(k,1)*ee(i,k,1)
     x          +co(k+3,1)*ee(i,k+3,1)*(zexp(om*xnu(k+3,1)*(-h2)))
          end do
c  check for the xtra terms associated with defective matrices
          if(idfct(1,1).ne.0) then
            ii=idfct(1,1)
            jj=idfct(2,1)
            ur(i)=ur(i)+co(jj,1)*adf(1,1)*ee(i,ii,1)*(-h1)
          endif
          if(idfct(3,1).ne.0) then
            ii=idfct(3,1)
            jj=idfct(4,1)
            ur(i)=ur(i)
     x +co(jj,1)*adf(2,1)*ee(i,ii,1)*(-h2)*(zexp(om*xnu(ii,1)*(-h2)))
          endif
        end do
c  copy the surface displacement into the response matrix
        do i=1,3
          resp(i,ips(iup))=ur(i)
        end do            
      end do     
      return
      end
      subroutine defective(i,j,n,adf,a,b,c,d,e,px)
c  kluge for dealing with nearly defective propagator matrices
c  in which the eigenvectors, 
c  which represent the particle motion of upgoing and downgoing waves
c  become nearly parallel.
c  in this case the solution for system of ODEs is
c  a_1 \bf_1 e^xnu*(z-z0) + a_2*(\bf_2 + adf*(z-z0)*\bf_1)e^xnu*(z-z0)
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      complex*16 pp,u0,ee,z1,z0,znu,xnu,e1,e2,zla,xl,u,pfac,eye
      complex*16 zq1,zq2,u1,u2,zq3,xee
      common/stfff/w(3,101),t(3,3),ttl(3,3),s(3,3),stl(3,3),
     x                     			 r(3,3),x(3),y(3)
      common/propag/xnu(6,101),xl(6,100),pfac(6,3),u(3,6)
      common/defect1/zq1(3,3),zq2(3,3),u1(3),u2(3),zq3(2,2),xee(3)
      common/defect2/edr(6),edi(6),qdr(5,5),qdi(5,5),ydr(6),ydi(6)
      common/defect3/q1r(3,3),q1i(3,3),q2r(3,3),q2i(3,3),fv2(3),fv3(3)
      common/mstff/qq(6,6),wr(6),wi(6),zr(6,6),zi(6,6),iv(6),fv(6)
      common/pstff/pp(3),u0(3),ee(6,6,101),e1(6,6),e2(6,6),zla(6)
      common/qstff/qi(6,6),xr(6),xi(6),yr(6),yi(6),ips(3)
      z1=dcmplx(1.d0,0.d0)
      z0=dcmplx(0.d0,0.d0)
      eye=dcmplx(0.d0,1.d0)
c  for the extravector, need to solve system of equations
c  based on original 6x6 Q matrix
c  the plane-wave solutions generalize to the form
c  u0*e^{i*nu*(z-z0)}  and  u1*e^{i*nu*(z-z0)} + adf* u0*(z-z0)*e^{i*nu*(z-z0)}
c  u1 is the solution to 
c  (\bTtil + nu*\bStil + nu^2*\bI).u1=i*adf*(\bStil + 2*nu*\bI).u0
c  in practice, we absorb the adf factor into u1, then normalize
c  (\bTtil + nu*\bStil + nu^2*\bI).(u1/adf)=i*(\bStil + 2*nu*\bI).u0
c  since nu is the known eigenvalue of u0, the solution is easier
c  form the matrices on either side
      znu=-eye*xnu(i,n)
      do ii=1,3
        do jj=1,3
          zq1(jj,ii)=dcmplx(ttl(jj,ii),0.d0)+znu*stl(jj,ii)
          zq2(jj,ii)=dcmplx(stl(jj,ii),0.d0)
        end do
        zq1(ii,ii)=zq1(ii,ii)+znu*znu
        zq2(ii,ii)=zq2(ii,ii)+2.d0*znu
      end do
c  we wish to find the eigenvector of the near-defective matrix 
c  in the region where its eigenvectors are numerically unstable
c   we explicitly calculate the eigenvector with smallest right-eigenvalue of
c  (\bTtil + nu*\bStil + nu^2*\bI)=zq1
c  copy into real, imag matrices
      do ii=1,3
        do jj=1,3
          q1r(jj,ii)=dreal(zq1(jj,ii))
          q1i(jj,ii)=dimag(zq1(jj,ii))
        end do
      end do      
c  into eispack
      call cbal(3,3,q1r,q1i,low,igh,fv)
      call corth(3,3,low,igh,q1r,q1i,fv2,fv3)
      call comqr2(3,3,low,igh,fv2,fv3,q1r,q1i,wr,wi,q2r,q2i,ierr)
      if(ierr.ne.0) go to 400
      call cbabk2(3,3,low,igh,fv,3,q2r,q2i)
      amn=wr(1)**2+wi(1)**2
      ij=1
      do ii=2,3
        amm=wr(ii)**2+wi(ii)**2
        if(amm.lt.amn) then
          ij=ii
          amn=amm
        endif
      end do
      sum=0.d0
      do ii=1,3
        u0(ii)=dcmplx(q2r(ii,ij),q2i(ii,ij))
        sum=sum+zabs(u0(ii))**2
      end do
      sum=dsqrt(sum)
      do ii=1,3
        u0(ii)=u0(ii)/sum
      end do
c  assemble the ith stress-displacement vector
c  calculate the traction components, with i removed
      pp(1)=dcmplx(px,0.d0)
      pp(2)=z0
      pp(3)=znu
      pu=z0
      pw=z0
      uw=z0
      abcde=a-b+c-2.d0*d+2.d0*e
      bce=b-4.d0*c-4.d0*e
      de=d-e
      do ii=1,3
        pu=pu+pp(ii)*u0(ii)
        pw=pw+pp(ii)*w(ii,n)
        uw=uw+u0(ii)*w(ii,n)
      end do
      do ii=1,3
        ee(ii,i,n)= u0(ii)
        ee(ii+3,i,n)=w(ii,n)*(pu*w(3,n)*bce+8.d0*pw*uw*w(3,n)*c
     x                                  +2.d0*(pw*u0(3)+uw*pp(3))*e)
        ee(ii+3,i,n)=ee(ii+3,i,n)+pp(ii)*(u0(3)*de+2.d0*uw*w(3,n)*e)
        ee(ii+3,i,n)=ee(ii+3,i,n)+u0(ii)*(pp(3)*de+2.d0*pw*w(3,n)*e)
      end do
      ee(6,i,n)=ee(6,i,n)+pu*abcde+pw*uw*bce
c  almost lastly, mult traction by i
      do ii=1,3
        ee(ii+3,i,n)=eye*ee(ii+3,i,n)
      end do
c  extract u0 from ee(*,i,n) use it to calculate the additional traction terms
c  and store in ee(*,j,n)
c  additional traction terms involve gradient of (z-z0)
c  so can be calculated from standard formulas with \bk=zhat
c  we dont multiply by i
      pp(1)=z0
      pp(2)=z0
      pp(3)=z1
      pu=z0
      pw=z0
      uw=z0
      abcde=a-b+c-2.d0*d+2.d0*e
      bce=b-4.d0*c-4.d0*e
      de=d-e
      do ii=1,3
        u0(ii)=ee(ii,i,n)
        pu=pu+pp(ii)*u0(ii)
        pw=pw+pp(ii)*w(ii,n)
        uw=uw+u0(ii)*w(ii,n)
      end do
      do ii=1,3
        xee(ii)=w(ii,n)*(pu*w(3,n)*bce+8.d0*pw*uw*w(3,n)*c
     x            		        +2.d0*(pw*u0(3)+uw*pp(3))*e)
        xee(ii)=xee(ii)+pp(ii)*(u0(3)*de+2.d0*uw*w(3,n)*e)
        xee(ii)=xee(ii)+u0(ii)*(pp(3)*de+2.d0*pw*w(3,n)*e)
      end do
      xee(3)=xee(3)+pu*abcde+pw*uw*bce
c  extract u0 from ee(*,i,n), mult by i*(\bStil + 2*nu*\bI), replace in u0
      do ii=1,3
        u0(ii)=z0
        do jj=1,3
          u0(ii)=u0(ii)+zq2(ii,jj)*ee(jj,i,n)
        end do
        u0(ii)=eye*u0(ii)
      end do
 1002 format(3(2g14.6,3x))
c  for znu NOT an eigenvalue, 
c  but rather the average of closely-space eigenvalues
c  in this case, zq1 is nonsingular, and we just solve for u1 
      do ii=1,3
        yr(ii)=dreal(u0(ii))
        yi(ii)=dimag(u0(ii))
        do jj=1,3
          q1r(jj,ii)=dreal(zq1(jj,ii))
          q1i(jj,ii)=dimag(zq1(jj,ii))
        end do
      end do
      call csolve(3,q1r,q1i,xr,xi,yr,yi)
      do ii=1,3
        u1(ii)=dcmplx(xr(ii),xi(ii))
      end do
c   End, different tactic
c
c  normalize  
      sum=0.d0
      do ii=1,3
        sum=sum+zabs(u1(ii))**2
      end do
      sum=dsqrt(sum)
      do ii=1,3
        u1(ii)=u1(ii)/sum
      end do
c  adf is the normalization constant
      adf=1.d0/sum
c  calculate the traction
c  and place the new stress-displacement vector in column j
c  pp is the wavenumber vector, and first two components are already in place
      pp(1)=dcmplx(px,0.d0)
      pp(2)=z0
      pp(3)=znu
      pu=z0
      pw=z0
      uw=z0
      abcde=a-b+c-2.d0*d+2.d0*e
      bce=b-4.d0*c-4.d0*e
      de=d-e
      do ii=1,3
        pu=pu+pp(ii)*u1(ii)
        pw=pw+pp(ii)*w(ii,n)
        uw=uw+u1(ii)*w(ii,n)
      end do
      do ii=1,3
        ee(ii,j,n)=u1(ii)
        ee(ii+3,j,n)=w(ii,n)*(pu*w(3,n)*bce+8.d0*pw*uw*w(3,n)*c
     x                                  +2.d0*(pw*u1(3)+uw*pp(3))*e)
        ee(ii+3,j,n)=ee(ii+3,j,n)+pp(ii)*(u1(3)*de+2.d0*uw*w(3,n)*e)
        ee(ii+3,j,n)=ee(ii+3,j,n)+u1(ii)*(pp(3)*de+2.d0*pw*w(3,n)*e)
      end do
      ee(6,j,n)=ee(6,j,n)+pu*abcde+pw*uw*bce
c  almost lastly, mult traction by i 
c  and add extra traction from (z-z0) term (not mult by i)
c  TEST - mult xee by zero, see if it is important --- it IS important
      do ii=1,3
        ee(ii+3,j,n)=eye*ee(ii+3,j,n)+adf*xee(ii)
      end do
      return
 400  print *,'eispack error'
      stop
      end
      subroutine read_ttimes
      implicit real*4(a-h,o-z)
      implicit integer*4(i-n)
      common/tt_p/d_p(200),d_pkp(200),s_p(200),s_pkp(200),n_p,n_pkp
c  read in the P- and PKP-slowness values from ttimes
c  the code interpolates across the upper-mantle triplication for DELTA=19-27
c  and chooses P or PKP based on the start time of the data-analysis window
c  for DELTA.ge.105
      do i=1,200
        d_p(i)=0.
      end do
c  fill in the head-wave slowness for DELTA.lt.14
c  13.2 sec/degree is 8.4+eps km/sec
      do i=1,13
        d_p(i)=float(i)
        s_p(i)=13.2
      end do
c      open(7,file='/park/backup/Vary/Research/Layers/ttimes_P.dat',
c     x  form='formatted')
      open(7,file='/Users/jjpark/Research/Layers/ttimes_P.dat',
     x  form='formatted')
      do i=14,200
        read(7,*,end=110) d_p(i),s_p(i)
      end do
  110 close(7)
      print *,i,' = counter at end of ttimes file'
      do i=1,200
        if(d_p(i).eq.0.) go to 120
      end do
  120 n_p=i-1
c  interpolates the triplication
c  linearly interpolate slowness between Delta=17 and p=11.06 sec/degree
c                                    and Delta=28 and p=8.92 sec/degree
      do i=18,27
        d_p(i)=float(i)
        s_p(i)=11.06-(i-17)*(11.06-8.92)/11.
      end do
      do i=45,n_p
        d_p(i-18)=d_p(i)
        s_p(i-18)=s_p(i)
      end do
      n_p=n_p-18
      do i=1,200
        d_pkp(i)=0.
      end do
c      open(7,file='/park/backup/Vary/Research/Layers/ttimes_PKP.dat',
c     x  form='formatted')
      open(7,file='/Users/jjpark/Research/Layers/ttimes_PKP.dat',
     x  form='formatted')
      do i=1,200
        read(7,*,end=130) d_pkp(i),s_pkp(i)
      end do
  130 close(7)
      print *,i,' = counter at end of ttimes file'
      do i=1,200
        if(d_pkp(i).eq.0.) go to 140
      end do
  140 n_pkp=i-1
      factor=180./(3.14159265358979*6371.)
      do i=1,n_p
        s_p(i)=s_p(i)*factor
      end do
      do i=1,n_pkp
        s_pkp(i)=s_pkp(i)*factor
      end do
      print *,'n_p,n_pkp',n_p,n_pkp
c      call plotit_axes(0.,0.,0.,0.)
c      call plotit(d_p,s_p,dum,n_p,' ',' ',' ',2,0,0,0,0)
c      call plotit(d_pkp,s_pkp,dum,n_pkp,' ',
c     x 'Epicentral Distance (degrees)','Slowness (sec/km)',2,0,0,0,1)
      return
      end