https://github.com/cran/ReacTran
Tip revision: e114557
tran.cylspher.R
``````##==============================================================================
## Diffusive transport in cylindrical coordinates (r, theta, z)
##==============================================================================

tran.cylindrical <- function(C, C.r.up = NULL, C.r.down = NULL,
C.theta.up = NULL, C.theta.down = NULL,
C.z.up = NULL, C.z.down = NULL,
flux.r.up = NULL, flux.r.down = NULL,
flux.theta.up = NULL, flux.theta.down = NULL,
flux.z.up = NULL, flux.z.down = NULL,
cyclicBnd = NULL,
D.r = NULL, D.theta = D.r, D.z = D.r,
r = NULL,  theta = NULL, z = NULL)
{
# ------------------------------------------------------------------------------
# Initialisation
# ------------------------------------------------------------------------------

dimens <- dim(C)
if (length(dimens) != 3)
stop("'C' should an array of dimensionality 3")

Nr     <- dimens[1]
Ntheta <- dimens[2]
Nz     <- dimens[3]

if (length(r) != Nr+1)
stop("Length of 'r' should equal first dimension of 'C' + 1")
if (length(theta) != Ntheta+1)
stop("Length of 'theta' should equal second dimension of 'C' + 1")
if (length(z) != Nz+1)
stop("Length of 'z' should equal third dimension of 'C' + 1")

if (max(theta) > 2 * pi)
stop("theta should be < 2pi")
if (min(theta) < 0)
stop("theta should be > 0 ")

# boundary conditions default = zero-gradient
Bc.r.Up       <- 3
Bc.r.Down     <- 3
Bc.theta.Up   <- 3
Bc.theta.Down <- 3
Bc.z.Up       <- 3
Bc.z.Down     <- 3

# boundaries in r-direction
if (! is.null(C.r.up))
Bc.r.Up <- 2
else
C.r.up <- 0.

if (! is.null(flux.r.up))
Bc.r.Up <- 1
else
flux.r.up <- 0.

if (! is.null(C.r.down))
Bc.r.Down <- 2
else
C.r.down <- 0.

if (! is.null(flux.r.down))
Bc.r.Down <- 1
else
flux.r.down <- 0.

if (1 %in% cyclicBnd) {
Bc.r.Up   <- 5
Bc.r.Down <- 5
}

# boundaries in theta-direction
if (! is.null(C.theta.up))
Bc.theta.Up <- 2
else
C.theta.up <- 0.

if (! is.null(flux.theta.up))
Bc.theta.Up <- 1
else
flux.theta.up <- 0.

if (! is.null(C.theta.down))
Bc.theta.Down <- 2
else
C.theta.down <- 0.

if (! is.null(flux.theta.down))
Bc.theta.Down <- 1
else
flux.theta.down <- 0.

if (2 %in% cyclicBnd) {
Bc.theta.Up   <- 5
Bc.theta.Down <- 5
}

# boundaries in z-direction
if (! is.null(C.z.up))
Bc.z.Up <- 2
else
C.z.up <- 0.

if (! is.null(flux.z.up))
Bc.z.Up <- 1
else
flux.z.up <- 0.

if (! is.null(C.z.down))
Bc.z.Down <- 2
else
C.z.down <- 0.

if (! is.null(flux.z.down))
Bc.z.Down <- 1
else
flux.z.down <- 0.

if (3 %in% cyclicBnd) {
Bc.phi.Up   <- 5
Bc.phi.Down <- 5
}

# Diffusion coefficients
D.r <- rep(D.r, length.out=Nr+1)
D.theta <- rep(D.theta, length.out=Ntheta+1)
D.z <- rep(D.z, length.out=Nz+1)

if (length(C.r.up) + length(flux.r.up) +
length(C.r.down) + length(flux.r.down) +
length(C.theta.up) + length(flux.theta.up)+
length(C.theta.down) + length(flux.theta.down)+
length(C.z.up) + length(flux.z.up) +
length(C.z.down)+length(flux.z.down)!=12)
stop("length of all 'boundary conditions' should be 1")

tr <- .Fortran("diffcyl", as.integer(Nr), as.integer(Ntheta),
as.integer(Nz), as.double(C),
as.double(C.r.up), as.double(C.r.down),
as.double(C.theta.up), as.double(C.theta.down),
as.double(C.z.up), as.double(C.z.down),
as.double(flux.r.up),  as.double(flux.r.down),
as.double(flux.theta.up),  as.double(flux.theta.down),
as.double(flux.z.up),  as.double(flux.z.down),
as.integer(Bc.r.Up), as.integer(Bc.r.Down),
as.integer(Bc.theta.Up), as.integer(Bc.theta.Down),
as.integer(Bc.z.Up), as.integer(Bc.z.Down),
as.double(D.r), as.double(D.theta), as.double(D.z),
as.double(r), as.double(theta), as.double(z),
Fluxrup= numeric(Ntheta*Nz),Fluxrdown= numeric(Ntheta*Nz),
Fluxthetaup= numeric(Nr*Nz),Fluxthetadown= numeric(Nr*Nz),
FluxZup= numeric(Nr*Ntheta),FluxZdown= numeric(Nr*Ntheta),
dC = numeric(Nr*Ntheta*Nz))

list(dC = array(dim=c(Nr, Ntheta, Nz), tr\$dC),
flux.r.up   = matrix(nr=Ntheta, nc = Nz, tr\$Fluxrup),
flux.r.down = matrix(nr=Ntheta, nc = Nz, tr\$Fluxrdown),
flux.theta.up   = matrix(nr=Nr, nc = Nz, tr\$Fluxthetaup),
flux.theta.down = matrix(nr=Nr, nc = Nz, tr\$Fluxthetadown),
flux.z.up   = matrix(nr=Nr, nc = Ntheta, tr\$FluxZup),
flux.z.down = matrix(nr=Nr, nc = Ntheta, tr\$FluxZdown)
)
}

##==============================================================================
## Diffusive transport in spherical coordinates (r, theta, phi)
##==============================================================================

tran.spherical <- function(C, C.r.up = NULL, C.r.down = NULL,
C.theta.up = NULL, C.theta.down = NULL,
C.phi.up = NULL, C.phi.down = NULL,
flux.r.up = NULL, flux.r.down = NULL,
flux.theta.up = NULL, flux.theta.down = NULL,
flux.phi.up = NULL, flux.phi.down = NULL,
cyclicBnd = NULL,
D.r = NULL, D.theta = D.r, D.phi = D.r,
r = NULL,  theta = NULL, phi = NULL)
{
# ------------------------------------------------------------------------------
# Initialisation
# ------------------------------------------------------------------------------

dimens <- dim(C)
if (length(dimens) != 3)
stop("'C' should an array of dimensionality 3")

Nr     <- dimens[1]
Ntheta <- dimens[2]
Nphi     <- dimens[3]

if (length(r) != Nr+1)
stop("Length of 'r' should equal first dimension of 'C' + 1")
if (length(theta) != Ntheta+1)
stop("Length of 'theta' should equal second dimension of 'C' + 1")
if (length(phi) != Nphi+1)
stop("Length of 'phi' should equal third dimension of 'C' + 1")

if (max(theta) > 2 * pi)
stop("theta should be <= 2*pi")
if (min(theta) < 0)
stop("theta should be >= 0 ")

if (max(phi) > 2 * pi)
stop("phi should be <= 2*pi")
if (min(phi) < 0)
stop("phi should be >= 0 ")

# boundary conditions default = zero-gradient
Bc.r.Up       <- 3
Bc.r.Down     <- 3
Bc.theta.Up   <- 3
Bc.theta.Down <- 3
Bc.phi.Up       <- 3
Bc.phi.Down     <- 3

# boundaries in r-direction
if (! is.null(C.r.up))
Bc.r.Up <- 2
else
C.r.up <- 0.

if (! is.null(flux.r.up))
Bc.r.Up <- 1
else
flux.r.up <- 0.

if (! is.null(C.r.down))
Bc.r.Down <- 2
else
C.r.down <- 0.

if (! is.null(flux.r.down))
Bc.r.Down <- 1
else
flux.r.down <- 0.

if (1 %in% cyclicBnd) {
Bc.r.Up   <- 5
Bc.r.Down <- 5
}

# boundaries in theta-direction
if (! is.null(C.theta.up))
Bc.theta.Up <- 2
else
C.theta.up <- 0.

if (! is.null(flux.theta.up))
Bc.theta.Up <- 1
else
flux.theta.up <- 0.

if (! is.null(C.theta.down))
Bc.theta.Down <- 2
else
C.theta.down <- 0.

if (! is.null(flux.theta.down))
Bc.theta.Down <- 1
else
flux.theta.down <- 0.

if (2 %in% cyclicBnd) {
Bc.theta.Up   <- 5
Bc.theta.Down <- 5
}

# boundaries in phi-direction
if (! is.null(C.phi.up))
Bc.phi.Up <- 2
else
C.phi.up <- 0.

if (! is.null(flux.phi.up))
Bc.phi.Up <- 1
else
flux.phi.up <- 0.

if (! is.null(C.phi.down))
Bc.phi.Down <- 2
else
C.phi.down <- 0.

if (! is.null(flux.phi.down))
Bc.phi.Down <- 1
else
flux.phi.down <- 0.

if (3 %in% cyclicBnd) {
Bc.phi.Up   <- 5
Bc.phi.Down <- 5
}

# Diffusion coefficients
D.r <- rep(D.r, length.out=Nr+1)
D.theta <- rep(D.theta, length.out=Ntheta+1)
D.phi <- rep(D.phi, length.out=Nphi+1)

if (length(C.r.up) + length(flux.r.up) +
length(C.r.down) + length(flux.r.down) +
length(C.theta.up) + length(flux.theta.up)+
length(C.theta.down) + length(flux.theta.down)+
length(C.phi.up) + length(flux.phi.up) +
length(C.phi.down)+length(flux.phi.down)!=12)
stop("length of all 'boundary conditions' should be 1")

tr <- .Fortran("diffspher", as.integer(Nr), as.integer(Ntheta),
as.integer(Nphi), as.double(C),
as.double(C.r.up), as.double(C.r.down),
as.double(C.theta.up), as.double(C.theta.down),
as.double(C.phi.up), as.double(C.phi.down),
as.double(flux.r.up),  as.double(flux.r.down),
as.double(flux.theta.up),  as.double(flux.theta.down),
as.double(flux.phi.up),  as.double(flux.phi.down),
as.integer(Bc.r.Up), as.integer(Bc.r.Down),
as.integer(Bc.theta.Up), as.integer(Bc.theta.Down),
as.integer(Bc.phi.Up), as.integer(Bc.phi.Down),
as.double(D.r), as.double(D.theta), as.double(D.phi),
as.double(r), as.double(theta), as.double(phi),
Fluxrup= numeric(Ntheta*Nphi),Fluxrdown= numeric(Ntheta*Nphi),
Fluxthetaup= numeric(Nr*Nphi),Fluxthetadown= numeric(Nr*Nphi),
Fluxphiup= numeric(Nr*Ntheta),Fluxphidown= numeric(Nr*Ntheta),
dC = numeric(Nr*Ntheta*Nphi))

list(dC = array(dim=c(Nr, Ntheta, Nphi), tr\$dC),
flux.r.up   = matrix(nr=Ntheta, nc = Nphi, tr\$Fluxrup),
flux.r.down = matrix(nr=Ntheta, nc = Nphi, tr\$Fluxrdown),
flux.theta.up   = matrix(nr=Nr, nc = Nphi, tr\$Fluxthetaup),
flux.theta.down = matrix(nr=Nr, nc = Nphi, tr\$Fluxthetadown),
flux.phi.up   = matrix(nr=Nr, nc = Ntheta, tr\$Fluxphiup),
flux.phi.down = matrix(nr=Nr, nc = Ntheta, tr\$Fluxphidown)
)
}
``````