##==============================================================================
## 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) 
       )
}