R/kappasMetal.R
In increasechief/chromastar: Pedagogical Stellar Atmposhere and Spectrum Modelling of Stars
Defines functions
masterMetal
opacC1
seaton
opacMg1
opacMg2
opacAl1
opacSi1
opacSi2
opacFe1
## 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")
masterMetal <- function(numDeps, numLams, temp, lambdaScale, stagePops) {
#Metal b-f opacity routines taken from Moog (moogjul2014/)
#Chris Sneden (University of Texas at Austin) and collaborators
#https://www.as.utexas.edu/~chris/moog.html
#From Moog source file Opacitymetals.f
#From how values such as aC1[] are used in Moog file Opacit.f to compute the total
#opacity and then the optical depth scale, I infer that they are extinct coefficients
#in cm^-1. There does not seem to be any correction for simulated emission
logE <- log10(exp(1))
masterBF <- matrix(0.0, numDeps, numLams)
logUC1<- rep("numeric", 5)
logUMg1<- rep("numeric", 5)
logUMg2<- rep("numeric", 5)
logUAl1<- rep("numeric", 5)
logUSi1<- rep("numeric", 5)
logUSi2<- rep("numeric", 5)
logUFe1<- rep("numeric", 5)
logStatWC1 <- logStatWMg1 <- logStatWMg2 <- logStatWAl1 <- logStatWSi1 <- logStatWSi2 <- logStatWFe1 <- 0.0
theta <- 1.0
species <- ""
logGroundPopsC1 <- logGroundPopsMg1 <- logGroundPopsMg2 <- logGroundPopsAl1 <- logGroundPopsSi1 <- logGroundPopsSi2 <- logGroundPopsFe1 <- rep("numeric", numDeps)
# C I: Z=6 --> iZ=5:
aC1 <- rep("numeric", numDeps)
# Mg I: Z=12 --> iZ=11:
aMg1 <- rep("numeric", numDeps)
# Mg II: Z=12 --> iZ=11:
aMg2 <- rep("numeric", numDeps)
# Al I: Z=13 --> iZ=12:
aAl1 <- rep("numeric", numDeps)
# Si I: Z=14 --> iZ=13:
aSi1 <- rep("numeric", numDeps)
# Si II: Z=14 --> iZ =13:
aSi2 <- rep("numeric", numDeps)
# Fe I: Z=26 --> iZ=25
aFe1 <- rep("numeric", numDeps)
species <- "CI"
logUC1 <- getPartFn2(species)
species <- "MgI"
logUMg1 <- getPartFn2(species)
species <- "MgII"
logUMg2 <- getPartFn2(species)
species <- "AlI"
logUAl1 <- getPartFn2(species)
species <- "SiI"
logUSi1 <- getPartFn2(species)
species <- "SiII"
logUSi2 <- getPartFn2(species)
species <- "FeI"
logUFe1 <- getPartFn2(species)
for (iD in 1:numDeps) {
#neutral stage assumes state stat weight, g_1 is 1.0
thisTemp <- temp[iD, 1]
if (thisTemp <= 130) {
logStatWC1 <- logUC1[1]
logStatWMg1 <- logUMg1[1]
logStatWMg2 <- logUMg2[1]
logStatWAl1 <- logUAl1[1]
logStatWSi1 <- logUSi1[1]
logStatWSi2 <- logUSi2[1]
logStatWFe1 <- logUFe1[1]
} else if (thisTemp > 130 && thisTemp <= 500) {
#// Add in interpolation here lburns
logStatWC1 <- logUC1[2] * (thisTemp - 130)/(500 - 130)
+ logUC1[1] * (500 - thisTemp)/(500 - 130)
logStatWMg1 <- logUMg1[2] * (thisTemp - 130)/(500 - 130)
+ logUMg1[1] * (500 - thisTemp)/(500 - 130)
logStatWMg2 <- logUMg2[2] * (thisTemp - 130)/(500 - 130)
+ logUMg2[1] * (500 - thisTemp)/(500 - 130)
logStatWAl1 <- logUAl1[2] * (thisTemp - 130)/(500 - 130)
+ logUAl1[1] * (500 - thisTemp)/(500 - 130)
logStatWSi1 <- logUSi1[2] * (thisTemp - 130)/(500 - 130)
+ logUSi1[1] * (500 - thisTemp)/(500 - 130)
logStatWSi2 <- logUSi2[2] * (thisTemp - 130)/(500 - 130)
+ logUSi2[1] * (500 - thisTemp)/(500 - 130)
logStatWFe1 <- logUFe1[2] * (thisTemp - 130)/(500 - 130)
+ logUFe1[1] * (500 - thisTemp)/(500 - 130)
} else if (thisTemp > 500 && thisTemp <= 3000) {
logStatWC1 <- logUC1[3] * (thisTemp - 500)/(3000 - 500)
+ logUC1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWMg1 <- logUMg1[3] * (thisTemp - 500)/(3000 - 500)
+ logUMg1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWMg2 <- logUMg2[3] * (thisTemp - 500)/(3000 - 500)
+ logUMg2[2] * (3000 - thisTemp)/(3000 - 500)
logStatWAl1 <- logUAl1[3] * (thisTemp - 500)/(3000 - 500)
+ logUAl1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWSi1 <- logUSi1[3] * (thisTemp - 500)/(3000 - 500)
+ logUSi1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWSi2 <- logUSi2[3] * (thisTemp - 500)/(3000 - 500)
+ logUSi2[2] * (3000 - thisTemp)/(3000 - 500)
logStatWFe1 <- logUFe1[3] * (thisTemp - 500)/(3000 - 500)
+ logUFe1[2] * (3000 - thisTemp)/(3000 - 500)
} else if (thisTemp > 3000 && thisTemp <= 8000) {
logStatWC1 <- logUC1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUC1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWMg1 <- logUMg1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUMg1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWMg2 <- logUMg2[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUMg2[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWAl1 <- logUAl1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUAl1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWSi1 <- logUSi1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUSi1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWSi2 <- logUSi2[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUSi2[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWFe1 <- logUFe1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUFe1[3] * (8000 - thisTemp)/(8000 - 3000)
} else if (thisTemp > 8000 && thisTemp < 10000) {
logStatWC1 <- logUC1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUC1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWMg1 <- logUMg1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUMg1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWMg2 <- logUMg2[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUMg2[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWAl1 <- logUAl1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUAl1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWSi1 <- logUSi1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUSi1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWSi2 <- logUSi2[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUSi2[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWFe1 <- logUFe1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUFe1[4] * (10000 - thisTemp)/(10000 - 8000)
} else {
# for temperatures greater than or equal to 10000
logStatWC1 <- logUC1[5]
logStatWMg1 <- logUMg1[5]
logStatWMg2 <- logUMg2[5]
logStatWAl1 <- logUAl1[5]
logStatWSi1 <- logUSi1[5]
logStatWSi2 <- logUSi2[5]
logStatWFe1 <- logUFe1[5]
}
logGroundPopsC1[iD] <- stagePops[1,6][iD] - logStatWC1
logGroundPopsMg1[iD] <- stagePops[1,12][iD] - logStatWMg1
logGroundPopsMg2[iD] <- stagePops[2,12][iD] - logStatWMg2
logGroundPopsAl1[iD] <- stagePops[1,13][iD] - logStatWAl1
logGroundPopsSi1[iD] <- stagePops[1,14][iD] - logStatWSi1
logGroundPopsSi2[iD] <- stagePops[2,14][iD] - logStatWSi2
logGroundPopsFe1[iD] <- stagePops[1,26][iD] - logStatWFe1
}
for (iL in 1:numLams) {
for (i in 1:numDeps) {
aC1[i] <- 0.0
aMg1[i] <- 0.0
aMg2[i] <- 0.0
aAl1[i] <- 0.0
aSi1[i] <- 0.0
aSi2[i] <- 0.0
aFe1[i] <- 0.0
}
waveno <- 1.0/lambdaScale[iL]
logFreq <- useful[,"logC"]-log(lambdaScale[i])
freq <- exp(logFreq)
stimEmLogExpHelp <- useful[,"logH"]+logFreq-useful[,"logK"]
if (freq >= 2.0761e15)
aC1 = opacC1(numDeps, temp, lambdaScale[iL], logGroundPopsC1)
if (freq >= 2.997925e+14)
aMg1 = opacMg1(numDeps, temp, lambdaScale[iL], logGroundPopsMg1)
if (freq >= 2.564306e15)
aMg2 = opacMg2(numDeps, temp, lambdaScale[iL], logGroundPopsMg2)
if (freq >= 1.443e15)
aAl1 = opacAl1(numDeps, temp, lambdaScale[iL], logGroundPopsAl1)
if (freq >= 2.997925e+14)
aSi1 = opacSi1(numDeps, temp, lambdaScale[iL], logGroundPopsSi1)
if (freq >= 7.6869872e14)
aSi2 = opacSi2(numDeps, temp, lambdaScale[iL], logGroundPopsSi2)
if (waveno >= 21000.0)
aFe1 = opacFe1(numDeps, temp, lambdaScale[iL], logGroundPopsFe1)
for (iD in 1:numDeps) {
kapBF <- 1.0e-99 #minimum safe value
stimEmLogExp <- stimEmLogExpHelp - temp[iD,2]
stimEmExp <- -1.0 * exp(stimEmLogExp)
stimEm <- 1.0 - exp(stimEmExp) #//LTE correction for stimulated emission
kapBF <- kapBF + aC1[iD] + aMg1[iD] + aMg2[iD] + aAl1[iD] + aSi1[iD] + aSi2[iD] + aFe1[iD]
}
}
masterBF
}
opacC1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to CI
waveno <- 1.0/lambda2
freq <- useful[,"c"]/lambda2
c1240 <- c1444 <- tkev <- aC1 <- rep("numeric", numDeps)
freq1 <- sigma <- 0.0
for (i in 1:numDeps) {
logTkev <- temp[i,2]+useful[,"logK"]-useful[,"logEv"]
tkev[i] <- exp(logTkev)
}
#initialize some quantities for each new model atmosphere
for (i in 1:numDeps) {
c1240[i] <- 5.0*exp(-1.264/tkev[i])
c1444[i] <- exp(-2.683/tkev[i])
}
#initialize some quantities for each new model atmosphere or new frequency
#Luo, D. and Pradhan, A.K. 1989, J.Phys. B, 22. 3377-3395.
#Burke, P.G. and Taylor, K.T. 1979, J.Phys. B, 12, 2971-2984.
ryd <- 109732.298
xs0 <- xs1 <- 0.0
xd0 <- xd1 <- xd2 <- 0.0
x1444 <- x1240 <- x1100 <- 0.0
if (freq >= 2.7254E15) {
arg <- -16.80-((waveno-90777.00)/3.0/ryd)
x1100 <- (10**arg)*seaton(2.7254E15, 1.219E-7, 2.0E0, 3.317E0, freq)
}
if (freq >= 2.4196E15) {
arg <- -16.80-((waveno-80627.76)/3.0/ryd)
xd0 <- 10**arg
eeps <- (waveno-93917.0)*2.0/9230.0
aa <- 22.0E-18
bb <- 26.0E-18
xd1 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
eeps <- (waveno-111130.0)*2.0/2743.0
aa <- -10.5E-18
bb <- 46.0E-18
xd2 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
x1240 <- xd0+xd1+xd2
}
if (freq >= 2.0761E15) {
arg <- -16.80-((waveno-69172.40)/3.0/ryd)
xs0 <- 10**arg
eeps <- (waveno-97700.0)*2.0/2743.0
aa <- 68.0E-18
bb <- 118.0E-18
xs1 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
x1444 <- xs0+xs1
}
for (i in 1:numDeps) {
if (freq >= 2.0761E15) {
sigma <- (x1100*9.0+xs1240*c1240[i]+x1444*c1444[i])
aC1[i] <- sigma*exp(logGroundPops[i])
}
}
aC1
}
seaton <- function(freq0, xsect, power, a, freq) {
#This function is a general representation for b-f absorption above
#a given ionization limit freq0, using cross-section xsect
freqratio <- freq0/freq
help <- floor((2*power)+0.01)
seaton <- xsect*(a+freqratio*(1.0-a))*sqrt(freqratio**help)
seaton
}
opacMg1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Mg I
sigma <- 0.0
aMg1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 7)
dt <- rep("numeric", 100)
nt <- rep("numeric", 100)
#4000K
peach0 = c(-42.474, -41.808, -41.273, -45.583, -44.324, -50.969,
-50.633, -53.028, -51.785, -52.285, -52.028, -52.384,
-52.363, -54.704, -54.359)
#5000 K
peach1 =c(-42.350, -41.735, -41.223, -44.008, -42.747, -48.388,
-48.026, -49.643, -48.352, -48.797, -48.540, -48.876,
-48.856, -50.772, -50.349)
#6000 K
peach2 =c(-42.109, -41.582, -41.114, -42.957, -41.694, -46.630,
-46.220, -47.367, -46.050, -46.453, -46.196, -46.513,
-46.493, -48.107, -47.643)
#7000 K
peach3 =c(-41.795, -41.363, -40.951, -42.205, -40.939, -45.344,
-44.859, -45.729, -44.393, -44.765, -44.507, -44.806,
-44.786, -46.176, -45.685)
#8000 K
peach4 =c(-41.467, -41.115, -40.755, -41.639, -40.370, -44.355,
-43.803, -44.491, -43.140, -43.486, -43.227, -43.509,
-43.489, -44.707, -44.198)
#9000 K
peach5 =c(-41.159, -40.866, -40.549, -41.198, -39.925, -43.568,
-42.957, -43.520, -42.157, -42.480, -42.222, -42.488,
-42.467, -43.549, -43.027)
#10000 K
peach6 =c(-40.883, -40.631, -40.347, -40.841, -39.566, -42.924,
-42.264, -42.736, -41.363, -41.668, -41.408, -41.660,
-41.639, -42.611, -42.418)
peach <- rbind(peach0, peach1, peach2, peach3, peach4, peach5, peach6)
freqMg <- c(1.9341452e15, 1.8488510e15, 1.1925797e15,
7.9804046e14, 4.5772110e14, 4.1440977e14,
4.1113514e14)
flog <- c(35.23123, 35.19844, 35.15334, 34.71490, 34.31318,
33.75728, 33.65788, 33.64994, 33.43947)
tlg <- c(8.29405, 8.51719, 8.69951, 8.85367, 8.98720, 9.10498,
9.21034)
freq1 <- 0.0
for (i in 1:numDeps) {
thelp <- floor((temp[i,1]/1000.0))-3
nn <- max(min(6, thelp), 1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
freq1 <- freq
nn <- 0
for (i in 1:7) {
if (freq > freqMg[n]) {
break
}
n <- n+1
}
if (freq <= freqMg[6]) {
nn <- 7
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
dd1 <- 1.0-dd
for (it in 1:numDeps) {
if (freq >= 2.997925E14) {
nn <- nt[i]
sigma <- exp(xx[nn]*(1.0E0-dt[i]))+(xx[nn+1]*dt[i])
aMg1[i] <- sigma*exp(logGroundPops[i])
}
}
aMg1
}
opacMg2 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Mg II
sigma <- 0.0
aMg2 <- tkev <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
c1169 <- rep("numeric", 100)
freq1 <- x824 <- x1169 <- 0.0
for (i in 1:numDeps) {
logTkev <- temp[i,2]+useful[,"logK"]-useful[,"logEv"]
tkev[i] <- exp(logTkev)
}
for (i in 1:numDeps) {
c1169[i] <- 6.0*exp(-4.43E0/tkev[i])
}
if (freq >= 3.635492E15) {
x824 <- seaton(3.635492E15, 1.40E-19, 4.0E0, 6.7E0, freq)
} else {
x824 <- 1.0E-99
}
if (freq >= 2.564306E15) {
x1169 <- 5.11E-19*((2.564306E15/freq)**3.0)
} else {
x1169 <- 1.0E-99
}
for (i in 1:numDeps) {
if (x1169 >= 1.0E-90) {
sigma <- (x824*2.0+x1169*c1169[i])
aMg2[i] <- sigma*exp(logGroundPops[i])
}
}
aMg2
}
opacAl1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Al I
sigma <- 0.0
aAl1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
for (i in 1:numDeps) {
if (freq >= 1.443E15) {
sigma <- 6.0*6.5E-17*((1.443E15/freq)**5.0)
aAl1[i] <- sigma*exp(logGroundPops[i])
}
}
aAl1
}
opacSi1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Si I
sigma <- 0.0
aSi1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 9)
dt <- nt <- rep("numeric", 100)
#4000
peach0 <- c(38.136, 37.834, 37.898, 40.737, 40.581, 45.521,
45.520, 55.068, 53.868, 54.133, 54.051, 54.442,
54.320, 55.691, 55.661, 55.973, 55.922, 56.828, 56.657)
#5000
peach1 <- c(38.138, 37.839, 37.898, 40.319, 40.164, 44.456,
44.455, 51.783, 50.369, 50.597, 50.514, 50.854,
50.722, 51.965, 51.933, 52.193, 52.141, 52.821, 52.653)
#6000
peach2 <- c(38.140, 37.843, 37.897, 40.047, 39.893, 43.753,
43.752, 49.553, 48.031, 48.233, 48.150, 48.455,
48.313, 49.444, 49.412, 49.630, 49.577, 50.110, 49.944)
#7000
peach3 <- c(38.141, 37.847, 37.897, 39.855, 39.702, 43.254,
43.251, 47.942, 46.355, 46.539, 46.454, 46.733,
46.583, 47.615, 47.582, 47.769, 47.715, 48.146, 47.983)
#8000
peach4 <- c(38.143, 37.850, 37.897, 39.714, 39.561, 42.878,
42.871, 46.723, 45.092, 45.261, 45.176, 45.433,
45.277, 46.221, 46.188, 46.349, 46.295, 46.654, 46.491)
#9000
peach5 <- c(38.144, 37.853, 37.896, 39.604, 39.452, 42.580,
42.569, 45.768, 44.104, 44.262, 44.175, 44.415,
44.251, 45.119, 45.085, 45.226, 45.172, 45.477, 45.315)
#10000
peach6 <- c(38.144, 37.855, 37.895, 39.517, 39.366, 42.332,
42.315, 44.997, 43.308, 43.456, 43.368, 43.592,
43.423, 44.223, 44.189, 44.314, 44.259, 44.522, 44.360)
#11000
peach7 <- c(38.145, 37.857, 37.895, 39.445, 39.295, 42.119,
42.094, 44.360, 42.652, 42.790, 42.702, 42.912,
42.738, 43.478, 43.445, 43.555, 43.500, 43.730, 43.569)
#12000
peach8 <- c(38.145, 37.858, 37.894, 39.385, 39.235, 41.930,
41.896, 43.823, 42.100, 42.230, 42.141, 42.340,
42.160, 42.848, 42.813, 42.913, 42.858, 43.061, 42.901)
peach <- rbind(peach0, peach1, peach2, peach3, peach4, peach5, peach6,
peach7, peach8)
freqSi <- c(2.1413750e15, 1.9723165e15, 1.7879689e15,
1.5152920e15, 5.5723927e14, 5.3295914e14,
4.7886458e14, 4.7216422e14, 4.6185133e14)
flog <- c(35.45438, 35.30022, 35.21799, 35.11986, 34.95438,
33.95402, 33.90947, 33.80244, 33.78835, 33.76626,
33.70518)
tlg <- c(8.29405, 8.51719, 8.69951, 8.85367, 8.98720, 9.10498,
9.21034, 9.30565, 9.39266)
freq1 <- 0.0
for (i in 1:numDeps) {
thelp <- floor(temp[i,1]/1000.0)-3
nn <- max(min(8,thelp),1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
freq1 <- freq
nn <- 0
for (n in 1:9) {
if (freq > freqSi[n]) {
break
}
nn <- nn+1
}
if (freq <= freqSi[8]) {
nn <- 9
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
dd1 <- 1.0-dd
for (it in 1:9) {
if (freq >= 2.997925E14) {
nn <- nt[i]
sigma <- (9.0*exp(-(xx[nn]*(1.-dt[i])+xx[nn+1]*dt[i])))
sSi1[i] <- sigma*exp(logGroundPops[i])
}
}
aSi1
}
opacSi2 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Si II
sigma <- 0.0
aSi2 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 6)
dt <- nt <- rep("numeric", 100)
#c 10000
peach0 <- c( -43.8941, -42.2444, -40.6054, -54.2389, -50.4108,
-52.0936, -51.9548, -54.2407, -52.7355, -53.5387,
-53.2417, -53.5097, -54.0561, -53.8469)
#c 12000
peach1 <- c( -43.8941, -42.2444, -40.6054, -52.2906, -48.4892,
-50.0741, -49.9371, -51.7319, -50.2218, -50.9189,
-50.6234, -50.8535, -51.2365, -51.0256)
#c 14000
peach2 <- c( -43.8941, -42.2444, -40.6054, -50.8799, -47.1090,
-48.5999, -48.4647, -49.9178, -48.4059, -49.0200,
-48.7252, -48.9263, -49.1980, -48.9860)
#c 16000
peach3 <- c( -43.8941, -42.2444, -40.6054, -49.8033, -46.0672,
-47.4676, -47.3340, -48.5395, -47.0267, -47.5750,
-47.2810, -47.4586, -47.6497, -47.4368)
#c 18000
peach4 <- c( -43.8941, -42.2444, -40.6054, -48.9485, -45.2510,
-46.5649, -46.4333, -47.4529, -45.9402, -46.4341,
-46.1410, -46.2994, -46.4302, -46.2162)
#c 20000
peach5 <- c( -43.8941, -42.2444, -40.6054, -48.2490, -44.5933,
-45.8246, -45.6947, -46.5709, -45.0592, -45.5082,
-45.2153, -45.3581, -45.4414, -45.2266)
peach <- c(peach0, peach1, peach2, peach3, peach4, peach5)
freqSi <- c(4.9965417e15, 3.9466738e15, 1.5736321e15,
1.5171539e15, 9.2378947e14, 8.3825004e14,
7.6869872e14)
flog <- c(36.32984, 36.14752, 35.91165, 34.99216, 34.95561,
34.45951, 34.36234, 34.27572, 34.20161)
tlg <- c(9.21034, 9.39266, 9.54681, 9.68034, 9.79813, 9.90349)
freq1 <- 0.0
#setup some data upon first entrance with a new model atmosphere
for (i in 1:numDeps) {
thelp <- floor(temp[i,1]/2000.0)-4.0
nn <- max(min(5, thelp), 1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
#initialize some quantities for each new model atmosphere or frequency
freq1 <- freq
nn <- 0
for (n in 1:7) {
if (freq > freqSi[n]) {
break
}
nn <- nn + 1
}
if (freq <= freqSi[6]) {
nn <- 7
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
if (nn == 13) {
nn <- 12
}
dd1 <- 1.0 - dd
for (it in 1:6) {
xx[it] <- peach[nn+1,it]*dd+peach[nn,it]*dd1
}
for (i in 1:numDeps) {
if (freq >= 7.6869872E14) {
nn <- nt[i]
sigma <- (6.0*exp(xx[nn]*(1.0-dt[i])+xx[nn+1]*dt[i]))
aSi2[i] <- sigma*exp(logGroundPops[i])
}
}
aSi2
}
opacFe1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Fe I
sigma <- 0.0
aFe1 <- rep("numeric", numDeps)
#cross-section is zero below threshold, so initialize
for (i in 1:numDeps) {
aFe1[i] <- 0.0
}
wavenoo <- 1.0/lambda2
freq <- useful[,"c"]/lambda2
bolt <- matrix(0.0, 100, 48)
xsect <- rep("numeric", 1:48)
gg <- c(25.0, 35.0, 21.0, 15.0, 9.0, 35.0, 33.0, 21.0, 27.0, 49.0,
9.0, 21.0, 27.0, 9.0, 9.0, 25.0, 33.0, 15.0, 35.0, 3.0,
5.0, 11.0, 15.0, 13.0, 15.0, 9.0, 21.0, 15.0, 21.0, 25.0, 35.0,
9.0, 5.0, 45.0, 27.0, 21.0, 15.0, 21.0, 15.0, 25.0, 21.0,
35.0, 5.0, 15.0, 45.0, 35.0, 55.0, 25.0)
ee<- c(500.0, 7500.0, 12500.0, 17500.0, 19000.0, 19500.0, 19500.0,
21000.0, 22000.0, 23000.0, 23000.0, 24000.0, 24000.0, 24500.0,
24500.0, 26000.0, 26500.0, 26500.0, 27000.0, 27500.0, 28500.0,
29000.0, 29500.0, 29500.0, 29500.0, 30000.0, 31500.0, 31500.0,
33500.0, 33500.0, 34000.0, 34500.0, 34500.0, 35000.0, 35500.0,
37000.0, 37000.0, 37000.0, 38500.0, 40000.0, 40000.0, 41000.0,
41000.0, 43000.0, 43000.0, 43000.0, 43000.0, 44000.0)
wno <- c(63500.0, 58500.0, 53500.0, 59500.0, 45000.0, 44500.0, 44500.0,
43000.0, 58000.0, 41000.0, 54000.0, 40000.0, 40000.0, 57500.0,
55500.0, 38000.0, 57500.0, 57500.0, 37000.0, 54500.0, 53500.0,
55000.0, 34500.0, 34500.0, 34500.0, 34000.0, 32500.0, 32500.0,
32500.0, 32500.0, 32000.0, 29500.0, 29500.0, 31000.0, 30500.0,
29000.0, 27000.0, 54000.0, 27500.0, 24000.0, 47000.0, 23000.0,
44000.0, 42000.0, 42000.0, 21000.0, 42000.0, 42000.0)
freq1 <- 0.0
#set up some data upon first entrance with a new model atmosphere
for (i in 1:numDeps) {
hkt <- 6.6256E-27/(1.38054E-16*temp[i,1])
for (k in 1:48) {
bolt[i,k] <- gg[k]*exp(-ee[k]*useful[,"c"]*hkt)
}
}
#initialize some quantities for each new model atmosphere or frequency
freq1 <- freq
if (waveno >= 21000.0) {
for (k in 1:48) {
xsect[k] <- 0.0
if (wno[k] < waveno) {
xsect[k] <- 3.0E-18/(1.0+((wno[k]+3000.0-waveno)/wno[k])/0.1)**4
}
}
}
for (i in 1:numDeps) {
aFe1[i] <- 0.0
if (waveno >= 21000.0) {
for (k in 1:48) {
aFe[i] <- 0.0
sigma <- aFe1[i]+xsect[k]*bolt[i,k]
aFe1[i] <- sigma*exp(logGroundPops[i])
}
}
}
aFe1
}
increasechief/chromastar documentation built on May 14, 2019, 5:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions masterMetal opacC1 seaton opacMg1 opacMg2 opacAl1 opacSi1 opacSi2 opacFe1
## 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")
masterMetal <- function(numDeps, numLams, temp, lambdaScale, stagePops) {
#Metal b-f opacity routines taken from Moog (moogjul2014/)
#Chris Sneden (University of Texas at Austin) and collaborators
#https://www.as.utexas.edu/~chris/moog.html
#From Moog source file Opacitymetals.f
#From how values such as aC1[] are used in Moog file Opacit.f to compute the total
#opacity and then the optical depth scale, I infer that they are extinct coefficients
#in cm^-1. There does not seem to be any correction for simulated emission
logE <- log10(exp(1))
masterBF <- matrix(0.0, numDeps, numLams)
logUC1<- rep("numeric", 5)
logUMg1<- rep("numeric", 5)
logUMg2<- rep("numeric", 5)
logUAl1<- rep("numeric", 5)
logUSi1<- rep("numeric", 5)
logUSi2<- rep("numeric", 5)
logUFe1<- rep("numeric", 5)
logStatWC1 <- logStatWMg1 <- logStatWMg2 <- logStatWAl1 <- logStatWSi1 <- logStatWSi2 <- logStatWFe1 <- 0.0
theta <- 1.0
species <- ""
logGroundPopsC1 <- logGroundPopsMg1 <- logGroundPopsMg2 <- logGroundPopsAl1 <- logGroundPopsSi1 <- logGroundPopsSi2 <- logGroundPopsFe1 <- rep("numeric", numDeps)
# C I: Z=6 --> iZ=5:
aC1 <- rep("numeric", numDeps)
# Mg I: Z=12 --> iZ=11:
aMg1 <- rep("numeric", numDeps)
# Mg II: Z=12 --> iZ=11:
aMg2 <- rep("numeric", numDeps)
# Al I: Z=13 --> iZ=12:
aAl1 <- rep("numeric", numDeps)
# Si I: Z=14 --> iZ=13:
aSi1 <- rep("numeric", numDeps)
# Si II: Z=14 --> iZ =13:
aSi2 <- rep("numeric", numDeps)
# Fe I: Z=26 --> iZ=25
aFe1 <- rep("numeric", numDeps)
species <- "CI"
logUC1 <- getPartFn2(species)
species <- "MgI"
logUMg1 <- getPartFn2(species)
species <- "MgII"
logUMg2 <- getPartFn2(species)
species <- "AlI"
logUAl1 <- getPartFn2(species)
species <- "SiI"
logUSi1 <- getPartFn2(species)
species <- "SiII"
logUSi2 <- getPartFn2(species)
species <- "FeI"
logUFe1 <- getPartFn2(species)
for (iD in 1:numDeps) {
#neutral stage assumes state stat weight, g_1 is 1.0
thisTemp <- temp[iD, 1]
if (thisTemp <= 130) {
logStatWC1 <- logUC1[1]
logStatWMg1 <- logUMg1[1]
logStatWMg2 <- logUMg2[1]
logStatWAl1 <- logUAl1[1]
logStatWSi1 <- logUSi1[1]
logStatWSi2 <- logUSi2[1]
logStatWFe1 <- logUFe1[1]
} else if (thisTemp > 130 && thisTemp <= 500) {
#// Add in interpolation here lburns
logStatWC1 <- logUC1[2] * (thisTemp - 130)/(500 - 130)
+ logUC1[1] * (500 - thisTemp)/(500 - 130)
logStatWMg1 <- logUMg1[2] * (thisTemp - 130)/(500 - 130)
+ logUMg1[1] * (500 - thisTemp)/(500 - 130)
logStatWMg2 <- logUMg2[2] * (thisTemp - 130)/(500 - 130)
+ logUMg2[1] * (500 - thisTemp)/(500 - 130)
logStatWAl1 <- logUAl1[2] * (thisTemp - 130)/(500 - 130)
+ logUAl1[1] * (500 - thisTemp)/(500 - 130)
logStatWSi1 <- logUSi1[2] * (thisTemp - 130)/(500 - 130)
+ logUSi1[1] * (500 - thisTemp)/(500 - 130)
logStatWSi2 <- logUSi2[2] * (thisTemp - 130)/(500 - 130)
+ logUSi2[1] * (500 - thisTemp)/(500 - 130)
logStatWFe1 <- logUFe1[2] * (thisTemp - 130)/(500 - 130)
+ logUFe1[1] * (500 - thisTemp)/(500 - 130)
} else if (thisTemp > 500 && thisTemp <= 3000) {
logStatWC1 <- logUC1[3] * (thisTemp - 500)/(3000 - 500)
+ logUC1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWMg1 <- logUMg1[3] * (thisTemp - 500)/(3000 - 500)
+ logUMg1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWMg2 <- logUMg2[3] * (thisTemp - 500)/(3000 - 500)
+ logUMg2[2] * (3000 - thisTemp)/(3000 - 500)
logStatWAl1 <- logUAl1[3] * (thisTemp - 500)/(3000 - 500)
+ logUAl1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWSi1 <- logUSi1[3] * (thisTemp - 500)/(3000 - 500)
+ logUSi1[2] * (3000 - thisTemp)/(3000 - 500)
logStatWSi2 <- logUSi2[3] * (thisTemp - 500)/(3000 - 500)
+ logUSi2[2] * (3000 - thisTemp)/(3000 - 500)
logStatWFe1 <- logUFe1[3] * (thisTemp - 500)/(3000 - 500)
+ logUFe1[2] * (3000 - thisTemp)/(3000 - 500)
} else if (thisTemp > 3000 && thisTemp <= 8000) {
logStatWC1 <- logUC1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUC1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWMg1 <- logUMg1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUMg1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWMg2 <- logUMg2[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUMg2[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWAl1 <- logUAl1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUAl1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWSi1 <- logUSi1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUSi1[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWSi2 <- logUSi2[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUSi2[3] * (8000 - thisTemp)/(8000 - 3000)
logStatWFe1 <- logUFe1[4] * (thisTemp - 3000)/(8000 - 3000)
+ logUFe1[3] * (8000 - thisTemp)/(8000 - 3000)
} else if (thisTemp > 8000 && thisTemp < 10000) {
logStatWC1 <- logUC1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUC1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWMg1 <- logUMg1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUMg1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWMg2 <- logUMg2[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUMg2[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWAl1 <- logUAl1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUAl1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWSi1 <- logUSi1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUSi1[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWSi2 <- logUSi2[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUSi2[4] * (10000 - thisTemp)/(10000 - 8000)
logStatWFe1 <- logUFe1[5] * (thisTemp - 8000)/(10000 - 8000)
+ logUFe1[4] * (10000 - thisTemp)/(10000 - 8000)
} else {
# for temperatures greater than or equal to 10000
logStatWC1 <- logUC1[5]
logStatWMg1 <- logUMg1[5]
logStatWMg2 <- logUMg2[5]
logStatWAl1 <- logUAl1[5]
logStatWSi1 <- logUSi1[5]
logStatWSi2 <- logUSi2[5]
logStatWFe1 <- logUFe1[5]
}
logGroundPopsC1[iD] <- stagePops[1,6][iD] - logStatWC1
logGroundPopsMg1[iD] <- stagePops[1,12][iD] - logStatWMg1
logGroundPopsMg2[iD] <- stagePops[2,12][iD] - logStatWMg2
logGroundPopsAl1[iD] <- stagePops[1,13][iD] - logStatWAl1
logGroundPopsSi1[iD] <- stagePops[1,14][iD] - logStatWSi1
logGroundPopsSi2[iD] <- stagePops[2,14][iD] - logStatWSi2
logGroundPopsFe1[iD] <- stagePops[1,26][iD] - logStatWFe1
}
for (iL in 1:numLams) {
for (i in 1:numDeps) {
aC1[i] <- 0.0
aMg1[i] <- 0.0
aMg2[i] <- 0.0
aAl1[i] <- 0.0
aSi1[i] <- 0.0
aSi2[i] <- 0.0
aFe1[i] <- 0.0
}
waveno <- 1.0/lambdaScale[iL]
logFreq <- useful[,"logC"]-log(lambdaScale[i])
freq <- exp(logFreq)
stimEmLogExpHelp <- useful[,"logH"]+logFreq-useful[,"logK"]
if (freq >= 2.0761e15)
aC1 = opacC1(numDeps, temp, lambdaScale[iL], logGroundPopsC1)
if (freq >= 2.997925e+14)
aMg1 = opacMg1(numDeps, temp, lambdaScale[iL], logGroundPopsMg1)
if (freq >= 2.564306e15)
aMg2 = opacMg2(numDeps, temp, lambdaScale[iL], logGroundPopsMg2)
if (freq >= 1.443e15)
aAl1 = opacAl1(numDeps, temp, lambdaScale[iL], logGroundPopsAl1)
if (freq >= 2.997925e+14)
aSi1 = opacSi1(numDeps, temp, lambdaScale[iL], logGroundPopsSi1)
if (freq >= 7.6869872e14)
aSi2 = opacSi2(numDeps, temp, lambdaScale[iL], logGroundPopsSi2)
if (waveno >= 21000.0)
aFe1 = opacFe1(numDeps, temp, lambdaScale[iL], logGroundPopsFe1)
for (iD in 1:numDeps) {
kapBF <- 1.0e-99 #minimum safe value
stimEmLogExp <- stimEmLogExpHelp - temp[iD,2]
stimEmExp <- -1.0 * exp(stimEmLogExp)
stimEm <- 1.0 - exp(stimEmExp) #//LTE correction for stimulated emission
kapBF <- kapBF + aC1[iD] + aMg1[iD] + aMg2[iD] + aAl1[iD] + aSi1[iD] + aSi2[iD] + aFe1[iD]
}
}
masterBF
}
opacC1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to CI
waveno <- 1.0/lambda2
freq <- useful[,"c"]/lambda2
c1240 <- c1444 <- tkev <- aC1 <- rep("numeric", numDeps)
freq1 <- sigma <- 0.0
for (i in 1:numDeps) {
logTkev <- temp[i,2]+useful[,"logK"]-useful[,"logEv"]
tkev[i] <- exp(logTkev)
}
#initialize some quantities for each new model atmosphere
for (i in 1:numDeps) {
c1240[i] <- 5.0*exp(-1.264/tkev[i])
c1444[i] <- exp(-2.683/tkev[i])
}
#initialize some quantities for each new model atmosphere or new frequency
#Luo, D. and Pradhan, A.K. 1989, J.Phys. B, 22. 3377-3395.
#Burke, P.G. and Taylor, K.T. 1979, J.Phys. B, 12, 2971-2984.
ryd <- 109732.298
xs0 <- xs1 <- 0.0
xd0 <- xd1 <- xd2 <- 0.0
x1444 <- x1240 <- x1100 <- 0.0
if (freq >= 2.7254E15) {
arg <- -16.80-((waveno-90777.00)/3.0/ryd)
x1100 <- (10**arg)*seaton(2.7254E15, 1.219E-7, 2.0E0, 3.317E0, freq)
}
if (freq >= 2.4196E15) {
arg <- -16.80-((waveno-80627.76)/3.0/ryd)
xd0 <- 10**arg
eeps <- (waveno-93917.0)*2.0/9230.0
aa <- 22.0E-18
bb <- 26.0E-18
xd1 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
eeps <- (waveno-111130.0)*2.0/2743.0
aa <- -10.5E-18
bb <- 46.0E-18
xd2 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
x1240 <- xd0+xd1+xd2
}
if (freq >= 2.0761E15) {
arg <- -16.80-((waveno-69172.40)/3.0/ryd)
xs0 <- 10**arg
eeps <- (waveno-97700.0)*2.0/2743.0
aa <- 68.0E-18
bb <- 118.0E-18
xs1 <- ((aa*eeps)+bb)/((eeps**2)+1.0)
x1444 <- xs0+xs1
}
for (i in 1:numDeps) {
if (freq >= 2.0761E15) {
sigma <- (x1100*9.0+xs1240*c1240[i]+x1444*c1444[i])
aC1[i] <- sigma*exp(logGroundPops[i])
}
}
aC1
}
seaton <- function(freq0, xsect, power, a, freq) {
#This function is a general representation for b-f absorption above
#a given ionization limit freq0, using cross-section xsect
freqratio <- freq0/freq
help <- floor((2*power)+0.01)
seaton <- xsect*(a+freqratio*(1.0-a))*sqrt(freqratio**help)
seaton
}
opacMg1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Mg I
sigma <- 0.0
aMg1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 7)
dt <- rep("numeric", 100)
nt <- rep("numeric", 100)
#4000K
peach0 = c(-42.474, -41.808, -41.273, -45.583, -44.324, -50.969,
-50.633, -53.028, -51.785, -52.285, -52.028, -52.384,
-52.363, -54.704, -54.359)
#5000 K
peach1 =c(-42.350, -41.735, -41.223, -44.008, -42.747, -48.388,
-48.026, -49.643, -48.352, -48.797, -48.540, -48.876,
-48.856, -50.772, -50.349)
#6000 K
peach2 =c(-42.109, -41.582, -41.114, -42.957, -41.694, -46.630,
-46.220, -47.367, -46.050, -46.453, -46.196, -46.513,
-46.493, -48.107, -47.643)
#7000 K
peach3 =c(-41.795, -41.363, -40.951, -42.205, -40.939, -45.344,
-44.859, -45.729, -44.393, -44.765, -44.507, -44.806,
-44.786, -46.176, -45.685)
#8000 K
peach4 =c(-41.467, -41.115, -40.755, -41.639, -40.370, -44.355,
-43.803, -44.491, -43.140, -43.486, -43.227, -43.509,
-43.489, -44.707, -44.198)
#9000 K
peach5 =c(-41.159, -40.866, -40.549, -41.198, -39.925, -43.568,
-42.957, -43.520, -42.157, -42.480, -42.222, -42.488,
-42.467, -43.549, -43.027)
#10000 K
peach6 =c(-40.883, -40.631, -40.347, -40.841, -39.566, -42.924,
-42.264, -42.736, -41.363, -41.668, -41.408, -41.660,
-41.639, -42.611, -42.418)
peach <- rbind(peach0, peach1, peach2, peach3, peach4, peach5, peach6)
freqMg <- c(1.9341452e15, 1.8488510e15, 1.1925797e15,
7.9804046e14, 4.5772110e14, 4.1440977e14,
4.1113514e14)
flog <- c(35.23123, 35.19844, 35.15334, 34.71490, 34.31318,
33.75728, 33.65788, 33.64994, 33.43947)
tlg <- c(8.29405, 8.51719, 8.69951, 8.85367, 8.98720, 9.10498,
9.21034)
freq1 <- 0.0
for (i in 1:numDeps) {
thelp <- floor((temp[i,1]/1000.0))-3
nn <- max(min(6, thelp), 1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
freq1 <- freq
nn <- 0
for (i in 1:7) {
if (freq > freqMg[n]) {
break
}
n <- n+1
}
if (freq <= freqMg[6]) {
nn <- 7
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
dd1 <- 1.0-dd
for (it in 1:numDeps) {
if (freq >= 2.997925E14) {
nn <- nt[i]
sigma <- exp(xx[nn]*(1.0E0-dt[i]))+(xx[nn+1]*dt[i])
aMg1[i] <- sigma*exp(logGroundPops[i])
}
}
aMg1
}
opacMg2 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Mg II
sigma <- 0.0
aMg2 <- tkev <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
c1169 <- rep("numeric", 100)
freq1 <- x824 <- x1169 <- 0.0
for (i in 1:numDeps) {
logTkev <- temp[i,2]+useful[,"logK"]-useful[,"logEv"]
tkev[i] <- exp(logTkev)
}
for (i in 1:numDeps) {
c1169[i] <- 6.0*exp(-4.43E0/tkev[i])
}
if (freq >= 3.635492E15) {
x824 <- seaton(3.635492E15, 1.40E-19, 4.0E0, 6.7E0, freq)
} else {
x824 <- 1.0E-99
}
if (freq >= 2.564306E15) {
x1169 <- 5.11E-19*((2.564306E15/freq)**3.0)
} else {
x1169 <- 1.0E-99
}
for (i in 1:numDeps) {
if (x1169 >= 1.0E-90) {
sigma <- (x824*2.0+x1169*c1169[i])
aMg2[i] <- sigma*exp(logGroundPops[i])
}
}
aMg2
}
opacAl1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Al I
sigma <- 0.0
aAl1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
for (i in 1:numDeps) {
if (freq >= 1.443E15) {
sigma <- 6.0*6.5E-17*((1.443E15/freq)**5.0)
aAl1[i] <- sigma*exp(logGroundPops[i])
}
}
aAl1
}
opacSi1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Si I
sigma <- 0.0
aSi1 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 9)
dt <- nt <- rep("numeric", 100)
#4000
peach0 <- c(38.136, 37.834, 37.898, 40.737, 40.581, 45.521,
45.520, 55.068, 53.868, 54.133, 54.051, 54.442,
54.320, 55.691, 55.661, 55.973, 55.922, 56.828, 56.657)
#5000
peach1 <- c(38.138, 37.839, 37.898, 40.319, 40.164, 44.456,
44.455, 51.783, 50.369, 50.597, 50.514, 50.854,
50.722, 51.965, 51.933, 52.193, 52.141, 52.821, 52.653)
#6000
peach2 <- c(38.140, 37.843, 37.897, 40.047, 39.893, 43.753,
43.752, 49.553, 48.031, 48.233, 48.150, 48.455,
48.313, 49.444, 49.412, 49.630, 49.577, 50.110, 49.944)
#7000
peach3 <- c(38.141, 37.847, 37.897, 39.855, 39.702, 43.254,
43.251, 47.942, 46.355, 46.539, 46.454, 46.733,
46.583, 47.615, 47.582, 47.769, 47.715, 48.146, 47.983)
#8000
peach4 <- c(38.143, 37.850, 37.897, 39.714, 39.561, 42.878,
42.871, 46.723, 45.092, 45.261, 45.176, 45.433,
45.277, 46.221, 46.188, 46.349, 46.295, 46.654, 46.491)
#9000
peach5 <- c(38.144, 37.853, 37.896, 39.604, 39.452, 42.580,
42.569, 45.768, 44.104, 44.262, 44.175, 44.415,
44.251, 45.119, 45.085, 45.226, 45.172, 45.477, 45.315)
#10000
peach6 <- c(38.144, 37.855, 37.895, 39.517, 39.366, 42.332,
42.315, 44.997, 43.308, 43.456, 43.368, 43.592,
43.423, 44.223, 44.189, 44.314, 44.259, 44.522, 44.360)
#11000
peach7 <- c(38.145, 37.857, 37.895, 39.445, 39.295, 42.119,
42.094, 44.360, 42.652, 42.790, 42.702, 42.912,
42.738, 43.478, 43.445, 43.555, 43.500, 43.730, 43.569)
#12000
peach8 <- c(38.145, 37.858, 37.894, 39.385, 39.235, 41.930,
41.896, 43.823, 42.100, 42.230, 42.141, 42.340,
42.160, 42.848, 42.813, 42.913, 42.858, 43.061, 42.901)
peach <- rbind(peach0, peach1, peach2, peach3, peach4, peach5, peach6,
peach7, peach8)
freqSi <- c(2.1413750e15, 1.9723165e15, 1.7879689e15,
1.5152920e15, 5.5723927e14, 5.3295914e14,
4.7886458e14, 4.7216422e14, 4.6185133e14)
flog <- c(35.45438, 35.30022, 35.21799, 35.11986, 34.95438,
33.95402, 33.90947, 33.80244, 33.78835, 33.76626,
33.70518)
tlg <- c(8.29405, 8.51719, 8.69951, 8.85367, 8.98720, 9.10498,
9.21034, 9.30565, 9.39266)
freq1 <- 0.0
for (i in 1:numDeps) {
thelp <- floor(temp[i,1]/1000.0)-3
nn <- max(min(8,thelp),1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
freq1 <- freq
nn <- 0
for (n in 1:9) {
if (freq > freqSi[n]) {
break
}
nn <- nn+1
}
if (freq <= freqSi[8]) {
nn <- 9
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
dd1 <- 1.0-dd
for (it in 1:9) {
if (freq >= 2.997925E14) {
nn <- nt[i]
sigma <- (9.0*exp(-(xx[nn]*(1.-dt[i])+xx[nn+1]*dt[i])))
sSi1[i] <- sigma*exp(logGroundPops[i])
}
}
aSi1
}
opacSi2 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Si II
sigma <- 0.0
aSi2 <- rep("numeric", numDeps)
freq <- useful[,"c"]/lambda2
freqlg <- log(freq)
xx <- rep("numeric", 6)
dt <- nt <- rep("numeric", 100)
#c 10000
peach0 <- c( -43.8941, -42.2444, -40.6054, -54.2389, -50.4108,
-52.0936, -51.9548, -54.2407, -52.7355, -53.5387,
-53.2417, -53.5097, -54.0561, -53.8469)
#c 12000
peach1 <- c( -43.8941, -42.2444, -40.6054, -52.2906, -48.4892,
-50.0741, -49.9371, -51.7319, -50.2218, -50.9189,
-50.6234, -50.8535, -51.2365, -51.0256)
#c 14000
peach2 <- c( -43.8941, -42.2444, -40.6054, -50.8799, -47.1090,
-48.5999, -48.4647, -49.9178, -48.4059, -49.0200,
-48.7252, -48.9263, -49.1980, -48.9860)
#c 16000
peach3 <- c( -43.8941, -42.2444, -40.6054, -49.8033, -46.0672,
-47.4676, -47.3340, -48.5395, -47.0267, -47.5750,
-47.2810, -47.4586, -47.6497, -47.4368)
#c 18000
peach4 <- c( -43.8941, -42.2444, -40.6054, -48.9485, -45.2510,
-46.5649, -46.4333, -47.4529, -45.9402, -46.4341,
-46.1410, -46.2994, -46.4302, -46.2162)
#c 20000
peach5 <- c( -43.8941, -42.2444, -40.6054, -48.2490, -44.5933,
-45.8246, -45.6947, -46.5709, -45.0592, -45.5082,
-45.2153, -45.3581, -45.4414, -45.2266)
peach <- c(peach0, peach1, peach2, peach3, peach4, peach5)
freqSi <- c(4.9965417e15, 3.9466738e15, 1.5736321e15,
1.5171539e15, 9.2378947e14, 8.3825004e14,
7.6869872e14)
flog <- c(36.32984, 36.14752, 35.91165, 34.99216, 34.95561,
34.45951, 34.36234, 34.27572, 34.20161)
tlg <- c(9.21034, 9.39266, 9.54681, 9.68034, 9.79813, 9.90349)
freq1 <- 0.0
#setup some data upon first entrance with a new model atmosphere
for (i in 1:numDeps) {
thelp <- floor(temp[i,1]/2000.0)-4.0
nn <- max(min(5, thelp), 1)-1
nt[i] <- nn
dt[i] <- (temp[i,2]-tlg[nn])/(tlg[nn+1]-tlg[nn])
}
#initialize some quantities for each new model atmosphere or frequency
freq1 <- freq
nn <- 0
for (n in 1:7) {
if (freq > freqSi[n]) {
break
}
nn <- nn + 1
}
if (freq <= freqSi[6]) {
nn <- 7
}
dd <- (freqlg-flog[nn])/(flog[nn+1]-flog[nn])
if (nn > 1) {
nn <- 2*nn-2
}
if (nn == 13) {
nn <- 12
}
dd1 <- 1.0 - dd
for (it in 1:6) {
xx[it] <- peach[nn+1,it]*dd+peach[nn,it]*dd1
}
for (i in 1:numDeps) {
if (freq >= 7.6869872E14) {
nn <- nt[i]
sigma <- (6.0*exp(xx[nn]*(1.0-dt[i])+xx[nn+1]*dt[i]))
aSi2[i] <- sigma*exp(logGroundPops[i])
}
}
aSi2
}
opacFe1 <- function(numDeps, temp, lambda2, logGroundPops) {
#This routine computes the bound-free absorption due to Fe I
sigma <- 0.0
aFe1 <- rep("numeric", numDeps)
#cross-section is zero below threshold, so initialize
for (i in 1:numDeps) {
aFe1[i] <- 0.0
}
wavenoo <- 1.0/lambda2
freq <- useful[,"c"]/lambda2
bolt <- matrix(0.0, 100, 48)
xsect <- rep("numeric", 1:48)
gg <- c(25.0, 35.0, 21.0, 15.0, 9.0, 35.0, 33.0, 21.0, 27.0, 49.0,
9.0, 21.0, 27.0, 9.0, 9.0, 25.0, 33.0, 15.0, 35.0, 3.0,
5.0, 11.0, 15.0, 13.0, 15.0, 9.0, 21.0, 15.0, 21.0, 25.0, 35.0,
9.0, 5.0, 45.0, 27.0, 21.0, 15.0, 21.0, 15.0, 25.0, 21.0,
35.0, 5.0, 15.0, 45.0, 35.0, 55.0, 25.0)
ee<- c(500.0, 7500.0, 12500.0, 17500.0, 19000.0, 19500.0, 19500.0,
21000.0, 22000.0, 23000.0, 23000.0, 24000.0, 24000.0, 24500.0,
24500.0, 26000.0, 26500.0, 26500.0, 27000.0, 27500.0, 28500.0,
29000.0, 29500.0, 29500.0, 29500.0, 30000.0, 31500.0, 31500.0,
33500.0, 33500.0, 34000.0, 34500.0, 34500.0, 35000.0, 35500.0,
37000.0, 37000.0, 37000.0, 38500.0, 40000.0, 40000.0, 41000.0,
41000.0, 43000.0, 43000.0, 43000.0, 43000.0, 44000.0)
wno <- c(63500.0, 58500.0, 53500.0, 59500.0, 45000.0, 44500.0, 44500.0,
43000.0, 58000.0, 41000.0, 54000.0, 40000.0, 40000.0, 57500.0,
55500.0, 38000.0, 57500.0, 57500.0, 37000.0, 54500.0, 53500.0,
55000.0, 34500.0, 34500.0, 34500.0, 34000.0, 32500.0, 32500.0,
32500.0, 32500.0, 32000.0, 29500.0, 29500.0, 31000.0, 30500.0,
29000.0, 27000.0, 54000.0, 27500.0, 24000.0, 47000.0, 23000.0,
44000.0, 42000.0, 42000.0, 21000.0, 42000.0, 42000.0)
freq1 <- 0.0
#set up some data upon first entrance with a new model atmosphere
for (i in 1:numDeps) {
hkt <- 6.6256E-27/(1.38054E-16*temp[i,1])
for (k in 1:48) {
bolt[i,k] <- gg[k]*exp(-ee[k]*useful[,"c"]*hkt)
}
}
#initialize some quantities for each new model atmosphere or frequency
freq1 <- freq
if (waveno >= 21000.0) {
for (k in 1:48) {
xsect[k] <- 0.0
if (wno[k] < waveno) {
xsect[k] <- 3.0E-18/(1.0+((wno[k]+3000.0-waveno)/wno[k])/0.1)**4
}
}
}
for (i in 1:numDeps) {
aFe1[i] <- 0.0
if (waveno >= 21000.0) {
for (k in 1:48) {
aFe[i] <- 0.0
sigma <- aFe1[i]+xsect[k]*bolt[i,k]
aFe1[i] <- sigma*exp(logGroundPops[i])
}
}
}
aFe1
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.