Revision 69b0f9dca8eb051f132725ecc679fe1997246e50 authored by Adrian Baddeley on 18 January 2006, 21:47:25 UTC, committed by cran-robot on 18 January 2006, 21:47:25 UTC
1 parent cb2215f
Raw File
methas.f
C Output from Public domain Ratfor, version 1.0
      subroutine methas(nmbr,aw,par,period,xprop,yprop,mprop,ntypes,ptyp
     *es,iseed, nrep,mrep,p,q,npmax,nverb,x,y,marks,aux,npts,fixall)
      implicit double precision(a-h,o-z)
      dimension par(1), ptypes(1), iseed(3)
      dimension x(1), y(1), marks(1), period(2)
      dimension xprop(1), yprop(1), mprop(1)
      integer aux(1)
      logical verb, marked, fixall
      eps = 2.22d-16
      one = 1.d0
      zero = 0.d0
      m1 = -1
      verb = .not.(nverb.eq.0)
      marked = ntypes .gt. 1
      irep = mrep
23000 if(irep .le. nrep)then
      if(verb .and. mod(irep,nverb).eq.0)then
      iprt = irep/nverb
      call intpr('irep/nverb=',-1,iprt,1)
      endif
      itype = 0
      call aru(1,zero,one,iseed,urp)
      if(urp.gt.p)then
      call aru(1,zero,one,iseed,urq)
      if(urq.gt.q)then
      u = xprop(irep)
      v = yprop(irep)
      if(marked)then
      mrk = mprop(irep)
      else
      mrk = -1
      endif
      call cif(nmbr,u,v,mrk,m1,x,y,marks,npts,ntypes, par,period,cifval,
     *aux)
      if(marked)then
      ffm = ptypes(mrk)
      else
      ffm = one
      endif
      a1 = q*aw*cifval/(ffm*(1-q)*(npts+1))
      a = min(one,a1)
      call aru(1,zero,one,iseed,bp)
      if(bp.lt.a)then
      npts = npts + 1
      if(npts .gt. npmax)then
      mrep = irep
      return
      endif
      itype = 1
      endif
      else
      call aru(1,zero,one,iseed,xi)
      ix = 1 + int(npts*xi)
      if(marked)then
      mrki = marks(ix)
      ffm = ptypes(mrki)
      else
      mrki = -1
      ffm = one
      endif
      call cif(nmbr,x(ix),y(ix),mrki,ix,x,y, marks,npts,ntypes,par,perio
     *d,cifval,aux)
      a1 = (one-q)*npts*ffm/(q*aw*cifval)
      a = min(one,a1)
      call aru(1,zero,one,iseed,dp)
      if(dp.lt.a)then
      itype = 2
      endif
      endif
      else
      call aru(1,zero,one,iseed,xi)
      ix = 1 + int(npts*xi)
      u = xprop(irep)
      v = yprop(irep)
      if(marked)then
      mrki = marks(ix)
      if(fixall)then
      mrk = mrki
      else
      mrk = mprop(irep)
      endif
      ffm = ptypes(mrki)/ptypes(mrk)
      else
      mrki = -1
      mrk = -1
      ffm = one
      endif
      call cif(nmbr,x(ix),y(ix),mrki,ix,x,y,marks,npts, ntypes,par,perio
     *d,cvd,aux)
      call cif(nmbr,u,v,mrk,ix,x,y,marks,npts, ntypes,par,period,cvn,aux
     *)
      if(cvd .lt. eps)then
      if(cvn .lt. eps)then
      goto 23000
      else
      a = one
      endif
      else
      a1 = ffm*cvn/cvd
      a = min(one,a1)
      endif
      call aru(1,zero,one,iseed,sp)
      if(sp.lt.a)then
      itype = 3
      endif
      endif
      if(itype .gt. 0)then
      if(nmbr .eq. 8)then
      call updaux(itype,x,y,u,v,npts,ix,par,period,aux)
      endif
      if(itype.eq.1)then
      ix = npts
      x(ix) = u
      y(ix) = v
      if(marked)then
      marks(ix) = mrk
      endif
      else
      if(itype.eq.2)then
      call death(x,y,marks,marked,npts,ix)
      else
      x(ix) = u
      y(ix) = v
      if(marked)then
      marks(ix) = mrk
      endif
      endif
      endif
      endif
      irep = irep+1
      goto 23000
      endif
23001 continue
      return
      end
back to top