Nothing
##==============================================================================
## Diffusive transport in polar coordinates (r, theta)
##==============================================================================
tran.polar <- function(C, C.r.up=NULL, C.r.down=NULL,
C.theta.up=NULL, C.theta.down=NULL,
flux.r.up=NULL, flux.r.down=NULL,
flux.theta.up=NULL, flux.theta.down=NULL,
cyclicBnd = NULL,
D.r=1, D.theta=D.r,
r=NULL, theta=NULL,
full.output = FALSE )
{
# ------------------------------------------------------------------------------
# Initialisation
# ------------------------------------------------------------------------------
# a function to check the dimensionality of the system
checkDim <- function (dr, dtheta, nx="dr", ntheta="dtheta") {
if (length(dr) != 1 )
stop (nx,", should have length 1")
if (length(dtheta) != 1)
stop (ntheta,", should have length 1")
}
checkDim2 <- function (dr,dtheta,nx="dr",ntheta="dtheta") {
if (length(dr) != dimC[1]+1)
stop (nx,", should have length equal to dim(C)[1]+1 = ",dimC[1]+1," it is ",length(dr))
if (length(dtheta) != dimC[2]+1)
stop (ntheta,", should have length equal to dim(C)[2]+1 = ",dimC[2]+1," it is ",length(dtheta))
}
dimC <- dim (C)
Type <- 2
if (max(theta)> 2*pi)
stop ("theta should be < 2pi")
if (min(theta)< 0)
stop ("theta should be > 0 ")
checkDim2(r, theta,"r, the grid in x-direction",
"theta, the grid in theta-direction")
# central values
Nr <- dimC[1]
r.c <- 0.5*(r[-1]+r[-(Nr+1)])
dr <- diff(r)
drint <- c(r.c[1]-r[1],diff(r.c),r[Nr+1]-r.c[Nr])
divr <- 1/r
divr[is.na(divr)] <- 0
Ntheta <- dimC[2]
theta.c <- 0.5*(theta[-1]+theta[-(Ntheta +1)])
dtheta <- diff(theta)
dthetaint <- c(theta.c[1]-theta[1],diff(theta.c),theta[Ntheta+1]-theta.c[Ntheta])
checkDim (D.r,D.theta,"D.r, the diffusion in x-direction ",
"D.theta, the diffusion in theta-direction ")
# check the dimensionality of the boundaries
if (! is.null(cyclicBnd)) {
# check if other boundaries not prescribed in this direction
}
if (! is.null(C.r.up)){
if( length(C.r.up) != 1 && length(C.r.up) != Ntheta)
stop ("'C.r.up' should have length 1 or equal to ", Ntheta)
} else if (! is.null(flux.r.up)) {
if( length(flux.r.up) != 1 && length(flux.r.up) != Ntheta)
stop ("'flux.r.up' should have length 1 or equal to ", Ntheta)
} else if (!1 %in% cyclicBnd)
C.r.up = C[1,]
if (! is.null(C.r.down)) {
if( length(C.r.down) != 1 && length(C.r.down) != Ntheta)
stop ("'C.r.down' should have length 1 or equal to ", Ntheta)
} else if (! is.null(flux.r.down)){
if(length(flux.r.down) != 1 && length(flux.r.down) != Ntheta)
stop ("'flux.r.down' should have length 1 or equal to ", Ntheta)
} else if (!1 %in% cyclicBnd)
C.r.down = C[Nr,]
if (! is.null(C.theta.up)) {
if ( length(C.theta.up) != 1 && length(C.theta.up) != Nr)
stop ("'C.theta.up' should have length 1 or equal to ", Nr)
} else if (! is.null(flux.theta.up)){
if (length(flux.theta.up) != 1 && length(flux.theta.up) != Nr)
stop ("'flux.theta.up' should have length 1 or equal to ", Nr)
} else if (!2 %in% cyclicBnd)
stop ("'flux.theta.up' OR 'C.theta.up' should be specified")
if (! is.null(C.theta.down)) {
if(length(C.theta.down) != 1 && length(C.theta.down) != Nr)
stop ("'C.theta.down' should have length 1 or equal to ", Nr)
} else if (! is.null(flux.theta.down)) {
if( length(flux.theta.down) != 1 && length(flux.theta.down) != Nr)
stop ("'flux.theta.down' should have length 1 or equal to ", Nr)
} else if (!2 %in% cyclicBnd)
stop ("'flux.theta.down' OR 'C.theta.down' should be specified")
# ------------------------------------------------------------------------------
# Function body: calculation
# ------------------------------------------------------------------------------
## Initialise 'fluxes' in all directions
Flux.r <- 0
Flux.theta <- 0
## First rewrite boundary values as "concentration"
## then perform diffusive transport
if (1 %in% cyclicBnd) {
C.r.up <- (C[1,]*drint[Nr+1] + C[Nr,]*drint[1]) /(drint[1]+drint[Nr+1])
C.r.down <- C.r.up
}
if (is.null(C.r.up))
C.r.up <- flux.r.up / D.r * drint[1]+ C[1,]
if (is.null(C.r.down))
C.r.down <- - flux.r.down / D.r * drint[Nr+1] + C[Nr,]
Flux.r <- D.r * (rbind(C.r.up,C, deparse.level = 0) -
rbind(C,C.r.down, deparse.level = 0))/drint
if (2 %in% cyclicBnd) {
C.theta.up <- (C[,1]*dthetaint[Ntheta+1] + C[,Ntheta]*dthetaint[1]) /
(dthetaint[1]+dthetaint[Ntheta+1])
C.theta.down <- C.theta.up
}
if (is.null(C.theta.up))
C.theta.up <- flux.theta.up/D.theta * dthetaint[1]*r.c + C[,1]
if (is.null(C.theta.down))
C.theta.down <- - flux.theta.down /D.theta * dthetaint[Ntheta+1]*r.c + C[,Ntheta]
Flux.theta <- D.theta * (cbind(C.theta.up,C,deparse.level = 0) -
cbind(C,C.theta.down, deparse.level = 0))/
matrix(data= dthetaint,nrow=Nr,ncol=(Ntheta+1),byrow=TRUE) /r.c
## Calculate rate of change = flux gradient
dFdtheta <- 0.
dFdr <- -diff(Flux.r * r)/dr/r.c
dFdtheta <- -t(diff(t(Flux.theta))/dtheta)/r.c
if (full.output)
return (list (dC = dFdr + dFdtheta , # Rate of change
C.r.up = C.r.up,
C.r.down = C.r.down,
C.theta.up = C.theta.up,
C.theta.down = C.theta.down,
r.flux = Flux.r,
theta.flux = Flux.theta,
flux.r.up = Flux.r[1,],
flux.r.down = Flux.r[Nr+1,],
flux.theta.up = Flux.theta[,1],
flux.theta.down = Flux.theta[,Ntheta+1]
))
else
return (list (dC = dFdr + dFdtheta ,
flux.r.up = Flux.r[1,],
flux.r.down = Flux.r[Nr+1,],
flux.theta.up = Flux.theta[,1],
flux.theta.down = Flux.theta[,Ntheta+1]
))
} # end tran.polar
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.