c  simulate magnitude for fixed locations 
!  simplified from ../test.f
      program main
      implicit real*8(a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
      dimension blog(npa),bv(npa)
      dimension adt(mdl),ldt(mdl,3),ndt(mdl,3)
      dimension alon(npa),alat(npa)
      dimension tt(nda),xx(nda),yy(nda),am(nda),bm(nda),cm(nda)
      dimension np(nda),xp(nda),yp(nda),tp(nda)
      real*4 r,u
      data pi/3.14159265358979d0/
c
      open(4,file='delo2d-bvalues.conf')
        read(4,*) tx0,ty0,am0
      close(4)
!     am0=6.0  !test
c
      open(1,file='./delone2.out')
        read(1,*) nev,npext,ndelo,tx,ty
        write(6,1) nev,ndelo
        do 21 i=1,npext
          read(1,*,end=25) idum,alon(i),alat(i),bm(i),tt(i),dd,ee,ff
   21   continue
        tend=int(tt(nev))+1
   25   continue
        do 509 i=1,ndelo
          read(1,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509   continue
      close(1)
c
      tx1=tx
      aspect=cos( pi/180*(ty0+ty/2) )
      tx=tx*aspect
c
      open(3,file='./delo2d-bvalues.omap')
      read(3,*) w1,abic,npex
      read(3,*) (blog(i),i=1,npex)
    9 format(4f20.0)
      close(3)
c
c
      open(2, file='simnhPoi.data')
      read(2,*) npo, tx, ty amc, tpend
      read(2,*) ((np(i),xp(i),yp(i),cm(i),tp(i)), i=1,npo)
      close(2)
c
      do i=1,npo
*       xi=alon(i)
*       yi=alat(i)
!       write(6,*) 'zz',xp(i),yp(i),tx0,ty0
        xi=xp(i)-tx0
        yi=yp(i)-ty0
!       write(6,*) 'before call', i,xi,yi,alon(i),alat(i)
        call interp(ndelo,alon,alat,blog,ldt,xi,yi,bpi)
!       write(6,3) i,xi,yi,bpi
        bv(i)=exp(bpi)
      enddo
c
c
!     do i=1,npex
!       bv(i)=exp(blog(i))
!     end do
c simulation of magnitude sequence
      do 40 i=1,nev
        call pseudo(u)
        am(i)=am0-log10(u)/bv(i)
   40 continue
c
      do 110 i=1,nev
*       xx(i)=tx0+alon(i)
*       yy(i)=ty0+alat(i)
c
        xx(i)=xp(i)
        yy(i)=yp(i)
  110 continue
      write(6,5) tend,am0
      write(2,2) nev,tx1,ty,am0,tend
      write(2,2) (i,xx(i),yy(i),am(i),tt(i), i=1,nev)
    1 format(3i10)
    2 format(i5,3f12.5,f14.5)
    3 format(8f10.5)
    4 format(2f10.1,i10)
    5 format(1h ,6f13.4)
      stop
      end
c
      subroutine pseudo(random)
c     wichmann+hill appl. statist (31) 188-? (1982)
      data ix,iy,iz /1993,1119,1046/    ! Seeds for different
!     data ix,iy,iz /1993,1119,1614/    ! uniform pseudo-random
c     data ix,iy,iz /1993,1119,1529/    ! number series.
      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+
     &           float(iz)/30323.0,1.0)
      return
      end
c-----------------------------------------------------------------------
      subroutine interp(nd,xx,yy,zxy,ldt,xi,yj,fij)
c  inputs:
c         nd: number of the delaunay triangulars.
c   xx(i)
c   yy(i)
c   zxy(i)
c   ldt(i,.): identification of the verteces of the i-th triangular.
c   xi 
c   yj
c  output
c     fij: values of the delaunay surface on the lattices (xi,yj)
c
      implicit real * 8 (a-h,o-z)
      parameter(npa=12345,nda=23456,mdl=33333,nr=030)
!     parameter(mdl=33333,npa=12345)
      dimension ldt(mdl,3)
      dimension xx(npa),yy(npa),zxy(npa)
      dimension xd(3),yd(3),zd(3)
      det(x1,x2,x3,y1,y2,y3)=x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2
      data eps/-0.1d-7/
      do ipt=1,nd
        do k=1,3
          xd(k)=xx(ldt(ipt,k))
          yd(k)=yy(ldt(ipt,k))
          zd(k)=zxy(ldt(ipt,k))
        end do
!       write(6,*) 'aa',nd,xd(1),xd(2),xd(3),yd(1),yd(2),yd(3)
!     stop 4444
        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 20
        aa1=det(xi,xd(2),xd(3),yj,yd(2),yd(3))/denom
        aa2=det(xd(1),xi,xd(3),yd(1),yj,yd(3))/denom
        aa3=det(xd(1),xd(2),xi,yd(1),yd(2),yj)/denom
c           if(a1.lt.0. .or. a2.lt.0. .or. a3.lt.0.) go to 20
      if(aa1.lt.eps.or. aa2.lt.eps.or. aa3.lt.eps) go to 20
        fij=aa1*zd(1)+aa2*zd(2)+aa3*zd(3)
c      write(6,2) ipt,aa1,aa2,aa3,fij,zd(1),zd(2),zd(3)
c      write(6,*) xi,yj
    2 format('ipt',i5,7f10.2,e20.10)
c      write(6,*) (ldt(ipt,k),k=1,3)
c      write(6,*) (xx(ldt(ipt,k)),k=1,3)
c      write(6,*) (yy(ldt(ipt,k)),k=1,3)
c      write(6,*) (zxy(ldt(ipt,k)),k=1,3)
        return
   20 continue
      end do
      write(6,*) 'error'
      stop 119
      end
