      implicit real*8(a-h,o-z)
      dimension x(1400000),am(1400000),xs(2,1400000),id(12),no(1400000)
      dimension ix(1400000),iy(1400000),idp(1400000),pm(2,1400000)
      dimension ia(1400000),ib(20,1400000),ic(1400000),title(20)
      dimension sxx(1400000),syy(1400000),iz(1400000),sxy(1400000)
      dimension ix1(1400000)
      character*80 hypodata
      data id/0,31,59,90,120,151,181,212,243,273,304,334/
c  tokoku-oki
c     n=4333
c     bmg=4.5
c  cwjap98
c     n=3007
c     bmg=4.0
c     n=7399
c     bmg=2.5
c     n=2635
c     bmg=2.5
!     n=10308
!     bmg=2.5
c     ****   read configuration file   ****
      open (7,file='etas2aniso.conf')
      read (7,112) hypodata
c     read (7,*) isw
      isw=1
      read (7,*) clsm,cutm
c     read (7,*) xxx
      read (7,*) xxx1
      xxx=xxx1
  112 format(a)
      xxmin=1000.
      yymin=1000.
      xxmax=-1000.
      yymax=-1000.
c     ****  specify output file names  ****
      open (3,file='etas2aniso.out3')
      if(isw.gt.0) then
      open (2,file='etas2aniso.out2')
      open (4,file='etas2aniso.out4')
      open (8,file='etas2aniso.out8')
      open (9,file='etas2aniso.out9')
      end if
c     *************************************
      open(1,file=hypodata)
c     open (1,file='nzgnsM00.etas')
c     open (1,file='jma2008all.etas')
      read(1,3) title
      jj=1
      nc=0
   11 continue
!     do 10 i=1,n
!     read(1,3) (ib(j,jj),j=1,20)
      read(1,3,end=12) (ib(j,jj),j=1,20)
    3 format(20a4)
      backspace 1
!     read(1,1) ino,e,an,amg,xt,ip
      read(1,*) ino,e,an,amg,xt,dep
    1 format(i5,4f12.5,i5)
      if(amg.ge.cutm-0.05) nc=nc+1
c     if(amg.lt.bmg) go to 10
      if(xxmin.gt.e) xxmin=e
      if(xxmax.lt.e) xxmax=e
      if(yymin.gt.an) yymin=an
      if(yymax.lt.an) yymax=an
      no(jj)=ino
      idp(jj)=dep
      x(jj)=xt
      xs(1,jj)=e
      xs(2,jj)=an
      am(jj)=amg
      jj=jj+1
c     if(mod(jj,1000).eq.0) write(6,*) jj
c
      go to 11
   12 continue
!  10 continue
c     write(6,*) 'no',(no(i),i=1,5)
      n=jj-1
      write(6,*) 'n',n
      do 20 i=1,n
      ix(i)=0
      iy(i)=0
      ia(i)=0
      ic(i)=0
      ix1(i)=0
   20 continue
! ordering according to magnitude & time; till 30 continue
c     do 30 i=1,n
c     bmax=0.0
c     do 40 j=1,n
c     if(iy(j).ne.0) go to 40
c     if(bmax.ge.am(j)) go to 40
c     ii=j
c     bmax=am(j)
c  40 continue
c     ix(i)=ii
c     iy(ii)=i
c     if(mod(i,1000).eq.0) write(6,*) 'i =', i
c  30 continue
      do i=1,n
      ix(i)=i
      iy(i)=i
      end do
      pi=3.141565
      i1=0
      i2=0
      i3=0
      do 50 i=1,n
      ii=ix(i)
c     if(ii.eq.23858) write(6,*) '23858',ii,ia(ii)
      if(ia(ii).ne.0) go to 50
c     cutm now read from configuration file
c     cutm=4.5
      if(am(ii).lt.clsm) go to 50
      x1=xs(1,ii)
      y1=xs(2,ii)
      if(xxx1.eq.0.0) then
      xxx=10**(0.5*am(ii)-1)
      if(xxx.lt.100.) xxx=100.
      end if
c the following 'xxx' is the cutt-off time after main shocks.
c for the cutt-off time before main shocks, see 'do 60' loop.
!     xxx=10**(0.5*am(ii)-1)
c utsu's criterion for jma data with m_j>=4.0.
c     xh=(0.015*10.**(0.5*am(ii)-2.)+0.3)/cos(2*pi*y1/360)
c     yh=(0.015*10.**(0.5*am(ii)-2.)+0.3)
c     if(xxx.lt.100.) xxx=100.
c utsu's criterion for isc m>=5.0 data
*     xh=(0.015*10.**(0.5*am(ii)-2.)+0.3)/cos(2*pi*(y1-90)/360)
*     yh=(0.015*10.**(0.5*am(ii)-2.)+0.3)
!     xh=(0.015*10.**(0.5*am(ii)-2.)+0.01)/cos(2*pi*(y1-90)/360)
      xh=(0.015*10.**(0.5*am(ii)-2.)+0.01)/cos(2*pi*(y1)/360)
      yh=(0.015*10.**(0.5*am(ii)-2.)+0.01)
*     if(xxx.lt.100.) xxx=100.
c     xxx=1.0
c     xxx=3.0/24.d0
c     xxx=1.0/24.d0
c criterion for jap univ data with m>=2.0  ---> use (meqmgcls)
c     xh= 0.015*10.**(0.5*am(ii)-2.)/cos(2*pi*y1/360)+0.001
c     yh= 0.015*10.**(0.5*am(ii)-2.)+0.001
c     if(xxx.lt.2.) xxx=1.
c     if(am(ii).gt.3.95) write(6,9) xh,yh,xxx,am(ii),x1,y1
    9 format(1h ,6f10.4)
c     ia(ii)=i
      i1=i1+1
      i2=0
c     i2=1
c     write(6,4) i1,i2,(ib(j,ii),j=1,20)
c   4 format(1h ,2i5,5x,20a4)
c
c
      if(am(ii).lt.clsm) go to 65
c
      do 60 j=1,n
c     if(ii.eq.j) go to 60
      if(ia(j).ne.0) go to 60
      if(am(ii).lt.am(j)) go to 60
c 2011.3.24
      if(am(ii).gt.am(j)+4.0) go to 60
c for the present objective, we are not intersted in taking foreshocks.
      if(x(j).lt.x(ii)) go to 60
c     if(x(j).lt.x(ii)-30.) go to 60
      if(x(j).gt.x(ii)+xxx) go to 50
      x12=abs(x1-xs(1,j))
      y12=abs(y1-xs(2,j))
      x12=min(x12,360.-x12)
      if(x12.gt.xh) go to 60
      if(y12.gt.yh) go to 60
      i2=i2+1
      if(ii.eq.j) go to 60
      ia(ii)=-i1
      ia(j)=i1
c     if(i1.eq.1547) write(6,*) '23858j',i,ii,i1,ia(j)
c     if(i1.eq.1547) write(6,9) x12,y12,xs(1,j),xs(2,j),x1,y1,x(j),x(ii
c    &)  ,xxx,am(ii)
c     if(i1.eq.1547) write(6,*) 'j',j,i2
c     if(am(ii).gt.6.95) write(6,9) x12,y12,xs(1,j),xs(2,j)
    4 format(i7,i5,2x,20a4,3i3,i6)
   60 continue
   65 continue
c     if(i2.eq.1) i1=i1-1
c     if(i2.eq.1) ia(ii)=0
      if(i2.eq.0) i1=i1-1
      if(i2.eq.0) ia(ii)=0
      write(6,*) i,ii,am(i)
   50 continue
      i1=0
      do 70 i=1,n
      if(ia(i).ge.0) go to 70
      i1=i1+1
      ic(i1)=i
      write(6,*) i1,i,ia(i)
   70 continue
c     write(6,*) (ia(i),i=1,n)
      if(isw.gt.0) write(2,7) title
    7 format('     ',20a4/)
      if(isw.gt.0) write(2,6) bmg
    6 format(' mag ',f5.1/)
      i3=0
      do 80 i=1,i1
      ii=abs(ia(ic(i)))
      i2=0
      xm=0.0
      ym=0.0
      do 90 j=1,n
      if(ii.ne.abs(ia(j))) go to 90
      i2=i2+1
      i3=i3+1
      call dat(x(j),in,im,ih)
      if(isw.gt.0) then
!     if(ic(i).eq.j) write(2,4) -i,i2,(ib(k,j),k=1,20),in,im,ih,i3
      if(ic(i).eq.j) write(2,4) -i,i2,(ib(k,j),k=1,20)
      if(ic(i).eq.j) write(8,44) (ib(k,j),k=1,20)
   44 format(20a4)
!     if(ic(i).ne.j) write(2,4) i,i2,(ib(k,j),k=1,20),in,im,ih,i3
      if(ic(i).ne.j) write(2,4) i,i2,(ib(k,j),k=1,20)
      end if
      xm=xm+xs(1,j)
      ym=ym+xs(2,j)
   90 continue
      ix1(ii)=i2
      xm=xm/i2
      ym=ym/i2
      pm(1,ii)=0.0
      pm(2,ii)=0.0
c     if(i2.ge.5) then
      pm(1,ii)=xm
      pm(2,ii)=ym
c     end if
      i2=0
      xm1=0.0
      ym1=0.0
      xm2=0.0
      ym2=0.0
      xy1=0.0
      xy2=0.0
      do 95 j=1,n
      if(ii.ne.abs(ia(j))) go to 95
      i2=i2+1
      if(i2.eq.1) xm0=xs(1,j)
      if(i2.eq.1) ym0=xs(2,j)
      xm2=xm2+(xs(1,j)-xm)**2
      ym2=ym2+(xs(2,j)-ym)**2
      xy2=xy2+(xs(1,j)-xm )*(xs(2,j)-ym )
      xm1=xm1+(xs(1,j)-xm0)**2
      ym1=ym1+(xs(2,j)-ym0)**2
      xy1=xy1+(xs(1,j)-xm0)*(xs(2,j)-ym0)
   95 continue
      sx=sqrt(xm2/(i2))
      sy=sqrt(ym2/(i2))
      sx0=sqrt(xm1/(i2))
      sy0=sqrt(ym1/(i2))
      sxy1=(xy1/(i2))
      sxy2=(xy2/(i2))
      sxy3=0.0
!     if(sx.ne.0.0.and.sy.ne.0.0) sxy3=xy2/(sx*sy)
      if(sx.ne.0.0.and.sy.ne.0.0) sxy3=xy2/sqrt(xm2)/sqrt(ym2)
      iz(ii)=0
      sxx(ii)=1.
      syy(ii)=1.
      sxy(ii)=0.
      aic0=1.0e6
      aic1=1.0e6
      aic2=1.0e6
      aic3=1.0e6
c     if(i2.ge.5) iz(ii)=3
c     go to 88
      if(sx.eq.0.0) sx=0.1d-5
      if(sx0.eq.0.0) sx0=0.1d-5
      if(sy.eq.0.0) sy=0.1d-5
      if(sy0.eq.0.0) sy0=0.1d-5
      if(i2.ge.3) then
      aic0=i2*log(sx0**2)+i2*log(sy0**2)+2
      aic1=i2*log(sx **2)+i2*log(sy **2)+6
      if(i2.gt.5) then
      aic2=i2*log((sx0*sy0)**2-sxy1**2) +6
      aic3=i2*log((sx *sy )**2-sxy2**2) +10
      end if
      aic=min(aic0,aic1,aic2,aic3)
      if(aic3.eq.aic) iz(ii)=3
      if(aic2.eq.aic) iz(ii)=2
      if(aic1.eq.aic) iz(ii)=1
      if(aic0.eq.aic) iz(ii)=0
   88 continue
      if(iz(ii).ge.2) then
      sxx(ii)=sx
      syy(ii)=sy
      xy2=xy2/(i2)
      sxy(ii)=0.0
      if(sx.ne.0.0.and.sy.ne.0.0) sxy(ii)=xy2/(sx*sy)
      end if
      end if
      if(sx.le.0.1e-5.and.sy.le.0.1e-5) then
      iz(ii)=0
      sxx(ii)=1.
      syy(ii)=1.
      sxy(ii)=0.
      end if
      if(isw.gt.0) then
      write(2,5) xm,ym,sx,sy,sxy3,iz(ii)
      write(2,66) aic0-aic0,aic1-aic0,aic2-aic0,aic3-aic0
      end if
!   5 format(/15x,5f10.4,i5/)
    5 format(15x,5f10.4,i5)
   66 format(' aic  ',4f12.1/)
   80 continue
      sx1=1.
      sy1=1.
      sxy1=0.
c     write(3,3) title
c     write(3,33) nc,xxmin,xxmax-xxmin,yymin,yymax-yymin
c  33 format(i8,4e15.7)
      write(3,33) nc,xxmin,xxmax,xxmax-xxmin,yymin,yymax,yymax-yymin
   33 format(i8,6e11.3)
      is=0
      do 190 j=1,n
      ii=abs(ia(j))
      if(ia(j).lt.0) then
      xss1=xs(1,j)
      xss2=xs(2,j)
      if(iz(ii).eq.1.or.iz(ii).eq.3) then
c     is=is+2
      xss1=pm(1,ii)
      xss2=pm(2,ii)
      end if
c     write(6,*) 'sxx',j,sxx(ii),syy(ii),sxy(ii)
      if(am(j).ge.cutm) then
      write(3,8) no(j),xss1,xss2,am(j),x(j),
     &           sxx(ii),syy(ii),sxy(ii)
      end if
      if(am(j).ge.clsm) then
      if(isw.gt.0) write(9,8) no(j),xss1,xss2,am(j),x(j),
     &           sxx(ii),syy(ii),sxy(ii)
      end if
c     if(is.ne.iss) write(4,8) no(j),xss1,xss2,am(j),x(j),
c    &           iz(ii),(ib(k,j),k=16,20),ia(j),xs(1,j),xs(2,j)
      if(iz(ii).eq.1) is=is+2
      if(iz(ii).eq.2) is=is+3
      if(iz(ii).eq.3) is=is+5
      if(isw.gt.0) then
      if(iz(ii).eq.1)
     &write(4,39)is,iz(ii),xss1,xss2,am(j),x(j),ix1(ii)
      if(iz(ii).eq.2)
     &write(4,29)is,iz(ii),sxx(ii),syy(ii),sxy(ii),am(j),x(j),ix1(ii)
      if(iz(ii).eq.3)
     &write(4,19)is,iz(ii),xss1,xss2,sxx(ii),syy(ii),sxy(ii),am(j),x(j)
     &           ,ix1(ii)
      end if
      else
      if(am(j).ge.cutm) then
      write(3,8) no(j),xs(1,j),xs(2,j),am(j),x(j),
     &           sx1,sy1,sxy1
      end if
      if(am(j).ge.clsm) then
      if(isw.gt.0) write(9,8) no(j),xs(1,j),xs(2,j),am(j),x(j),
     &           sx1,sy1,sxy1
      end if
      end if
      iss=is
  190 continue
c     write(6,*)'is',is
    8 format(i8,2f11.5,f6.2,f13.5,3f10.5)
   19 format(2i5,2f8.3,3f6.2,f4.1,f12.5,i6)
   29 format(2i5,16x,3f6.2,f4.1,f12.5,i6)
   39 format(2i5,2f8.3,18x,f4.1,f12.5,i6)
      stop
      end
      subroutine dat(x,iy,im,id)
      implicit real*8(a-h,o-z)
      dimension imm(12)
      data imm/0,31,59,90,120,151,181,212,243,273,304,334/
      iy=26+x/365.25
      do 10 i=1,12
      id1=x-(iy-26)*365.25
      if(id1.lt.imm(i)) go to 10
      im=i
   10 continue
      id=id1-imm(im)
      return
      end
