c  simulate HIST-ETAS models with given magnitude sequence
      program main
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /domain/ amulog(nda),aklog(nda),alflog(nda),plog(nda),
     &                qlog(nda),ainteg,bgrmax,mx,my,nmain
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      common /basepa/ amu0,ak0,c,alf0,p0,d,q0
      common /orgepi/ alon(npa),alat(npa),aspect
      common /ancestor/ j1
      common /nsimevents/ nn
      common /pseudorandom/ ix,iy,iz
      dimension ak(nda),alf(nda),p(nda),q(nda)
      dimension tt(nda),xx(nda),yy(nda),am(nda),sam(nda),bm(nda)
*     real*4 r,u
      data pi/3.14159265358979d0/
Choose either of the two model.
      model=7  ! hist-etas-mk
      model=5  ! hist-etas5pa
c Parameters of the b-value and space-time volume, where b-value is not
c used in this particular program: you can modify below if you simulate
c magnitude sequence.
      data bmg,cm0,tx0,tx,ty0,ty,tend/ 
     &     1.1, 6.0, 128.0, 21.0, 30.0, 17.0, 14012.0/ !alljapanMc6
c    &     0.9, 6.0, 128.0, 21.0, 30.0, 17.0, 14012.0/ !alljapanMc6
c
      open(2,file='histetasim.conf')
      read(2,*) model
      write(6,*) 'model=', model
      read(2,*) bmg,cm0,tx0,tx,ty0,ty,tend
      write(6,8) bmg,cm0,tx0,tx,ty0,ty,tend
    8 format(7f8.1)
      read(2,7) ix,iy,iz
      write(6,7) ix,iy,iz
    7 format(3i5)
      close(2) 
!     stop 7777
c
      beta=exp(bmg)
      tx1=tx
      tx=tx*cos( pi/180*(ty0+ty/2) )
      aspect=cos( pi/180*(ty0+ty/2) )
c     aspect=1.0d0
      area=tx*ty
c
c option of models, either hist-etas-mk or hist-etas5pa
      if(model.eq.7) then
        open(3,file='./hist-etas-mk.upda')
        read(3,*) w1,w2,abic,npa2
        read(3,*) amu0,ak0,c,alf0,p0,d,q0 ! baseline parameters
        read(3,*) (amulog(i),i=1,npa2/2),(aklog(i),i=1,npa2/2)
        close(3)
        do i=1,npa2/2 ! paramters are const. except for amu(.) and ak(.)
          ak(i)=ak0*exp(aklog(i))
          alflog(i)=0.0
          plog(i)=0.0
          qlog(i)=0.0
          alf(i)=alf0
          p(i)=p0
          q(i)=q0
        end do
      endif
c 
      if(model.eq.5) then
        open(3,file='./hist-etas5pa.upda')
        read(3,*) w1,w2,w3,w4,w5,abic,npa5,w00
        read(3,*) amu0,ak0,c,alf0,p0,d,q0    ! baseline param.
        read(3,*) (amulog(i),i=1,npa5/5),(aklog(i),i=1,npa5/5),
     & (alflog(i),i=1,npa5/5),(plog(i),i=1,npa5/5),(qlog(i),i=1,npa5/5)
        close(3)
      endif
c
      open(1,file='./delone2.out')
      read(1,*) nev,npext,ndelo,tx,ty
      nn=nev
*     nn=308
      write(6,1) nn,nev
      do 21 i=1,npext
        read(1,*,end=25) idum,alon(i),alat(i),bm(i),otime,dd,ee,ff
   21 continue
      i=nev+1
   25 nev=i-1
      do 509 i=1,ndelo
        read(1,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509 continue
      close(1)
c
c Integration of the background intensity, and maximum b.g. intensity
      call bgrint(tx,ty,bmg)
      write(6,*) 'ainteg,bgrmax',ainteg,bgrmax
c
c Magnitude sequence
      do 40 i=1,nn
        am(i)=bm(i)
*       call pseudo(u)
*       am(i)=cm0-log10(u)/bmg
*       am(i)=cm0-log(1+u*(exp(-beta*(7.5-cm0))-1.0))/beta  !trunc.G-R
   40 continue
      write(6,5) t,cm0
      write(6,3) amu0,ak0,c,alf0,p0,d,q0,bmg
      do 30 j=1,nn
        exalf=exp(alf0*(am(j)-cm0))
        sam(j)=exalf
   30 continue
c
      t=0.0
      aity=ainteg
      uity=ainteg
      sity=nn/tend
      usty=nn/tend
      i=1
      j0=0
      write(7,*) nn
   10 continue
      call pseudo(r)
      e=-log(r+0.1d-8)/usty
      t=t+e
c     if(t.gt.tend) go to 20
      if(i.gt.1) then
        call ft(i-1,t,ak,alf,p,q,cm0,area,am,sam,tt,sity)
      endif
      probx=sity/usty
      if(probx.gt.1.d0) write(6,*) 'prob',probx
      if(probx.gt.1.d0) stop 119
      call pseudo(r)
      if(r.le.probx) then
        call fs(i-1,x,y,t,ak,alf,p,q,cm0,tx,ty,xx,yy,am,tt,sam,sity)
        if(x.ge..0 .and. x.le.tx .and. y.ge..0 .and. y.le.ty) then
          tt(i)=t
          xx(i)=x
          yy(i)=y
          if(j1.eq.0) j0=j0+1
          write(6,6) i,tx0+xx(i),ty0+yy(i),am(i),tt(i),j1,am(j1),j0
          write(7,6) i,tx0+xx(i),ty0+yy(i),am(i),tt(i),j1,am(j1),j0
    6 format(i6,f10.3,f9.3,f6.2,f14.5,i6,f5.2,i6)
c         if(mod(i,1).eq.0) write(6,2) i,am(i),tt(i),tt(i)-tt(i-1)
          if(i.ge.nn) go to 20
          call ft(i,t,ak,alf,p,q,cm0,area,am,sam,tt,usty)
          i=i+1
          go to 10
        else
          write(6,*) 'out',t,x,y
          go to 50
        end if
      else
        go to 50
      end if
   50 continue
      call ft(i-1,t,ak,alf,p,q,cm0,area,am,sam,tt,usty)
      go to 10
   20 continue
      tend=int(tt(i)+1)
      nn=i
      do 110 i=1,nn
c       xx(i)=tx0+xx(i)*tx1/tx
        xx(i)=tx0+xx(i)
        yy(i)=ty0+yy(i)
  110 continue
      write(6,1) nn
      write(6,5) tend,cm0
      write(2,2) nn,tx1,ty,cm0,tend
      write(2,2) (i,xx(i),yy(i),am(i),tt(i), i=1,nn)
c     write(2,2) nn,tx1,ty,cm0,tend/10**bmg
c     write(2,2) (i,xx(i),yy(i),am(i),tt(i)/10**bmg, i=1,nn)
    1 format(2i10)
    2 format(i5,4f12.5)
    3 format(8f10.5)
    4 format(2f10.1,i10)
    5 format(1h ,6f13.4)
      stop
      end
c
      subroutine ft(i,t,ak,alf,p,q,cm0,area,am,sam,tt,sity)
      implicit real*8(a-h,o-z)
      parameter(nda=23456)
      common /domain/ amulog(nda),aklog(nda),alflog(nda),plog(nda),
     &                qlog(nda),ainteg,bgrmax,mx,my,nmain
      common /basepa/ amu0,ak0,c,alf0,p0,d,q0
      common /nsimevents/ nn
      dimension tt(nda),xx(nda),yy(nda),am(nda),sam(nda),ak(nda),p(nda)
      dimension alf(nda),q(nda)
      data pi/3.14159265358979d0/
      sity=ainteg
      do 10 j=1,i
        sity=sity+ak(j)/(t-tt(j)+c)**p(j)*sam(j)
   10 continue
      return
      end
c
      subroutine fs(i,x,y,t,ak,alf,p,q,
     &              cm0,tx,ty,xx,yy,am,tt,sam,sity)
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /orgepi/ alon(npa),alat(npa),aspect
      common /nsimevents/ nn
      dimension tt(npa),xx(npa),yy(npa),am(npa),sam(npa)
      dimension alf(nda),p(nda),q(nda)
      dimension xd(3),yd(3),zmu(3),zk(3),za(3),zp(3),zq(3)
      common /domain/ amulog(nda),aklog(nda),alflog(nda),plog(nda),
     &                qlog(nda),ainteg,bgrmax,mx,my,nmain
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      common /basepa/ amu0,ak0,c,alf0,p0,d,q0
      common /ancestor/ j1
      data ic,icreplace/ 0,0/
      dimension ak(nda)
      save ic,icreplace
*     real*4 r,r1,r2
      data pi/3.14159265358979d0/
      det(x1,x2,x3,y1,y2,y3)=x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2
      upbgr=bgrmax
      call pseudo(r)
      sr=sity*r
      psum=ainteg
      do 10 j=1,i
        if(psum.gt.sr) go to 20
        psum=psum+ak(j)/(t-tt(j)+c)**p(j)*sam(j)
   10 continue
   20 continue
      j1=j-1
  100 continue
      call pseudo(r1)
      call pseudo(r2)
      if(j1.eq.0) then     ! This means background events.
        x=tx*r1
        y=ty*r2
        do 120 ipt=1,ndelo
          xmin=1.d69
          ymin=1.d69
          xmax=-1.d69
          ymax=-1.d69
          do 130 k=1,3
            xd(k)=alon(ldt(ipt,k))
            yd(k)=alat(ldt(ipt,k))
            zmu(k)=amulog(ldt(ipt,k))
            if(xd(k).lt.xmin) xmin=xd(k)
            if(xd(k).gt.xmax) xmax=xd(k)
            if(yd(k).lt.ymin) ymin=yd(k)
            if(yd(k).gt.ymax) ymax=yd(k)
  130     continue
          if(x.lt.xmin) go to 120
          if(x.gt.xmax) go to 120
          if(y.lt.ymin) go to 120
          if(y.gt.ymax) go to 120
          denom=det(xd(1),xd(2),xd(3),yd(1),yd(2),yd(3))
          if(denom.eq.0.0) write(6,*) 'denominator=0.0 for ipt =', ipt
          if(denom.eq.0.0) go to 120
          aa1=det(x    ,xd(2),xd(3),y    ,yd(2),yd(3))/denom
          aa2=det(xd(1),x    ,xd(3),yd(1),y    ,yd(3))/denom
          aa3=det(xd(1),xd(2),x    ,yd(1),yd(2),y    )/denom
          if(aa1.lt.eps.or. aa2.lt.eps.or. aa3.lt.eps) go to 120
          sfk0=aa1*zmu(1)+aa2*zmu(2)+aa3*zmu(3)
          go to 140
  120   continue
  140   continue
        bgr=amu0*exp(sfk0) 
        call pseudo(r)
        ic=ic+1
        if(bgr/upbgr.gt.1.d0) write(6,*) 'check', r, bgr,upbgr
        if(r.gt.bgr/upbgr) go to 100
        main=main+1
        nmain=nmain+1
c       write(6,*) 'main',main,nmain,ic
      else                               ! This means decendents.
c       exalf=exp(alf*(am(j1)-cm0))
        exalf=exp(alf(j1)*(am(j1)-cm0))
!       dist=sqrt( d*exalf*(r1**(1.d0/(1-q))-1.d0) )
        dist=sqrt( d*exalf*(r1**(1.d0/(1-q(j1)))-1.d0) )
        x=xx(j1)+dist*cos(2*pi*r2)
        y=yy(j1)+dist*sin(2*pi*r2)
      end if
      do 220 ipt=1,ndelo
        xmin=1.d69
        ymin=1.d69
        xmax=-1.d69
        ymax=-1.d69
        do 230 k=1,3
          xd(k)=alon(ldt(ipt,k))
          yd(k)=alat(ldt(ipt,k))
          zk(k)=aklog(ldt(ipt,k))
          za(k)=alflog(ldt(ipt,k))
          zp(k)=plog(ldt(ipt,k))
          zq(k)=qlog(ldt(ipt,k))
          if(xd(k).lt.xmin) xmin=xd(k)
          if(xd(k).gt.xmax) xmax=xd(k)
          if(yd(k).lt.ymin) ymin=yd(k)
          if(yd(k).gt.ymax) ymax=yd(k)
  230   continue
        if(x.lt.xmin) go to 220
        if(x.gt.xmax) go to 220
        if(y.lt.ymin) go to 220
        if(y.gt.ymax) go to 220
        denom=det(xd(1),xd(2),xd(3),yd(1),yd(2),yd(3))
        if(denom.eq.0.0) write(6,*) 'denominator=0.0 for ipt =', ipt
        if(denom.eq.0.0) go to 220
        aa1=det(x    ,xd(2),xd(3),y    ,yd(2),yd(3))/denom
        aa2=det(xd(1),x    ,xd(3),yd(1),y    ,yd(3))/denom
        aa3=det(xd(1),xd(2),x    ,yd(1),yd(2),y    )/denom
        if(aa1.lt.eps.or. aa2.lt.eps.or. aa3.lt.eps) go to 220
        sfkk=aa1*zk(1)+aa2*zk(2)+aa3*zk(3)
        sfka=aa1*za(1)+aa2*za(2)+aa3*za(3)
        sfkp=aa1*zp(1)+aa2*zp(2)+aa3*zp(3)
        sfkq=aa1*zq(1)+aa2*zq(2)+aa3*zq(3)
c       write(6,*) 'check', ipt,sfkk,aa1,zk(1),aa2,zk(2),aa3,zk(3)
        go to 240
  220 continue
  240 continue
      ak(i+1)=ak0*exp(sfkk)
      alf(i+1)=alf0*exp(sfka) 
      p(i+1)=p0*exp(sfkp) 
      q(i+1)=q0*exp(sfkq) 
      exalf=exp(alf(i+1)*(am(i+1)-cm0))
      sam(i+1)=pi*exalf/(1-q(i+1))*( 0 - d**(1-q(i+1)) )
    9 format(i5,f5.1,e20.8,2f20.8)
    8 format(i5,f5.1,e20.8,4f12.5)
c     write(6,*) 'parents',j1,x,y
      return
      end
c
      subroutine pseudo(random)
c     wichmann+hill appl. statist (31) 188-? (1982)
c The followings are seeds for uniform random numbers.
      implicit real*8(a-h,o-z)
      common /pseudorandom/ ix,iy,iz
c     data ix,iy,iz /1993,1119,1112/
!     data ix,iy,iz /1993,1119,1046/
c     data ix,iy,iz /1993,1119,1614/
c     data ix,iy,iz /1993,1119,1529/
c     data ix,iy,iz /1993,1119,1443/
c     data ix,iy,iz /1993,1119,1011/
c     data ix,iy,iz /1993,1119,1116/
c     data ix,iy,iz /1996,1118,1127/
!     save ix,iy,iz
      ix=171*mod(ix,177)-2*(ix/177)
      iy=172*mod(iy,176)-35*(iy/176)
      iz=170*mod(iz,178)-63*(iz/178)
      if (ix.lt.0) ix=ix+30269
      if (iy.lt.0) iy=iy+30307
      if (iz.lt.0) iz=iz+30323
      random=mod(float(ix)/30269.0+float(iy)/30307.0+
     1           float(iz)/30323.0,1.0)
      return
      end
c*********************************************************************
      subroutine bgrint(tx,ty,bmg)
c Integral (ainteg) and maximum (bgrmax) of the background intensity
c.....................................................................
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /orgepi/ alon(npa),alat(npa),aspect
      common /domain/ amulog(nda),aklog(nda),alflog(nda),plog(nda),
     &                qlog(nda),ainteg,bgrmax,mx,my,nmain
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      common /basepa/ amu0,ak0,c,alf0,p0,d,q0
      common /nsimevents/ nn
      dimension amu(npa)
      bgrmax=0.0
      sum=0.0
      do 70 j=1,npext
        amu(j)=amu0*exp(amulog(j))
        if(amu(j).ge.bgrmax) bgrmax=amu(j)
        sum=sum+amu(j)
   70 continue
      ainteg=0.0
      do 180 id=1,ndelo
        do 190 k=1,3
          idk=ldt(id,k)
          ainteg=ainteg+amu(idk)/3*adt(id)
  190   continue
  180 continue
      ainteg=ainteg*aspect*10.0**bmg
c     ainteg=ainteg*aspect
      write(6,*) 'ainteg,bgrmax',ainteg,bgrmax
      return
      end
