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