CHNOSZ: R/util.units.R

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
# CHNOSZ/util.units.R
# set units and convert values between units

P.units <- function(units=NULL) {
  ## change units of pressure or list the current one
  # show the current units, if none is specified
  if(is.null(units)) return(get("thermo")$opt$P.units)
  # argument handling
  units <- tolower(units)
  if(!units %in% c("bar","mpa")) stop("units of pressure must be either bar or MPa")
  # set the units and return them
  if(units=="bar") with(as.environment("CHNOSZ"), thermo$opt$P.units <- "bar")
  if(units=="mpa") with(as.environment("CHNOSZ"), thermo$opt$P.units <- "MPa")
  message("changed pressure units to ", get("thermo")$opt$P.units)
}

T.units <- function(units=NULL) {
  ## change units of temperature or list the current one
  # show the current units, if none is specified
  if(is.null(units)) return(get("thermo")$opt$T.units)
  # argument handling
  units <- tolower(units)
  if(!units %in% c("c","k")) stop("units of temperature must be either C or K")
  # set the units and return them
  if(units=="c") with(as.environment("CHNOSZ"), thermo$opt$T.units <- "C")
  if(units=="k") with(as.environment("CHNOSZ"), thermo$opt$T.units <- "K")
  message("changed temperature units to ", get("thermo")$opt$T.units)
}

E.units <- function(units=NULL) {
  ## change units of energy or list the current one
  # show the current units, if none is specified
  if(is.null(units)) return(get("thermo")$opt$E.units)
  # argument handling
  units <- tolower(units)
  if(!units %in% c("cal","j")) stop("units of energy must be either cal or J")
  # set the units and return them
  if(units=="cal") with(as.environment("CHNOSZ"), thermo$opt$E.units <- "cal")
  if(units=="j") with(as.environment("CHNOSZ"), thermo$opt$E.units <- "J")
  message("changed energy units to ", get("thermo")$opt$E.units)
}

convert <- function(value, units, T=get("thermo")$opt$Tr,
  P=get("thermo")$opt$Pr, pH=7, logaH2O=0) {
  # converts value(s) to the specified units

  if(is.null(value)) return(NULL)
  ### argument handling
  if(!is.character(units)) stop(paste('convert: please specify',
    'a character argument for the destination units.\n',
    'possibilities include (G or logK) (C or K) (J or cal) (cm3bar or calories) (Eh or pe)\n',
    'or their lowercase equivalents.\n'),call.=FALSE)
  Units <- units # for the possible message to user
  units <- tolower(units)

  # tests and calculations for the specified units
  if(units %in% c('c','k')) {
    CK <- 273.15
    if(units=='k') value <- value + CK
    if(units=='c') value <- value - CK 
  }
  else if(units[1] %in% c('j','cal')) {
    Jcal <- 4.184
    if(units=='j') value <- value * Jcal
    if(units=='cal') value <- value / Jcal
  }
  else if(units %in% c('g','logk')) {
    R <- get("thermo")$opt$R
    if(units=='logk') value <- value / (-log(10) * R * T)
    if(units=='g') value <- value * (-log(10) * R * T)
  }
  else if(units %in% c('cm3bar','calories')) {
    if(units=='cm3bar') value <- convert(value,'J') * 10
    if(units=='calories') value <- convert(value,'cal') / 10
  }
  else if(units %in% c('eh','pe')) {
    R <- 0.00831470
    F <- 96.4935
    if(units=='pe') value <- value * F / ( log(10) * R * T )
    if(units=='eh') value <- value * ( log(10) * R * T ) / F
  }
  else if(units %in% c('bar','mpa')) {
    barmpa <- 10
    if(units=='mpa') value <- value / barmpa
    if(units=='bar') value <- value * barmpa
  }
  else if(units %in% c('e0','logfo2')) {
    # convert between Eh and logfO2
    supcrt.out <- subcrt(c("H2O", "oxygen", "H+", "e-"), c(-1, 0.5, 2, 2), T=T, P=P, convert=FALSE)
    if(units=='logfo2') value <- 2*(supcrt.out$out$logK + logaH2O + 2*pH + 2*(convert(value,'pe',T=T)))
    if(units=='e0') value <- convert(( -supcrt.out$out$logK - 2*pH + value/2 - logaH2O )/2, 'Eh',T=T)
  }
  else cat(paste('convert: no conversion to ',Units,' found.\n',sep=''))
  return(value)
}

### unexported functions ###

outvert <- function(value,units) {
  # converts the given value from the given units to
  # those specified in thermo$opt
  units <- tolower(units)
  opt <- get("thermo")$opt
  if(units %in% c('c','k')) {
    if(units=='c' & opt$T.units=='K') return(convert(value,'k'))
    if(units=='k' & opt$T.units=='C') return(convert(value,'c'))
  }
  if(units %in% c('j','cal')) {
    if(units=='j' & opt$E.units=='Cal') return(convert(value,'cal'))
    if(units=='cal' & opt$E.units=='J') return(convert(value,'j'))
  }
  if(units %in% c('bar','mpa')) {
    if(units=='mpa' & opt$P.units=='bar') return(convert(value,'bar'))
    if(units=='bar' & opt$P.units=='MPa') return(convert(value,'mpa'))
  }
  return(value)
}

envert <- function(value,units) {
  # convert values to the specified units
  # from those given in thermo$opt
  if(!is.numeric(value[1])) return(value)
  units <- tolower(units)
  opt <- get("thermo")$opt
  if(units %in% c('c','k','t.units')) {
    if(units=='c' & opt$T.units=='K') return(convert(value,'c'))
    if(units=='k' & opt$T.units=='C') return(convert(value,'k'))
  }
  if(units %in% c('j','cal','e.units')) {
    if(units=='j' & opt$T.units=='Cal') return(convert(value,'j'))
    if(units=='cal' & opt$T.units=='J') return(convert(value,'cal'))
  }
  if(units %in% c('bar','mpa','p.units')) {
    if(units=='mpa' & opt$P.units=='bar') return(convert(value,'mpa'))
    if(units=='bar' & opt$P.units=='MPa') return(convert(value,'bar'))
  }
  return(value)
}

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.