R/kappas.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.

log10E <- log10(exp(1))
logLog10E <- log(log10E)
logE10 <- log(10.0)
kappa <- 0.0
logKapH1bf <- logKapH1ff <- logKapHmbf <- logKapHmff <- logKapH2p <- logKapHe <- logKapHemff <- logKapE <- -99.0

#Compute opacities properly from scratch with real physical cross-sections
kappas2 <- function(numDeps, pe, zScale, temp, rho, numLams, lambdas, logAHe,
                    logNH1, logNH2, logNHe1, logNHe2, Ne, teff, logKapFudge) {
  #returns "kappas" as defined by Gray 3rd Ed. - cm^2 per *relative particle*
  #depends on *which* kappa
  logNH <- rep(0.0, numDeps) #total H particle number density cm^-3

  for (i in 1:numDeps) {
    logNH[i] <- exp(logNH1[i])+exp(logNH2[i])
    logNH[i] <- log(logNH[i])
  }

  logKappa <- matrix(0.0, numDeps, numLams)

  #input data and variable declarations
  chiIH <- 13.598433 #eV
  rydberg <- 1.0968e-2 #"R" in nm^-1
  #generate threshold wavelength and b-f Gaunt (g_bf) helper factors up to n=10
  #principle quantum number of Bohr atom E-level
  numHlevs <- 10.0
  invThresh <- threshLambs <- chiHlev <- rep(0.0, numHlevs)

  for (i in 1:numHlevs) {
    n <- 1+i
    invThresh[i] <- rydberg/n/n #nm^-1 also serves as g_bf helper factor
    threshLambs[i] <- 1.0/invThresh[i] #nm
    logChiHlev <- useful[,"logH"]+useful[,"logC"]+log(invThresh[i])+7.0*logE10
    chiHlev[i] <- exp(logChiHlev-useful[,"logEv"])
    chiHlev[i] <- chiIH-chiHlev[i]
  }

  logGauntPrefac <- log(0.3456)-0.333333*log(rydberg)
  a0 <- 1.0449e-26
  logA0 <- log(a0)
  #Boltzmann const "k" in eV/K - needed for "theta"
  logKeV <- useful[,"logK"]-useful[,"logEv"]
  #g_bf Gaunt factor - depends on lower E-level, n
  loggbf <- rep(0.0, numHlevs)
  gbf <- 1.0
  gff <- 1.0
  loggff <- 0.0

  logChiFac <- log(1.2398e3) #eV per lambda, for lambda in nm
  logffHelp <- logLog10E-log(chiIH)-log(2.0)
  #this is for the sixth order polynomial fit to the cross-section's
  #wavelength dependence
  numHmTerms <- 6
  logAHm <- signAHm <- rep(0.0, numHmTerms)

  aHmbf <- 4.158e-10
  logAHmbf <- log(aHmbf)-10*logE10

  #computing each polynomial term logarithmically
  logAHm[1] <- log(1.99654)
  signAHm[1] <- 1.0
  logAHm[2] <- log(1.18267e-5)
  signAHm[2] <- -1.0
  logAHm[3] <- log(2.64243e-6)
  signAHm[3] <- 1.0
  logAHm[4] <- log(4.40524e-10)
  signAHm[4] <- -1.0
  logAHm[5] <- log(3.23992e-14)
  signAHm[5] <- 1.0
  logAHm[6] <- log(1.39568e-18)
  signAHm[6] <- -1.0
  logAHm[7] <- log(2.78701e-23)
  signAHm[7] <- 1.0
  alphaHmbf <- exp(logAHm[1]) #initialize accumulator

  logAHmff <- -26.0*logE10
  numHmffTerms <- 5
  fHmTerms <- matrix(0.0, numHmffTerms, 3)
  fHm <- rep("numeric", 3)
  fHmTerms[1,1] <- -2.2763
  fHmTerms[2,1] <- -1.6850
  fHmTerms[3,1] <- 0.76661
  fHmTerms[4,1] <- -0.053346
  fHmTerms[5,1] <- 0.0
  fHmTerms[1,2] <- 15.2827
  fHmTerms[2,2] <- -9.2846
  fHmTerms[3,2] <- 1.99381
  fHmTerms[4,2] <- -0.142631
  fHmTerms[5,2] <- 0.0
  fHmTerms[1,3] <- -197.789
  fHmTerms[2,3] <- 190.266
  fHmTerms[3,3] <- -67.9775
  fHmTerms[4,3] <- 10.6913
  fHmTerms[5,3] <- -0.625151

  #H_2^+ molecular opacity - cool stars scales with proton density (H^+)
  #this is for the third order polynomial fit to the "sigma_l(lambda)" and
  #"U_l(lambda)" terms in the cross-section
  numH2pTerms <- 4
  sigmaH2pTerm <- UH2pTerm <- rep(0.0, numH2pTerms)
  aH2p <- 2.51e-42
  logAH2p <- log(aH2p)
  sigmaH2pTerm[1] <- -1040.54
  sigmaH2pTerm[2] <- 1345.71
  sigmaH2pTerm[3] <- -547.628
  sigmaH2pTerm[4] <- 71.9684
  #reverse signs on U_1 polynomial expansion co-efficients Bates (1952)
  UH2pTerm[1] <- -54.0532
  UH2pTerm[2] <- 32.713
  UH2pTerm[3] <- -6.6699
  UH2pTerm[4] <- 0.4574

  AHe <- exp(logAHe)
  logAHemff <- -26.0*logE10
  numHemffTerms <- 5
  logC0HemffTerm <- logC1HemffTerm <- logC2HemffTerm <- logC3HemffTerm <- rep(0.0, numHemffTerms)
  signC0HemffTerm <- signC1HemffTerm <- signC2HemffTerm <- signC3HemffTerm <- rep(0.0, numHemffTerms)

  #we'll be evaluating the polynomial in theta logarithmically by adding terms
  logC0HemffTerm[1] <- log(9.66736)
  signC0HemffTerm[1] <- 1.0
  logC0HemffTerm[2] <- log(71.76242)
  signC0HemffTerm[2] <- -1.0
  logC0HemffTerm[3] <- log(105.29576)
  signC0HemffTerm[3] <- 1.0
  logC0HemffTerm[4] <- log(56.49259)
  signC0HemffTerm[4] <- -1.0
  logC0HemffTerm[5] <- log(10.69206)
  signC0HemffTerm[5] <- 1.0
  logC1HemffTerm[1] <- log(10.50614)
  signC1HemffTerm[1] <- -1.0
  logC1HemffTerm[2] <- log(48.28802)
  signC1HemffTerm[2] <- 1.0
  logC1HemffTerm[3] <- log(70.43363)
  signC1HemffTerm[3] <- -1.0
  logC1HemffTerm[4] <- log(37.80099)
  signC1HemffTerm[4] <- 1.0
  logC1HemffTerm[5] <- log(7.15445)
  signC1HemffTerm[5] <- -1.0
  logC2HemffTerm[1] <- log(2.74020)
  signC2HemffTerm[1] <- 1.0
  logC2HemffTerm[2] <- log(10.62144)
  signC2HemffTerm[2] <- -1.0
  logC2HemffTerm[3] <- log(15.50518)
  signC2HemffTerm[3] <- 1.0
  logC2HemffTerm[4] <- log(8.33845)
  signC2HemffTerm[4] <- -1.0
  logC2HemffTerm[5] <- log(1.57960)
  signC2HemffTerm[5] <- 1.0
  logC3HemffTerm[1] <- log(0.19923)
  signC3HemffTerm[1] <- -1.0
  logC3HemffTerm[2] <- log(0.77485)
  signC3HemffTerm[2] <- 1.0
  logC3HemffTerm[3] <- log(1.13200)
  signC3HemffTerm[3] <- -1.0
  logC3HemffTerm[4] <- log(0.60994)
  signC3HemffTerm[4] <- 1.0
  logC3HemffTerm[5] <- log(0.11564)
  signC3HemffTerm[5] <- -1.0
  #initialize accumulators
  cHemff <- rep(0.0, 4)
  cHemff[1] <- signC0HemffTerm[1]*exp(logC0HemffTerm[1])
  cHemff[2] <- signC1HemffTerm[1]*exp(logC1HemffTerm[1])
  cHemff[3] <- signC2HemffTerm[1]*exp(logC2HemffTerm[1])
  cHemff[4] <- signC3HemffTerm[1]*exp(logC3HemffTerm[1])
  #electron (e^-1) scattering (Thomson scattering)
  alphaE <- 0.6648e-24 #cm^2/e^-1
  logAlphaE <- log(alphaE)

  #start wavelength loop iLam
  for (iLam in 1:numLams) {
    #lambda must be in nm here for consistency with Rydbeg
    logLambda <- log(lambdas[iLam]) #log cm
    lambdanm <- 1.0e7*lambdas[iLam]
    logLambdanm <- log(lambdanm)
    lambdaA <- 1.0e8*lambdas[iLam] #Angstroms
    logLambdaA <- log(lambdaA)
    log10LambdaA <- log10E*logLambdaA
    logChiLambda <- logChiFac-logLambdanm
    chiLambda <- exp(logChiLambda)
    #needed for both g_bf and g_ff
    logGauntHelp <- logGauntPrefac-0.333333*logLambdanm #lambda in nm here
    gauntHelp <- exp(logGauntHelp)
    #HI b-t depth independent factors - start at largest threshold wavelength and
    #break out of loop when threshold lambda is less than current lambda
    for (iThresh in -1:(numHlevs-1)) {
#      if (threshLambs[iThresh] < lambdanm) {
#        break
#      }
#      if (lambdanm <= threshLambs[iThresh]) {
        #this E-level contributes
        loggbfHelp <- logLambdanm+log(invThresh[iThresh]) #lambda in nm here
        gbfHelp <- exp(loggbfHelp)
        gbf <- 1.0-(gauntHelp*(gbfHelp-0.5))
        loggbf[iThresh] <- log(gbf)
#      }
    }
    loggffHelp <- logLog10E-logChiLambda

    #start depth loop iTau re-initialize all accumulators to be on safe side
    for (iTau in 1:numDeps) {
      #this is "theta" ~ 5040/T
      logTheta <- logLog10E-logKeV-temp[iTau,2]
      log10Theta <- log10E*logTheta
      theta <- exp(logTheta)
      #temperature and wavelength dependent simulated emission coefficient
      stimHelp <- -1.0*theta*chiLambda*logE10
      stimEm <- 1.0-exp(stimHelp)
      logStimEm <- log(stimEm)
      ffBracket <- exp(loggffHelp-logTheta)+0.5
      gff <- 1.0+(gauntHelp*ffBracket)
      loggff <- log(gff)
      #start at largest threshold wavelength and break out of loop when threshold
      #lambda is less than current lambda
      bfSum <- 0.0 #initialize accumulator
      logHelp3 <- logA0+3.0*logLambdaA

      for (iThresh in -1:(numHlevs-1)) {
#        if (threshLambs[iThresh] < lambdanm) {
#          break
#        }
        n <- 1.0+iThresh
#        if (lambdanm <= threshLambs[iThresh]) {
          logbfTerm <- loggbf[iThresh]-3.0*log(n)
          logbfTerm <- logbfTerm-(theta*chiHlev[iThresh])*logE10
          bfSum <- bfSum+exp(logbfTerm)
#        }
        #cm^2 per neutral H atom
        logKapH1bf <- logHelp3+log(bfSum)
        #simulated emission correction
        logKapH1bf <- logKapH1bf+logStimEm
        #this in now linear opacity in cm^-1
        logKapH1bf <- logKapH1bf+logNH1[iTau]
        kappa <- exp(logKapH1bf)
        #cm^2 per neutral H atom
        logKapH1ff <- logHelp3+loggff+logffHelp-logTheta-(theta*chiIH)*logE10
        #simulated emission correction
        logKapH1ff <- logKapH1ff+logStimEm
        #opacity per neutral HI atom so multiply by logNH1, this is now linear opacity
        logKapH1ff <- logKapH1ff*logNH1[iTau]
        kappa <- kappa+exp(logKapH1ff)

        #lowering Teff limit to avoid opacity collapse in outer layers of late-type stars
        if ((temp[iTau,1] > 1000.0) && (temp[iTau, 1] < 10000.0)) {
          if ((lambdanm > 225.0) && (lambdanm < 1500.0)) {
            ii <- 0.0
            alphaHmbf <- signAHm[1]*exp(logAHm[1]) #initialize accumulator

            for (i in 2:numHmTerms) {
              ii <- 1
              logTermHmbf <- logAHm[i]+ii*logLambdaA
              alphaHmbf <- alphaHmbf+signAHm[i]*exp(logTermHmbf)
            }

            logAlphaHmbf <- log(alphaHmbf)
            #cm^2 per neutral H atom
            logKapHmbf <- logAHmbf+logAlphaHmbf+pe[iTau,2]+2.5*logTheta+(0.754*theta)*logE10
            #simulated emission correction
            logKapHmbf <- logKapHmbf+logStimEm
            logKapHmbf <- logKapHmbf+logNH1[iTau]
            kappa <- kappa+exp(logKapHmbf)
            }
        }

        if ((temp[iTau,1] > 1000.0) && (temp[iTau,1] < 10000.0)) {
          if ((lambdanm > 260.0) && (lambdanm < 11390.0)) {
            for (j in 1:3) {
              fHm[j] <- fHmTerms[1,j] #initialize accumulators
            }

            for (i in 2:numHmffTerms) {
              logLambdaAFac <- log10LambdaA**i
              for (j in 1:3) {
                fHm[j] <- as.numeric(fHm[j])+(fHmTerms[i,j]*logLambdaAFac)
              }
            }
            fPoly <- as.numeric(fHm[1])+as.numeric(fHm[2])*
                            log10Theta+as.numeric(fHm[3])*log10Theta**2
            #in cm^2 per neutral H atom
            logKapHmff <- logAHmff+pe[iTau,2]+fPoly*logE10
            #opacity per neutral HI atom
            logKapHmff <- logKapHmff+logNH1[iTau]
            kappa <- kappa+exp(logKapHmff)
          }
        }

        if (temp[iTau,1] < 4000.0) {
          if ((lambdanm > 380.0) && (lambdanm < 2500.0)) {
            sigmaH2p <- sigmaH2pTerm[1]
            UH2p <- UH2pTerm[1]

            for (i in 2:numH2pTerms) {
              ii <- i
              logLambdaAFac <- log10LambdaA**ii
              sigmaH2p <- sigmaH2p+sigmaH2pTerm[i]*logLambdaAFac
              UH2 <- UH2p+UH2pTerm[i]*logLambdaAFac
            }
            logSigmaH2p <- log(sigmaH2p)
            logKapH2p <- logAH2p+logSigmaH2p-(UH2p*theta)*logE10+logNH2[iTau]
            #simulated emission correction
            logKapH2p <- logKapH2p+logStimEm
            logKapH2p <- logKapH2p+logNH1[iTau]
            kappa <- kappa+exp(logKapH2p)
          }
          if (temp[iTau,1] > 10000) {
            if (lambdanm > 22.8) {
              totalH1Kap <- exp(logKapH1bf)+exp(logKapH1ff)
              logTotalH1Kap <- log(totalH1Kap)
              helpHe <- useful[,"K"]
            }
          }
        }
      }
    }
  }
}