      program etaspa
c-----------------------------------------------------------------------
c  space-time version of the etas model.
c  copied from #ogatake.blgmagnt.fort; func2 only.
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /kkxy/kkx,kky
      common /trend/kkt
      common     / ccc /  isw , ipr,igrd
      common     / ddd /  r , f , aic , sd
      dimension r(31,31)
      common /sv/bx(4096),by(4096),index,mbx,mby
      index=0
c     call xclock(time1,3)
c     call clock(time1)
c     time1=dclock()
      call redata(xx,yy,zz,zmg,tx,ty,tz,nn,nfunct,mx,my,a0)
      call dav(nfunct)
c     call xclock(time2,4)
c     call clock(time2)
c     time2=dclock()
c     write(6,*) 'time',time,'sec'
c     write(6,*) 'time',time2-time1,'sec'
      stop
      end
c***********************************************************************
      subroutine redata(xx,yy,zz,zmg,tx,ty,tz,nn,nfunct,mx,my,a0)
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /jfun/jxy(npa,6000),jnn(npa)
      dimension xx(1),yy(1),zz(1),zmg(1)
      dimension a0(1)
      character*8 fmt(10),ftl(6)
      character*80 hypodata,splcoef
      common /approx/ xmg1
      common /range/ tstar,nstar
      open(1,file='st-etas.conf')
c     open(1,file='init.data')
c     open(1,file='etaspac.init')
c     open(1,file='inland.init')
c     open(1,file='tohok95.init')
c     open(1,file='alljp95.init')
      read(1,112) hypodata   !./etas2aniso.out3
  112 format(a)
c     read(1,2) fmt
c     write(6,2) fmt
      write(6,112) hypodata
      read(1,1) nfunct    !nfunct = 7
      write(6,*) nfunct
    1 format(i10)
c     read(1,2) fmt
c     write(6,2) fmt
    2 format(10a8)
      read(1,*) tx,ty,tz,nn
      write(6,3) tx,ty,tz,nn
    3 format(3f10.0,i10)
      read(1,*) xmin,ymin,xmg0,zmin,xmg1,tstar,bi2
      write(6,8) xmin,ymin,xmg0,zmin,xmg1,tstar,bi2
    8 format(8f10.1)
      in=10
      open(in,file=hypodata)
      read(in,2) ftl
      write(6,7) ftl
    7 format(1h ,'data set  ',10a8)
      nnc=0
      do 20 i=1,nn
      if(nfunct.ge.15) then
c     read(in,fmt,end=25) xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
      read(in,*,end=25) i1,xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
c     dd=1.d0
c     ee=1.d0
c     ff=0.d0
      aa(i)=ee/dd/sqrt(1-ff**2)
      bb(i)=ff   /sqrt(1-ff**2)
      cc(i)=dd/ee/sqrt(1-ff**2)
c     write(6,11) i,aa(i),bb(i),cc(i),aa(i)*cc(i)-bb(i)**2
   11 format(i5,4f8.4)
      else
c     read(in,fmt,end=25) xx(i),yy(i),zmg(i),zz(i)
      read(in,*,end=25) i1,xx(i),yy(i),zmg(i),zz(i)
c     write(6,*) i,i1,xx(i),yy(i),zmg(i),zz(i)
      end if
   20 continue
      i=nn+1
   25 nn=i-1
      write(6,4) in
    4 format(1h ,'input device',i10)
c setection of time interval (xmg1,xmg2) in iday from jmacatlg 1926-90.5
      ii=0
      nstar=0
      write(6,*) 'nn=',nn
      do 10 i=1,nn
        if(zmg(i).lt.xmg0) go to 30
c       if(xx(i).gt.tx+xmin) go to 10
c       if(xx(i).lt.   xmin) go to 10
c       if(yy(i).gt.ty+ymin) go to 10
c       if(yy(i).lt.   ymin) go to 10
        if(yy(i).eq.   ymin) write(6,*) i,yy(i)
        if(yy(i).eq.   ymin) go to 10
        if(zz(i).gt.tz+zmin) go to 30
        if(zz(i).lt.   zmin) go to 30
c       write(6,44) i,xx(i),yy(i),zz(i),zmg(i)
   44   format(i8,2f11.5,f15.5,f8.2)
        ii=ii+1
c       xx(ii)=xx(i)-xmin
c       yy(ii)=yy(i)-ymin
        xx(ii)=xx(i)
        yy(ii)=yy(i)
        zmg(ii)=zmg(i)-xmg0
        zz(ii)=zz(i)-zmin
        aa(ii)=aa(i)
        bb(ii)=bb(i)
        cc(ii)=cc(i)
c       if(aa(ii).ne.1.d0.and.nfunct.ge.15)
c    &       write(6,*) ii,aa(ii),bb(ii),cc(ii)
        if(zz(ii).gt.tstar) nnc=nnc+1
        if(zz(ii).lt.tstar) nstar=nstar+1
c transformation from degree to km
c       xx(ii)=xx(ii)*cos(ymin*3.14159/180.)
        xx(ii)=xx(ii)*cos((ymin+ty/2)*3.14159/180.)
c       xx(ii)=xx(ii)*111.11
c       yy(ii)=yy(ii)*111.11
        go to 10
   30   continue
        write(6,44) i,xx(i),yy(i),zz(i),zmg(i)
   10 continue
c transformation from degree to km
      tx=tx*cos((ymin+ty/2)*3.14159/180.)
      tz=tz-zmin
c     tx=tx*111.11
c     ty=ty*111.11
c
      write(6,*) 'nfunct=',nfunct
      nn=ii
      do i=1,nn-1
        if(zz(i+1).le.zz(i)) then
          write(6,*) 'Data error! interval-length =< 0',
     &                i,zz(i+1)-zz(i),zz(i),zz(i+1)
        end if
      end do
      write(6,5) tx,ty,tz,xmin,ymin,xmg0,xmg1,zmin,tstart,nn,nnc
    5 format(1h0,'tx,ty,tz,xmin,ymin,xmg0,xmg1,zmin,tsta'/9f10.3,
     &/' nn =',i5,' nnc =',i5)
c     stop 119
      xmg1=xmg1-xmg0
c     write(6,6) ((i,xx(i),yy(i),zmg(i),zz(i),zz(i)-zz(i-1)),i=1,nn)
    6 format(///1h ,i5,5f13.5/(1h ,i5,5f13.5))
c
c     if(nfunct.ge.15) read(1,112) splcoef
c     if(nfunct.ge.15) then
c     open (2, file=splcoef)
c
c     read(2,*) mx,my
      mx=1
      my=1
c     read(2,9) dum1,dum2,dum3
c     read(2,9) (a0(i),i=1,(mx+3)*(my+3))
      do i=1,(mx+3)*(my+3)
        a0(i)=0.0
      end do
c
    9 format(4f20.0)
c     write(6,*) 'mx,my',mx,my
c     end if
      do i=1,nn
      do j=1,6000
      jxy(i,j)=0
      end do
      jnn(i)=0
      end do
      jmax=0
c     bi2=2.0
      write(6,*) 'bi2',bi2
      do 500 i=2,nn
      do 510 j=1,i-1
      xxx=max(10**(0.5*(zmg(j)+xmg0)-1),100.d0)*bi2
      yh=(0.015*10.**(0.5*(zmg(j)+xmg0)-2.)+0.3)*bi2
c     if(zz(i)-zz(j).gt.xxx) go to 510
      if(abs(xx(i)-xx(j)).gt.yh) go to 510
      if(abs(yy(i)-yy(j)).gt.yh) go to 510
      jnn(i)=jnn(i)+1
      ji=jnn(i)
      if(ji.gt.6000) then
      write(6,*) 'jmax over',ji
      stop 119
      end if
      jxy(i,ji)=j
  510 continue
      ji=jnn(i)
      if(jmax.lt.ji) jmax=ji
c     write(6,511) i,jnn(i),(jxy(i,j),j=1,ji)
  511 format(10i6/(10x,9i6))
  500 continue
      write(6,*) 'jmax',jmax
      write(6,*) 'tstar,nstar', tstar,nstar
c     stop
c
      return
      end
c***********************************************************************
      subroutine dav(nfunct)
      implicit real * 8 (a-h,o-z)
      external func5,func6,func7,func8,func15,func16,func17,func18
      external func19
      common /ddd/r,fff,aic,sd
      common /kkxy/kkx,kky
      common /trend/kkt
      common     / ccc /  isw , ipr,igrd       
      dimension x(50),r(31,31),xo(50)
   10 read(1,*) n
      write(6,1020) n
      read(1,*) (x(i),i=1,n)
      write(6,1030)  (x(i),i=1,n)
      read(1,*) ipr,igrd
      write(6,*) 'linear ipr,igrd',ipr,igrd
c     write(6,*) 'HERE'
c     stop 1020
      do 70 i=1,n
   70 x(i)=sqrt(x(i))
      do 30 ii=1,1
c     do 30 ii=1,2
c     do 30 ii=1,5
      if(nfunct.eq.5) call davidn(x,n,0,func5)
      if(nfunct.eq.6) call davidn(x,n,0,func6)
      if(nfunct.eq.7) call davidn(x,n,0,func7)
      if(nfunct.eq.8) call davidn(x,n,0,func8)
      if(nfunct.eq.15) call davidn(x,n,0,func15)
      if(nfunct.eq.16) call davidn(x,n,0,func16)
      if(nfunct.eq.17) call davidn(x,n,0,func17)
      if(nfunct.eq.18) call davidn(x,n,0,func18)
      if(nfunct.eq.19) call davidn(x,n,0,func19)
c     call davidn(x,n,0,func5)
   30 continue
      do 80 i=1,n
   80 x(i)=x(i)**2
      write(6,1040) (x(i),i=1,n)
      return
 1000 format(3i10,2f15.6)
 1010 format(8f10.4)
 1020 format(1h0,3x,'input data'/1h ,5x,'n=',i3,3x,'itr=',i5,3x,
     2          'ier=',i3,3x,'eps=',1pe16.7)
 1030 format(1h0,'x=',5e14.5)
 1040 format(
     3      /' mle = ',5e12.5/('       ',9e12.5))
 1050 format(4d20.13)
 1060 format(e25.15)
 1070 format(1h0,'  c = ',e25.15)
 1080 format(1h0,'          ',10f12.5)
 1090 format(1h ,'          ',6d20.10)
 1100 format(1h ,'          ',6d20.10)
      end
c***********************************************************************
      subroutine linear( x,h,ram,ee,k,ig,funct )
c
c     this subroutine performs the linear search along the direction spe
c     by the vector h
c       ----------------------------------------------------------------
c       the following subroutine is directly called by this subroutine:
c             funct
c       ----------------------------------------------------------------
c
c     inputs:
c        x:       vector of position
c        h:       search direction
c        k:       dimension of vector x
c
c     outputs:
c        ram:     optimal step width
c        e2:      minimum function value
c        ig:      error code
c
      implicit  real  *8 ( a-h,o-z )
      parameter(nda=22222,npa=22222)
      integer  return,sub
      dimension  x(npa) , h(npa) , x1(npa)
      dimension  g(npa)
      common     / ccc /  isw , ipr
      external funct
c
c     ipr = 0
c     ipr = 7
      if( ram .le. 1.0d-30 )  ram = 0.01d0
      const2 = 1.0d-60
      hnorm = 0.d0
      do 10  i=1,k
   10 hnorm = hnorm + h(i)**2
      hnorm = sqrt( hnorm )
c
      ram2 = ram
      e1 =ee
      ram1 = 0.d0
c
      do 20  i=1,k
   20 x1(i) = x(i) + ram2*h(i)
      call  funct( k,x1,e2,g,ig )
      if(ipr.ge.7)  write(6,2)  ram2,e2
c
      if( ig .eq. 1 )  go to  50
      if( e2 .gt. e1 )  go to 50
   30 ram3 = ram2*2.d0
      do 40  i=1,k
   40 x1(i) = x(i) + ram3*h(i)
      call  funct( k,x1,e3,g,ig )
      if( ig.eq.1 )  go to  500
      if( ipr.ge.7 )  write(6,3)  ram3,e3
      if( e3 .gt. e2 )  go to 70
      ram1 = ram2
      ram2 = ram3
      e1 = e2
      e2 = e3
      go to 30
c
   50 ram3 = ram2
      e3 = e2
      ram2 = ram3*0.1d0
      if( ram2*hnorm .lt. const2 )  go to  400
      do 60  i=1,k
   60 x1(i) = x(i) + ram2*h(i)
      call  funct( k,x1,e2,g,ig )
      if(ipr.ge.7)  write(6,4)  ram2,e2
c     write(6,8) e1,e2,e3,ee
    8 format(1h ,5d20.10)
      if( e2.gt.e1 )  go to 50
c
   70 assign 80 to return
      go to 200
c
   80 do 90  i=1,k
   90 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,ee,g,ig )
      if(ipr.ge.7)  write(6,5)  ram,ee
c
      ifg = 0
      assign  300 to  sub
      assign 200 to sub
   95 assign 130 to return
      if( ram .gt. ram2 )  go to 110
      if( ee .ge. e2 )  go to 100
      ram3 = ram2
      ram2 = ram
      e3 =e2
      e2 =ee
      go to  sub,( 200,300 )
c
  100 ram1 = ram
      e1 = ee
      go to  sub,( 200,300 )
c
  110 if( ee .le. e2 )  go to 120
      ram3 = ram
      e3 = ee
      go to  sub,( 200,300 )
c
  120 ram1 = ram2
      ram2 = ram
      e1 = e2
      e2 = ee
      go to  sub,( 200,300 )
c
  130 do 140  i=1,k
  140 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,ee,g,ig )
      if( ipr.ge.7 )  write(6,6)  ram,ee
      assign 200 to sub
      ifg = ifg+1
      ifg = 0
      if( ifg .eq. 1 )  go to 95
c
      if( e2 .lt. ee )  ram = ram2
      if( e2 .lt. ee )  ee = e2
      return
c
c      -------  internal subroutine sub1  -------
  200 a1 = (ram3-ram2)*e1
      a2 = (ram1-ram3)*e2
      a3 = (ram2-ram1)*e3
      b2 = (a1+a2+a3)*2.d0
      b1 = a1*(ram3+ram2) + a2*(ram1+ram3) + a3*(ram2+ram1)
      if( b2 .eq. 0.d0 )  go to 210
      ram = b1 /b2
      go to return ,( 80,130 )
c
  210 ig = 1
      ram = ram2
      return
c
c      -------  internal subroutine sub2  -------
c
  300 if( ram3-ram2 .gt. ram2-ram1 )  go to 310
      ram = (ram1+ram2)*0.5d0
      go to return ,( 80,130 )
c
  310 ram = (ram2+ram3)*0.5d0
      go to return ,( 80,130 )
c ------------------------------------------------------------
c
  400 ram = 0.d0
      return
c ------------------------------------------------------------
c
  500 ram = (ram2+ram3)*0.5d0
  510 do 520  i=1,k
  520 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,e3,g,ig )
      if( ipr.ge.7 )  write(6,7)  ram,e3
      if( ig.eq.1 )  go to 540
      if( e3.gt.e2 )  go to 530
      ram1 = ram2
      ram2 = ram
      e1 = e2
      e2 = e3
      go to 500
c
  530 ram3 = ram
      go to 70
c
  540 ram = (ram2+ram)*0.5d0
      go to 510
c
c ------------------------------------------------------------
    1 format( 1h ,'lambd1 =',d18.10, 10x,'e1 =',d25.17 )
    2 format( 1h ,'lambd2 =',d18.10, 10x,'e2 =',d25.17 )
    3 format( 1h ,'lambd3 =',d18.10, 10x,'e3 =',d25.17 )
    4 format( 1h ,'lambd4 =',d18.10, 10x,'e4 =',d25.17 )
    5 format( 1h ,'lambd5 =',d18.10, 10x,'e5 =',d25.17 )
    6 format( 1h ,'lambd6 =',d18.10, 10x,'e6 =',d25.17 )
    7 format( 1h ,'lambd7 =',d18.10, 10x,'e7 =',d25.17 )
      end
c***********************************************************************
      subroutine  davidn( x,n,ihes,funct )
c
c          minimization by davidon-fletcher-powell procedure
c
c       ----------------------------------------------------------------
c       the following subroutines are directly called by this subroutine
c             funct
c             hesian
c             linear
c       ----------------------------------------------------------------
c          inputs:
c             x:       vector of initial values
c             k:       dimension of the vector x
c             ihes:    =0   inverse of hessian matrix is not available
c                      =1   inverse of hessian matrix is available
c
c          output:
c             x:       vector of minimizing solution
c
      implicit  real * 8  ( a-h , o-z )
      external funct
      dimension  x(50) , dx(50) , g(50) , g0(50) , y(50)
      dimension  h(50,50) , wrk(50) , s(50)
      dimension  r(31,31)
      common     / ccc /  isw
      common     / ddd /  r , f , aic , sd
      data  tau1 , tau2  /  1.0d-5 , 1.0d-5  /
      data  eps1 , eps2  / 1.0d-5 , 1.0d-5  /
      iter=0
      ramda = 0.5d0
      const1 = 1.0d-70
c
c          initial estimate of inverse of hessian
c
      do  20   i=1,n
      do  10   j=1,n
   10 h(i,j) = 0.0d00
      s(i) = 0.0d00
      dx(i) = 0.0d00
   20 h(i,i) = 1.0d00
      isw = 0
c
      call  funct( n,x,xm,g,ig )
c
      write( 6,340 )     xm , 2*(xm+n)
c     write( 9,340 )     xm , 2*(xm+n)
c
c          inverse of hessian computation (if available)
c
      if( ihes .eq. 1 )   call  hesian( x,n,h )
c
      icc = 0
c      iteration
11110 continue
      icc = icc + 1
      do  11111   ic=1,n
      if( ic .eq. 1 .and. icc .eq. 1 )     go to 120
c
      do  40   i=1,n
   40 y(i) = g(i) - g0(i)
      do  60   i=1,n
      sum = 0.0d00
      do  50   j=1,n
   50 sum = sum + y(j) * h(i,j)
   60 wrk(i) = sum
      s1 = 0.0d00
      s2 = 0.0d00
      do  70   i=1,n
      s1 = s1 + wrk(i) * y(i)
   70 s2 = s2 + dx(i) * y(i)
      if( s1.le.const1 .or. s2.le.const1 )  go to 900
      if( s1 .le. s2 )     go to 100
c
c          update the inverse of hessian matrix
c
c               ---  davidon-fletcher-powell type correction  ---
c
      do  90   i=1,n
      do  90   j=i,n
      h(i,j) = h(i,j) + dx(i)*dx(j)/s2 - wrk(i)*wrk(j)/s1
   90 h(j,i) = h(i,j)
      go to  120
c
c               ---  fletcher type correction  ---
c
  100 continue
      stem = s1 / s2 + 1.0d00
      do  110   i=1,n
      do  110   j=i,n
      h(i,j) = h(i,j)- (dx(i)*wrk(j)+wrk(i)*dx(j)-dx(i)*dx(j)*stem)/s2
  110 h(j,i) = h(i,j)
c
c
c
  120 continue
      ss = 0.0d00
      do  150   i=1,n
      sum = 0.0d00
      do  140   j=1,n
  140 sum = sum + h(i,j)*g(j)
      ss = ss + sum * sum
  150 s(i) = -sum
c
c
      s1 = 0.0d00
      s2 = 0.0d00
      do  170   i=1,n
      s1 = s1 + s(i)*g(i)
  170 s2 = s2 + g(i)*g(i)
      ds2 = dsqrt(s2)
      gtem = dabs(s1) / ds2
c     write(6,610)gtem,ds2
      if( gtem .le. tau1  .and.  ds2 .le. tau2 )     go to  900
      if( s1 .lt. 0.0d00 )     go to  200
      do  190   i=1,n
      do  180   j=1,n
  180 h(i,j) = 0.0d00
      h(i,i) = 1.0d00
  190 s(i) = -s(i)
  200 continue
c
      ed = xm
c
c          linear  search
c
      call  linear( x,s,ramda,ed,n,ig,funct )
c
      write( 6,330 )     ramda , ed , s1 , s2
c     write( 9,330 )     ramda , ed , s1 , s2
c
      iter=iter+1
c     if(mod(iter,10).eq.0) call hes4(n,x,h)
c     if(mod(iter,10).eq.0) call invdet(h,tdet,n,50)
c
      s1 = 0.0d00
      do  210   i=1,n
      dx(i) = s(i) * ramda
      s1 = s1 + dx(i) * dx(i)
      g0(i) = g(i)
  210 x(i) = x(i) + dx(i)
      xmb = xm
      isw = 0
c
      call  funct( n,x,xm,g,ig )
c
      s2 = 0.d0
      do  220     i=1,n
  220 s2 = s2 + g(i)*g(i)
      if( dsqrt(s2) .gt. tau2 )   go to  11111
      if( xmb/xm-1.d0 .lt. eps1  .and.  dsqrt(s1) .lt. eps2 )  go to 900
11111 continue
      if( icc .ge. 5 )     go to 900
      go to 11110
  900 continue
c
      write( 6,340 )     xm , 2*(xm+n)
c     write( 9,340 )     xm , 2*(xm+n)
c
      write( 6,600 )
      write( 6,610 )     (x(i),i=1,n)
      write( 6,601 )
      write( 6,610 )     (g(i),i=1,n)
c     write( 9,600 )
c     write( 9,610 )     (x(i),i=1,n)
c     write( 9,601 )
c     write( 9,610 )     (g(i),i=1,n)
c@@@
c     call hes4(n,x,h)
c     do 653 i=1,n
c 653 write(6,604) (h(i,j),j=1,n)
c     call invdet(h,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (h(i,j),j=1,n)
  602 format(1h0,'*** estimated inverse hessian ***')
  604 format(1h ,8d15.5/(1h ,8d15.5))
c
      return
  330 format( 1h ,' lmbd = ',d14.7,2x,'-ll = ',d22.15,2x,d9.2,x,d9.2)
  340 format( 1h ,'- log likelihood =',d23.15,5x,'aic =',f10.1)
  600 format( 1h0,'-----  x  -----' )
  601 format( 1h0,'***  gradient  ***' )
  610 format( 1h ,10d13.5 )
      end
c***********************************************************************
      subroutine invdet(x,tdet,mm,mj)
c-----------------------------------------------------------------------
c       the inverse and determinant of x computation
c----------------------------------------------------------------------
c       inputs:
c          x:     mm*mm square matrix
c          mm:    dimension of x
c          mj:    absolute dimension of x in the main program
c
c       outputs:
c          x:     inverse of x
c          xdet:  determinant of x
c
      implicit  real * 8  ( a-h , o-z )
      parameter(ldat=1999, npa=234, mpa=npa*7)
      dimension x(mj,mj)
      dimension  ids(mpa)
      xdet = 1.0d00
      tdet=0.0d0
      gtdet=0.0d0
      do 10 l=1,mm
c     pivoting at l-th stage
      xmaxp=0.10000d-10
      maxi=0
      do 110 i=l,mm
    1 if( abs(xmaxp) .ge. abs(x(i,l)) )     go to 110
      xmaxp=x(i,l)
      maxi=i
  110 continue
      ids(l)=maxi
      if(maxi.eq.l) go to 120
      if(maxi.gt.0) go to 121
      xdet = 0.0d00
      go to 140
c     row interchange
  121 do 14 j=1,mm
      xc=x(maxi,j)
      x(maxi,j)=x(l,j)
   14 x(l,j)=xc
      xdet=-xdet
  120 tdet=tdet+log(abs(xmaxp))
c 120 xdet=xdet*xmaxp
c 120 continue
      xc = 1.0d00 / xmaxp
      x(l,l)=1.0d00
      do 11 j=1,mm
   11 x(l,j)=x(l,j)*xc
      do 12 i=1,mm
      if(i.eq.l) go to 12
      xc=x(i,l)
      x(i,l) = 0.0d00
      do 13 j=1,mm
   13 x(i,j)=x(i,j)-xc*x(l,j)
   12 continue
   10 continue
      if(mm.gt.1) go to 123
      go to 140
c     column interchange
  123 mm1=mm-1
      do 130 j=1,mm1
      mmj=mm-j
      jj=ids(mmj)
      if(jj.eq.mmj) go to 130
      do 131 i=1,mm
      xc=x(i,jj)
      x(i,jj)=x(i,mmj)
  131 x(i,mmj)=xc
  130 continue
  140 return
      end
c***********************************************************************
      subroutine hesian(x,n,h)
      implicit real * 8 (a-h,o-z)
      return
      end
c***********************************************************************
      subroutine func5(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized).
c     the space-time version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c               *\exp{ ( (x-x_j)^2+(y-y_j)^2 )/2/\exp{2\alpha m_j} }
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
      a4=b(4)**2
      a5=b(5)**2
      a6=b(6)**2
c     write(6,1) a1,a2,a3,a4,a5,a6
      if(a4.gt.9.99) go to 50
      if(a5.gt.3.d0) go to 50
      do 207 k=1,kkt
  207 at(k)=b(k+5)
      do 208 k=1,kkt
  208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
      do 420 k=1,kkt
  420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
      ramdai=a1+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
      d1ff=1.d0/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
c
      do 390 k=1,kkt
      ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      do 10 i=2,nn
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      do 20 j=1,i-1
      space=((xx(i)-xx(j))**2+(yy(i)-yy(j))**2)/2/exp(a4*zmg(j))/a6
      d4spa=-zmg(j)*space
      d6spa=-space/a6
      rmdi=rmdi+exp(-space)/(zz(i)-zz(j)+a3)**a5
      rmd3i=rmd3i+exp(-space)/(zz(i)-zz(j)+a3)**(a5+1)
      rmd4i=rmd4i-exp(-space)*d4spa/(zz(i)-zz(j)+a3)**a5
      rmd5i=rmd5i-exp(-space)/(zz(i)-zz(j)+a3)**a5*log(zz(i)-zz(j)+a3)
      rmd6i=rmd6i-exp(-space)*d6spa/(zz(i)-zz(j)+a3)**a5
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
      do 410 k=1,kkt
  410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
      ramdai=a1+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
      d1ff=d1ff+1.d0/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
c
      do 380 k=1,kkt
      ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
  380 continue
c
   10 continue
c
      ft=0.0
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      area=tx*ty
      do 70 i=1,nn
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rr1=xx(i)**2+( by(jy)  -yy(i) )**2
      rr2=xx(i)**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rr1=( bx(jx)  -xx(i) )**2+( ty-yy(i) )**2
      rr2=( bx(jx+1)-xx(i) )**2+( ty-yy(i) )**2
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rr1=( tx-xx(i) )**2+( by(jy)  -yy(i) )**2
      rr2=( tx-xx(i) )**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rr1=( bx(jx)  -xx(i) )**2+yy(i)**2
      rr2=( bx(jx+1)-xx(i) )**2+yy(i)**2
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rr1+rr2)/2
      
      exi=exp(-a4*zmg(i))
      spac2=2*pi*a6*( 1.d0-exp(-rr*exi/2/a6) )
      zexi=-zmg(i)*exp(-a4*zmg(i))
      d4sp= 2*pi*a6*( -exp(-rr*exi/2/a6) )*(-rr*zexi/2/a6 )
      d6sp= 2*pi*( 1.d0-exp(-rr*exi/2/a6) )
     &     +2*pi*a6 *( -exp(-rr*exi/2/a6)*rr*exi/2/a6**2 )
      ft = ft  +1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &       *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d3ft=d3ft+( (tz-zz(i)+a3)**(-a5)-a3**(-a5) )
     &         *exp(a4*zmg(i))*spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *(  zmg(i)*exp( a4*zmg(i) )*spac2+exp( a4*zmg(i) )*d4sp  )
     &     *delta/(2*pi)
      d5ft=d5ft + 1.d0/(1-a5)**2*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &                         *exp( a4*zmg(i) )*spac2*delta/(2*pi)
     &     + 1.d0/(1-a5)*( -(tz-zz(i)+a3)**(1-a5)*log(tz-zz(i)+a3)
     &      +a3**(1-a5)*log(a3) )*exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &       *exp( a4*zmg(i) )*d6sp*delta/(2*pi)
   80 continue
   70 continue
      area=tx*ty
      fs=a1*tz*area+a2*ft
      d1fs=tz     *area
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=a2*d5ft
      d6fs=a2*d6ft
c
  170 continue
c
      stsum=0.0
      do 360 k=1,kkt
  360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
c
c
      do 350 k=1,kkt
  350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
      h(4)=-g4*2.0d0*b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
c     h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
      do 267 k=1,kkt
  267 h(k+5)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,f
c     write(6,1)(h(i),i=1,n)
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func15(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized).
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c *\exp{ ( (x-x_j)^2+(y-y_j)^2 )/2/\exp{2\alpham_j} }
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      dimension rnx(4),rny(4)
      common /approx/ xmg1
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
      
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
      a4=b(4)**2
      a5=b(5)**2
      a6=b(6)**2
c     write(6,1) a1,a2,a3,a4,a5,a6
      if(a5.gt.3.d0) go to 50
c
c
      hx=tx/mx
      hy=ty/my
      ip=int(xx(1)/hx)+1
      jp=int(yy(1)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(1)/hx-(ip-1)
      ry=yy(1)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 45 i=0,3
      do 45 j=0,3
        ik=(ip+i)+(mx+3)*(jp+j-1)
        sfk0=sfk0+a0(ik)*rnx(4-i)*rny(4-j)
   45 continue
c
c
      do 207 k=1,kkt
  207 at(k)=b(k+5)
      do 208 k=1,kkt
  208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
      do 420 k=1,kkt
  420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
c     ramdai=a1+rmdti1
      ramdai=a1*exp(sfk0)+rmdti
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
c     d1ff=1.d0/ramdai
      d1ff=exp(sfk0)/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
c
      do 390 k=1,kkt
      ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      do 10 i=2,nn
c
      ip=int(xx(i)/hx)+1
      jp=int(yy(i)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(i)/hx-(ip-1)
      ry=yy(i)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 41 k=0,3
      do 41 l=0,3
        ik=(ip+k)+(mx+3)*(jp+l-1)
        sfk0=sfk0+a0(ik)*rnx(4-k)*rny(4-l)
   41 continue
c
c
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      do 20 j=1,i-1
      space=( aa(j)*(xx(i)-xx(j))**2-2*bb(j)*(xx(i)-xx(j))*(yy(i)-yy(j))
     &       +cc(j)*(yy(i)-yy(j))**2 ) /2/exp(a4*zmg(j))/a6
      d4spa=-zmg(j)*space
      d6spa=-space/a6
      rmdi=rmdi+exp(-space)/(zz(i)-zz(j)+a3)**a5
      rmd3i=rmd3i+exp(-space)/(zz(i)-zz(j)+a3)**(a5+1)
      rmd4i=rmd4i-exp(-space)*d4spa/(zz(i)-zz(j)+a3)**a5
      rmd5i=rmd5i-exp(-space)/(zz(i)-zz(j)+a3)**a5*log(zz(i)-zz(j)+a3)
      rmd6i=rmd6i-exp(-space)*d6spa/(zz(i)-zz(j)+a3)**a5
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
      do 410 k=1,kkt
  410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
c     ramdai=a1+a2*rmdi+rmdti1
      ramdai=a1*exp(sfk0)+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
c     d1ff=d1ff+1.d0/ramdai
      d1ff=d1ff+exp(sfk0)/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
c
      do 380 k=1,kkt
      ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
  380 continue
c
   10 continue
c
      ft=0.0
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
c
      area=tx*ty
      do 70 i=1,nn
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rrd1=q( aa(i), bb(i), cc(i), xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), xx(i), by(jy+1)-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy)-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), ty-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), ty-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), ty-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), ty-yy(i) )
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rrd1=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy+1)-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy)-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), yy(i) )
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rrd1+rrd2)/2
c
      exi=exp(-a4*zmg(i))
      spac2=2*pi*a6*( 1.d0-exp(-rr*exi/2/a6) )
      zexi=-zmg(i)*exp(-a4*zmg(i))
      d4sp= 2*pi*a6*( -exp(-rr*exi/2/a6) )*(-rr*zexi/2/a6 )
      d6sp= 2*pi*( 1.d0-exp(-rr*exi/2/a6) )
     &     +2*pi*a6 *( -exp(-rr*exi/2/a6)*rr*exi/2/a6**2 )
      ft = ft  +1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &       *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d3ft=d3ft+( (tz-zz(i)+a3)**(-a5)-a3**(-a5) )
     &         *exp(a4*zmg(i))*spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *(  zmg(i)*exp( a4*zmg(i) )*spac2+exp( a4*zmg(i) )*d4sp  )
     &     *delta/(2*pi)
      d5ft=d5ft + 1.d0/(1-a5)**2*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &                         *exp( a4*zmg(i) )*spac2*delta/(2*pi)
     &     + 1.d0/(1-a5)*( -(tz-zz(i)+a3)**(1-a5)*log(tz-zz(i)+a3)
     &      +a3**(1-a5)*log(a3) )*exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &       *exp( a4*zmg(i) )*d6sp*delta/(2*pi)
   80 continue
   70 continue
c
      f2=0.0
      nox=101
      noy=101
      dx=tx/(nox-1)
      dy=ty/(noy-1)
      do 150 kx=1,nox
      do 150 ky=1,noy
      xk=tx*(kx-1)/(nox-1)
      yk=ty*(ky-1)/(noy-1)
      ip=int(xk/hx)+1
      jp=int(yk/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xk/hx-(ip-1)
      ry=yk/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 120 j=0,3
      do 120 i=0,3
      ij=(ip+i)+(mx+3)*(jp+j-1)
  120 sfk0=sfk0+a0(ij)*rnx(4-i)*rny(4-j)
c................................................
      if(sfk0.gt.170.0) go to 50
      ram=exp(sfk0)
      rmd=ram
      if(kx.eq.1 .or. kx.eq.nox) rmd=ram/2
      if(ky.eq.1 .or. ky.eq.noy) rmd=ram/2
      if(kx.eq.1 .and. ky.eq.1) rmd=ram/3
      if(kx.eq.nox .and. ky.eq.noy) rmd=ram/3
      if(kx.eq.1 .and. ky.eq.noy) rmd=ram/6
      if(kx.eq.nox .and. ky.eq.1) rmd=ram/6
      f2=f2+rmd
  150 continue
      volume=f2*dx*dy
c     volume=area
c
c
      area=tx*ty
c     fs=a1*tz*area+a2*ft
c     d1fs= tz*area
      fs=a1*tz*volume+a2*ft
      d1fs=tz*volume
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=a2*d5ft
      d6fs=a2*d6ft
c
  170 continue
c
      stsum=0.0
      do 360 k=1,kkt
  360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
c
c
      do 350 k=1,kkt
  350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
      h(4)=-g4*2.0d0*b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
c     h(1)=0.0
c     h(2)=0.0
c     if(a3.le..1d-4) h(3)=0.0
c     if(a3.le..1d-4) b(3)=0.1d-3
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
      do 267 k=1,kkt
  267 h(k+5)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,f
c     write(6,1)(h(i),i=1,n)
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func6(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time version.
c \lambda(t,x,y)=\mu+\sum k\exp(\alpham_j)/(t-t_j+c)^p
c                         /( (x-x_j)^2+(y-y_j)^2+d )^q
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
      common /approx/ xmg1
c     dimension bx(4096),by(4096)
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
      a4=b(4)**2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
      if(a4.gt.3.d1) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
c     do 207 k=1,kkt
c 207 at(k)=b(k+7)
c     do 208 k=1,kkt
c 208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
c     do 420 k=1,kkt
c 420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
      ramdai=a1+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
      d1ff=1.d0/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
c
c     do 390 k=1,kkt
c     ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      do 10 i=2,nn
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      do 20 j=1,i-1
c@@@4 if(zmg(j).lt.xmg1) go to 20
      space=((xx(i)-xx(j))**2+(yy(i)-yy(j))**2+a6)**(-a7)
      d6spa=-a7*((xx(i)-xx(j))**2+(yy(i)-yy(j))**2+a6)**(-a7-1)
      d7spa=-space*log((xx(i)-xx(j))**2+(yy(i)-yy(j))**2+a6)
      rmdi=rmdi+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+zmg(j)*exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5*space
      rmd5i=rmd5i-exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5
     &                              *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d7spa
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
c     do 410 k=1,kkt
c 410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
      ramdai=a1+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
      d1ff=d1ff+1.d0/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
c
c     do 380 k=1,kkt
c     ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
c 380 continue
c
   10 continue
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      ft=0.0
      do 70 i=1,nn
c@@@4 if(zmg(i).lt.xmg1) go to 70
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rr1=xx(i)**2+( by(jy)  -yy(i) )**2
      rr2=xx(i)**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rr1=( bx(jx)  -xx(i) )**2+( ty-yy(i) )**2
      rr2=( bx(jx+1)-xx(i) )**2+( ty-yy(i) )**2
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rr1=( tx-xx(i) )**2+( by(jy)  -yy(i) )**2
      rr2=( tx-xx(i) )**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rr1=( bx(jx)  -xx(i) )**2+yy(i)**2
      rr2=( bx(jx+1)-xx(i) )**2+yy(i)**2
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rr1+rr2)/2
c
      spac2=pi/(1-a7)*( (rr+a6)**(1-a7)-a6**(1-a7) )
      ft=ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d3ft=d3ft+( (tz-zz(i)+a3)**(-a5)-a3**(-a5) )
     &         *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1-a5)
     &*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &*zmg(i)*exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d5ft=d5ft + 1.d0/(1-a5)**2
     &           *( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &           *exp( a4*zmg(i) )*spac2*delta/(2*pi)
     &+ 1.d0/(1-a5)*( -( tz-zz(i)+a3 )**(1-a5)*log( tz-zz(i)+a3 )
     &                +a3**(1-a5)*log(a3) )
     &             *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d6fta= pi*( (rr+a6)**( -a7)-a6**( -a7) )
      d6ft=d6ft+1.d0/(1-a5)*((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
     &     *exp(a4*zmg(i))*d6fta*delta/(2*pi)
      d7fta= pi/(1-a7)**2*( (rr+a6)**(1-a7)-a6**(1-a7) )
     &      +pi/(1-a7)   *(-(rr+a6)**(1-a7)*log(rr+a6)
     &                         +a6 **(1-a7)*log(   a6) )
      d7ft=d7ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *exp(a4*zmg(i))*d7fta*delta/(2*pi)
   80 continue
   70 continue
      area=tx*ty
      fs=a1*tz*area+a2*ft
      d1fs=tz     *area
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*              d5ft
      d6fs=a2*     d6ft
      d7fs=a2*     d7ft
      stsum=0.0
c     do 360 k=1,kkt
c 360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
c
c
c     do 350 k=1,kkt
c 350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
      h(4)=-g4*2.0d0*b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
c     h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
c     do 267 k=1,kkt
c 267 h(k+7)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func16(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum k\exp(\alpham_j)/(t-t_j+c)^p
c                         /( (x-x_j)^2+(y-y_j)^2+d )^q
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
      common /approx/ xmg1
c     dimension bx(4096),by(4096)
      dimension rnx(4),rny(4)
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
      
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
      a4=b(4)**2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
c
c
      hx=tx/mx
      hy=ty/my
      ip=int(xx(1)/hx)+1
      jp=int(yy(1)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(1)/hx-(ip-1)
      ry=yy(1)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 45 i=0,3
      do 45 j=0,3
        ik=(ip+i)+(mx+3)*(jp+j-1)
        sfk0=sfk0+a0(ik)*rnx(4-i)*rny(4-j)
   45 continue
c
c
c     do 207 k=1,kkt
c 207 at(k)=b(k+7)
c     do 208 k=1,kkt
c 208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
c     do 420 k=1,kkt
c 420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
c     ramdai=a1+rmdti1
      ramdai=a1*exp(sfk0)+rmdti
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
c     d1ff=1.d0/ramdai
      d1ff=exp(sfk0)/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
c
c     do 390 k=1,kkt
c     ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      do 10 i=2,nn
c
      ip=int(xx(i)/hx)+1
      jp=int(yy(i)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(i)/hx-(ip-1)
      ry=yy(i)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 41 k=0,3
      do 41 l=0,3
        ik=(ip+k)+(mx+3)*(jp+l-1)
        sfk0=sfk0+a0(ik)*rnx(4-k)*rny(4-l)
   41 continue
c
c
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      do 20 j=1,i-1
c@@@4 if(zmg(j).lt.xmg1) go to 20
      space=( aa(j)*(xx(i)-xx(j))**2-2*bb(j)*(xx(i)-xx(j))*(yy(i)-yy(j))
     &       +cc(j)*(yy(i)-yy(j))**2+a6)**(-a7)
      d6spa=-a7*(
     &        aa(j)*(xx(i)-xx(j))**2-2*bb(j)*(xx(i)-xx(j))*(yy(i)-yy(j))
     &       +cc(j)*(yy(i)-yy(j))**2+a6 )**(-a7-1)
      d7spa=-space*log(
     &        aa(j)*(xx(i)-xx(j))**2-2*bb(j)*(xx(i)-xx(j))*(yy(i)-yy(j))
     &       +cc(j)*(yy(i)-yy(j))**2+a6 )
      rmdi=rmdi+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+zmg(j)*exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5*space
      rmd5i=rmd5i-exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5
     &                              *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+exp(a4*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d7spa
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
c     do 410 k=1,kkt
c 410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
c     ramdai=a1+a2*rmdi+rmdti1
      ramdai=a1*exp(sfk0)+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
c     d1ff=d1ff+1.d0/ramdai
      d1ff=d1ff+exp(sfk0)/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
c
c     do 380 k=1,kkt
c     ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
c 380 continue
c
   10 continue
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      ft=0.0
c
      do 70 i=1,nn
c@@@4 if(zmg(i).lt.xmg1) go to 70
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rrd1=q( aa(i), bb(i), cc(i), xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), xx(i), by(jy+1)-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy)-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), ty-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), ty-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), ty-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), ty-yy(i) )
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rrd1=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy+1)-yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy)-yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), yy(i) )
      rr1 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), yy(i) )
      rr2 =q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), yy(i) )
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rrd1+rrd2)/2
c
      spac2=pi/(1-a7)*( (rr+a6)**(1-a7)-a6**(1-a7) )
      ft=ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d3ft=d3ft+( (tz-zz(i)+a3)**(-a5)-a3**(-a5) )
     &         *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1-a5)
     &*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &*zmg(i)*exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d5ft=d5ft + 1.d0/(1-a5)**2
     &           *( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &           *exp( a4*zmg(i) )*spac2*delta/(2*pi)
     &+ 1.d0/(1-a5)*( -( tz-zz(i)+a3 )**(1-a5)*log( tz-zz(i)+a3 )
     &                +a3**(1-a5)*log(a3) )
     &             *exp( a4*zmg(i) )*spac2*delta/(2*pi)
      d6fta= pi*( (rr+a6)**( -a7)-a6**( -a7) )
      d6ft=d6ft+1.d0/(1-a5)*((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
     &     *exp(a4*zmg(i))*d6fta*delta/(2*pi)
      d7fta= pi/(1-a7)**2*( (rr+a6)**(1-a7)-a6**(1-a7) )
     &      +pi/(1-a7)   *(-(rr+a6)**(1-a7)*log(rr+a6)
     &                         +a6 **(1-a7)*log(   a6) )
      d7ft=d7ft+1.d0/(1-a5)*( (tz-zz(i)+a3)**(1-a5)-a3**(1-a5) )
     &     *exp(a4*zmg(i))*d7fta*delta/(2*pi)
   80 continue
   70 continue
c
      f2=0.0
      nox=101
      noy=101
      dx=tx/(nox-1)
      dy=ty/(noy-1)
      do 150 kx=1,nox
      do 150 ky=1,noy
      xk=tx*(kx-1)/(nox-1)
      yk=ty*(ky-1)/(noy-1)
      ip=int(xk/hx)+1
      jp=int(yk/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xk/hx-(ip-1)
      ry=yk/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 120 j=0,3
      do 120 i=0,3
      ij=(ip+i)+(mx+3)*(jp+j-1)
  120 sfk0=sfk0+a0(ij)*rnx(4-i)*rny(4-j)
c................................................
      if(sfk0.gt.170.0) go to 50
      ram=exp(sfk0)
      rmd=ram
      if(kx.eq.1 .or. kx.eq.nox) rmd=ram/2
      if(ky.eq.1 .or. ky.eq.noy) rmd=ram/2
      if(kx.eq.1 .and. ky.eq.1) rmd=ram/3
      if(kx.eq.nox .and. ky.eq.noy) rmd=ram/3
      if(kx.eq.1 .and. ky.eq.noy) rmd=ram/6
      if(kx.eq.nox .and. ky.eq.1) rmd=ram/6
      f2=f2+rmd
  150 continue
      volume=f2*dx*dy
c     volume=area
c
c
      area=tx*ty
c     fs=a1*tz*area+a2*ft
c     d1fs=tz     *area
      fs=a1*tz*volume+a2*ft
      d1fs=tz*volume
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*d5ft
      d6fs=a2*     d6ft
      d7fs=a2*     d7ft
      stsum=0.0
c     do 360 k=1,kkt
c 360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
c
c     do 350 k=1,kkt
c 350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
      h(4)=-g4*2.0d0*b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
c     h(1)=0.0
c     h(2)=0.0
c     if(a3.le..1d-4) h(3)=0.0
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
c     do 267 k=1,kkt
c 267 h(k+7)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func7(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c boundary condition neatly considered.
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common /range/ tstar, nstar
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      common /approx/ xmg1
      common /jfun/jxy(npa,6000),jnn(npa)
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
c     a4=b(4)**2
      a4=b(4)**2/2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
      if(a4.gt.5.d0) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
c     do 207 k=1,kkt
c 207 at(k)=b(k+7)
c     do 208 k=1,kkt
c 208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
c     do 420 k=1,kkt
c 420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
      ramdai=a1+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
      d1ff=1.d0/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
c
c     do 390 k=1,kkt
c     ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      ic=0
      do 10 i=2,nn
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      ji=jnn(i)
      do 20 j1=1,ji
c     do 20 j=1,i-1
c
c     xxx=max(10**(0.5*(zmg(j)+4.5)-1),100.d0)*2
c     yh=(0.015*10.**(0.5*(zmg(j)+4.5)-2.)+0.3)*2
c     if(zz(i)-zz(j).gt.xxx) go to 25
c     if(abs(xx(i)-xx(j)).gt.yh) go to 25
c     if(abs(yy(i)-yy(j)).gt.yh) go to 25
c
c@@@6 if(zmg(j).lt.xmg1) go to 20
c     if(zmg(j)+4.5.lt. 5.0) go to 20
      j=jxy(i,j1)
      xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
      dxxij4=(xx(i)-xx(j))*exp(-a4*zmg(j))*(-zmg(j))
      yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
      dyyij4=(yy(i)-yy(j))*exp(-a4*zmg(j))*(-zmg(j))
      space=(xxij**2+yyij**2+a6)**(-a7)
      d6spa=-a7*(xxij**2+yyij**2+a6)**(-a7-1)
      d4spa=d6spa*(2*xxij*dxxij4+2*yyij*dyyij4)
      d7spa=-space*log(xxij**2+yyij**2+a6)
      rmdi=rmdi+1.d0/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+1.d0/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+1.d0       /(zz(i)-zz(j)+a3)**a5*d4spa
      rmd5i=rmd5i-1.d0       /(zz(i)-zz(j)+a3)**a5
     &                    *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+1.d0/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+1.d0/(zz(i)-zz(j)+a3)**a5       *d7spa
      ic=ic+1
      go to 20
   25 continue
c     write(6,*) i,j,xxx,yh,zmg(j)
c     write(6,*) zz(i)-zz(j),abs(xx(i)-xx(j)),abs(yy(i)-yy(j))
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
c     do 410 k=1,kkt
c 410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
chistory
            if(zz(i).lt.tstar) go to 10
c
      ramdai=a1+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
      d1ff=d1ff+1.d0/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
c
c     do 380 k=1,kkt
c     ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
c 380 continue
c
   10 continue
c     write(6,*) ic, nn*(nn-1)/2
c     stop 'I am tired.'
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      ft=0.0
      area=tx*ty
      do 70 i=1,nn
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rr1=xx(i)**2+( by(jy)  -yy(i) )**2
      rr2=xx(i)**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rr1=( bx(jx)  -xx(i) )**2+( ty-yy(i) )**2
      rr2=( bx(jx+1)-xx(i) )**2+( ty-yy(i) )**2
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rr1=( tx-xx(i) )**2+( by(jy)  -yy(i) )**2
      rr2=( tx-xx(i) )**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rr1=( bx(jx)  -xx(i) )**2+yy(i)**2
      rr2=( bx(jx+1)-xx(i) )**2+yy(i)**2
      rr3=(tx/mbx)**2
      end if
c
c     if(rr1.eq.0.0 .or. rr2.eq.0.0) go to 80
c     if( rr1+rr2-rr3.gt.2*sqrt(rr1)*sqrt(rr2) ) go to 80
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rr1+rr2)/2
      
      scl=exp(2*a4*zmg(i))
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
      spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
     &                                  +a6 **(1-a7)*log(       a6) )
c
      if(a5.ne.1.d0) then
      ft=ft+1.d0/(1.d0-a5)
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))     *spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)
     &                 -a3**(1.d0-a5))               * d4s2*delta/(2*pi)
      d5ft=d5ft+1.d0/(1.d0-a5)**2
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
     &    +1.d0/(1.d0-a5)*(-(tz-zz(i)+a3)**(1.d0-a5)
     &      *log(tz-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                                               * d6s2*delta/(2*pi)
      d7ft=d7ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                                               * d7s2*delta/(2*pi)
      else
      ft=ft+(log(tz-zz(i)+a3)-log(a3))*spac2*delta/(2*pi)
      d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))     *spac2*delta/(2*pi)
      d4ft=d4ft+(log(tz-zz(i)+a3)-log(a3))           * d4s2*delta/(2*pi)
      d6ft=d6ft+(log(tz-zz(i)+a3)-log(a3))           * d6s2*delta/(2*pi)
      d7ft=d7ft+(log(tz-zz(i)+a3)-log(a3))           * d7s2*delta/(2*pi)
      endif
   80 continue
   70 continue
      d3ftstar=0.0
      d4ftstar=0.0
      d5ftstar=0.0
      d6ftstar=0.0
      d7ftstar=0.0
      ftstar=0.0
      do 170 i=1,nstar
cn
      do 180 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 180
      jy=j
      rr1=xx(i)**2+( by(jy)  -yy(i) )**2
      rr2=xx(i)**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 180
      jx=j-mby
      rr1=( bx(jx)  -xx(i) )**2+( ty-yy(i) )**2
      rr2=( bx(jx+1)-xx(i) )**2+( ty-yy(i) )**2
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 180
      jy=j-mbx-mby
      rr1=( tx-xx(i) )**2+( by(jy)  -yy(i) )**2
      rr2=( tx-xx(i) )**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 180
      jx=j-mbx-2*mby
      rr1=( bx(jx)  -xx(i) )**2+yy(i)**2
      rr2=( bx(jx+1)-xx(i) )**2+yy(i)**2
      rr3=(tx/mbx)**2
      end if
c
c     if(rr1.eq.0.0 .or. rr2.eq.0.0) go to 80
c     if( rr1+rr2-rr3.gt.2*sqrt(rr1)*sqrt(rr2) ) go to 80
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rr1+rr2)/2

      scl=exp(2*a4*zmg(i))
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
      spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
     &                                  +a6 **(1-a7)*log(       a6) )
c
      if(a5.ne.1.d0) then
      ftstar=ftstar+1.d0/(1.d0-a5)
     &  *((tstar-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      d3ftstar=d3ftstar+((tstar-zz(i)+a3)**(-a5)-a3**(-a5))   
     &  *spac2*delta/(2*pi)
      d4ftstar=d4ftstar+1.d0/(1.d0-a5)*((tstar-zz(i)+a3)**(1.d0-a5)
     &                -a3**(1.d0-a5))               * d4s2*delta/(2*pi)
      d5ftstar=d5ftstar+1.d0/(1.d0-a5)**2
     &  *((tstar-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
     &    +1.d0/(1.d0-a5)*(-(tstar-zz(i)+a3)**(1.d0-a5)
     &  *log(tstar-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d6ftstar=d6ftstar+1.d0/(1.d0-a5)*((tstar-zz(i)+a3)
     & **(1.d0-a5)-a3**(1.d0-a5))* d6s2*delta/(2*pi)
      d7ftstar=d7ftstar+1.d0/(1.d0-a5)*((tstar-zz(i)+a3)
     &  **(1.d0-a5)-a3**(1.d0-a5)) * d7s2*delta/(2*pi)
      else
      ftstar=ftstar+(log(tstar-zz(i)+a3)-log(a3))*spac2*delta/(2*pi)
      d3ftstar=d3ftstar+((tstar-zz(i)+a3)**(-a5)-a3**(-a5))
     &    *spac2*delta/(2*pi)
      d4ftstar=d4ftstar+(log(tstar-zz(i)+a3)-log(a3))*d4s2*delta/(2*pi)
      d6ftstar=d6ftstar+(log(tstar-zz(i)+a3)-log(a3))*d6s2*delta/(2*pi)
      d7ftstar=d7ftstar+(log(tstar-zz(i)+a3)-log(a3))*d7s2*delta/(2*pi)
      endif
  180 continue
  170 continue
c     write(6,*) 'FT',ft,tstar,nstar
      fs=a1*(tz-tstar)*area+a2*(ft-ftstar)
      d1fs=(tz-tstar)*area
      d2fs=ft-ftstar
      d3fs=a2*(d3ft-d3ftstar)
      d4fs=a2*(d4ft-d4ftstar)
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*(d5ft-d5ftstar)
      d6fs=a2*(d6ft-d6ftstar)
      d7fs=a2*(d7ft-d7ftstar)
c
      stsum=0.0
c     do 360 k=1,kkt
c 360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
c
c     do 350 k=1,kkt
c 350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
c     h(4)=-g4*2.0d0*b(4)
      h(4)=-g4      *b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
      if(igrd.eq.0) h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
c     do 267 k=1,kkt
c 267 h(k+7)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
c     write(6,*) f,ff,fs
c     stop 'Hi, its me.'
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func8(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c boundary condition neatly considered.
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      common /approx/ xmg1
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
c     a4=b(4)**2
      a4=b(4)**2/2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
      a8=b(8)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
      if(a4.gt.5.d0) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
      if(a8.gt.5.d0) go to 50
c     do 207 k=1,kkt
c 207 at(k)=b(k+7)
c     do 208 k=1,kkt
c 208 ggt(k)=0.0
      rmdti1=0.0
      if(kkt.eq.0) go to 83
c     do 420 k=1,kkt
c 420 rmdti1=rmdti1+at(k)*(zz(1)/tz)**k
   83 continue
      ramdai=a1+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
      d1ff=1.d0/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
      d8ff=0.0
c
c     do 390 k=1,kkt
c     ggt(k)=ggt(k)+(zz(1)/tz)**k/ramdai
  390 continue
c
      do 10 i=2,nn
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      rmd8i=0.0
      do 20 j=1,i-1
c@@@6 if(zmg(j).lt.xmg1) go to 20
      xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
      dxxij4=(xx(i)-xx(j))*exp(-a4*zmg(j))*(-zmg(j))
      yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
      dyyij4=(yy(i)-yy(j))*exp(-a4*zmg(j))*(-zmg(j))
      space=(xxij**2+yyij**2+a6)**(-a7)
      d6spa=-a7*(xxij**2+yyij**2+a6)**(-a7-1)
      d4spa=d6spa*(2*xxij*dxxij4+2*yyij*dyyij4)
      d7spa=-space*log(xxij**2+yyij**2+a6)
      rmdi=rmdi+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5*d4spa
      rmd5i=rmd5i-exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5
     &                    *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d7spa
      rmd8i=rmd8i+zmg(j)*exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5*space
   20 continue
c
      rmdti1=0.0
      if(kkt.eq.0) go to 84
c     do 410 k=1,kkt
c 410 rmdti1=rmdti1+at(k)*(zz(i)/tz)**k
   84 continue
c     write(6,*) rmdi
      ramdai=a1+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
      d1ff=d1ff+1.d0/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
      d8ff=d8ff+a2*rmd8i/ramdai
c
c     do 380 k=1,kkt
c     ggt(k)=ggt(k)+(zz(i)/tz)**k/ramdai
c 380 continue
c
   10 continue
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      d8ft=0.0
      ft=0.0
      area=tx*ty
      do 70 i=1,nn
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rr1=xx(i)**2+( by(jy)  -yy(i) )**2
      rr2=xx(i)**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rr1=( bx(jx)  -xx(i) )**2+( ty-yy(i) )**2
      rr2=( bx(jx+1)-xx(i) )**2+( ty-yy(i) )**2
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rr1=( tx-xx(i) )**2+( by(jy)  -yy(i) )**2
      rr2=( tx-xx(i) )**2+( by(jy+1)-yy(i) )**2
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rr1=( bx(jx)  -xx(i) )**2+yy(i)**2
      rr2=( bx(jx+1)-xx(i) )**2+yy(i)**2
      rr3=(tx/mbx)**2
      end if
c
c     if(rr1.eq.0.0 .or. rr2.eq.0.0) go to 80
c     if( rr1+rr2-rr3.gt.2*sqrt(rr1)*sqrt(rr2) ) go to 80
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rr1+rr2)/2
      
      scl=exp(2*a4*zmg(i))
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
      spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
     &                                  +a6 **(1-a7)*log(       a6) )
c
      ft=ft+1.d0/(1.d0-a5)*exp(a8*zmg(i))
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))*spac2*delta/(2*pi)
     &          *exp(a8*zmg(i))
      d4ft=d4ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &          *exp(a8*zmg(i)) * d4s2*delta/(2*pi)
      d5ft=d5ft+1.d0/(1.d0-a5)**2 *exp(a8*zmg(i))
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
     &     +1.d0/(1.d0-a5) *exp(a8*zmg(i)) *(-(tz-zz(i)+a3)**(1.d0-a5)
     &      *log(tz-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                        * exp(a8*zmg(i)) * d6s2*delta/(2*pi)
      d7ft=d7ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &          * exp(a8*zmg(i)) * d7s2*delta/(2*pi)
      d8ft=d8ft+1.d0/(1.d0-a5)*exp(a8*zmg(i))*zmg(i)
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
   80 continue
   70 continue
      fs=a1*tz*area+a2*ft
      d1fs=tz*area
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*d5ft
      d6fs=a2*d6ft
      d7fs=a2*d7ft
      d8fs=a2*d8ft
c
      stsum=0.0
c     do 360 k=1,kkt
c 360 stsum=stsum+at(k)*tz*(tz/tz)**(k+1)/(k+1)
c
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
      g8=d8ff-d8fs
c
c     do 350 k=1,kkt
c 350 gt(k)=ggt(k)-tz*(tz/tz)**(k+1)/(k+1)
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
c     h(4)=-g4*2.0d0*b(4)
      h(4)=-g4      *b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
      h(8)=-g8*2.0d0*b(8)
c     h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
c     do 267 k=1,kkt
c 267 h(k+7)=-gt(k)
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func17(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      common /range/ tstart, ntstar
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      dimension rnx(4),rny(4)
      common /approx/ xmg1
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
      
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
c     a4=b(4)**2
      a4=b(4)**2/2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
      if(a4.gt.5.d0) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
c
c
      hx=tx/mx
      hy=ty/my
      ip=int(xx(1)/hx)+1
      jp=int(yy(1)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(1)/hx-(ip-1)
      ry=yy(1)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 45 i=0,3
      do 45 j=0,3
        ik=(ip+i)+(mx+3)*(jp+j-1)
        sfk0=sfk0+a0(ik)*rnx(4-i)*rny(4-j)
   45 continue
c
c
      rmdti1=0.0
c     ramdai=a1+rmdti1
      ramdai=a1*exp(sfk0)+rmdti
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
c     d1ff=1.d0/ramdai
      d1ff=exp(sfk0)/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
c
      do 10 i=2,nn
c
      ip=int(xx(i)/hx)+1
      jp=int(yy(i)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(i)/hx-(ip-1)
      ry=yy(i)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 41 k=0,3
      do 41 l=0,3
        ik=(ip+k)+(mx+3)*(jp+l-1)
        sfk0=sfk0+a0(ik)*rnx(4-k)*rny(4-l)
   41 continue
c
c
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      do 20 j=1,i-1
c@@@6 if(zmg(j).lt.xmg1) go to 20
      xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
      dxxij4=(xx(i)-xx(j))*exp(-a4*zmg(j))*(-zmg(j))
      yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
      dyyij4=(yy(i)-yy(j))*exp(-a4*zmg(j))*(-zmg(j))
      space=(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7)
      d6spa=
     &-a7*(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7-1)
      d4spa=d6spa*( 2*aa(j)*xxij*dxxij4-
     &  4*bb(j)*xxij*yyij*(-zmg(j))+2*cc(j)*yyij*dyyij4 )
      d7spa=-space*
     &  log(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)
      rmdi=rmdi+1.d0/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+1.d0/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+1.d0       /(zz(i)-zz(j)+a3)**a5*d4spa
      rmd5i=rmd5i-1.d0       /(zz(i)-zz(j)+a3)**a5
     &                    *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+1.d0/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+1.d0/(zz(i)-zz(j)+a3)**a5       *d7spa
   20 continue
c
c2012.7.24
      if(zz(i).lt.tstart) go to 10
c2012.7.24
      rmdti1=0.0
c     write(6,*) rmdi
c     ramdai=a1+a2*rmdi+rmdti1
      ramdai=a1*exp(sfk0)+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
c     d1ff=d1ff+1.d0/ramdai
      d1ff=d1ff+exp(sfk0)/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
c
   10 continue
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      ft=0.0
      area=tx*ty
      do 70 i=1,nn
c@@@6 if(zmg(i).lt.xmg1) go to 70
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rrd1=q( aa(i), bb(i), cc(i), xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), ty-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), ty-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), ty-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), ty-yy(i) )
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rrd1=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), yy(i) )
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rrd1+rrd2)/2
      
      scl=exp(2*a4*zmg(i))
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
      spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
     &                                  +a6 **(1-a7)*log(       a6) )
c
c2021.7.24
      if(zz(i).lt.tstart) go to 80
c2021.7.24
      ft=ft+1.d0/(1.d0-a5)
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))     *spac2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)
     &                 -a3**(1.d0-a5))               * d4s2*delta/(2*pi)
      d5ft=d5ft+1.d0/(1.d0-a5)**2
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
     &    +1.d0/(1.d0-a5)*(-(tz-zz(i)+a3)**(1.d0-a5)
     &      *log(tz-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                                               * d6s2*delta/(2*pi)
      d7ft=d7ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                                               * d7s2*delta/(2*pi)
   80 continue
   70 continue
c
      f2=0.0
      nox=101
      noy=101
      dx=tx/(nox-1)
      dy=ty/(noy-1)
      do 150 kx=1,nox
      do 150 ky=1,noy
      xk=tx*(kx-1)/(nox-1)
      yk=ty*(ky-1)/(noy-1)
      ip=int(xk/hx)+1
      jp=int(yk/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xk/hx-(ip-1)
      ry=yk/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 120 j=0,3
      do 120 i=0,3
      ij=(ip+i)+(mx+3)*(jp+j-1)
  120 sfk0=sfk0+a0(ij)*rnx(4-i)*rny(4-j)
c................................................
      if(sfk0.gt.170.0) go to 50
      ram=exp(sfk0)
      rmd=ram
      if(kx.eq.1 .or. kx.eq.nox) rmd=ram/2
      if(ky.eq.1 .or. ky.eq.noy) rmd=ram/2
      if(kx.eq.1 .and. ky.eq.1) rmd=ram/3
      if(kx.eq.nox .and. ky.eq.noy) rmd=ram/3
      if(kx.eq.1 .and. ky.eq.noy) rmd=ram/6
      if(kx.eq.nox .and. ky.eq.1) rmd=ram/6
      f2=f2+rmd
  150 continue
      volume=f2*dx*dy
c     volume=area
c
c
c     fs=a1*tz*area+a2*ft
c     d1fs=tz*area
      fs=a1*tz*volume+a2*ft
      d1fs=tz*volume
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*d5ft
      d6fs=a2*d6ft
      d7fs=a2*d7ft
c
      stsum=0.0
      f=ff-fs-stsum
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
c
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
c     h(4)=-g4*2.0d0*b(4)
      h(4)=-g4      *b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
      if(igrd.eq.0) h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
      if(a3.le..1d-4) h(3)=0.0
c     if(a3.le..1d-4) b(3)=0.1d-3
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine func18(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified Omori type
c     etas point process (with trend: not realized)
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum K*\exp(\beta_j)/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      dimension rnx(4),rny(4)
      common /approx/ xmg1
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
      
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
c     a4=b(4)**2
      a4=b(4)**2/2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
      a8=b(8)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
      if(a4.gt.5.d0) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
      if(a8.gt.5.d0) go to 50
c
c
      hx=tx/mx
      hy=ty/my
      ip=int(xx(1)/hx)+1
      jp=int(yy(1)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(1)/hx-(ip-1)
      ry=yy(1)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 45 i=0,3
      do 45 j=0,3
        ik=(ip+i)+(mx+3)*(jp+j-1)
        sfk0=sfk0+a0(ik)*rnx(4-i)*rny(4-j)
   45 continue
c
c
      rmdti1=0.0
c     ramdai=a1+rmdti1
      ramdai=a1*exp(sfk0)+rmdti
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
c     d1ff=1.d0/ramdai
      d1ff=exp(sfk0)/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
      d8ff=0.0
c
      do 10 i=2,nn
c
      ip=int(xx(i)/hx)+1
      jp=int(yy(i)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(i)/hx-(ip-1)
      ry=yy(i)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 41 k=0,3
      do 41 l=0,3
        ik=(ip+k)+(mx+3)*(jp+l-1)
        sfk0=sfk0+a0(ik)*rnx(4-k)*rny(4-l)
   41 continue
c
c
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      rmd8i=0.0
      do 20 j=1,i-1
c@@@6 if(zmg(j).lt.xmg1) go to 20
      xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
      dxxij4=(xx(i)-xx(j))*exp(-a4*zmg(j))*(-zmg(j))
      yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
      dyyij4=(yy(i)-yy(j))*exp(-a4*zmg(j))*(-zmg(j))
      space=(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7)
      d6spa=
     &-a7*(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7-1)
      d4spa=d6spa*( 2*aa(j)*xxij*dxxij4-
     &  4*bb(j)*xxij*yyij*(-zmg(j))+2*cc(j)*yyij*dyyij4 )
      d7spa=-space*
     &  log(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)
      rmdi=rmdi+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5         *space
      rmd3i=rmd3i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd4i=rmd4i+exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5*d4spa
      rmd5i=rmd5i-exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5
     &                    *log(zz(i)-zz(j)+a3)    *space
      rmd6i=rmd6i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd7i=rmd7i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d7spa
      rmd8i=rmd8i+zmg(j)*exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5*space
   20 continue
c
      rmdti1=0.0
c     write(6,*) rmdi
c     ramdai=a1+a2*rmdi+rmdti1
      ramdai=a1*exp(sfk0)+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
c     write(6,*) 'ff',i,ff
c     d1ff=d1ff+1.d0/ramdai
      d1ff=d1ff+exp(sfk0)/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
      d8ff=d8ff+a2*rmd8i/ramdai
c
   10 continue
c     write(6,*) ff
c     stop
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      d8ft=0.0
      ft=0.0
      area=tx*ty
      do 70 i=1,nn
c@@@6 if(zmg(i).lt.xmg1) go to 70
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rrd1=q( aa(i), bb(i), cc(i), xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), ty-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), ty-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), ty-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), ty-yy(i) )
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rrd1=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), yy(i) )
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rrd1+rrd2)/2
      
      scl=exp(2*a4*zmg(i))
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
      spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
     &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
     &                                  +a6 **(1-a7)*log(       a6) )
c
      ft=ft+1.d0/(1.d0-a5)*exp(a8*zmg(i))
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))*spac2*delta/(2*pi)
     &          *exp(a8*zmg(i))
      d4ft=d4ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &          *exp(a8*zmg(i)) * d4s2*delta/(2*pi)
      d5ft=d5ft+1.d0/(1.d0-a5)**2 *exp(a8*zmg(i))
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
     &     +1.d0/(1.d0-a5) *exp(a8*zmg(i)) *(-(tz-zz(i)+a3)**(1.d0-a5)
     &      *log(tz-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &                        * exp(a8*zmg(i)) * d6s2*delta/(2*pi)
      d7ft=d7ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
     &          * exp(a8*zmg(i)) * d7s2*delta/(2*pi)
      d8ft=d8ft+1.d0/(1.d0-a5)*exp(a8*zmg(i))*zmg(i)
     &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
   80 continue
   70 continue
c
      f2=0.0
      nox=101
      noy=101
      dx=tx/(nox-1)
      dy=ty/(noy-1)
      do 150 kx=1,nox
      do 150 ky=1,noy
      xk=tx*(kx-1)/(nox-1)
      yk=ty*(ky-1)/(noy-1)
      ip=int(xk/hx)+1
      jp=int(yk/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xk/hx-(ip-1)
      ry=yk/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 120 j=0,3
      do 120 i=0,3
      ij=(ip+i)+(mx+3)*(jp+j-1)
  120 sfk0=sfk0+a0(ij)*rnx(4-i)*rny(4-j)
c................................................
      if(sfk0.gt.170.0) go to 50
      ram=exp(sfk0)
      rmd=ram
      if(kx.eq.1 .or. kx.eq.nox) rmd=ram/2
      if(ky.eq.1 .or. ky.eq.noy) rmd=ram/2
      if(kx.eq.1 .and. ky.eq.1) rmd=ram/3
      if(kx.eq.nox .and. ky.eq.noy) rmd=ram/3
      if(kx.eq.1 .and. ky.eq.noy) rmd=ram/6
      if(kx.eq.nox .and. ky.eq.1) rmd=ram/6
      f2=f2+rmd
  150 continue
      volume=f2*dx*dy
c     volume=area
c
c
c     fs=a1*tz*area+a2*ft
c     d1fs=tz*area
      fs=a1*tz*volume+a2*ft
      d1fs=tz*volume
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*d5ft
      d6fs=a2*d6ft
      d7fs=a2*d7ft
      d8fs=a2*d8ft
c
      stsum=0.0
      f=ff-fs-stsum
c     write(6,*) 'ff',f,ff,fs,stsum
c     stop
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
      g8=d8ff-d8fs
c
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
c     h(4)=-g4*2.0d0*b(4)
      h(4)=-g4      *b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
      h(8)=-g8*2.0d0*b(8)
c@@@@
c     h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
      if(a3.le..1d-4) h(3)=0.0
c     if(a3.le..1d-4) b(3)=0.1d-3
c     h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c*********************************************************************
      subroutine bsp3(r,h,idef,bn)
c
c     calculation of 4 non-zero b-splines
c.....................................................................
c input
c      r    : normalized coordinate
c      h    : knot interval
c      idef : order of differenciation
c output
c      bn(4): four b-splines
c.....................................................................
      implicit real*8(a-h,o-z)
      dimension bn(4)
c
      if(idef.eq.2) then
        a=1.d0/(h*h)
        bn(1)=r*a
        bn(2)=(-3*r+1)*a
        bn(3)=(3*r-2)*a
        bn(4)=(-r+1)*a
      endif
      r2=r*r
      if(idef.eq.1) then
        a=0.5d0/h
        bn(1)=r2*a
        bn(2)=(-3*r2+2*r+1)*a
        bn(3)=(3*r2-4*r)*a
        bn(4)=(-r2+2*r-1)*a
      endif
      r3=r2*r
      if(idef.eq.0) then
        a=1.d0/6.d0
        bn(1)=r3*a
        bn(2)=(3*(-r3+r2+r)+1)*a
        bn(3)=(3*r3-6*r2+4)*a
        bn(4)=(-r3+3*(r2-r)+1)*a
      endif
c
      return
      end
c***********************************************************************
      subroutine func19(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified Omori type
c     etas point process (with trend: not realized)
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum K*\exp(\beta_j)/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222)
      common /domain/ tx,ty,tz,a0(999),mx,my
      common /data/ xx(npa),yy(npa),zz(npa),zmg(npa),nn,mm
      common /coef/ aa(npa),bb(npa),cc(npa)
      common /ddd/r,fff,aic,sd
      common /sv/bx(4096),by(4096),index,mbx,mby
      common     / ccc /  isw , ipr,igrd
      data pi /3.14159265358979d0/
      dimension r(31,31),he(50,50)
      dimension b(50),h(50),ggt(50),gt(50),at(50)
c     dimension bx(4096),by(4096)
      dimension rnx(4),rny(4)
      common /approx/ xmg1
!sgi  --- start ---
c     save mbx,mby,bx,by
      real*8 tmp1,logtmp1,tmp2,logtmp2,exptmp2
      real*8 tmp3,tmp4,tmp5,tmp6,tmp30,tmp40,tmp50,tmp60
!sgi  ---- end ----
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
      
cn
c     data index /0/
      if(index.eq.0) then
      mbx=2**7
      mby=2**7
      do 30 i=1,mbx+1
      bx(i)=(i-1)*tx/mbx
   30 continue
      do 40 i=1,mby+1
      by(i)=(i-1)*ty/mby
   40 continue
      end if
      index=1
cn
      kkt=0
      ifg=0
      a1=b(1)**2
      a2=b(2)**2
      a3=b(3)**2
c     a4=b(4)**2
      a4=b(4)**2/2
      a5=b(5)**2
      a6=b(6)**2
      a7=b(7)**2
      a8=b(8)**2
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
      if(a4.gt.5.d0) go to 50
      if(a5.gt.3.d0) go to 50
      if(a7.gt.5.d0) go to 50
      if(a8.gt.5.d0) go to 50
c
c
      hx=tx/mx
      hy=ty/my
      ip=int(xx(1)/hx)+1
      jp=int(yy(1)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(1)/hx-(ip-1)
      ry=yy(1)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 45 i=0,3
      do 45 j=0,3
        ik=(ip+i)+(mx+3)*(jp+j-1)
        sfk0=sfk0+a0(ik)*rnx(4-i)*rny(4-j)
   45 continue
c
c
      rmdti1=0.0
c     ramdai=a1+rmdti1
      ramdai=a1*exp(sfk0)+rmdti
      if(ramdai.le.0.0) go to 50
      ff=log(ramdai)
c     d1ff=1.d0/ramdai
      d1ff=exp(sfk0)/ramdai
      d2ff=0.0
      d3ff=0.0
      d4ff=0.0
      d5ff=0.0
      d6ff=0.0
      d7ff=0.0
      d8ff=0.0
c
      do 10 i=2,nn
c
      ip=int(xx(i)/hx)+1
      jp=int(yy(i)/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xx(i)/hx-(ip-1)
      ry=yy(i)/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 41 k=0,3
      do 41 l=0,3
        ik=(ip+k)+(mx+3)*(jp+l-1)
        sfk0=sfk0+a0(ik)*rnx(4-k)*rny(4-l)
   41 continue
c
c
      rmdi=0.0
      rmd3i=0.0
      rmd4i=0.0
      rmd5i=0.0
      rmd6i=0.0
      rmd7i=0.0
      rmd8i=0.0
      do 20 j=1,i-1
c@@@6 if(zmg(j).lt.xmg1) go to 20
      xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
      dxxij4=(xx(i)-xx(j))*exp(-a4*zmg(j))*(-zmg(j))
      yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
      dyyij4=(yy(i)-yy(j))*exp(-a4*zmg(j))*(-zmg(j))
c sgi
      tmp1=(aa(j)*xxij*xxij-2*bb(j)*xxij*yyij+cc(j)*yyij*yyij+a6)
      logtmp1=log(tmp1)
      space=exp(-a7*logtmp1)
      d6spa=
     &-a7*space/tmp1
      d4spa=d6spa*( 2*aa(j)*xxij*dxxij4-
     &  4*bb(j)*xxij*yyij*(-zmg(j))+2*cc(j)*yyij*dyyij4 )
      d7spa=-space*
     &  logtmp1
      tmp2=zz(i)-zz(j)+a3
      logtmp2=log(tmp2)
      exptmp2=exp(-a5*logtmp2)
      exptmp3=exp(a8*zmg(j))
c     space=(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7)
c     d6spa=
c    &-a7*(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)**(-a7-1)
c     d4spa=d6spa*( 2*aa(j)*xxij*dxxij4-
c    &  4*bb(j)*xxij*yyij*(-zmg(j))+2*cc(j)*yyij*dyyij4 )
c     d7spa=-space*
c    &  log(aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6)
c     rmdi=rmdi+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5         *space
      rmdi=rmdi+exptmp3*exptmp2         *space 
c     rmd3i=rmd3i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**(a5+1)   *space
      rmd3i=rmd3i+exptmp3*exptmp2/tmp2   *space
c     rmd4i=rmd4i+exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5*d4spa
      rmd4i=rmd4i+exptmp3       *exptmp2*d4spa
c     rmd5i=rmd5i-exp(a8*zmg(j))       /(zz(i)-zz(j)+a3)**a5
c    &                    *log(zz(i)-zz(j)+a3)    *space
      rmd5i=rmd5i-exptmp3       *exptmp2
     &                    *logtmp2    *space
c     rmd6i=rmd6i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d6spa
      rmd6i=rmd6i+exptmp3*exptmp2       *d6spa
c     rmd7i=rmd7i+exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5       *d7spa
      rmd7i=rmd7i+exptmp3*exptmp2       *d7spa
c     rmd8i=rmd8i+zmg(j)*exp(a8*zmg(j))/(zz(i)-zz(j)+a3)**a5*space
      rmd8i=rmd8i+zmg(j)*exptmp3*exptmp2*space
csgi
   20 continue
c
      rmdti1=0.0
c     write(6,*) rmdi
c     ramdai=a1+a2*rmdi+rmdti1
      ramdai=a1*exp(sfk0)+a2*rmdi+rmdti1
      if(ramdai.le.0.0) go to 50
      ff=ff+log(ramdai)
c     write(6,*) 'ff',i,ff,rmdi,rmdii
c     d1ff=d1ff+1.d0/ramdai
      d1ff=d1ff+exp(sfk0)/ramdai
      d2ff=d2ff+rmdi/ramdai
      d3ff=d3ff-a2*a5*rmd3i/ramdai
      d4ff=d4ff+a2*rmd4i/ramdai
      if(a5.ne.1.d0) d5ff=d5ff+a2*rmd5i/ramdai
      d6ff=d6ff+a2*rmd6i/ramdai
      d7ff=d7ff+a2*rmd7i/ramdai
      d8ff=d8ff+a2*rmd8i/ramdai
c
   10 continue
c     write(6,*) ff
c     stop
c
      d3ft=0.0
      d4ft=0.0
      d5ft=0.0
      d6ft=0.0
      d7ft=0.0
      d8ft=0.0
      ft=0.0
      area=tx*ty
csgi  --- start ---
      tmp5=a6**(-a7)
      tmp6=a6*tmp5
      tmp50=a3**(-a5)
      tmp60=a3*tmp50
csgi  ---- end ----
      do 70 i=1,nn
c@@@6 if(zmg(i).lt.xmg1) go to 70
csgi  --- start ---
      tmp30=(tz-zz(i)+a3)**(-a5)
      tmp40=(tz-zz(i)+a3)*tmp30
      exptmp3=exp(a8*zmg(i))
csgi  ---- end ----
cn
      do 80 j=1,2*(mbx+mby)
      if(j.le.mby) then
      if(xx(i).le.0.0) go to 80
      jy=j
      rrd1=q( aa(i), bb(i), cc(i), xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mby .and. j.le. mbx+mby) then
      if(yy(i).ge.ty) go to 80
      jx=j-mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), ty-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), ty-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), ty-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), ty-yy(i) )
      rr3=(tx/mbx)**2
      end if
      if(j.gt.mbx+mby .and. j.le. mbx+2*mby) then
      if(xx(i).ge.tx) go to 80
      jy=j-mbx-mby
      rrd1=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy)-yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), tx-xx(i), by(jy+1)-yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy)-yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, tx-xx(i), by(jy+1)-yy(i) )
      rr3=(ty/mby)**2
      end if
      if(j.gt.mbx+2*mby .and. j.le. 2*mbx+2*mby) then
      if(yy(i).le.0.0) go to 80
      jx=j-mbx-2*mby
      rrd1=q( aa(i), bb(i), cc(i), bx(jx)-xx(i), yy(i) )
      rrd2=q( aa(i), bb(i), cc(i), bx(jx+1)-xx(i), yy(i) )
      rr1=q( 1.0d0, 0.0d0, 1.0d0, bx(jx)-xx(i), yy(i) )
      rr2=q( 1.0d0, 0.0d0, 1.0d0, bx(jx+1)-xx(i), yy(i) )
      rr3=(tx/mbx)**2
      end if
c
      delta=dacos(  (rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )  )
      rr=(rrd1+rrd2)/2
      
      scl=exp(2*a4*zmg(i))
!sgi  --- start ---
      tmp3=(rr/scl+a6)**(-a7)
      tmp4=(rr/scl+a6)*tmp3
      dscl4=exp(2*a4*zmg(i))*2*zmg(i)
c     spac2=pi/(1-a7)   *scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
      spac2=pi/(1-a7)   *scl  *( tmp4-tmp6 )
c     d4s2= pi/(1-a7)   *dscl4*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
c    &     +pi          *scl  *  (rr/scl+a6)**( -a7)*(-rr/scl**2*dscl4)
      d4s2= pi/(1-a7)   *dscl4*( tmp4-tmp6 )
     &     +pi          *scl  *  tmp3*(-rr/scl**2*dscl4)
c     d6s2= pi          *scl  *( (rr/scl+a6)**( -a7)-a6**( -a7) )
      d6s2= pi          *scl  *( tmp3-tmp5 )
c     d7s2= pi/(1-a7)**2*scl  *( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
c    &     +pi/(1-a7)   *scl  *(-(rr/scl+a6)**(1-a7)*log(rr/scl+a6)
c    &                                  +a6 **(1-a7)*log(       a6) )
      d7s2= pi/(1-a7)**2*scl  *( tmp4-tmp6 )
     &     +pi/(1-a7)   *scl  *(-tmp4*log(rr/scl+a6)
     &                                  +tmp6*log(       a6) )
c
c     ft=ft+1.d0/(1.d0-a5)*exp(a8*zmg(i))
c    &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
      ft=ft+1.d0/(1.d0-a5)*exptmp3
     &      *(tmp40-tmp60)*spac2*delta/(2*pi)
c     d3ft=d3ft+((tz-zz(i)+a3)**(-a5)-a3**(-a5))*spac2*delta/(2*pi)
c    &          *exp(a8*zmg(i))
      d3ft=d3ft+(tmp30-tmp50)*spac2*delta/(2*pi)
     &          *exptmp3
c     d4ft=d4ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
c    &          *exp(a8*zmg(i)) * d4s2*delta/(2*pi)
      d4ft=d4ft+1.d0/(1.d0-a5)*(tmp40-tmp60)
     &          *exptmp3 * d4s2*delta/(2*pi)
c     d5ft=d5ft+1.d0/(1.d0-a5)**2 *exp(a8*zmg(i))
c    &      *((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))*spac2*delta/(2*pi)
c    &     +1.d0/(1.d0-a5) *exp(a8*zmg(i)) *(-(tz-zz(i)+a3)**(1.d0-a5)
c    &      *log(tz-zz(i)+a3)+a3**(1.d0-a5)*log(a3)) *spac2*delta/(2*pi)
      d5ft=d5ft+1.d0/(1.d0-a5)**2 *exptmp3
     &      *(tmp40-tmp60)*spac2*delta/(2*pi)
     &     +1.d0/(1.d0-a5) *exptmp3 *(-tmp40
     &      *log(tz-zz(i)+a3)+tmp60*log(a3)) *spac2*delta/(2*pi)
c     d6ft=d6ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
c    &                        * exp(a8*zmg(i)) * d6s2*delta/(2*pi)
      d6ft=d6ft+1.d0/(1.d0-a5)*(tmp40-tmp60)
     &                        * exptmp3 * d6s2*delta/(2*pi)
c     d7ft=d7ft+1.d0/(1.d0-a5)*((tz-zz(i)+a3)**(1.d0-a5)-a3**(1.d0-a5))
c    &          * exp(a8*zmg(i)) * d7s2*delta/(2*pi)
      d7ft=d7ft+1.d0/(1.d0-a5)*(tmp40-tmp60)
     &          * exptmp3 * d7s2*delta/(2*pi)
      d8ft=d8ft+1.d0/(1.d0-a5)*exptmp3*zmg(i)
     &      *(tmp40-tmp60)*spac2*delta/(2*pi)
   80 continue
   70 continue
c
      f2=0.0
      nox=101
      noy=101
      dx=tx/(nox-1)
      dy=ty/(noy-1)
      do 150 kx=1,nox
      do 150 ky=1,noy
      xk=tx*(kx-1)/(nox-1)
      yk=ty*(ky-1)/(noy-1)
      ip=int(xk/hx)+1
      jp=int(yk/hy)+1
      if(ip.gt.mx) ip=mx
      if(jp.gt.my) jp=my
      rx=xk/hx-(ip-1)
      ry=yk/hy-(jp-1)
      call bsp3(rx,hx,0,rnx)
      call bsp3(ry,hy,0,rny)
      sfk0=0.0
      do 120 j=0,3
      do 120 i=0,3
      ij=(ip+i)+(mx+3)*(jp+j-1)
  120 sfk0=sfk0+a0(ij)*rnx(4-i)*rny(4-j)
c................................................
      if(sfk0.gt.170.0) go to 50
      ram=exp(sfk0)
      rmd=ram
      if(kx.eq.1 .or. kx.eq.nox) rmd=ram/2
      if(ky.eq.1 .or. ky.eq.noy) rmd=ram/2
      if(kx.eq.1 .and. ky.eq.1) rmd=ram/3
      if(kx.eq.nox .and. ky.eq.noy) rmd=ram/3
      if(kx.eq.1 .and. ky.eq.noy) rmd=ram/6
      if(kx.eq.nox .and. ky.eq.1) rmd=ram/6
      f2=f2+rmd
  150 continue
      volume=f2*dx*dy
c     volume=area
c
c
c     fs=a1*tz*area+a2*ft
c     d1fs=tz*area
      fs=a1*tz*volume+a2*ft
      d1fs=tz*volume
      d2fs=ft
      d3fs=a2*d3ft
      d4fs=a2*d4ft
      d5fs=0.0
      if(a5.ne.1.d0) d5fs=a2*d5ft
      d6fs=a2*d6ft
      d7fs=a2*d7ft
      d8fs=a2*d8ft
c
      stsum=0.0
      f=ff-fs-stsum
c     write(6,*) 'ff',f,ff,fs,stsum
c     stop
      g1=d1ff-d1fs
      g2=d2ff-d2fs
      g3=d3ff-d3fs
      g4=d4ff-d4fs
      g5=d5ff-d5fs
      g6=d6ff-d6fs
      g7=d7ff-d7fs
      g8=d8ff-d8fs
c
      f=-f
      h(1)=-g1*2.0d0*b(1)
      h(2)=-g2*2.0d0*b(2)
      h(3)=-g3*2.0d0*b(3)
c     h(4)=-g4*2.0d0*b(4)
      h(4)=-g4      *b(4)
      h(5)=-g5*2.0d0*b(5)
      h(6)=-g6*2.0d0*b(6)
      h(7)=-g7*2.0d0*b(7)
      h(8)=-g8*2.0d0*b(8)
c@@@@
c     h(1)=0.0
c     h(2)=0.0
c     h(3)=0.0
      if(a3.le..1d-4) h(3)=0.0
c     if(a3.le..1d-4) b(3)=0.1d-3
      h(4)=0.0
c     h(5)=0.0
c     h(6)=0.0
c     h(7)=0.0
      if(a5.eq.1.d0) h(5)=0.0
c     write(6,1) a1,a2,a3,a4,a5,a6,a7,a8
c     write(6,1) f,ff,fs
c     write(6,1)(h(i),i=1,n)
c
c@@@  n=7
c     call hes4(n,b,he)
c     do 653 i=1,n
c 653 write(6,604) (he(i,j),j=1,n)
c     call invdet(he,hdet,n,50)
c     write(6,602)
c     do 652 i=1,n
c 652 write(6,604) (he(i,j),j=1,n)
c 602 format(1h0,'*** estimated inverse hessian ***')
c 604 format(1h ,8d15.5/(1h ,8d15.5))
c     stop 666
      fff=f
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
      return
   50 continue
      ifg=1
      f=1.0d30
      return
      end
c*********************************************************************
