
      subroutine st_ilc
c
c ...... calculate initial location
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      character*4 stn
      character*72 id
      character*1 iwtp,iwts,iwt,polar
      real*8 pt,st,syf,sxf,psaf,ssaf
      integer*4 nsreq
      common/files/ ncr,nsp,ns1,ns2,ns3,ns4
      common/obs/pt(1000),st(1000),nsreq
      common/obs1/stn(1000),polar(1000)
      common/epid/id
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common/stn/syf(1000),sxf(1000),psaf(1000),ssaf(1000),nstns
      common/time/eqhour,isyear,isdoy,ishour,ismin,fssec
      common/wght/iwtp(1000),iwts(1000),iwt(1000)
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
c
c ...... set initial depth to 8 km
c
      edepth=8.0d0
      ez=8.0d0
      call eq_ilc(nsreq,pt,st,sxf,syf,ex,ey,ez,et,
     1            seex,seey,seez,seet)
c
      return
      end

      subroutine eq_ilc(nsreq,pt,st,sxf,syf,ex,ey,ez,et,
     1                  seex,seey,seez,seet)
c
c ...... calculate initial location
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 pt(1000),st(1000),sxf(1000),syf(1000)
      real*8 ex,ey,ez,et
c
      real*8 tsmp(1000),tpmo(1000),t0est(1000)
      real*8 xx(1000),yy(1000),dd(1000)
      real*8 kpd
      data kpd /111.195d0/
c
c ...... standard error is set to -1 if it is not computed
c
c
c ...... assume initial depth of 8 km
c
      ez=8.d0
      nsmp=0
c
c ...... check dispersion of recording station locations and 
c        phase onset times
c
      call stn_dspr(pt,st,sxf,syf,nsreq,d_dspr,t_dspr)
c
      do 1 i=1,nsreq
      if((pt(i).gt.0.d0).and.(st(i).gt.0.d0)) then
        nsmp=nsmp+1
        tsmp(nsmp)=st(i)-pt(i)
        call smp2pt(tsmp(nsmp),tpmo(nsmp))
        t0est(nsmp)=pt(i)-tpmo(nsmp)
      endif
    1 continue
      if(nsmp.lt.1) go to 10
      if(d_dspr.lt.50.d0) go to 10
      if(t_dspr.lt.1.d0) go to 10
      t0sum=0.d0
      do 2 i=1,nsmp
      t0sum=t0sum+t0est(i)
    2 continue
      et=t0sum/dble(nsmp)
c      write(*,*) '1 - st_ilc - et = ',et
c
      if(nsmp.gt.1) then
        t0var=0.d0
        do 25 i=1,nsmp
        t0var=t0var+(t0est(i)-et)*(t0est(i)-et)
   25   continue
        nnn=1
        seet=dsqrtf(nnn,t0var/dble(nsmp-1))
      else
        seet=-1.d0
      endif
      seez=-1.d0
c
      do 3 i=1,nsreq
      if(pt(i).gt.0.d0) then
        tpmot=pt(i)-et
        call pmodel(tpmot,dd(i))
      else
        tsmot=st(i)-et
        call smodel(tsmot,dd(i))
      endif
      xx(i)=sxf(i)
      yy(i)=syf(i)
    3 continue
      xxsum=0.d0
      yysum=0.d0
      do 4 i=1,nsreq
      xxsum=xxsum+xx(i)
      yysum=yysum+yy(i)
    4 continue
      x0=xxsum/dble(nsreq)
      y0=yysum/dble(nsreq)
      do 5 i=1,nsreq
      yyi=(yy(i)-y0)*kpd
      xx(i)=(xx(i)-x0)*kpd*dcosd((yy(i)+y0)/2.d0)
      yy(i)=yyi
    5 continue
      call eiloc(nsreq,xx,yy,dd,exkm,eykm,seex,seey,corxy)
c
c ...... initial location in geocentric degrees
c        origin time in seconds from beginning of day
c        depth set to 8 km
c
      ey=y0+eykm/kpd
      ex=x0+(exkm/kpd)*dcosd((2.d0*y0+eykm/kpd)/2.d0)
      ez=8.d0
      et=et
c      write(*,*) '2 - st_ilc - et = ',et
      go to 20
c
c ...... chose location of station with earliest P onset time
c        if no stations have both P- & S-wave onset times given
c
   10 continue
      imin=1
      ptmin=100000.d0
      do 11 i=1,nsreq
      if(pt(i).gt.0.d0) then
        if(pt(i).lt.ptmin) then
          imin=i
          ptmin=pt(i)
        endif
      endif
   11 continue
      ex=sxf(imin)
      ey=syf(imin)
      ez=8.d0
      et=pt(imin)-2.d0
c      write(*,*) '3 - st_ilc - et = ',et
      seex=-1.d0
      seey=-1.d0
      seez=-1.d0
      seet=-1.d0
      corxy=0.d0
c
   20 continue
c      write(*,*) 'st_ilc - ex, ey, ez & et = ',ex,ey,ez,et
c      
      return
      end

      subroutine stn_dspr(pt,st,sxf,syf,nsreq,d_dspr,t_dspr)
c
c ...... calculate station median distance and onset time
c
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      real*8 pt(1000),st(1000),sxf(1000),syf(1000)
      real*8 stn_dist(10000),stn_time(10000)
c
      if(nsreq.gt.30) then
        d_dspr=100.d0
        t_dspr=10.d0
        return
      endif
      k=0
      do 2 i=1,nsreq-1
      do 1 j=i+1,nsreq
      yy=111.195d0*(syf(j)-syf(i))
      xx=111.195d0*dcosd((syf(i)+syf(j))/2.d0)*(sxf(j)-sxf(i))
      arg=xx*xx+yy*yy
      if(arg.gt.0.d0) then
        k=k+1
        stn_dist(k)=dsqrt(arg)
        if(pt(i).gt.0.d0) then
          timei=pt(i)
        else
          timei=st(i)
        endif
        if(pt(j).gt.0.d0) then
          timej=pt(j)
        else
          timej=st(j)
        endif
        d_time=timej-timei
        if(d_time.lt.-600.d0) then
          stn_time(k)=d_time+86400.d0
        elseif(d_time.lt.600.d0) then
          stn_time(k)=d_time
        else
          stn_time(k)=d_time-86400.d0
        endif
      endif
    1 continue
    2 continue
c
      do 4 i=1,k-1
      do 3 j=i+1,k
      if(stn_time(j).lt.stn_time(i)) then
        t_temp=stn_time(i)
        stn_time(i)=stn_time(j)
        stn_time(j)=t_temp
      endif
      if(stn_dist(j).lt.stn_dist(i)) then
        d_temp=stn_dist(i)
        stn_dist(i)=stn_dist(j)
        stn_dist(j)=d_temp
      endif
    3 continue
    4 continue
	if(k/2.ge.1) then
	  d_dspr=stn_dist(k/2)
	  t_dspr=stn_time(k/2)
	else
	  d_dspr=stn_dist(1)
	  t_dspr=stn_time(1)
	endif
c
c      do 5 i=1,k
c      write(*,*) 'stn_dspr - i, stn_dist & stn_time = ',
c     1           i,stn_dist(i),stn_time(i)
c    5 continue
c      write(*,*) 'stn_dspr - d_dspr & t_dspr = ',d_dspr,t_dspr
c
      return
      end

      subroutine eiloc(n,x,y,d,ex,ey,seex,seey,corxy)
c
c ...... calculate initial epicenter
c
c        finds epicenter when given station locations and inferred
c        distance to epicenter
c
c            1 - computes distance from travel time using the input model
c
c            2 - linearize n equations by squaring and subtraction
c                  (x1-x0)^2 + (y1-y0)^2 =   d1^2
c                - (x2-x0)^2 - (y2-y0)^2 = - d2^2
c                --------------------------------
c                 2*(x2-x1)*x0 + 2*(y2-y1)*y0 = (x2^2-x1^2) + (y2^2-y1^2)
c                                               - (d2^2-d1^2)
c
c            3 - solve resulting n-1 equations by least-squares
c
      implicit real*8 (a-h,o-z)
      real*8 x(1000),y(1000),d(1000)
      real*8 x2(1000),y2(1000),d2(1000)
      real*8 a(1000,2),bb(1000),ata(2,2),atb(2)
      real*8 atai(2,2)
      real*8 ex,ey,seex,seey,corxy
c
      nm1=n-1
      do 2 i=1,nm1
      ip1=i+1
      do 1 j=ip1,n
      if(d(j).lt.d(i)) then
        dt=d(i)
        xt=x(i)
        yt=y(i)
        d(i)=d(j)
        x(i)=x(j)
        y(i)=y(j)
        d(j)=dt
        x(j)=xt
        y(j)=yt
      endif
    1 continue
    2 continue
      x0=0.d0
      y0=0.d0
      do 3 i=1,n
      x0=x0+x(i)
      y0=y0+y(i)
    3 continue
      x0=x0/dble(n)
      y0=y0/dble(n)
      do 4 i=1,n
      x(i)=x(i)-x0
      y(i)=y(i)-y0
      x2(i)=x(i)*x(i)
      y2(i)=y(i)*y(i)
      d2(i)=d(i)*d(i)
    4 continue
      do 5 i=1,nm1
      j=i+1
      a(i,1)=2.d0*(x(j)-x(i))
      a(i,2)=2.d0*(y(j)-y(i))
      bb(i)=(d2(i)-d2(j))-(x2(i)-x2(j))-(y2(i)-y2(j))
    5 continue
      do 9 i=1,2
      do 7 k=1,2
      sum=0.d0
      do 6 j=1,nm1
      sum=sum+a(j,i)*a(j,k)
    6 continue
      ata(i,k)=sum
    7 continue
      sum=0.d0
      do 8 j=1,nm1
      sum=sum+a(j,i)*bb(j)
    8 continue
      atb(i)=sum
    9 continue
      det=ata(1,1)*ata(2,2)-ata(1,2)*ata(2,1)
      atai(1,1)=ata(2,2)/det
      atai(1,2)=-ata(1,2)/det
      atai(2,1)=-ata(2,1)/det
      atai(2,2)=ata(1,1)/det
      delx=atai(1,1)*atb(1)+atai(1,2)*atb(2)
      dely=atai(2,1)*atb(1)+atai(2,2)*atb(2)
      ex=x0+delx
      ey=y0+dely
c
      if(nm1.lt.3) then
        seex=-1.d0
        seey=-1.d0
        covxy=0.d0
      else
        varb=0.d0
        do 11 i=1,nm1
        eta=bb(i)-(a(i,1)*delx+a(i,2)*dely)
        varb=varb+eta*eta
   11   continue
        varb=varb/dble(nm1-2)
        varx=atai(1,1)*varb
        nnn=2
        seex=dsqrtf(nnn,varx)
        vary=atai(2,2)*varb
        nnn=3
        seey=dsqrtf(nnn,vary)
        nnn=4
        corxy=atai(1,2)*varb/dsqrtf(nnn,varx*vary)
      endif
c
      return
      end

      subroutine pmodel(tpmo,delta)
c
c ...... calculate distance from P-wave travel time
c
      implicit real*8 (a-h,o-z)
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      real*8 tpmo,delta
      real*8 cd(21),ct(21)
      data va,vb,vc,vt/0.d0,0.d0,0.d0,0.d0/
      data cd /  0.d0,  1.d0,  2.d0,  5.d0, 10.d0, 20.d0, 40.d0,
     1          60.d0, 80.d0,100.d0,150.d0,200.d0,250.d0,300.d0,
     2         400.d0,500.d0,600.d0,700.d0,800.d0,1000.d0,5000.d0/
      if(vpa.ne.va) then
        va=vpa
        vb=vpb
        vc=vpc
        vt=vd
        err=1.e-3
        e0=0.d0
        do 1 i=1,21
        call ptcalc(va,vb,vc,vt,e0,ez,err,cd(i),p,ct(i),h,a,dtz)
    1   continue
      endif
      do 2 j=1,21
      if(ct(j).ge.tpmo) then
        if(j.eq.1) then
          delta=0.d0
          go to 3
        else
          delta=cd(j-1)+((cd(j)-cd(j-1))/
     1         (ct(j)-ct(j-1)))*(tpmo-ct(j-1))
          go to 3
        endif
      else
        if(j.eq.21) then
          delta=cd(j)+((cd(j)-cd(j-1))/(ct(j)-ct(j-1)))*(tpmo-ct(j))
        endif
      endif
    2 continue
    3 continue
c
      return
      end

      subroutine smodel(tsmo,delta)
c
c ...... calculate distance from S-wave travel time
c
      implicit real*8 (a-h,o-z)
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      real*8 tsmo,delta
      real*8 cd(21),ct(21)
      data va,vb,vc,vt/0.d0,0.d0,0.d0,0.d0/
      data cd /  0.d0,  1.d0,  2.d0,  5.d0, 10.d0, 20.d0, 40.d0,
     1          60.d0, 80.d0,100.d0,150.d0,200.d0,250.d0,300.d0,
     2         400.d0,500.d0,600.d0,700.d0,800.d0,1000.d0,5000.d0/
      if(vsa.ne.va) then
        va=vsa
        vb=vsb
        vc=vsc
        vt=vd
        err=1.e-3
        e0=0.d0
        do 1 i=1,21
        call ptcalc(va,vb,vc,vt,e0,ez,err,cd(i),p,ct(i),h,a,dtz)
    1   continue
      endif
      do 2 j=1,21
      if(ct(j).ge.tsmo) then
        if(j.eq.1) then
          delta=0.d0
          go to 3
        else
          delta=cd(j-1)+((cd(j)-cd(j-1))/
     1         (ct(j)-ct(j-1)))*(tsmo-ct(j-1))
          go to 3
        endif
      else
        if(j.eq.21) then
          delta=cd(j)+((cd(j)-cd(j-1))/(ct(j)-ct(j-1)))*(tsmo-ct(j))
        endif
      endif
    2 continue
    3 continue
c
      return
      end

      subroutine smp2pt(tsmp,tpmo)
c
c ...... calculate P-wave travel time from S-P observed time
c
      implicit real*8 (a-h,o-z)
      common /epi/elatd,elongd,edepth,eth,etm,ets,ex,ey,ez,et
      common /model/vpa,vpb,vpc,vsa,vsb,vsc,vd
      real*8 tsmp,tpmo
      real*8 cd(21),csmp(21),cpmo(21)
      data vap,vbp,vcp,vas,vbs,vcs,vt/7*0.d0/
      data cd /  0.d0,  1.d0,  2.d0,  5.d0, 10.d0, 20.d0, 40.d0,
     1          60.d0, 80.d0,100.d0,150.d0,200.d0,250.d0,300.d0,
     2         400.d0,500.d0,600.d0,700.d0,800.d0,1000.d0,5000.d0/
      if(vpa.ne.va) then
        vap=vpa
        vbp=vpb
        vcp=vpc
        vas=vsa
        vbs=vsb
        vcs=vsc
        vt=vd
        err=1.e-3
        e0=0.d0
        do 1 i=1,21
        call ptcalc(vap,vbp,vcp,vt,e0,ez,err,cd(i),p,ptime,h,a,dtz)
        cpmo(i)=ptime
        call ptcalc(vas,vbs,vcs,vt,e0,ez,err,cd(i),p,stime,h,a,dtz)
        csmp(i)=stime-ptime
    1   continue
      endif
      do 2 j=1,21
      if(csmp(j).ge.tsmp) then
        if(j.eq.1) then
          tpmo=0.d0
          go to 3
        else
          tpmo=cpmo(j-1)+((cpmo(j)-cpmo(j-1))/
     1        (csmp(j)-csmp(j-1)))*(tsmp-csmp(j-1))
          go to 3
        endif
      else
        if(j.eq.21) then
          tpmo=cpmo(j)+((cpmo(j)-cpmo(j-1))/
     1        (csmp(j)-csmp(j-1)))*(tpmo-cpmo(j))
        endif
      endif
    2 continue
    3 continue
c
      return
      end
