R/state.R

In increasechief/chromastar: Pedagogical Stellar Atmposhere and Spectrum Modelling of Stars

## MIT License
##
## Copyright (c) 2018 Oliver Dechant
##
## Permission is hereby granted, free of charge, to any person obtaining a copy
## of this software and associated documentation files (the "Software"), to deal
## in the Software without restriction, including without limitation the rights
## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
## copies of the Software, and to permit persons to whom the Software is
## furnished to do so, subject to the following conditions {
##
## The above copyright notice and this permission notice shall be included in all
## copies or substantial portions of the Software.
##
## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
## OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
## SOFTWARE.

useful <- readRDS("./data/useful.rds")

massDensity <- function(numDeps, temp, press, mmw, zScale) {
  #press is a 4x numDeps array
  #rows 0-1 are linear and log *gas* pressure
  #rows 2-3 are linear and log *radiation* pressure
  # row 0 of mmwNe is mean molecular weight in amu
  logE <- log10(exp(1))
  k <- useful[,"K"]
  logK <- useful[,"logK"]
  amu <- useful[,"amu"]
  logAmu <- useful[,"logAmu"]

  rho <- matrix(0.0, numDeps, 2)
  #declare scatch variables
  for (i in 1:numDeps) {
    logMuAmu <- log(mmw[i]) + logAmu
    #compute LTE bolometric radiation contribution to total HSE pressure
    rho[i][2] <- press[i][2]-temp[i][2]+(logMuAmu-logK)
    rho[i][1] <- exp(rho[i][2])
  }

  rho
}

mmwFn <- function(numDeps, temp, zScale) {
  #mean molecular weight in amu
  #Carrol and Ostlie 2nd Ed., p. 293: mu_N = 1.3, mu_I = 0.62
  mmw <- rep(0.0, numDeps)
  logMuN <- log(1.3)
  logMuI <- log(0.62)
  logTempN <- log(4000.0) #teff in K for fully neutral gas
  logTempI <- log(10000.0) #teff in K for *Hydrogen* fully ionized

  for (i in 1:numDeps) {
    if (is.na(temp[i,2]) || is.na(logTempN)) {
      stop("Na/NAN incorporated")
    }
    if (temp[i,2] < logTempN) {
      mmw[i] <- exp(logMuN)
    } else if ((temp[i,2] > logTempN) && (temp[i,2] < logTempI)) {
      logMu <- logMuN+((temp[i,2]-logTempN)/(logTempI-logTempN))*(logMuI-logMuN)
      #mean molecular weight in amu
      mmw[i] <- exp(logMu)
    } else {
      mmw[i] <- exp(logMuI)
    }
  }

  mmw
}

getNz <- function(numDeps, temp, pGas, guessPe, ATot, nelemAbnd, logAz) {
  logNz <- matrix(0.0, numDeps, nelemAbnd)
  logATot <- log(ATot)

  for (i in 1:numDeps) {
    #initial safety check to avoid negative logNz as Pg and Pe each converge
    #maximum physical Pe is about 0.5*Pgas (complete ionization of pure H)
    guessPe[is.na(guessPe)] <- 0
    if (guessPe[i,1] > 0.5*pGas[i,1]) {
      guessPe[i,1] <- 0.5*pGas[i,1]
      guessPe[i,2] <- log(guessPe[i,1])
    }
    #H (Z=1) is a special case
    logHelp <- guessPe[i,2]-pGas[i,2]
    help <- 1.0-exp(logHelp)
    logHelp <- log(help)
    logNumerator <- pGas[i,2]+logHelp
    logNz[i,1] <- logNumerator-useful[,"logK"]-temp[i,2]-logATot
    #remaining elements
    for (j in 1:nelemAbnd) {
      logNz[i,j] <- logAz[j]+logNz[i,1]
    }
  }

  logNz
}

massDensity2 <- function(numDeps, nelemAbnd, logNz, cname) {
  rho <- matrix(0.0, numDeps, 2)

  lAmu <- useful[,"logAmu"]
  #prepare log atomic masses once for each element
  logAMass <- rep(0.0,nelemAbnd)

  for (i in 1:nelemAbnd) {
    logAMass[i] <- log(getMass(cname[i]))

    for (i in 1:numDeps) {
      rho[i,1] <- 0.0
      for (j in 1:nelemAbnd) {
        logAddend <- logNz[i][j]+lAmu+logAMass[j]
        rho[i,1] <- rho[i][1]+exp(logAddend)
        rho[i,2] <- log(rho[i][1])
      }
    }
  }

  rho
}