[Stpp-commits] r57 - in pkg: man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 15 20:58:21 CEST 2012
Author: gabriele
Date: 2012-08-15 20:58:20 +0200 (Wed, 15 Aug 2012)
New Revision: 57
Added:
pkg/man/stpp-package.Rd
Removed:
pkg/man/stpp.Rd
Modified:
pkg/man/PCFhat.Rd
pkg/man/STIKhat.Rd
pkg/man/rinfec.Rd
pkg/man/rpp.Rd
pkg/src/pcffunction.f
Log:
Modified: pkg/man/PCFhat.Rd
===================================================================
--- pkg/man/PCFhat.Rd 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/man/PCFhat.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -5,7 +5,7 @@
correlation function.}
\usage{PCFhat(xyt, s.region, t.region, dist, times, lambda,
-ks="box", hs, kt="box", ht, lambda, correction = TRUE }
+ks="box", hs, kt="box", ht, correction = TRUE) }
\arguments{
\item{xyt}{coordinates and times (x,y,t) of the point pattern.}
@@ -94,21 +94,21 @@
\examples{\dontrun{
data(fmd)
data(northcumbria)
-FMD=as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
+FMD<-as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
Northcumbria=northcumbria/1000
# estimation of the temporal intensity
-Mt=density(FMD[,3],n=1000)
-mut=Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]
+Mt<-density(FMD[,3],n=1000)
+mut<-Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]
# estimation of the spatial intensity
-h = mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
-h = h$h[which.min(h$mse)]
-Ms = kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=5000, ny=5000)
-atx=findInterval(x=FMD[,1],vec=Ms$x)
-aty=findInterval(x=FMD[,2],vec=Ms$y)
-mhat=NULL
-for(i in 1:length(atx)) mhat=c(mhat,Ms$z[atx[i],aty[i]])
+h<-mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
+h<-h$h[which.min(h$mse)]
+Ms<-kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=5000, ny=5000)
+atx<-findInterval(x=FMD[,1],vec=Ms$x)
+aty<-findInterval(x=FMD[,2],vec=Ms$y)
+mhat<-NULL
+for(i in 1:length(atx)) mhat<-c(mhat,Ms$z[atx[i],aty[i]])
# estimation of the pair correlation function
g <- PCFhat(xyt=FMD, dist=1:20, times=1:20, lambda=lambda=mhat*mut/dim(FMD)[1], s.region=northcumbria/1000,t.region=c(1,200))
Modified: pkg/man/STIKhat.Rd
===================================================================
--- pkg/man/STIKhat.Rd 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/man/STIKhat.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -103,21 +103,21 @@
\examples{\dontrun{
data(fmd)
data(northcumbria)
-FMD=as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
+FMD<-as.3dpoints(fmd[,1]/1000,fmd[,2]/1000,fmd[,3])
Northcumbria=northcumbria/1000
# estimation of the temporal intensity
-Mt=density(FMD[,3],n=1000)
-mut=Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]
+Mt<-density(FMD[,3],n=1000)
+mut<-Mt$y[findInterval(FMD[,3],Mt$x)]*dim(FMD)[1]
# estimation of the spatial intensity
-h = mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
-h = h$h[which.min(h$mse)]
-Ms = kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=5000, ny=5000)
-atx=findInterval(x=FMD[,1],vec=Ms$x)
-aty=findInterval(x=FMD[,2],vec=Ms$y)
-mhat=NULL
-for(i in 1:length(atx)) mhat=c(mhat,Ms$z[atx[i],aty[i]])
+h<-mse2d(as.points(FMD[,1:2]), Northcumbria, nsmse=50, range=4)
+h<-h$h[which.min(h$mse)]
+Ms<-kernel2d(as.points(FMD[,1:2]), Northcumbria, h, nx=5000, ny=5000)
+atx<-findInterval(x=FMD[,1],vec=Ms$x)
+aty<-findInterval(x=FMD[,2],vec=Ms$y)
+mhat<-NULL
+for(i in 1:length(atx)) mhat<-c(mhat,Ms$z[atx[i],aty[i]])
# estimation of the STIK function
u <- seq(0,10,by=1)
Modified: pkg/man/rinfec.Rd
===================================================================
--- pkg/man/rinfec.Rd 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/man/rinfec.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -69,6 +69,7 @@
}
\examples{
+\dontrun{
# inhibition; spatial distribution: uniform
inf1 = rinfec(npoints=100, alpha=0.2, beta=0.6, gamma=0.5, maxrad=c(0.075,0.5), t.region=c(0,50), s.distr="uniform", t.distr="uniform", h="step", p="min", recent="all", inhibition=TRUE)
\dontrun{animation(inf1$xyt, cex=0.8, runtime=10)}
@@ -81,7 +82,7 @@
Ls = kernel2d(as.points(fmd[,1:2]), northcumbria, h, nx=50, ny=50)
inf2 = rinfec(npoints=100, alpha=4, beta=0.6, gamma=20, maxrad=c(12000,20), s.region=northcumbria, t.region=c(1,2000), s.distr="poisson", t.distr="uniform", h="step", p="min", recent=1, lambda=Ls$z, inhibition=FALSE)
-\dontrun{
+
image(Ls$x, Ls$y, Ls$z, col=grey((1000:1)/1000)); polygon(northcumbria,lwd=2)
animation(inf2$xyt, add=TRUE, cex=0.7, runtime=15)
}
Modified: pkg/man/rpp.Rd
===================================================================
--- pkg/man/rpp.Rd 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/man/rpp.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -8,8 +8,7 @@
\usage{
rpp(lambda, s.region, t.region, npoints=NULL, nsim=1, replace=TRUE,
- discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL,
- Lambda=NULL, ...)
+ discrete.time=FALSE, nx=100, ny=100, nt=100, lmax=NULL, ...)
}
\arguments{
Added: pkg/man/stpp-package.Rd
===================================================================
--- pkg/man/stpp-package.Rd (rev 0)
+++ pkg/man/stpp-package.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -0,0 +1,99 @@
+\name{stpp}
+\docType{package}
+\alias{stpp-package}
+
+\title{Space-Time Point Pattern simulation, visualisation and analysis}
+
+\description{ This package provides models of spatio-temporal
+point processes in a region S x T and statistical tools for
+analysing second-order properties of such processes. It also
+includes static and dynamic (2D and 3D) plots. \code{stpp} is
+the first dedicated unified computational environment in the
+area of spatio-temporal point processes.
+
+The \code{stpp} package depends upon some other packages:
+
+\code{splancs}: spatial and space-time point pattern analysis
+
+\code{rgl}: interactive 3D plotting of densities and surfaces
+
+\code{rpanel}: simple interactive controls for R using
+\code{tcltk} package
+}
+
+\details{
+\code{stpp} is a package for simulating, analysing and visualising
+patterns of points in space and time.
+
+Following is a summary of the main functions and the dataset in the
+\code{stpp} package.
+
+\emph{To visualise a spatio-temporal point pattern}
+
+\itemize{
+\item \code{\link{animation}}: space-time data animation.
+\item \code{\link{as.3dpoints}}: create data in spatio-temporal point
+ format.
+\item \code{\link{plot.stpp}}: plot spatio-temporal point object.
+Either a two-panels plot showing spatial locations and cumulative times,
+or
+a one-panel plot showing spatial locations with times treated as a
+quantitative mark
+attached to each location.
+\item \code{\link{stan}}: 3D space-time animation.
+}
+
+\emph{To simulate spatio-temporal point patterns}
+
+\itemize{
+\item \code{\link{rinfec}}: simulate an infection point process,
+\item \code{\link{rinter}}: simulate an interaction (inhibition or
+ contagious) point process,
+\item \code{\link{rlgcp}}: simulate a log-Gaussian Cox point process,
+\item \code{\link{rpcp}}: simulate a Poisson cluster point process,
+\item \code{\link{rpp}}: simulate a Poisson point process.
+}
+
+\emph{To analyse spatio-temporal point patterns}
+
+\itemize{
+\item \code{\link{PCFhat}}: space-time inhomogeneous pair correlation
+ function,
+\item \code{\link{STIKhat}}: space-time inhomogeneous K-function.
+}
+
+\emph{Dataset}
+
+\code{\link{fmd}}: 2001 food-and-mouth epidemic in north Cumbria (UK).
+}
+
+\author{
+Edith Gabriel <edith.gabriel at univ-avignon.fr>, Barry Rowlingson
+and Peter J Diggle. }
+
+\references{ Baddeley A., Moller J. and Waagepetersen R.
+(2000). Non- and semi-parametric estimation of interaction in
+inhomogeneous point patterns. Statistica Neerlandica, 54,
+329--350.
+
+Chan, G. and Wood A. (1997). An algorithm for simulating
+stationary Gaussian random fields. Applied Statistics,
+Algorithm Section, 46, 171--181.
+
+Chan, G. and Wood A. (1999). Simulation of stationary Gaussian
+vector fields. Statistics and Computing, 9, 265--268.
+
+Diggle P. , Chedwynd A., Haggkvist R. and Morris S. (1995).
+Second-order analysis of space-time clustering. Statistical
+Methods in Medical Research, 4, 124--136.
+
+Gabriel E., Diggle P. (2009) Second-order analysis of
+inhomogeneous spatio-temporal point process data. Statistica
+Neerlandica, 63, 43--51.
+
+Gneiting T. (2002). Nonseparable, stationary covariance functions
+for space-time data. Journal of the American Statistical Association,
+97, 590--600.
+}
+
+
Deleted: pkg/man/stpp.Rd
===================================================================
--- pkg/man/stpp.Rd 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/man/stpp.Rd 2012-08-15 18:58:20 UTC (rev 57)
@@ -1,85 +0,0 @@
-\name{stpp-package}
-\alias{stpp}
-\title{Space-Time Point Pattern simulation, visualisation and analysis}
-
-\description{
-This package provides models of spatio-temporal point processes in a region S x T
-and statistical tools for analysing second-order properties of such processe.
-It also inludes static and dynamic (2D and 3D) plots. \code{stpp} is the first dedicated
-unified computational environment in the area of spatio-temporal point processes.
-}
-
-\details{
-\code{stpp} is a package for simulating, analysing and visualising patterns of points in space and time.
-
-Following is a summary of the main functions and the dataset in the stpp package.
-
-\emph{To visualise a spatio-temporal point pattern}
-
-\itemize{
-\item \code{\link{animation}}: space-time data animation,
-\item \code{\link{as.3dpoints}}: create data in spatio-temporal point format,
-\item \code{\link{plot.stpp}}: plot spatio-temporal point object.,
-Either a two-panels plot showing spatial locations and cumulative times, or
-a one-panel plot shwoing spatial locations with times treated as a quantitative mark
-attached to each location.
-\item \code{\link{stan}}: 3D space-time animation.
-}
-
-\emph{To simulate spatio-temporal point patterns}
-
-\itemize{
-\item \code{\link{rinfec}}: simulate an infection point process,
-\item \code{\link{rinter}}: simulate an interaction (inhibition or contagious) point process,
-\item \code{\link{rlgcp}}: simulate a log-Gaussian Cox point process,
-\item \code{\link{rpcp}}: simulate a Poisson cluster point process,
-\item \code{\link{rpp}}: simulate a Poisson point process.
-}
-
-\emph{To analyse spatio-temporal point patterns}
-
-\itemize{
-\item \code{\link{PCFhat}}: pair correlation function,
-\item \code{\link{STIKhat}}: space-time inhomogeneous K-function.
-}
-
-\emph{Dataset}
-
-\code{\link{fmd}}: 2001 food-and-mouth epidemic in north Cumbria (UK).
-
-}
-
-\dependencies{
-The \code{stpp} package depends upon some other packages:
-
-\code{splancs} - Spatial and space-time point pattern analysis,
-
-\code{rgl} - Interactive 3D plotting of densities and surfaces,
-
-\code{rpanel} - Simple interactive controls for R using the \code{tcltk} package.
-}
-
-\author{
-Edith Gabriel <edith.gabriel at univ-avignon.fr>, Barry Rowlingson and Peter J Diggle.
-}
-
-\references{
-Baddeley A., Moller J. and Waagepetersen R. (2000). Non- and semi-parametric estimation
-of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54, 329--350.
-
-Chan, G. and Wood A. (1997). An algorithm for simulating stationary
-Gaussian random fields. Applied Statistics, Algorithm Section, 46, 171--181.
-
-Chan, G. and Wood A. (1999). Simulation of stationary Gaussian vector fields.
-Statistics and Computing, 9, 265--268.
-
-Diggle P. , Chedwynd A., Haggkvist R. and Morris S. (1995). Second-order analysis
-of space-time clustering. Statistical Methods in Medical Research, 4, 124–-136.
-
-Gabriel E., Diggle P. (2009) Second-order analysis of inhomogeneous spatio-temporal point process data.
-Statistica Neerlandica, 63, 43--51.
-
-Gneiting T. (2002). Nonseparable, stationary covariance functions
-for space-time data. Journal of the American Statistical Association,
-97, 590--600.
-}
Modified: pkg/src/pcffunction.f
===================================================================
--- pkg/src/pcffunction.f 2012-08-15 13:40:25 UTC (rev 56)
+++ pkg/src/pcffunction.f 2012-08-15 18:58:20 UTC (rev 57)
@@ -92,18 +92,6 @@
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c functions called by pcffunction:
-c -----------------------------------------
-c
-c * boxkernel
-c * iplace
-c * weight
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-
c--------------------------------------------------------------------
c
c boxkernel
@@ -188,388 +176,3 @@
return
end
-
-
-c--------------------------------------------------------------------
-c
-c iplace
-c
-c--------------------------------------------------------------------
-
- function iplace(s,ns,t)
-c
-c which of the variable width bins s is t in?
-c
- implicit real*8 (a-h,o-z)
-
- dimension s(ns)
-
- do ib=1,ns
- if(s(ib).ge.t)then
- iplace=ib
- return
- end if
- end do
-c
-c if it is outside the range of s
-c
- iplace=ns+1
-
- return
- end
-
-c--------------------------------------------------------------------
-c
-c weight
-c
-c--------------------------------------------------------------------
-
- function weight(x,y,r,xp,yp,np)
-c
-c find the weight for the point at x,y, radius r
-c
- implicit real*8 (a-h,o-z)
-
-c include 'bounds.cmn'
-
- dimension xp(np+1),yp(np+1)
-
- weight=cncvwt(x,y,r,xp,yp,np)
-
- return
- end
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c function called by weight:
-c --------------------------
-c
-c * cncvwt
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-c--------------------------------------------------------------------
-c
-c cncvwt
-c
-c--------------------------------------------------------------------
-
- function cncvwt(x,y,r,xp,yp,np)
-c
-c compute the weight given to a point at x,y according to how much
-c of a circle of radius r is inside the bounding polygon
-c
-c
- implicit real*8 (a-h,o-z)
-c include 'bounds.cmn'
- dimension xp(np+1),yp(np+1)
- parameter(pi=3.141592654d0)
-c store circle/poly intersections here
- parameter(maxcrs=40)
- dimension cross(maxcrs+1)
- parameter(tiny=1.0e-7)
-c set count of crossing points to zero
- ncross = 0
-
-c first loop over the boundary and find the crossing points
- do ic=1,np
-c work with the trial point at origin
- x1=xp(ic)-x
- y1=yp(ic)-y
- x2=xp(ic+1)-x
- y2=yp(ic+1)-y
-
- cx=x2-x1
- cy=y2-y1
-
-c these are the coefficients of the quadratic giving the
-c intercept of line and circle.
- a=cx*cx+cy*cy
- b=2*(x1*cx+y1*cy)
- c=x1*x1+y1*y1-r*r
-
-c find out if real solutions exist...
- b2m4ac=b*b-4*a*c
-
-c ... and if they do, find them.
- if (b2m4ac.ge.0) then
- t1=(-b+sqrt(b2m4ac))/(2*a)
- t2=(-b-sqrt(b2m4ac))/(2*a)
-
-c see if the solutions lie in the line segments
- if ((t1.gt.tiny).and.(t1-1.0.le.tiny)) then
- ncross=ncross+1
-c find the angle to this point on the circle
- ctemp=atan2(y1+t1*cy,x1+t1*cx)
- if(ctemp.lt.0)ctemp=2*pi+ctemp
- cross(ncross)=ctemp
-c check crossing of circle with vertex
- else if (abs(t1).le.tiny) then
-c compare this polygon segment's direction with that of the
-c previous one
- nprev = (mod((ic+ (np-2)),np)+1)
- x0 = xp(nprev) - x
- y0 = yp(nprev) - y
- idp1 = isig8((x2-x1)*x1+ (y2-y1)*y1,tiny)
- idp2 = isig8((x1-x0)*x1+ (y1-y0)*y1,tiny)
-c see if the polygon passes through the circle here
- if ((idp1-idp2).ne.1 .and.
- + abs(idp1+idp2).ne.2) then
- ncross = ncross + 1
- ctemp = atan2(y1+t1*cy,x1+t1*cx)
- if (ctemp.lt.0.0) ctemp = 2*pi + ctemp
- cross(ncross) = ctemp
- end if
- end if
-
- if ((t2.gt.tiny).and.(t2-1.0.lt.tiny)) then
- ncross=ncross+1
- ctemp=atan2(y1+t2*cy,x1+t2*cx)
- if(ctemp.lt.0)ctemp=2*pi+ctemp
- cross(ncross)=ctemp
-c check crossing of circle with vertex
- else if (abs(t2).le.tiny)then
-c compare this polygon segment's direction with that of the
-c previous one
- nprev = (mod((ic+ (np-2)),np)+1)
- x0 = xp(nprev) - x
- y0 = yp(nprev) - y
- idp1 = isig8((x2-x1)*x1+ (y2-y1)*y1,tiny)
- idp2 = isig8((x1-x0)*x1+ (y1-y0)*y1,tiny)
-c see if the polygon passes through the circle here
- if ((idp1-idp2).ne.1 .and.
- + abs(idp1+idp2).ne.2) then
- ncross = ncross + 1
- ctemp = atan2(y1+t2*cy,x1+t2*cx)
- if (ctemp.lt.0.0) ctemp = 2*pi + ctemp
- cross(ncross) = ctemp
- end if
- end if
- end if
- end do
-
-c now we have all the crossing point angles stored in
-c cross(1:ncross)
-
-c if ncross = 0 then the total angle within the poly is 2*pi
-c unless the circle is large and spans the polygon. this should
-c be checked beforehand so it's okay to assume 2*pi here.
-
- if (ncross.eq.0) then
- totang=2*pi
- else
-
-c sort into ascending order
- call sort2(cross,ncross)
-
-c fix the ncross+1'th element to be the first plus 2pi so that
-c the list is circular...
- cross(ncross+1)=cross(1)+2*pi
-
-c check that the number of crossings is even - if not then error.
- if (mod(ncross,2).ne.0) then
- cncvwt=-1
- return
- end if
-c now find a nice spot to do the point-in-poly search
- sepmax=0.0
- icm=0
-
- do ic=1,ncross
- if (cross(ic+1)-cross(ic).gt.sepmax) then
- sepmax=cross(ic+1)-cross(ic)
- icm=ic
- end if
- end do
-
-c icm is now the index of the crossing with the largest gap
-c between it and the next crossing point
-
-c test for point in poly of the point on the circle between these
-c points angtes=(cross(icm)+cross(icm+1))/2.
-
- xtest=x+r*cos(angtes)
- ytest=y+r*sin(angtes)
-
-c find out if test point is in the polygon boundary
- linpol=ipippa(xtest,ytest,xp,yp,np)
-
-c find the total angle between (odd-even) crossings
-c (i.e. 1-2 + 3-4 + ...)
- totang = 0.
- do ic=1,ncross-1,2
- totang = totang + (cross(ic+1)-cross(ic))
- end do
-
-c If the point we tested for p-i-p was on an odd-even
-c section and was in the poly, then totang is the amount of circle
-c inside the polygon. if the point was outside the polygon, then
-c we need to subtract totang from 2*pi radians to get the angle
-c inside the polygon. conversely, if the point tested was between
-c even-odd crossings and outside the polygon, then totang is the
-c angle we want, and if inside the polygon then again we have to
-c do 2*pi-totang
-
- if ( (((mod(icm,2).eq.1).and.(linpol.eq.0)) .or.
- & ((mod(icm,2).eq.0).and.(linpol.eq.1)) ) ) then
- totang = 2*pi-totang
- end if
-
- end if
-c now totang is the angle contained in the polygon
-
-c weight is proportion of total angle in the poly
- cncvwt = (2*pi)/(totang)
- return
- end
-
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c
-c functions called by cncvwt:
-c ---------------------------
-c
-c * isig8
-c * sort2
-c * ipippa
-c
-cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-c--------------------------------------------------------------------
-c
-c isig8
-c
-c--------------------------------------------------------------------
-
- integer function isig8(value,tiny)
-
-c return the sign (+1,0,-1) of a value
-
- real*8 tiny,value
- if (value.gt.tiny) then
- isig8 = 1
- else if (value.lt.-tiny) then
- isig8 = -1
- else
- isig8 = 0
- end if
- return
- end
-
-c--------------------------------------------------------------------
-c
-c sort2
-c
-c--------------------------------------------------------------------
-
- subroutine sort2(x,n)
-c
-c shellsort algorithm
-c n : number of elements to be sorted
-c x : on enter an array of dimension at least n containing
-c real numbers
-c on output first n elements of x are sorted from smallest
-c to largest
-c
- implicit real*8 (a-h,o-z)
- dimension x(n)
- i=1
- 1 i=i+1
- if (i.le.n) goto 1
- m=i-1
- 2 m=m/2
- if (m.eq.0) return
- k=n-m
- do 4 j=1,k
- kk=j
- 3 if (kk.lt.1) goto 4
- if (x(kk+m).ge.x(kk)) goto 4
- w=x(kk+m)
- x(kk+m)=x(kk)
- x(kk)=w
- kk=kk-m
- goto 3
- 4 continue
- goto 2
- end
-
-
-c--------------------------------------------------------------------
-c
-c ipippa
-c
-c--------------------------------------------------------------------
-
- function ipippa(x,y,xc,yc,nc)
-c
-c point in polygon routine.
-c
-c returns 0 if point x,y not in the bound polygon defined by xc,yc
-c
-c fortran version of C routine by Ken McElvain
-c
-
- implicit real*8 (a-h,o-z)
-c include 'bounds.cmn'
-
-
- dimension xc(nc+1),yc(nc+1)
-
- iwind = 0
- xlastp = xc(nc)
- ylastp = yc(nc)
- ioldq = iquad(xlastp,ylastp,x,y)
- do i=1,nc
-c for each point in the polygon
- xthisp=xc(i)
- ythisp=yc(i)
- inewq = iquad(xthisp,ythisp,x,y)
- if(ioldq.ne.inewq) then
- if(mod(ioldq+1,4).eq.inewq) then
- iwind=iwind+1
- else if(mod(inewq+1,4).eq.ioldq) then
- iwind = iwind - 1
- else
- a = (ylastp-ythisp)*(x-xlastp)
- b = xlastp-xthisp
- a = a + ylastp * b
- b=b*y
- if (a.gt.b) then
- iwind=iwind+2
- else
- iwind=iwind-2
- end if
- end if
- end if
- xlastp=xthisp
- ylastp=ythisp
- ioldq=inewq
- end do
-c
-c quadrant winding is either -4,0,+4 so divide down and take abs.
-c
- ipippa = abs(iwind/4)
-
- end
-
- function iquad(xp,yp,xo,yo)
-c
-c determine which quadrant xp,yp is in relative to xo,yo as origin
-c
- implicit real*8 (a-h,o-z)
-
- if(xp.lt.xo)then
- if(yp.lt.yo) then
- iquad=2
- else
- iquad=1
- end if
- else
- if(yp.lt.yo)then
- iquad = 3
- else
- iquad = 0
- end if
- end if
-
- return
- end
More information about the Stpp-commits
mailing list