swh:1:snp:bbc0b94b1c0463480b615c86f1ae2df3fcc700c0
Tip revision: 72984f1
polar2cart.R
``````## =============================================================================
## Maps a matrix u from (x,y) to (x2,y2) by 2-D linear interpolation
## =============================================================================
mapxy <- function(x, y, x2, y2, u) {
Nx <- length(x)
Ny <- length(y)

if (min(x2) < min(x)) stop("'x' should embrace 'x2'")
if (min(y2) < min(y)) stop("'y' should embrace 'y2'")

dx   <- c(diff(x),1)  # 1= for last value
dy <- c(diff(y),1)

Du <- dim(u)
if (Du[1] != Nx)
stop("u and x not compatible")
if (Du[2] != Ny)
stop("u and y not compatible")

Transf <- function (x2,y2,u) {

# find embracing values : first interval
ix <- findInterval(x2, x )
iy <- findInterval(y2, y)

#    ix <- pmax(ix, 1)
#    iy <- pmax(iy, 1)

# next interval
ixp1 <- pmin(ix+1,Nx)
iyp1 <- pmin(iy+1,Ny)

# interpolation factor
xfac <- (x2-x[ix])/dx[ix]
yfac <- (y2-y[iy])/dy[iy]

# interpolate
(1-yfac)*((1-xfac)*u[cbind(ix,iy)]+xfac*u[cbind(ixp1,iy)]) +
yfac*((1-xfac)*u[cbind(ix,iyp1)]+xfac*u[cbind(ixp1,iyp1)])

} # end Transf

outer(x2, y2, FUN = Transf, u = u)
}

#mapxy(1:87, 1:61, seq(1,87,0.25), seq(1,61,0.25),volcano)->VOLCANO
#IM <- matrix(nrow=3,ncol=3,c(1,1,1,1,2,1,1,1,1))
#mapxy(1:3,1:3,seq(0,3,0.1),seq(0,3,0.1),IM)->im2

## =============================================================================
## From polar to cartesian coordinates
## =============================================================================

polar2cart <- function(out, r, theta, x = NULL, y = NULL)
{
if (is.null(theta)) {   # one-D mapping to 2-D
if (!is.vector(out))
stop ("'out' should be a vector if 'theta' = NULL")
afun <- approxfun(r, out, yleft = out[1])
r2 <- matrix(nrow = length(x), ncol = length(y),
afun(outer (FUN = function(x, y) sqrt(x^2 + y^2), X = x, Y = y)))
return(r2)

}

classout <- class(out)[1]
if (!classout %in% c("steady2D", "deSolve", "matrix"))
stop ("Class of 'out' should be one of steady2D, deSolve or matrix")

ATTR <- attributes(out)
if (length(ATTR\$dimens)!= 2)
stop ("input should be 2-dimensional")

Nr   <- length(r)  -1
Np   <- length(theta)-1
maxr   <- max(r)
maxtheta <- max(theta)
minr   <- min(r)
mintheta <- min(theta)

r      <- c(r[1],   0.5*(r[-1]+r[-(Nr+1)]), r[Nr+1])
theta  <- c(theta[1], 0.5*(theta[-1]+theta[-(Np+1)]), theta[Np+1])
# check dimensions...
dr     <- c(diff(r),1)       # 1= for last value: no interpolation
dtheta <- c(diff(theta),1)

if (is.null (x))
if (! is.null(y)) x<-y
if (is.null (y)) {
mr <- max(abs(r))
x <- y <- seq(-mr, mr, len=Nr)
}

Transf <- function(x, y,  u) {

Du <- dim(u)
if (Du[1] != Nr)
stop("u and r not compatible")
if (Du[2] != Np)
stop("u and theta not compatible")
# augment u with boundary values
u    <- rbind(c(u[1,1],u[1,], u[1,Np]) ,
cbind(u[,1],u, u[,Np]),
c(u[Nr,1],u[Nr,], u[Nr,Np]) )

R   <- sqrt(x^2+y^2)
Theta <- atan(y/x)
Theta[x<0]           <- Theta[x < 0] + pi
Theta[x>0 & y < 0]   <- Theta[x > 0 & y < 0] + 2*pi
Theta[x==0 & y > 0]  <- pi/2
Theta[x==0 & y < 0]  <- 3*pi/2
Theta[x==0 & y == 0] <- 0

# find embracing values : interval location
iR <- findInterval(R, r )
iP <- findInterval(Theta,theta)
iR [iR == 0]   <- NA
iR [iR > Nr+1] <- NA
iP [iP == 0]   <- NA
iP [iP > Np+1] <- NA

# next interval
iPp1 <- pmin(iP+1,Np+1)
iRp1 <- pmin(iR+1,Nr+1)

# interpolation factor
iRfac <- (R-r[iR])/dr[iR]
iPfac <- (Theta-theta[iP])/dtheta[iP]

# interpolate inbetween 4 values
(1-iPfac)*((1-iRfac)*u[cbind(iR,iP)]+iRfac*u[cbind(iRp1,iP)]) +
iPfac*((1-iRfac)*u[cbind(iR,iPp1)]+iRfac*u[cbind(iRp1,iPp1)])

} # end Transf

Num   <- prod(ATTR\$dimens)
nspec <- ATTR\$nspec
MAP   <- list()
MAP\$y <- NULL
for (i in 1:nspec)  {
ii <- (1+(i-1)*Num):(i*Num)
U <- matrix(out\$y[ii], nrow = ATTR\$dimens[1], ncol= ATTR\$dimens[2])
MAP\$y <- c(MAP\$y, as.vector(outer(x, y, FUN = Transf, u = U)))
}
attributes(MAP) <- ATTR
} else {               # Dynamic output
MAP <- NULL
nr = nrow(out)
for ( j in 1: nr) {
mm <- out[j,1]
for (i in 1:nspec)  {
ii <- ((1+(i-1)*Num):(i*Num))+1
U <- matrix(out[j,ii], nrow = ATTR\$dimens[1], ncol= ATTR\$dimens[2])
mm <- c(mm, as.vector(outer(x, y, FUN = Transf, u = U)))
}
MAP <- rbind(MAP, mm)
}
attributes(MAP)[c("istate","rstate","class","type","nspec","dimens")] <-
ATTR[c("istate","rstate","class","type","nspec","dimens")]
}
attributes(MAP)\$dimens     <- c(length(x), length(y))
attributes(MAP)\$xgrid      <- x
attributes(MAP)\$ygrid      <- y

MAP
}
``````