c  simulate nonhomogeneous spatial Poisson pattern with stationary Poisson in time
      program main
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /domain/ ramdalog(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      common /orgepi/ alon(npa),alat(npa),aspect
      common /pseudorandom/ ix,iy,iz
      dimension tt(nda),xx(nda),yy(nda),am(nda),sam(nda)
      real*4 r,u
      data pi/3.14159265358979d0/
*     data tx0,tx,ty0,ty/128.0,21.0,30.0,17.0/ !main Japan area
      open(2,file='poisson.conf')
      read(2,*) tx0,tx,ty0,ty
      read(2,*) ix,iy,iz
      close(2)
      write(6,*) ix,iy,iz
!     stop 2222
c
      tx1=tx
      aspect=cos( pi/180*(ty0+ty/2) )
      tx=tx*aspect
      area=tx*ty
      open(3,file='./delo2d-poisson.omap')
      read(3,*) w1,abic,npex
c     write(6,*) w1,abic,npex
      read(3,*) (ramdalog(i),i=1,npex)
c     write(6,*) (ramdalog(i),i=1,npex)
      close(3)
c     do j=1,npex
c       write(6,*) j,ramdalog(j),exp(ramdalog(j))
c     enddo
      open(1,file='./delone2.out')
      read(1,*) nev,npext,ndelo,tx,ty
      asi=nev/area  !average spatial intensity: an event/deg^2
      write(6,*) '#shocks, #shocks/deg^2', nev,asi
      do 21 i=1,npext
        read(1,*,end=25) idum,alon(i),alat(i),am(i),otime,dd,ee,ff
        if(i.eq.nev) otend=otime
   21 continue
      cm0=minval(am(1:nev))
   25 continue
      do 509 i=1,ndelo
        read(1,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509 continue
      close(1)
c option
      call nhpint(tx,ty,aspect,ainteg,rmax) ! Integral of the intensity
      write(6,*) 'ainteg,rmax',ainteg,rmax
c option-end
      time=0.0
      i=1
      tend=int(otend)+1
      write(6,9) tend,otend,tx0,ty0
    9 format('tend,otend,tx0,ty0=',4f10.3)
   10 continue
      call pseudo(u)
      d=-log(u+0.1d-8)/(nev/tend)
c     write(6,*) 'time',i, d, u
      time=time+d
      tt(i)=time
   20 continue
      call fs(x,y,tx,ty,xx,yy)
      if(x.ge..0 .and. x.le.tx .and. y.ge..0 .and. y.le.ty) then
        xx(i)=x
        yy(i)=y
        i=i+1
        if(i.gt.nev) go to 50
        go to 10
      else
        write(6,*) 'out',t,x,y
        go to 20
      endif
    6 format(i6,f10.3,f9.3,f6.2,f14.5,i6,f5.2,i6)
      i=i+1
      go to 10
   50 continue
      tend=int(tt(i-1))+1.0
      nev=i-1
      do 110 i=1,nev
c       xx(i)=tx0+xx(i)*tx1/tx
        xx(i)=tx0+xx(i)
        yy(i)=ty0+yy(i)
  110 continue
      write(6,1) nev
      write(6,5) tend,cm0
      write(2,2) nev,tx1,ty,cm0,tend
      write(2,2) (i,xx(i),yy(i),am(i),tt(i), i=1,nev)
    1 format('#eqs=',2i10)
    2 format(i5,3f12.5,f14.5)
    3 format(8f10.5)
    4 format(2f10.1,i10)
    5 format('simtend,cutmag=',6f13.4)
      stop
      end
c
      subroutine fs(x,y,tx,ty,xx,yy)
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /orgepi/ alon(npa),alat(npa),aspect
      dimension tt(npa),xx(npa),yy(npa),am(npa),sam(npa)
      dimension xd(3),yd(3),zmu(3),zk(3),za(3),zp(3),zq(3)
      common /domain/ ramdalog(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      dimension ramda(nda)
      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
      do 70 j=1,npext
        ramda(j)=exp(ramdalog(j))
   70 continue
      rmax=maxval(ramda(1:npext))
  100 continue
      call pseudo(r1)
      call pseudo(r2)
      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)=ramdalog(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
      ram=exp(sfk0) 
      call pseudo(r)
      if(ram/rmax.gt.1.d0) write(6,*) 'check', r, ram,rmax
      if(r.gt.ram/rmax) go to 100
    9 format(i5,f5.1,e20.8,2f20.8)
    8 format(i5,f5.1,e20.8,4f12.5)
      return
      end
c
      subroutine pseudo(random)
c     wichmann+hill appl. statist. (31) 188-? (1982)
c the following triplet integeasr define random number seeds.
      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/
!     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 nhpint(tx,ty,aspect,ainteg,rmax)
c Integral of intensity and maximum of intensity
c.....................................................................
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      common /domain/ ramdalog(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),nev,npext,ndelo
      dimension ramda(npa)
      rmax=0.0
      do 70 j=1,npext
        ramda(j)=exp(ramdalog(j))
c       if(ramda(j).ge.rmax) rmax=ramda(j)
   70 continue
      rmax=maxval(ramda(1:npext))
      ainteg=0.0
      do 180 id=1,ndelo
        do 190 k=1,3
          idk=ldt(id,k)
          ainteg=ainteg+ramda(idk)/3*adt(id)
  190   continue
  180 continue
      ainteg=ainteg*aspect
      write(6,*) 'ainteg,rmax',ainteg,rmax
      return
      end
