## 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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.