How to elegantly determine the area of a hysteresis loop (inside/outside problem)? I have measured two parameters (Dissolved organic carbon DOC = y, and discharge = x). When these two variables are plotted against each other, we get a hysteresis loop  (see code example and picture).
Now, for further analysis, I would like to determine the area of this hysteretic loop. I figured out that this can be done using the Monte Carlo darting method. This method says that the area of an unknown area is proportional to the area of a known rectangular times the hits in the inside field (the loop). 
My problem now is, how to solve the inside / outside problem using R. How can I draw a rectangular with a known area and how can I excel the random hits inside and outside the hysteretic loop? 
Please note, that I am open to any other method...
I googled, and searched various statistical sites but could not find an answer. 
Any direct help or linkage to other websites/posts is greatly appreciated.

Data <- read.table("http://dl.dropbox.com/u/2108381/DOC_Q_hystersis.txt", sep = ";",
               header = T)

head(Data)
plot(Data$Q, Data$DOC, type = "o", xlab = "Discharge (m3 s-1)", ylab = "DOC (mg C l-1)",
 main = "Hystersis loop of the C/Q relationship")

 A: A completely different way would be to directly calculate the area of your polygon:
library(geometry)
polyarea(x=Data$Q, y=Data$DOC)

This yields 0.606.
A: One possibility would be this: it looks to me like the hysteresis loop should be convex, right? So one could generate points and test for each point whether it is part of the convex hull of the union of your dataset and the random point - if yes, it lies outside the original dataset. To speed up things, we can work with a subset of the original dataset, namely the points comprising its convex hull itself.
Data.subset <- Data[chull(Data$Q, Data$DOC),c("Q","DOC")]

x.min <- min(Data.subset$Q)
x.max <- max(Data.subset$Q)
y.min <- min(Data.subset$DOC)
y.max <- max(Data.subset$DOC)

n.sims <- 1000
random.points <- data.frame(Q=runif(n=n.sims,x.min,x.max),
  DOC=runif(n=n.sims,y.min,y.max))
hit <- rep(NA,n.sims)
for ( ii in 1:n.sims ) {
  hit[ii] <- !((nrow(Data.subset)+1) %in%
    chull(c(Data.subset$Q,random.points$Q[ii]),
      c(Data.subset$DOC,random.points$DOC[ii])))
}

points(random.points$Q[hit],random.points$DOC[hit],pch=21,bg="black",cex=0.6)
points(random.points$Q[!hit],random.points$DOC[!hit],pch=21,bg="red",col="red",cex=0.6)

estimated.area <- (y.max-y.min)*(x.max-x.min)*sum(hit)/n.sims


Of course, the R way of doing things would not use my for loop, but this is easy to understand and not too slow. I get an estimated area of 0.703.
EDIT: of course, this relies on the presumed convexity of the relationship. For instance, there appears to be a nonconvex part in the lower right. In principle, we could Monte-Carlo-estimate the area of such an area in the same way and subtract this from the original area estimate.
