[R-SIG-Mac] Contour lines not plotting in correct position

[David Winsemius]

On Oct 24, 2012, at 7:42 PM, Thomas Veltman wrote:

> Hey gang!
> 
> This is my first post to the list, but I will try to be as precise as possible.  I have googled around and can't find an easy explanation for why this is happening, so I figured I had better get some help from the experts, so thanks in advance!  My code is supplied below.  I am running R 2.15.1 on Mountain Lion, 64-bit.  OSX is 10.8.2
> 
> I want to plot iso-efficiency lines for a hypothetical device.  I have a fully-parameterized expression, and have checked its sanity by hand, and it appears to operate as it ought to.  However, when I use the contour() function, the contour lines it produces don't lie in the correct place.  I generally believe their shape, but they appear to be translated by some amount.  I know that the analytical expression produces a given value, but when I look for that point on the plot, the iso-efficiency line doesn't cut through the point.  This is illustrated in the following code by the blue circle.  That point should fall on the red line, but clearly it does not.  Nevertheless, the governing equation (as defined in the R code) produces the correct result (if you check the array "Electrolyzer", the numbers match what you would get from simply applying the appropriate values to the expression, so it isn't a problem in the iterative population of the array).  I tried varying the mesh !
> size (i.e. the number of points in the vectors etaELEC and FE), but of course this didn't have any effect besides smoothing the curves out a bit more.  I also thought that it might be a problem with the png output, but I tried generating the figure in the quartz window, and the error persists.  I also tried closing R, not saving the workspace, reopening, making a completely new working directory, and then re-entering the functions.  In all cases, the plot still ends up funky.
> 
> Is there something fundamental that I am missing about the contour plot function?

It seems unlikely that the contour functions is performing improperly. 

> str(Electrolyzer)
 num [1:96, 1:96] 0.0035 0.0042 0.00489 0.00559 0.00629 ...

My wild guess (in the absence of any specifics about why you thought this was the wrong answer) is that this is due to the fact that matrix indices in R at 1-based, not 0-based as in other languages common among electrical engineers.

And why is this being posted to the Mac-SIG group, anyway? I see no Mac-specific content.

>  Also, for those playing along at home, the code includes an "easy" way to add subscripts after special characters in a plot title (which is typically a useful feature).  Relevant code is pasted below:
> 
> 
> Ecap <- function(etaELEC, FE, etaFC, etaRP, ENH3, etaCOMP, etaST) {
>             (etaELEC)/((1+ ((1-FE)/(FE*1.43) )*(1-etaST* etaELEC+ ((1.2*1.43*etaELEC)/(etaRP*ENH3))) + ((7.3*etaELEC)/(2*etaCOMP*ENH3)) + ((25*etaELEC)/(etaRP*ENH3))))}
> 
> 
> etaELEC <- seq(from=0.05,to=1, by=0.01)
> FE <- seq(from=0.05,to=1, by=0.01)
> 
> 
> Electrolyzer <- array(0, dim=c(length(etaELEC), length(FE)))
> for (i in 1:length(etaELEC)){
>      for (j in 1:length(FE)){
>         Electrolyzer[i,j] <- Ecap(FE=FE[j], etaST=0, etaFC=0, etaELEC=etaELEC[i], etaRP=0.6, etaCOMP=0.95 ,ENH3=340)}} 

> png(file="Standalone.png",bg="transparent",height=5400,width=5400,res=864)
> par(family="serif",tck=0.01)
> plot(0, type="n", xlim=c(0,1), ylim=c(0,1), cex.lab=1.2,
>      main="Energy Capture Coefficient for Standalone Device",
>      xlab="Electrolyzer Voltage Efficiency", ylab="Faradaic Efficiency", cex.sub=0.75,
>      sub=expression(paste(eta,""[FC],"=0, ",eta,""[Comp],"=0.95, ",eta,""[Sterling],"=0, ",eta,""[RP],"=0.60, ")))
> #make blank plot box with correct axes, etc...
> contour(x=seq(0, 1, length.out=(nrow(Electrolyzer))),
>         y=seq(0, 1, length.out=(ncol(Electrolyzer))), 
>         Electrolyzer, axes=FALSE,   
>         levels=c(-0.1, 0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1))
> contour(x=seq(0,1, length.out=(nrow(Electrolyzer))),
>         y=seq(0,1,length.out=(ncol(Electrolyzer))),
>         Electrolyzer,axes=FALSE,levels=c(0.102),col="red")
> points(0.5,0.1547,col="blue")
> dev.off()

-- 
David Winsemius, MD
Alameda, CA, USA