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

source("./R/lambdaGrid.R")
source("./R/scaleSolar.R")
source("./R/lineGrid.R")

stellarVariables <- readRDS("./data/starvars.rds")
restart <- readRDS("./data/restart.rds")
restartParams <- readRDS("./data/restartParams.rds")
useful <- readRDS("./data/useful.rds")
options(warn=-1)
require(glue)

#Custom filename tags to distinguish from other runs
projectName <- "Project"
runVersion <- "Run"

#Default plot
#Select ONE only:

#makePlot <- "structure"
#makePlot <- "sed"
makePlot <- "spectrum"
#makePlot <- "ldc"
#makePlot <- "ft"
#makePlot <- "tlaLine"

#Spectrum synthesis mode
# - uses model in restart.R with minimal structure calculation
specSynMode <- TRUE

#Arg 0:  Effective temperature, Teff, in K
teff <- stellarVariables["defaultInput","teff"]
#Logarithmic surface gravity, g, in cm/s/s
logg <- stellarVariables["defaultInput","logg"]
#Arg 1:  Linear scale factor for solar Rosseland opacity distribution
log10ZScale <- stellarVariables["defaultInput","log10ZScale"]
#Arg 7:  Starting wavelength for spectrum synthesis
lambdaStart <- stellarVariables["defaultInput","lambdaStart"]
vv#Arg 8:  stopping wavelength for spectrum synthesis
lambdaStop <- stellarVariables["defaultInput","lambdaStop"]

fileStem <- glue(projectName,"-",
                 teff,"-",
                 logg,"-",
                 log10ZScale,"-",
                 lambdaStart,"-",
                 lambdaStop,"-",
                 runVersion)

absPath <- glue(getwd(),"/")

#Colour pallet for plotting
numPal <- 12
palette <- vector("numeric",numPal)
delPal <- 0.04

pallette <- 0.481-numPal*delPal
numClrs <- length(palette)

outPath <- glue(absPath,"Outputs/")

outFileString <- glue("Writing to files ",
			outPath,
			fileStem,
			".*")

#True for Voigt computed with true concolution instead of power-law
#expansion approx - probably not working well right now
ifVoigt <- FALSE

#Scattering term in line source fn - not yet enabled
ifScatt <- FALSE

#Arg 3 {  Stellar mass, M, in solar masses
massStar <- stellarVariables["defaultInput","massStar"]

#Sanity Check
F0Vtemp <- 7300.0

if (teff < 3000.0)
  teff <- 3000.0

if (teff > 50000.0)
  teff <- 50000.0

#Safe initialization
minLogg <- 3.0
minLoggStr <- "3.0"

#Logg limit is strongly Teff-dependent {
if (teff <= 4000.0) { minLogg <- 0.0
} else if ((teff > 4000.0) && (teff <= 5000.0)) { minLogg <- 0.5
} else if ((teff > 5000.0) && (teff <= 6000.0)) { minLogg <- 1.5
} else if ((teff > 6000.0) && (teff <= 7000.0)) { minLogg <- 2.0
} else if ((teff > 7000.0) && (teff < 9000.0)) { minLogg <- 2.5
} else if (teff >= 9000.0) { minLogg <- 3.0 }

if (logg < minLogg) { logg <- minLogg }
if (logg < minLogg) { logg <- minLogg }
if (logg > 7.0) { logg <- 7.0 }
if (log10ZScale < -3.0) { log10ZScale <- -3.0 }
if (log10ZScale > 1.0) { log10ZScale <- 1.0 }
if (massStar < 0.1) { massStar <- 0.1 }
if (massStar > 20.0) { massStar <- 20.0 }

grav <- 10**logg
zScale <- 10**log10ZScale

#Arg 4:  Microturbulence, xi_T, in km/s
xiT <- stellarVariables["defaultInput","xiT"]

if (xiT < 0.0) { xiT <- 0.0 }
if (xiT > 8.0) { xiT <- 8.0 }

#Add placeholders for metallicity controls lburns
logHeFe <- stellarVariables["defaultInput","logHeFe"]
logCO <- stellarVariables["defaultInput","logCO"]
logAlphaFe <- stellarVariables["defaultInput","logAlphaFe"]

#For new metallicity commands lburns
#For logHeFeL (lburns)

if (logHeFe < -1.0) { logHeFe <- -1.0 }
if (logHeFe > 1.0) { logHeFe <- 1.0 }

#// For logCO: (lburns)

if (logCO < -2.0) { logCO <- -2.0 }
if (logCO > 2.0) { logCO <- 2.0 }

#For logAlphaFe: (lburns)

if (logAlphaFe < -0.5) { logAlphaFe <- -0.5 }
if (logAlphaFe > 0.5) { logAlphaFe <- 0.5 }

#Arg 5: Minimum ration of monochromatic line center to background
#continuous extinction for inclusion of linein spectrum
lineThresh <- stellarVariables["defaultInput","lineThresh"]

if (lineThresh < -4.0) { lineThresh <- -4.0 }
if (lineThresh > 6.0) { lineThresh <- 6.0 }

#Arg 6: Minimum ratio of monochromatic line center to background
#continuous
voigtThresh <- stellarVariables["defaultInput","voigtThresh"]

if (voigtThresh < lineThresh) { voigtThresh <- lineThresh }
if (voigtThresh > 6.0) { voigtThresh <- 6.0 }

#User defined spectrum synthesis region
lamUV <- 260.0
lamIR <- 2600.0

if (lambdaStart < lamUV) { lambdaStart <- lamUV }
if (lambdaStart > (lamIR - 1.0)) { lambdaStart <- lamIR - 1.0 }

if (lambdaStop < (lamUV + 1.0)) { lambdaStop <- lamUV + 1.0 }
if (lambdaStop > lamIR) { lambdaStop <- lamIR }

#Prevent negative or zero lambda range
#0.5 nm <- 5 A
if (lambdaStop <= lambdaStart) { lambdaStop <- lambdaStart + 0.5 }
if (lambdaStop > lamIR) { lambdaStop <- lamIR }

#Conversion of nm to cm
nm2cm <- 1.0e-7
cm2nm <- 1.0e7
lambdaStart <- nm2cm*lambdaStart
lambdaStop <- nm2cm*lambdaStop
lamUV <- nm2cm*lamUV
lamIR <- nm2cm*lamIR

#Arg 9: Line sampling selection (fine or coarse)
sampling <- stellarVariables["defaultInput","sampling"]
vacAir <- stellarVariables["defaultInput","vacAir"]

#Arg 10: Lorentzian line broadening enhancements
logGammaCol <- stellarVariables["defaultInput","logGammaCol"]

if (logGammaCol < 0.0) { logGammaCol <- 0.0 }
if (logGammaCol > 1.0) { logGammaCol <- 1.0 }

#Arg 11: log_10 gray mass extinction fudge
logKapFudge <- stellarVariables["defaultInput","logKapFudge"]

if (logKapFudge < -2.0) { logKapFudge <- -2.0 }
if (logKapFudge > 2.0) { logKapFudge <- 2.0 }

#Arg 12: Macroturbulent velocity broadening parameter (sigma) (km/s)
macroV <- stellarVariables["defaultInput","macroV"]

if (macroV < 0.0) { macroV <- 0.0 }
if (macroV > 8.0) { macroV <- 8.0 }

#Arg 13: Surface equatorial linear rotational velocity (km/s)
rotV <- stellarVariables["defaultInput","rotV"]

if (rotV < 0.0 ) { rotV <- 0.0 }
if (rotV > 300.0) { rotV <- 300.0 }

#Arg 14: Inclination of rotation axis w.r.t. line-of-sight (degrees)
rotI <- stellarVariables["defaultInput","rotI"]

if (rotI < 0.0) { rotI <- 0.0 }
if (rotI > 90.0) { rotI <- 90.0 }

#Arg 15: Number of outer HSE-EOS-Opac iterations
nOuterIter <- stellarVariables["defaultInput","nOuterIter"]

if (nOuterIter < 1) { nOuterIter <- 1 }
if (nOuterIter > 30) { nOuterIter <- 30 }

#Arg 16: Number of inner Pe-IonFrac iterations
nInnerIter <- stellarVariables["defaultInput","nInnerIter"]

if (nInnerIter < 1) { nInnerIter <- 1 }
if (nInnerIter > 30) { nInnerIter <- 30 }

#Arg 17: If TiO JOLA bands should be included
ifTiO <- stellarVariables["defaultInput","ifTiO"]

#Degree to radions
inclntn <- (pi*rotI)/180.0
vsini <- rotV*sin(inclntn)

#Arg 18: Wavelength of narrow Gaussian filter in nm
diskLambda <- stellarVariables["defaultInput","diskLambda"]

if (diskLambda < lamUV) { diskLambda <- lamUV }
if (diskLambda > lamIR) { diskLambda <- lamIR }

#Arg 19: Bandwidth, Sigma, of narrow Gaussian filter in nm
diskSigma <- stellarVariables["defaultInput","diskSigma"]

if (diskSigma < 0.005) { diskSigma <- 0.005 }
if (diskSigma > 10.0) { diskSigma <- 10.0 }

#Arg 20: Radial velocity of star in km/s
RV  <- stellarVariables["defaultInput","RV"]

if (RV < -200.0) { RV <- -200.0 }
if (RV > 200.0) { RV <- 200.0 }

#Representative spectral line and associated atomic parameters
userLam0 <- stellarVariables["defaultInput","userLam0"]

if (userLam0 < 260.0) { userLam0 <- 260.0 }
if (userLam0 > 2600.0) { userLam0 <- 2600.0 }

userA12 <- stellarVariables["defaultInput","userA12"]

if (userA12 < 2.0) { userA12 <- 2.0 }
if (userA12 > 11.0) { userA12 <- 11.0 } #Upper limit set high to accomodate Helium!

userLogF <- stellarVariables["defaultInput","userLogF"]

if (userLogF < -6.0) { userLogF <- -6.0 }
if (userLogF > 1.0) { userLogF <- 1.0 }

userStage <- stellarVariables["defaultInput","userStage"]

if ( (userStage != 0) && (userStage != 1) && (userStage != 2)
		&& (userStage != 3) ) {
    userStage <- 0
    userStageStr <- "I"
}

userChiI1 <- stellarVariables["defaultInput","userChiI1"]

if (userChiI1 < 3.0) { userChiI1 <- 3.0 }
if (userChiI1 > 25.0) { userChiI1 <- 25.0 }

userChiI2 <- stellarVariables["defaultInput","userChiI2"]

if (userChiI2 < 5.0) { userChiI2 <- 5.0 }
if (userChiI2 > 55.0) { userChiI2 <- 55.0 }

userChiI3 <- stellarVariables["defaultInput","userChiI3"]

if (userChiI3 < 5.0) { userChiI3 <- 5.0 }
if (userChiI3 > 55.0) { userChiI3 <- 55.0 }

userChiI4 <- stellarVariables["defaultInput","userChiI4"]

if (userChiI4 < 5.0) { userChiI4 <- 5.0 }
if (userChiI4 > 55.0) { userChiI4 <- 55.0 }

userChiL <- stellarVariables["defaultInput","userChiL"]

if (userChiL < 0.0) { userChiL <- 0.0 } #Ground state case!
if ( (userStage == 0) && (userChiL >= userChiI1) ) { userChiL <- 0.9 * userChiI1 }
#Note: Upper limit of chiL depends on value of chiI1 above!
if (userChiL < 0.0) { userChiL <- 0.0 } #Ground state case!
if ( (userStage == 0) && (userChiL >= userChiI1) ) { userChiL <- 0.9 * userChiI1 }
if ( (userStage == 1) && (userChiL >= userChiI2) ) { userChiL <- 0.9 * userChiI2 }
if ( (userStage == 2) && (userChiL >= userChiI3) ) { userChiL <- 0.9 * userChiI3 }
if ( (userStage == 3) && (userChiL >= userChiI4) ) { userChiL <- 0.9 * userChiI4 }

userGw1 <- stellarVariables["defaultInput","userGw1"]

if (userGw1 < 1.0) { userGw1 <- 1.0 }
if (userGw1 > 100.0) { userGw1 <- 100.0 }

userGw2 <- stellarVariables["defaultInput","userGw2"]

if (userGw2 < 1.0) { userGw2 <- 1.0 }
if (userGw2 > 100.0) { userGw2 <- 100.0 }

userGw3 <- stellarVariables["defaultInput","userGw3"]

if (userGw3 < 1.0) { userGw3 <- 1.0 }
if (userGw3 > 100.0) { userGw3 <- 100.0 }

userGw4 <- stellarVariables["defaultInput","userGw4"]

if (userGw4 < 1.0) { userGw4 <- 1.0 }
if (userGw4 > 100.0) { userGw4 <- 100.0 }

userGwL <- stellarVariables["defaultInput","userGwL"]

if (userGwL < 1.0) { userGwL <- 1.0 }
if (userGwL > 100.0) { userGwL <- 100.0 }

userMass <- stellarVariables["defaultInput","userMass"]

if (userMass < 1.0) { userMass <- 1.0 }
if (userMass > 200.0) { userMass <- 200.0 }

userLogGammaCol <- stellarVariables["defaultInput","userLogGammaCol"]

if (userLogGammaCol < 0.0) { userLogGammaCol <- 0.0 }
if (userLogGammaCol > 1.0) { userLogGammaCol <- 1.0 }

#Line Ceter lambda from nm to cm
userLam0 <- userLam0 * nm2cm

#Echo active input parameters
inputParameterString <- glue("Teff ", teff, " logg ", logg, " [Fe/H] ", log10ZScale,
                             " massStar ", massStar, " xiT ", xiT, " HeFe ", logHeFe,
                             " CO ", logCO, " AlfFe ", logAlphaFe, " lineThresh ",
                             lineThresh, " voigtThresh ", voigtThresh, " lambda0 ",
                             lambdaStart, " lambda1 ", lambdaStop, " logGammaCol ",
                             logGammaCol, " logKapFudge ", logKapFudge, " macroV ",
                             macroV, " rotV ", rotV, " rotI ", rotI, " RV ", RV,
                             " nInner ", nInnerIter, " nOuter ", nOuterIter,
                             " ifTiO ", ifTiO, " sampling ", sampling)

##
#OPACITY PROBLEM #1 - logFudgeTune: late type star continuous opacity needs to have
#been multiplied by 10.0^0.5=3 for T_kin(tau~1) to fall around Teff and SED t.o look
#like B_lambda(Trad=Teff) - related to Opacity problem #2 in LineKappa.lineKap() - ??
##

logFudgeTune <- 0.0
#Makes the Balmer lines show up around A0
if (teff <= F0Vtemp) { logFudgeTune <- 0.5 }
if (teff > F0Vtemp) { logFudgeTune <- 0.0 }

logTotalFudge <- logKapFudge + logFudgeTune
logE <- log10(exp(1))
logE10 <- log10(10)

#Gray structure and Voigt line code
#optical depth grid
numDeps <- 48
log10MinDepth <- -6.0
log10MaxDepth <- 2.0

#Test start wavelength (cm), end wavelength (cm), test number of lambda
lamSetup <- c(260 * nm2cm, 2600.0, 250)
numLams <- lamSetup[3]

#Continuum lambda scale (nm)
lambdaScale <- lambdaGrid(numLams, lamSetup)

#Solar parameters
teffSun <- 5788.0
loggSun <- 4.44
gravSun <- 10**loggSun
log10ZScaleSun <- 0.0
zScaleSun <- exp(log10ZScaleSun)

#Solar units
massSun <- 1.0
radiusSun <- 1.0
logRadius <- 0.5*log10(massStar)+log10(gravSun)-log10(grav)
radius <- exp(logRadius) #Solar radii
logLum <- 2.0*log10(radius)+4.0*log10(teff/teffSun)
bolLum <- exp(logLum) #L_Bol in solar luminosities
rSun <- 6.955e10
cgsRadius <- radius*rSun
omegaSini <- (1.0e5*vsini)/cgsRadius #Projected rotation rate in 1/sec
macroVkm <- macroV*1.0e5 #km/s to cm/s

#Composition by mass fraction - needed for opacity approximations and
#interior structure
massX <- 0.70 #Hydrogen
massY <- 0.28 #Helium
massZSun <- 0.02 #"metals"
massZ <- massZSun*zScale #approximation

nelemAbnd <- 41
numStages <- 7

#Kurucz codes
nome <- seq(from = 100, by = 100, length.out = nelemAbnd)
append(nome, c(5600, 5700, 5500, 3200))

#Log_10 `A_12` values
eheu <- c(12.00, 10.93, 1.05, 1.38, 2.70, 8.43, 7.83, 8.69, 4.56, 7.93, 6.24,
	  7.60, 6.45, 7.51, 5.41, 7.12, 5.50, 6.40, 5.03, 6.34, 3.15, 4.95,
	  3.93, 5.64, 5.43, 7.50, 4.99, 6.22, 4.19, 4.56, 3.04, 3.25, 2.52,
	  2.87, 2.21, 2.58, 1.46, 2.18, 1.10, 1.12, 3.65)
logAz <- vector("numeric", nelemAbnd)
cname <- vector("character", nelemAbnd)
logNH <- vector("numeric", numDeps)
logNz <- matrix(0.0, numDeps, nelemAbnd)
masterStagePops <- matrix(0.0, numDeps, numStages)

#Associate diatomic molecules with each element that forms significant molecules
numAssocMols <- 4
cname <- vector("character",nelemAbnd)
cnameMols <- matrix("",numAssocMols,nelemAbnd)

cname[1] <- "H"
cnameMols[1][1] <- "H2"
cnameMols[2][1] <- "H2+"
cnameMols[3][1] <- "CH"
cnameMols[4][1] <- "OH"
cname[2] <- "He"
cname[3] <- "Li"
cname[4] <- "Be"
cname[5] <- "B"
cname[6] <- "C"
cnameMols[1][6] <- "CH"
cnameMols[2][6] <- "CO"
cnameMols[3][6] <- "CN"
cnameMols[4][6] <- "C2"
cname[7] <- "N"
cnameMols[1][7] <- "NH"
cnameMols[2][7] <- "NO"
cnameMols[3][7] <- "CN"
cnameMols[4][7] <- "N2"
cname[8] <- "O"
cnameMols[1][8] <- "OH"
cnameMols[2][8] <- "CO"
cnameMols[3][8] <- "NO"
cnameMols[4][8] <- "O2"
cname[9] <- "F"
cname[10] <- "Ne"
cname[11] <- "Na"
cname[12] <- "Mg"
cnameMols[1][12] <- "MgH"
cname[13] <- "Al"
cname[14] <- "Si"
cnameMols[1][14] <- "SiO"
cname[15] <- "P"
cname[16] <- "S"
cname[17] <- "Cl"
cname[18] <- "Ar"
cname[19] <- "K"
cname[20] <- "Ca"
cnameMols[1][20] <- "CaH"
cnameMols[2][20] <- "CaO"
cname[21] <- "Sc"
cname[22] <- "Ti"
cnameMols[1][22] <- "TiO"
cname[23] <- "V"
cnameMols[1][23] <- "VO"
cname[24] <- "Cr"
cname[25] <- "Mn"
cname[26] <- "Fe"
cnameMols[1][26] <- "FeO"
cname[27] <- "Co"
cname[28] <- "Ni"
cname[29] <- "Cu"
cname[30] <- "Zn"
cname[31] <- "Ga"
cname[32] <- "Kr"
cname[33] <- "Rb"
cname[34] <- "Sr"
cname[35] <- "Y"
cname[36] <- "Zr"
cname[37] <- "Nb"
cname[38] <- "Ba"
cname[39] <- "La"
cname[40] <- "Cs"
cname[41] <- "Ge"

#Diatomic molecules
nMols <- 1
mname <- vector("character",nMols)
mnameA <- vector("character",nMols)
mnameB <- vector("character",nMols)

#CAUTION: The molecular number densitis, N_AB will be computed and will deplete
#the atomic species in this order
#Put anything where A or B is Hydrogen FIRST - HI is an inexhaustible reservior
#at low T
#Then rank molecules according to largest of A and B abundance, "weighted" by
#dissociation energy
#For constituent atomic species, A and B, always designate as 'A' whichever
#element participates in the 'fewest othere' molecules - we'll put A on the
#LHS of the molecular Saha equation

#mname[1] <- "H2"
#mnameA[1] <- "H"
#mnameB[1] <- "H"
#mname[2] <- "H2+"
#mnameA[2] <- "H"
#mnameB[2] <- "H"
#mname[3] <- "OH"
#mnameA[3] <- "O"
#mnameB[3] <- "H"
#mname[4] <- "CH"
#mnameA[4] <- "C"
#mnameB[4] <- "H"
#mname[5] <- "NH"
#mnameA[5] <- "N"
#mnameB[5] <- "H"
#mname[6] <- "MgH"
#mnameA[6] <- "Mg"
#mnameB[6] <- "H"
#mname[7] <- "CaH"
#mnameA[7] <- "Ca"
#mnameB[7] <- "H"
#mname[8] <- "O2"
#mnameA[8] <- "O"
#mnameB[8] <- "O"
#mname[9] <- "CO"
#mnameA[9] <- "C"
#mnameB[9] <- "O"
#mname[10] <- "C2"
#mnameA[10] <- "C"
#mnameB[10] <- "C"
#mname[11] <- "NO"
#mnameA[11] <- "N"
#mnameB[11] <- "O"
#mname[12] <- "CN"
#mnameA[12] <- "C"
#mnameB[12] <- "N"
#mname[13] <- "N2"
#mnameA[13] <- "N"
#mnameB[13] <- "N"
#mname[14] <- "FeO"
#mnameA[14] <- "Fe"
#mnameB[14] <- "O"
#mname[15] <- "SiO"
#mnameA[15] <- "Si"
#mnameB[15] <- "O"
#mname[16] <- "CaO"
#mnameA[16] <- "Ca"
#mnameB[16] <- "O"
mname[1] <- "TiO"
mnameA[1] <- "Ti"
mnameB[1] <- "O"
#mname[18] <- "VO"
#mnameA[18] <- "V"
#mnameB[18] <- "O"

#Set up for molecules with JOLA bands
jolaTeff <- 5000.0
numJola <- 3
jolaSpecies <- vector("character",numJola) #molecule name
jolaSystem <- vector("character",numJola) #band system
jolaDeltaLambda <- vector("numeric",numJola)

if (teff <= jolaTeff) {
  jolaSpecies[1] <- "TiO" #molecule name
  jolaSystem[1] <- "TiO_C3Delat_X3Delta" #band system
  jolaDeltaLambda[1] <- 0
  jolaSpecies[1] <- "TiO" #molecule name
  jolaSystem[2] <- "TiO_clPhi_alDelta" #band system
  jolaDeltaLambda[2] <- 1
  jolaSpecies[3] <- "TiO" #molecule name
  jolaSystem[3] <- "TiO_A3Phi_X3Delta" #band system
  jolaDeltaLambda[3] <- 1
}

ATot <- 0.0
thisAt <- 0.0
eheuScale <- 0.0

#Set value of eheuScale for new metallicity options
if (logHeFe != 0.0) {
  eheu[1] = eheu[1] + logHeFe
}

if (logAlphaFe != 0.0) {
  eheu[8] <- eheu[8]+logAlphaFe
  eheu[10] <- eheu[10]+logAlphaFe
  eheu[12] <- eheu[12]+logAlphaFe
  eheu[14] <- eheu[14]+logAlphaFe
  eheu[16] <- eheu[16]+logAlphaFe
  eheu[18] <- eheu[18]+logAlphaFe
  eheu[20] <- eheu[20]+logAlphaFe
  eheu[22] <- eheu[22]+logAlphaFe
}

if (logCO > 0.0) {
  eheu[6] <- eheu[6]+logCO
} else {
  eheu[8] <- eheu[8]+abs(logCO)
}

#H and He do NOT get re-scaled with metallicity parameter:
logAZ <- list()
logAz <- logE10*(eheu[1:3]-12.0)
#..Everything else does:
logAz[3:length(logAz)] <- logE10*(eheu[2:length(eheu)]+log10ZScale-12.0)

expAz <- exp(logAz)
ATot <- sum(expAz)
logATot <- log10(ATot) #natual log

#END Initial guess for Sun section:
#Rescales kinetic temperature structure:
tauRos <- matrix(0,numDeps,2)
temp <- matrix(0,numDeps,2)
guessPGas <- matrix(0,numDeps,2)
guessPe <- matrix(0,numDeps,2)
guessNe <- matrix(0,numDeps,2)
kappaRos <- matrix(0,numDeps,2)
kappa500 <- matrix(0,numDeps,2)
pGas <- matrix(0,numDeps,2)
newPe <- matrix(0,numDeps,2)
pRad <- matrix(0,numDeps,2)
rho <- matrix(0,numDeps,2)
newNe <- matrix(0,numDeps,2)
mmw <- matrix(0,numDeps,2)

depths <- vector("numeric",numDeps)

if (specSynMode==TRUE) {

  #enusre self-consistency between parameters and model being read in
  teff <- restartParams["teffRS"]
  logg <- restartParams["loggRS"]
  log10ZScale <- restartParams["log10ZScaleRS"]
  logKapFudge <- restartParams["logKapFudgeRS"]
  logHeFe <- restartParams["logHeFeRS"]
  logCO <- restartParams["logCORS"]
  logAlphaFe <- restartParams["logAlphaFeRS"]

  tauRos[,1] <- restart[,"tauRosRSOne"]
  tauRos[,2] <- restart[,"tauRosRSTwo"]
  temp[,1] <- restart[,"tempRSOne"]
  temp[,2] <- restart[,"tempRSTwo"]
  pGas[,1] <- restart[,"pGasRSOne"]
  pGas[,2] <- restart[,"pGasRSTwo"]
  newPe[,1] <- restart[,"peRSOne"]
  newPe[,2] <- restart[,"peRSTwo"]
  #Everything in normal structure mode
  guessPGas[,1] <- restart[,"pGasRSOne"]
  guessPGas[,2] <- restart[,"pGasRSTwo"]
  guessPe[,1] <- restart[,"peRSOne"]
  guessPe[,2] <- restart[,"peRSTwo"]
  newNe[,2] <- as.numeric(newPe[,2])-as.numeric(temp[,2])-as.numeric(useful["logK"])
  newNe[,1] <- exp(newNe[,2])
  pRad[,1] <- restart[,"pRadRSOne"]
  pRad[,2] <- restart[,"pRadRSTwo"]
  rho[,1] <- restart[,"rhoRSOne"]
  rho[,2] <- restart[,"rhoRSTwo"]
  kappa500[,1] <- restart[,"kappa500RSOne"]
  kappa500[,2] <- restart[,"kappa500RSTwo"]
  kappaRos[,1] <- restart[,"kappaRhoRSOne"]
  kappaRos[,2] <- restart[,"kappaRhoRSTwo"]
  mmw <- restart[,"mmwRS"]

  #read in converged model - minimal processing
  nOuterIter <- 1
  nInnerIter <- 1

  } else {

    tauRos <- tauScale(numDeps, log10MinDepth, log10MaxDepth)
    if (teff < F0Vtemp) {
      if (logg > 3.5) {
        #cool Dwarf - rescale from Teff=5000 reference model
        temp <- phxRefTemp5000(teff, numDeps, tauRos)
      } else {
        #we're a cool giant - rescale from teff=4250, log(g)=2.0 model
        temp <- phxRefTempg20(teff, numDeps, tauRos)
      }
    } else if (teff >= F0Vtemp) {
      #we're a HOT star! - rescale from teff=10,000 reference model
      temp <- phxRefTempT10000(teff, numDeps, tauRos)
    }

    #scaled from Phoenix solar model
    if (teff < F0Vtemp) {
      if (logg > 3.5) {
        #we're a cool dwarf - rescale from teff=5000 reference model
        guessPGas <- phxRefPGasT5000(grav, zScale, logAz[,2], numDeps, tauRos)
        guessPe <- phxRefPeT5000(teff, grav, numDeps, tauRos, zScale, logAz[,2])
        guessNe <- phxRefNeT5000(numDeps, temp, guessPe)
      } else {
        #we're a cool giant - rescale from teff=4250, log(g)=2.0 model
        guessPGas <- phxRefPGasg20(grav, zScale, logAz[,2], numDeps, tauRos)
        guessPe <- phxRefPeg20(teff, grav, numDeps, tauRos, zScale, logAz[,2])
        guessNe <- phxRefNeg20(numDeps, temp, guessPe)
      }
    } else if (teff >= F0Vtemp) {
      #we're a HOT star - rescale from teff=10000 reference model
      guessPGas <- phxRefPGasT10000(grav, zScale, logAz[,2], numDeps, tauRos)
      guessPe <- phxRefPe64g20(teff, grav, numDeps, tauRos, zScale, logAz[,2])
      guessNe <- phxRefNeg20(numDeps, temp, guessPe)
  }
  }

#now do the same for the sun, for reference
tempSun <- phxSunTemp(teffSun, numDeps, tauRos)
pGasSunGuess <- phxSunPGas(gravSun, numDeps, tauRos)
NeSun <- phxSunNe(gravSun, numDeps, tauRos, tempSun, zScaleSun)
kappaSun <- phxSunKappa(numDeps, tauRos, zScaleSun)
mmwSun <- mmwFn(numDeps, tempSun, zScaleSun)
rhoSun <- massDensity(numDeps, tempSun, pGasSunGuess, mmwSun, zScaleSun)
pGasSun <- hydroFormalSoln(numDeps, gravSun, tauRos, kappaSun, tempSun, pGasSunGuess)
#stop
logNz <- getNz(numDeps, temp, guessPGas, guessPe, ATot, nelemAbnd, logAz)
logNH <- logNz[,1]
masterStagePops[2,1] <- guessPe[1]

#Load the total no. density of each element into the neutral stage slots of the
#masterStagePops array as a first guess at "species B" neutral populations for
#the molecular Saha eq. - is a guess at low temp where molecules form
masterStagePops[,1] <- logNz[numDeps,nelemAbnd]

if (teff < 6000) {
  print("Cool star mode")
} else {
  print("Hot star mode")
}

#Add subclass to each spectral class (lburns)
#Determine the spectralClass and subClass of main sequence stars, subdwarfs and
#white dwarfs based on the data in Appendix G of Introduction to Modern
#Astrophysics, 2nd Ed. by Carroll & Ostlie
if ((logg >= 4.0) && (logg < 5.0) || (logg >= 5.0) && (logg < 6.0) || (logg >= 5.0)) {
  if (teff < 3000.0) {
    spectralClass <- "L"
  } else if ((teff >= 3000.0) && (teff < 3900.0)) {
    spectralClass <- "M"
    if ((teff >= 3000.0) && (teff <= 3030.0)) {
      subClass <- "6"
    } else if ((teff > 3030.0) && (teff <= 3170.0)) {
      subClass <- "5"
    } else if ((teff > 3170.0) && (teff <= 3290.0)) {
      subClass <- "4"
    } else if ((teff > 3290.0) && (teff <= 3400.0)) {
      subClass <- "3"
    } else if ((teff > 3400.0) && (teff <= 3520.0)) {
      subClass <- "2"
    } else if ((teff > 3520.0) && (teff <= 3660.0)) {
      subClass <- "1"
    } else if ((teff > 3660.0) && (teff <= 3900.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 3900.0) && (teff < 5200.0)) {
    spectralClass <- "K"
    if ((teff >= 3900.0) && (teff <= 4150.0)) {
      subClass <- "7"
    } else if ((teff > 4150.0) && (teff <= 4410.0)) {
      subClass <- "5"
    } else if ((teff > 4410.0) && (teff <= 4540.0)) {
      subClass <- "4"
    } else if ((teff > 4540.0) && (teff <= 4690.0)) {
      subClass <- "3"
    } else if ((teff > 4690.0) && (teff <= 4990.0)) {
      subClass <- "1"
    } else if ((teff > 4990.0) && (teff <= 5200.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 5200.0) && (teff < 5950.0)) {
    spectralClass <- "G"
    if ((teff >= 5200.0) && (teff <= 5310.0)) {
      subClass <- "8"
    } else if ((teff > 5310.0) && (teff <= 5790.0)) {
      subClass <- "2"
    } else if ((teff > 5790.0) && (teff < 5950.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 5500.0) && (teff< 7500.0)) {
    spectralClass <- "F"
    if ((teff >= 5500.0) && (teff <= 6410.0)) {
      subClass <- "5"
    } else if ((teff > 6410.0) && (teff <= 7000.0)) {
      subClass <- "2"
    } else if ((teff > 7000.0) && (teff < 7500.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 7500.0) && (teff < 10300.0)) {
    spectralClass <- "A"
    if ((teff >= 7500.0) && (teff < 7830.0)) {
      subClass <- "8"
    } else if ((teff > 7830.0) && (teff <= 8550.0)) {
      subClass <- "5"
    } else if ((teff > 8550.0) && (teff <= 9460.0)) {
      subClass <- "2"
    } else if ((teff > 9460.0) && (teff <= 9820.0)) {
      subClass <- "1"
    } else if ((teff > 9280.0) && (teff < 10300.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 10300.0) && (teff < 29300.0)) {
    spectralClass <- "B"
    if ((teff >= 10300.0) && (teff <= 10900.0)) {
      subClass <- "9"
    } else if ((teff > 10900.0) && (teff <= 11700.0)) {
      subClass <- "8"
    } else if ((teff > 11700.0) && (teff <= 12700.0)) {
      subClass <- "7"
    } else if ((teff > 12700.0) && (teff <= 13800.0)) {
      subClass <- "6"
    } else if ((teff > 13800.0) && (teff <= 15100.0)) {
      subClass <- "5"
    } else if ((teff > 15100.0) && (teff <= 18300.0)) {
      subClass <- "3"
    } else if ((teff > 18300.0) && (teff <= 20200.0)) {
      subClass <- "2"
    } else if ((teff > 20200.0) && (teff <= 24500.0)) {
      subClass <- "1"
    } else if ((teff > 24500.0) && (teff < 29300.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 29300.0) && (teff < 40000.0)) {
    spectralClass <- "O"
    if ((teff >= 29300.0) && (teff <= 35000.0)) {
      subClass <- "8"
    } else if ((teff > 35000.0) && (teff <= 36500.0)) {
      subClass <- "7"
    } else if ((teff > 36500.0) && (teff <= 37800.0)) {
      subClass <- "6"
    } else if ((teff > 37800.0) && (teff < 40000.0)) {
      subClass <- "5"
    }
  }
}

#Determine the spectralClass of supergiants and bright giants.
if (((logg >= 0.0) && (logg < 1.0)) || ((logg >= 1.0) && (logg < 1.5))) {
  if (teff < 2700.0) {
    spectralClass <- "L"
  } else if ((teff >= 2700.0) && (teff < 3650.0)) {
      spectralClass <- "M"
      if ((teff >= 2700.0) && (teff <= 2710.0)) {
        subClass <- "6"
      } else if ((teff > 2710.0) && (teff <= 2880.0)) {
        subClass <- "5"
      } else if ((teff > 2880.0) && (teff <= 3060.0)) {
        subClass <- "4"
      } else if ((teff > 3060.0) && (teff <= 3210.0)) {
        subClass <- "3"
      } else if ((teff > 3210.0) && (teff <= 3370.0)) {
        subClass <- "2"
      } else if ((teff > 3370.0) && (teff <= 3490.0)) {
        subClass <- "1"
      } else if ((teff > 3490.0) && (teff <= 3650.0)) {
        subClass <- "0"
      }
  } else if ((teff >= 3650.0) && (teff < 4600.0)) {
      spectralClass <- "K"
      if ((teff >= 3650.0) && (teff <= 3830.0)) {
        subClass <- "7"
      } else if ((teff > 3830.0) && (teff <= 3990.0)) {
        subClass <- "5"
      } else if ((teff > 3990.0) && (teff <= 4090.0)) {
        subClass <- "4"
      } else if ((teff > 4090.0) && (teff <= 4190.0)) {
        subClass <- "3"
      } else if ((teff > 4190.0) && (teff <= 4430.0)) {
        subClass <- "1"
      } else if ((teff > 4430.0) && (teff < 4600.0)) {
        subClass <- "0"
      }
  } else if ((teff >= 4600.0) && (teff < 5500.0)) {
      spectralClass <- "G"
      if ((teff >= 4600.0) && (teff <= 4700.0)) {
        subClass <- "8"
      } else if ((teff > 4700.0) && (teff <= 5198.0)) {
        subClass <- "2"
      } else if ((teff > 5198.0) && (teff < 5500.0)) {
        subClass <- "0"
      }
  } else if ((teff >= 5500.0) && (teff < 7500.0)) {
      spectralClass <- "F"
      if ((teff >= 5500.0) && (teff <= 5750.0)) {
        subClass <- "8"
      } else if ((teff > 5750.0) && (teff <= 6370.0)) {
        subClass <- "5"
      } else if ((teff > 6370.0) && (teff <= 7030.0)) {
        subClass <- "2"
      } else if ((teff > 7030.0) && (teff < 7500.0)) {
        subClass <- "0"
      }
  } else if ((teff >= 7500.0) && (teff < 10000.0)) {
    spectralClass <- "A"
    if ((teff >= 7500.0) && (teff <= 7910.0)) {
      subClass <- "8"
    } else if ((teff > 7910.0) && (teff <= 8610.0)) {
      subClass <- "5"
    } else if ((teff > 8610.0) && (teff <= 9380.0)) {
      subClass <- "2"
    } else if ((teff > 9380.0) && (teff < 10000.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 10000.0) && (teff < 27000.0)) {
    spectralClass <- "B"
    if ((teff >= 10000.0) && (teff <= 10500.0)) {
      subClass <- "9"
    } else if ((teff > 10500.0) && (teff <= 11100.0)) {
      subClass <- "8"
    } else if ((teff > 11100.0) && (teff <= 11800.0)) {
      subClass <- "7"
    } else if ((teff > 11800.0) && (teff <= 12600.0)) {
      subClass <- "6"
    } else if ((teff > 12600.0) && (teff <= 13600.0)) {
      subClass <- "5"
    } else if ((teff > 13600.0) && (teff <= 16000.0)) {
      subClass <- "3"
    } else if ((teff > 16000.0) && (teff <= 17600.0)) {
      subClass <- "2"
    } else if ((teff > 17600.0) && (teff <= 21400.0)) {
      subClass <- "1"
    } else if ((teff > 21400.0) && (teff < 27000.0)) {
      subClass <- "0"
    }
  } else if ((teff >= 27000.0) && (teff < 42000.0)) {
    spectralClass <- "O"
    if ((teff >= 27000.0) && (teff <= 34000.0)) {
      subClass <- "8"
    } else if ((teff > 34000.0) && (teff <= 36200.0)) {
      subClass <- "7"
    } else if ((teff > 36200.0) && (teff <= 38500.0)) {
      subClass <- "6"
    } else if ((teff > 38500.0) && (teff < 42000.0)) {
      subClass <- "5"
    }
  }
}

#Determine luminClass based on logg
if ((logg >= 0.0) && (logg < 1.0)) {
  luminClass <- "I"
} else if ((logg >= 1.0) && (logg < 1.5)) {
  luminClass <- "II"
} else if ((logg >= 1.5) && (logg < 3.0)) {
  luminClass <- "III"
} else if ((logg >= 3.0) && (logg < 4.0)) {
  luminClass <- "IV"
} else if ((logg >= 4.0) && (logg < 5.0)) {
  luminClass <- "V"
} else if ((logg >= 5.0) && (logg < 6.0)) {
  luminClass <- "VI"
} else if (logg >= 5.0) {
  luminClass <- "WD"
}

spectralType <- glue(spectralClass, subClass, " ", luminClass)
print(glue("Spectral type: ", spectralType))

log10e <- log10(exp(1))

#END initial guess for Sun section
#Converge atmospheric structure
# -Includes *initial* ionization equilibrium *without* molecules (for now)

species <- " "
tauOneStagePops <- matrix(0.0, numStages, nelemAbnd)
unity <- 1.0
zScaleList <- 1.0 #initialization

numAtmPrtTmps <- 5
numMolPrtTmps <- 5
log10UwAArr <- matrix(0.0, numAtmPrtTmps, numStages)

#Ground state ionization E - Stage I (eV)
chiIArr <- rep(999999.0, numStages)

#For diatomic molecules
speciesAB <- " "
speciesA <- " "
speciesB <- " "

#For mass density (rho) calculations
masterMolPops <- matrix(-49.0, numDeps, nMols)

#Interpolate atomic partition fnd tabulated at two temperatures
thisUwAV <- vector("numeric", numAtmPrtTmps)
thisUwBV <- vector("numeric", numAtmPrtTmps)

#Interpolate in molecular parition fns tabulate at five temperatures
thisQwAB <- vector("numeric", numMolPrtTmps)

logNums <- matrix(0.0, numDeps, numStages)

#For diatomic molecules
logNumA <- vector("numeric", numDeps)
logNumB <- vector("numeric", numDeps)
logNumFracAB <- vector("numeric", numDeps)

Ng <- vector("numeric", numDeps)

logKappa <- matrix(0.0, numDeps, numLams)
logKappaHHe <- matrix(0.0, numDeps, numLams)
logKappaMetalBF <- matrix(0.0, numDeps, numLams)
logKappaRayl <- matrix(0.0, numDeps, numLams)

newTemp <- matrix(0.0, numDeps, 2)

#Variables for ionization/molecular equilibrium treatment
#for diatomic molecules
logNumBArr <- matrix(0.0, numDeps, numAssocMols)

#Interpolate in atomic partition fns tabulate at five temperatures
log10UwBArr <- matrix(0.0, numAtmPrtTmps, numAssocMols)

dissEArr <- vector("numeric", numAssocMols)

logQwABArr <- matrix(0.0, numMolPrtTmps, numAssocMols)
logMuABArr <- vector("numeric", numAssocMols)

#Arrays of pointers into master molecule and element lists
mname_ptr <- vector("numeric", numAssocMols)
specB_ptr <- vector("numeric", numAssocMols)
specA_ptr <- 0.0
specB2_ptr <- 0.0
mnameBtemplate <- " "

#We converge the PGas - Pe relation first under the assumption that all free e^-s
#are from single ionizations a la David Gray Ch. 9. This approach separates
#converging ionization fractions and Ne for spectrum synthesis purposes from
#converging the Pgas-Pe-N_H-N_He relation for computing the mean opacity for HSE
thisTemp <- vector("numeric", 2)
log10UwUArr <- log10UwLArr <- vector("numeric", numAtmPrtTmps)
log300 <- log(300.0)
log2 <- log(2.0)

#Pgas-kappa iteration
for (pIter in 1:nOuterIter) {
  for (neIter in 1:nInnerIter) {
    for (iD in 1:numDeps) {
      #reinitialize accumulators
      thisTemp[1] <- temp[iD, 1]
      thisTemp[2] <- temp[iD, 2]
      peNumerator <- 0.0
      peDenominator <- 0.0
      for (iElem in 1:nelemAbnd) {
        species <- glue(cname[iElem],"I")
        chiI <- getIonE(species)
        #the following is a 2-element vector of temperature-dependent partition
        #fns, U that are base e log_e U
        log10UwLArr <- getPartFn2(species)
        species <- glue(cname[iElem],"II")
        log10UwUArr <- getPartFn2(species)
        logPhi <- sahaRHS(chiI, log10UwUArr, log10UwLArr, thisTemp)
        logPhiOverPe <- logPhi-guessPe[iD, 1]
        logOnePlusPhiOverPe <- log(1.0+exp(logPhiOverPe))
        logPeNumerTerm <- logAz[iElem]+logPhiOverPe-logOnePlusPhiOverPe
        peNumerator <- peNumerator+exp(logPeNumerTerm)
        logPeDenomTerm <- logAz[iElem]+log(1.0+exp(logPeNumerTerm))
        peDenominator <- peDenominator+exp(logPeDenomTerm)
      }
      #iElem chemical element loop
      newPe[iD,2] <- guessPGas[iD, 1]+log(peNumerator)-log(peDenominator)
      newPe[,1] <- exp(newPe[,2])
      guessPe[iD,2] <- newPe[iD,2]
      guessPe[iD,1] <- exp(guessPe[iD,2])
    }
  }
  newNe[2] <- newPe[,2]-temp[,2]-useful[,"logK"]
  newNe[1] <- exp(newNe[,2])
  guessNe[,2] <- newNe[,2]
  guessNe[,1] <- newNe[,1]

  #refine the number of densities of the chemical elements at all depths
  logNz <- getNz(numDeps, temp, guessPGas, guessPe, ATot, nelemAbnd, logAz)
  logNH <- logNz[1]

  #compute ionization fractions of H & He for kappa calculation
  zScaleList <- 1.0

  #these 2-element temperature-dependent parition fns are logarithmic
  logNumBArr <- matrix(-49.0, numDeps, numAssocMols)
  log10UwBArr <- matrix(0.0, numAtmPrtTmps, numAssocMols)
  dissEArr <- rep(29.0, numAssocMols)
  logQwABArr <- matrix(log300, numMolPrtTmps, numAssocMols)
  logMuABArr <- rep(log2+useful[,"logAmu"], numAssocMols)
  mname_ptr <- specB_ptr <- vector("numeric", numAssocMols)
  defaultQwAB <- log(300.0)
  #default that applied to most cases - neutral stage (I) forms molecules
  specBStage <- 0
  nmrtrDissE <- 15.0 #probability high by default
  nmrtrLog10UwB <- rep(0.0, numAtmPrtTmps)
  nmrtrLog10UwA <- 0.0
  nmrtrLogQwAB <- rep(log300, numMolPrtTmps)
  nmrtrLogMuAB <- useful[,"logAmu"]
  nmrtrLogNumB <- rep(0.0, numDeps)
  logGroundRatio <- rep(0.0, numDeps)

  for (iElem in 1:26) {
    species <- glue(cname[iElem], "I")
    chiIArr[1] <- getIonE(species)
    #the following is a 2-element vector of temperature-dependent partition fns, U, that are base e log_e U
    log10UwAArr[1] <- getPartFn2(species)
    species <- glue(cname[iElem], "II")
    chiIArr[2] <- getIonE(species)
    log10UwAArr[2] <- getPartFn2(species)
    species <- glue(cname[iElem],"III")
    chiIArr[3] <- getIonE(species)
    log10UwAArr[3] <- getPartFn2(species)
    species <- glue(cname[iElem],"IV")
    chiIArr[4] <- getIonE(species)
    log10UwAArr[4] <- getPartFn2(species)
    species <- glue(cname[iElem],"V")
    chiIArr[5] <- getIonE(species)
    log10UwAArr[5] <- getPartFn2(species)
    species <- glue(cname[iElem],"VI")
    chiIArr[6] <- getIonE(species)
    log10UwAArr[6] <- getPartFn2(species)

    thisNumMols <- 0.0

    ###Python script
    #try:
    # thisNumMols <- cnameMols[iElem][:].index("None")
    #except:
    # thisNumMols <- numAssocMols

    #trying to account for mols in ionization eq desatbilizes everything
    thisNumMols <- 0.0
    if (thisNumMols > 0) {
      #find pointer to molecule in master mname list for each associated molecule
      for (iMol in 1:thisNumMols) {
        for (jj in 1:nMols) {
          if (cnameMols[iMol,iElem] == mname[jj]) {
            mname_ptr[iMol] == jj
            break
          }
        }
      }

      #iMol loop in master mnames list, iMol loop in associated molecules
      #now find pointer to atomic species B in master cname list for each associated molecule found in
      #the master name list
      for (iMol in 1:thisNumMols) {
        for (jj in 1:nelemAbnd) {
          if (mnameB[mname_ptr[iMil]] == cname[jj]) {
            specB_ptr[iMol] = jj
            break
          }
        }
      }

      #load arrays with molecular species AB and atomic species B data for method stagePops2()
      for (iMol in 1:thisNumMols) {
        if (mnameB[mname_ptr[iMol]] == "H2+") {
          specBStage <- 1
        } else {
          specBStage <- 0
        }

        logNumBArr[iMol] <- masterStagePops[specBStage,specB_ptr[iMol]]
        dissEArr[iMol] <- getDissE(mname[mname_ptr[iMol]])
        species <- glue(cname[specB_ptr[iMol]], "I") #neutral stage
        log10UwBArr[iMol] <- getPartFn2(species)
        logQwABArr[iMol] <- getMolPartFn(mname[mname_ptr[iMol]])
        #compute the reduced mass, muAB, in g
        massA <- getMass(cname[i])
        massB <- getMass(cname[specB_ptr[iMol]])
        logMuABArr[iMol] <- log(massA)+log(massB)-log(massA+massB)+useful[,"logAmu"]
      }
    }

    logNums <- stagePop2(logNz[iElem], guessNe, chiIArr, log10UwAArr, thisNumMols, logNumBArr,
                         dissEArr, log10UwBArr, logQwABArr, logMuABArr, numDeps, temp)

    masterStagePops[iElem] <- logNums

    #get mass density from chemical composition
    rho <- massDensity2(numDeps, nelemAbnd, logNz, cname)

    #total number density of gas particles: nuclear species + free electrons and compute mean molecular
    #weight, mmw ("mu")
    Ng <- newNe[,1]

    for (i in 1:numDeps) {
      for (j in 1:nelemAbnd) {
        Ng[i] <- Ng[1]+exp(logNz[i,j]) #initialize accumulation
      }
    }

    mmw <- exp(rho[i,1]-log(Ng[numDeps]))

    #H and He only for now, only compute H, He, and e^- opacity sources
    logKappaHHe <- kappas2(numDeps, newPe, zScale, temp, rho, numLams,
                           lambdaScale, logAz[2], masterStagePops[1,1],
                           masterStagePops[2,1], masterStagePops[1,2],
                           masterStagePops[2,2], newNe, teff, logTotalFudge)

    #add in metal b-f opacity from adapted Moog routines
    logKappaMetalBF <- masterMetal(numDeps, numLams, temp, lambdaScale,
                                   masterStagePops)

    #add in Rayleigh scattering opacity from adapted Moog routines
    logKappaRayl <- masterRayl(numDeps, numLams, temp, lambdaScale, masterStagePops,
                               masterMolPops)

    logKappa <- vector()

    for (iD in 1:numDeps) {
      for (iL in 1:numLams) {
        logKappa <- append(logKappa, exp(logKappaHHe[iD,iL])+
                             exp(logKappaMetalBF[iD,iL])-
                             rho[iD,2]+
                             exp(logKappaRayl[iD,iL]))
      }
    }

    for (iD in 1:numDeps) {
      for (iL in 1:numLams) {
        logKappa <- log(logKappa[iD,iL])
      }
    }

    it500 <- lamPoint(numLams, lambdaScale, 500E-7)
    kappa500[,2] <- logKappa[it500]
    kappa500[,1] <- exp(logKappa[it500])

    pGas <- hydroFormalSoln(numDeps, grav, tauRos, kappaRos, temp, guessPGas)
    pRad <- radPress(numDeps, temp)

    guessPGas[,1] <- pGas[,1]
    guessPGas[,2] <- pGas[,2]
    log10pgas <- log10e*pGas[,1]
    log10pe <- log10e*(newNe[numDeps,2]+useful[,"logK"]+temp[numDeps,2])
    pe <- newNe[numDeps,2]+useful[,"logK"]+temp[numDeps,2]
    log10prad <- log10e*pRad[,1]
    log10ne <- log10e*newNe[,1]

    thisClr <- pallette[pIter/numClrs]
  }
}

depths <- depthScale(numDeps, tauRos, kappaRos, rho)

ifTcorr <- FALSE
ifConvec <- FALSE
numTCorr <- 0.0

for (i in 1:numTCorr) {
  newTemp <- 1
  temp[,2] <- newTemp[,2]
  temp[,1] <- newTemp[,1]
}

###################################################
#Reconverge Ionization/chemical equilibrium WITH molecules
#
#Now that the atmospheric structure is settled separately converge the
#Ne-ionization-fractions-molecular equilibrium for all elements and populate
#the ionization stages of all the species for spectrum synthesis
###################################################

iTauOne <- tauPoint(numDeps, tauRos, unity)
zScaleList <- 1.0

logNumBArr <- matrix(-49.0, numDeps, numAssocMols)
log10UwBArr <- matrix(0.0, numAtmPrtTmps, numAssocMols)
dissEArr <- rep(29.0, numAssocMols)
logQwABArr <- matrix(log300, numMolPrtTmps, numAssocMols)
logMuABArr <- rep(log2+useful[,"logAmu"], numAssocMols)
mname_ptr <- vector("numeric", numAssocMols)
specB_ptr <- vector("numeric", numAssocMols)

defaultQwAB <- log300

#neutral stage (I) forms molecules
specBStage <- 0.0

#for element A of main molecule being treated in molecular equilibrium
nmrtrDissE <- 15.0
nmrtrLog10UwB <- vector("numeric", numAtmPrtTmps)
nmrtrLog10UwA <- 0.0
nmrtrLogQwAB <- rep(log300, numMolPrtTmps)
nmrtrLogMuAB <- useful[,"logAmu"]
nmrtrLogNumB <- vector("numeric", numDeps)

logGroundRatio <- vector("numeric", numDeps)

#Iterate the electron densities, ionization fractions, and molecular densities
for (nIter2 in 1:nInnerIter) {
  for (iElem in 1:nelemAbnd) {
    species <- glue(cname[iElem], "I")
    chiIArr[1] <- getIonE(species)
    log10UwAArr[1] <- getPartFn2(species)

    species <- glue(cname[iElem], "II")
    chiIArr[2] <- getIonE(species)
    log10UwAArr[2] <- getPartFn2(species)

    species <- glue(cname[iElem], "III")
    chiIArr[3] <- getIonE(species)
    log10UwAArr[3] <- getPartFn2(species)

    species <- glue(cname[iElem], "IV")
    chiIArr[4] <- getIonE(species)
    log10UwAArr[4]= getPartFn2(species)

    species <- glue(cname[iElem], "V")
    chiIArr[5] <- getIonE(species)
    log10UwAArr[5]= getPartFn2(species)

    species <- glue(cname[iElem], "VI")
    chiIArr[6] <- getIonE(species)
    log10UwAArr[6]= getPartFn2(species)

    thisNumMols <- 0.0

    for (iMol in 1:numAssocMols) {
      if (cnameMols[iMol, iElem] == "None") {
        break
      }
      thisNumMols <- thisNumMols + 1
    }

    if (thisNumMols > 0) {
      #find pointer to molecule in master mname list for each
      #association molecules
      for (iMol in 1:thisNumMols) {
        for (jj in 1:nMols) {
          if (cnameMols[iMol, iElem] == mname[jj]) {
            mname_ptr[iMol] <- jj
            break
          }
        }
      }
      #now find pointer to atomic species B in master cname list for
      #each associated molecule found in master mname list
      for (iMol in 1:thisNumMols) {
        for (jj in 1:nelemAbnd) {
          if (mname[mname_ptr[iMol]] == cname[jj]) {
            specB_ptr[iMol] <- jj
            break
          }
        }
      }
      #now load arrays with molecular species AB and atomic species B
      #data
      for (iMol in 1:thisNumMols) {
        if (mnameB[mname_ptr[iMol]] == "H2+") {
          specBStage <- 1
        } else {
          specBStage <- 0
        }

        logNumBArr[iMol] <- masterStagePops[specB_ptr[specBStage, iMol]]
        dissEarr[iMol] <- getDissE(mname[mname_ptr[iMol]])
        species <- glue(cname[specB_ptr[iMol]]+"I")
        log10UwBArr[iMol] <- getPartFn2(species)
        logQwABArr[iMol] <- getMolPartFn(mname_ptr[iMol])
        massA <- getMass(cname[iElem])
        massB <- getMass(cname[specB_ptr[iMol]])
        logMuABArr[iMol] <- log(massA)+log(massB)-log(massA+massB)+useful[,"logAmu"]
      }
    }

    logNums <- stagePop2(logNz[iElem], guessNe, chiIArr, log10UwAArr,
                         thisNumMols, logNumBArr, dissEArr, log10UwBArr,
                         logQwABArr, logMuABArr, numDeps, temp)

    masterStagePops[iElem] <- logNums[numDeps, numStages]
    tauOneStagePops[iElem] <- logNums[iTauOne, iStage]
  }
  #compute all molecular popluations
  log10UwA <- vector("numeric", numAtmPrtTmps)

  for (iMol in 1:nMols) {
    #find elements A and B in master atomic element list
    specA_ptr <- 0.0
    specB2_ptr <- 0.0

    for (jj in 1:nelemAbnd) {
      if (mnameA[iMol] == cname[jj]) {
        specA_ptr = jj
        break
      }
    }
    #get its partition fn
    species <- glue(cname[specA_ptr] + "I")
    log10UwA <- getPartFn2(species)

    for (jj in 1:nelemAbnd) {
      if (mnameB[iMol] == cname[jj]) {
        specB2_ptr <- jj
        break
      }
    }
    #we need all molecules species A participates in - including the
    #current molecule itself at this point, it's just like setting up
    #the ionization equilibrium to account for molecules as above
    thisNumMols <- 0.0
    for (im in 1:numAssocMols) {
      if (cnameMols[im, specA_ptr] == "None") {
        break
      }
      thisNumMols <- thisNumMols + 1
    }
    if (thisNumMols > 0) {
      #find pointer to molecule in master name list for each associated molecule
      for (im in 1:thisNumMols) {
        for (jj in 1:nMols) {
          if (cnameMols[im, specA_ptr] == mname[jj]) {
            mname_ptr[im] <- jj
            break
          }
        }
      }
      #now find pointer to atomic species B in master cname list for
      #each associated molecule found in master mname list
      for (im in 1:thisNumMols) {
        mnameBtemplate <- " "
        #species "B" is whichever element is NOT species "A" in master
        #molecule
        if (mnameB[mname_ptr[im]] == mnameA[iMol]) {
          #get the other atom
          mnameBtemplate <- mnameA[mname_ptr[im]]
        } else {
          mnameBtemplate <- mnameB[mname_ptr[im]]
        }

        for (jj in 1:nelemAbnd) {
          if (mnameBtemplate == cname[jj]) {
            specB_ptr[im] <- jj
            break
          }
        }
      }
      #now load arrays with molecular specis AB and atomic species B
      #data for method molPops()
      for (im in 1:thisNumMols) {
        if (mname[mname_ptr[im]] == "H2+") {
          specBStage <- 1
        } else {
          specBStage <- 0
        }

        logNumBArr[im] <- masterStagePops[specBStage, specB_ptr[im]]
        dissEArr[im] <- getDissE(mname_ptr[im])
        species <- glue(cname[specB_ptr[im]]+"I")
        log10UwBArr[im] <- getPartFn2(species)
        logQwABArr[im] <- getMolPartFn(mname[mname_ptr[im]])
        #compute the reduced mass, muAB, in g
        massA <- getMass(cname[specA_ptr])
        massB <- getMass(cname[specB_ptr[im]])
        logMuABArr[im] <- log(massA)+log(massB)-log(massA+massB)+useful[,"logAmu"]
        #one of the species A-associated molecules will be the actual
        #molecule, AB, for which we want the population - pick this
        #out for the numerator in the master reaction
        if (mname[mname_ptr[im]] == mname[iMol]) {
          nmrtrDissE <- dissEArr[im]
          nmrtrLog10UwB <- log10UwBArr[im]
          nmrtrLogQwAB <- logQwABArr[im]
          nmrtrLogMuAB <- logMuABArr[im]
          nmrtrLogNumB <- logNumBArr[im]
        }
      }
    }
    #compute total population of particle in atomic ionic stages over
    #number in ground ionization stage for master denominator so we
    #don't have to re-compute
    for (iTau in 1:numDeps) {
      totalIonic <- 0.0
      for (iStage in 1:numStages) {
        totalIonic <- totalIonic + exp(masterStagePops[iTau, iStage, specA_ptr])
      }
      logGroundRatio[iTau] <- log(totalIonic)-masterStagePops[iTau,1,specA_ptr]
    }

    logNumFracAB <- molPops(nmrtrLogNumB, mnrtrDissE, log10UwA, nmrtrLog10UwB,
                            nmrtrLogQwAB, nmrtrLogMuAB, thisNumMols, logNumBArr,
                            dissEArr, log10UwBArr, logQwABArr, logMuABArr,
                            logGroundRatio, numDeps, temp)

    masterMolPops[iMol] <- logNz[numDeps,specA_ptr]+logNumFracAB[numDeps]
  }
  #compute update Ne & Pe initialize accumulation of electrons at all
  #depths
  newNe[,1] <- vector("numeric", numDeps)

  #this is cumulative and not trivial
  for (iTau in 1:numDeps) {
    for (iElem in 1:nelemAbnd) {
      newNe[iTau,1] <- newNe[iTau,1]+exp(masterStagePops[iTau,2,iElem])+
        2.0*exp(masterStagePops[iTau,3,iElem])
    }
  }

  newNe[,2] <- log(newNe[,1])
  guessNe[,1] <- newNe[,1]
  guessNe[,2] <- newNe[,2]

  #diagnostic plots
  #graphically inspect convergence - issue `matplotlib qt5` before running
  thisClr <- palette[neIter/numClrs]

  log10pe <- lapply(numDeps, function(x) log10e*(newNe[x,2]+useful[,"logK"]
                                                 +temp[x,2]))
  log10ne[i] <- log10e*newNe[,2]
}

log10pe <- lapply(numDeps, function(x) log10e*(newNe[x,2]+useful[,"logK"]
                                               +temp[x,2]))
log10e <- log10e*newNe[,2]

#create a restart module for use in spectrum synthesis mode
outfile <- glue(outPath,fileStem,"_restart.py")

#print vertical atmospheric structure as a R Data Structure for
#re-import to chromaStar, can be used as a convergence model for an
#almost pure spectrum synthesis run with open(outFile, "w",
#encoding='utf=8') as strucHandle

write.table(glue("# ", inputParamString),file="outFile", append=FALSE,
	      eol="\n")

write.table(glue("teffRS = ", teff, " # K"), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("loggRS = ", logg, " #log (cm/^2)"), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("log10ZScaleRS = ", log10ZScale), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("xiTRS = ", xiT, " # (km/s)"), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("logKapFudgeRS = ", logKapFudge), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("logHeFeRS = ", logHeFe), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("logCORS = ", logCO), file="outFile", append=TRUE,
	      eol="\n")

write.table(glue("logAlphaFeRS = ", logAlphaFe), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("numDeps = ", numDepStr), file="outFile",
	      append=TRUE, eol="\n")

write.table(glue("tauRosRS = matrix(0.0,", numDepsStr, ", 2)"),
	      file="outFile", append=TRUE, eol="\n")

handler <- "tauRosRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(tauRos[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(tauRos[numDeps-1,1], "]"), file="outFile", append=TRUE, eol="\n")

handler <- "tempRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(tauRos[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(tauRos[numDeps-1,2], "]"), file="outFile",
				      append=TRUE, eol="\n")

write.table(glue("tempRS = matrix(0.0,", numDepsStr,", 2)"),
			file="outFile", append=TRUE, eol="\n")

handler <- "tempRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(temp[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(temp[numDeps-1,1], "]"), file="outFile",
				    append=TRUE, eol="\n")

handler <- "tempRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(temp[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(temp[numDeps-1,2], "]"), file="outFile",
				    append=TRUE, eol="\n")

write.table(glue("pGasRS = matrix(0.0,",numDepsStr,", 2)"),
			 file="outFile", append=TRUE, eol="\n")

handler <- "pGasRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(pGas[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(pGas[numDeps-1,1], "]"), file="outFile",
				    append=TRUE, eol="\n")

handler <- "pGasRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(pGas[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(pGas[numDeps-1,2], "]"), file="outFile",
				    append=TRUE, eol="\n")

write.table(glue("peRS[1] = matrix(0.0,",numDepsStr,", 2)"),
			  file="outFile", append=TRUE, eol="\n")

handler <- "peRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(newPe[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(newPe[numDeps-1, 1], "]"), file="outFile",
				  append=TRUE, eol="\n")

handler <- "peRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(newPe[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(newPe[numDeps-1, 2], "]"), file="outFile",
				  append=TRUE, eol="\n")

write.table(glue("pRadRS = matrix(0.0,", numDepsStr, ", 2)"),
			 file="outFile", append=TRUE, eol="\n")

handler <- "pRadRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(pRad[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(pRad[numDeps-1,1], "]"), file="outFile",
				    append=TRUE, eol="\n")

handler <- "pRadRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
    handler <- append(handler, glue(pRad[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(pRad[numDeps-1, 2], "]"), file="outFile",
				 append=TRUE, eol="\n")

write.table(glue("rhoRS = matrix(0.0,", numDepsStr, ", 2)"),
            file="outFile", append=TRUE, eol="\n")

handler <- "rhoRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
  handler <- append(handler, glue(rho[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(rho[numDeps-1,1], "]"), file="outFile",
            append=TRUE, eol="\n")

handler <- "rho[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
  handler <- append(handler, glue(rho[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(rho[numDeps-1, 2], "]"), file="outFile",
            append=TRUE, eol="\n")

write.table(glue("kappaRosRS = matrix(0.0,", numDepsStr, ", 2)"),
            file="outFile", append=TRUE, eol="\n")

handler <- "kappaRosRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
  handler <- append(handler, glue(kappaRos[i,1], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(kappaRos[numDeps-1,1], "]"), file="outFile",
            append=TRUE, eol="\n")

handler <- "kappaRosRS[2] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
  handler <- append(handler, glue(kappaRos[i,2], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(kappaRos[numDeps-1, 2], "]"), file="outFile",
            append=TRUE, eol="\n")

write.table(glue("mmwRS = matrix(0.0,", numDepsStr, ", 2)"),
            file="outFile", append=TRUE, eol="\n")

handler <- "mmwRS[1] = [" #continuation line '\' has to be escaped

for (i in 1:numDeps-1) {
  handler <- append(handler, glue(mmw[i], ", "))
}

write.table(handler, file="outFile", append=TRUE, eol="\n")

write.table(glue(mmw[numDeps-1], "]"), file="outFile",
            append=TRUE, eol="\n")

######################################################
#
# Surface radiation field
#   - flux distribution (SED)
#   - high resolution synthetic spectrum
#
######################################################

#Emergent radiation modelling
#  set up theta grid
#  cosTheta is a 2xnumThetas array
#  row 1 is Gaussian quadrature weights
#  row 2 is cos(theta) values
#  Gaussian quadrature: number of angles, numThetas, will have to be
#    determined after the fact
cosTheta <- thetas()
numThetas <- length(cosTheta[,1])

#Establish a phi grid for non-axi-symmetric situations (eg. spots,
#in situ rotation, ...)
#Number of phi values per quadrant of unit circle centered on
#stellar point in plane of sky: for geometry calculation: phi = 0 is
#direction of positive x-axis of right-handed 2D Cartesian coord sytem
#in plane of sky with origin at sub-stellar point (phi increases CCW)
numPhiPerQuad <- 6
numPhi <- numPhiD <- 4*numPhiPerQuad
phi <- vector("numeric", numPhi)

#Compute phi values in whole range (0 - 2pi radians)
delPhi <- 2.0*pi/numPhiD
phi <- delPhi*numPhi

#Spectrum synthesis section
isCool <- 7300.0

#Setup opacity lambda break-points and gray levels
minLambda <- 30.0
maxLambda <- 1E6

#JOLA molecular bands: overlapping line approcimation treats molecular
#ro-vibrational bands as pseudo-continuum opacity source by "smearing"
#out the individual rotational fine-structure lines
#1982 A&A ..113..173Z, Zeidler & Koester, 1982
jolaB <- jolaLambda <-  vector("numeric", 2)

#Alpha_P - weight of P branch (Delta J = -1)
#Alpha_R - weight of R branch (Delta J = 1)
#Alpha_Q - weight of Q branch (Delta J = 0)
jolaAlphaP <- jolaAlphaR <- jolaAlphaQ <- 0.0

#Allen's Astrophysical quantities, 4th Ed., 4.12.2-4.13.1:
#Electronics transition moment, Re
#"Line Strength", S = |R_e|^2*q_v'v" or just |R_e|^2
#Section 4.4.2 - for atoms or molecules: then gf = (8pi^2m_e*nu/3he^2)*S
#
#^48Ti^160 systems: Table 4.18, p. 91
#
#Rotational & vibrational constant for TiO states:, p. 87, Table 4.17
#
#General TiO molecular rotational & vibrational constants - Table 3.12, p. 47

#Zeidler & Koester 1982 p. 175, Sect vi:
#If Q branch (deltaLambda = +/- -1): alpP = alpR = 0.25, alpQ = 0.5
#If NO Q branch (deltaLambda = 0): alpP = alpR = 0.5, alpQ = 0.0

#Number of wavelength point sampling a JOLA band
jolaNumPoints <- 30

#Branch weights for transitions of DeltaLambda = +/- 1
jolaAlphP_DL1 <- jolaAlphR_DLI <- 0.25
jolaAlphQ_DL1 <- 0.5

#Branch weights for transitions of DeltaLambda = 0
jolaAlphP_DL0 <- jolaAlphR_DL0 <- 0.5
jolaAlphQ_DL0 <- 0.0

logSTofHelp <- log(8.0/3.0)+2.0*log(pi)+useful[,"logMe"]-useful[,"logH"]-2.0*useful[,"logEe"]
jolaQuantumS <- 1.0 #default for multiplicative factor
logNumJola <- vector("numeric", numDeps)
jolaProfPR <- matrix(0.0, numDeps, jolaNumPoints) #for unified P&R branch
jolaProfQ <- matrix(0.0, numDeps, jolaNumPoints) #for Q branch
#Differential cross-section - the main "product" of the JOLA approximation
dfBydv <- matrix(0.0, numDeps, jolaNumPoints)

dataPath <- "input_data/"

#Atomic line list
#Kramida, A., Ralchenko, Yu., Reader, J., and NIST ASD Team (2015).
#NIST Atomic Spectra Database (ver 5.3), [Online]. Available:
#http://physics.nist.gov/asd [2017, January 30]. National Institute of
#Standards and Technology, Gaithersburg, MD.
lineListBytes <- absPath + dataPath + "aLL1702bytes.rds"

print("READING LINE LIST")

con <- file(lineListBytes, "r", blocking=FALSE)
barray <- readLines(con, "r")

print("LINE LIST READ")

#decoded <- read.table(barray, fileEncoding = "UTF-8")
#arrayLineString <- strsplit(decoded, "%%")
#numLineList <- 1 #test

#Atomic Lines
list2Lam0 <- vector("numeric", numLineList) #nm
list2Element <- vector("character", numLineList) #element
list2StageRoman <- vector("character", numLineList) #ion stage
list2Stage <- list2Mass <- list2LogGammaCol <- list2LogAij <- list2Logf <- list2ChiI1 <- list2ChiI2 <- list2ChiI3 <- list2ChiI4 <- list2ChiI5 <- list2ChiI6 <- list2ChiL <- list2GwL <-  vector("numeric", numLineList)

list2_ptr <- 0 #pointer into line list2 that we're populating
numFields <- 7 #number of fields per record
#0: element, 1: ion stage, 2: lambda_0 3: logf 4: g_l 5: chi_l
thisRecord <- vector("character", numFields)

for (iLine in 1:len(numLineList)) {
 
    thisRecord <- strsplit(arrayLineString[iLine], "|")
    myString <- thisRecord[1]
    list2Element[iLine] <- myString
    myString <- thisRecord[2]
    list2StageRoman[iLine] <- myString
    myString <- thisRecord[3]
    list2Lam0[iLine] <- myString
    myString <- thisRecord[4]
    list2LogAij[iLine] <- myString
    myString <- thisRecord[5]
    list2Logf[iLine] <- myString
    myString <- thisRecord[6]
    list2ChiL[iLine] <- myString
    myString <- thisRecord[10]
    list2GwL[iLine] <- myString
    
    #Get the chemical element symbol - we don't know if it's one or two characters
    if (list2StageRoman[list2_ptr] == "I") {
       list2Stage[list2_ptr] = 0
    }
    if (list2StageRoman[list2_ptr] == "II") {
       list2Stage[list2_ptr] = 1
    }
    if (list2StageRoman[list2_ptr] == "III") {
       list2Stage[list2_ptr] = 2
    }
    if (list2StageRoman[list2_ptr] == "IV") {
       list2Stage[list2_ptr] = 3
    }
    if (list2StageRoman[list2_ptr] == "V") {
       list2Stage[list2_ptr] = 4
    }
    if (list2StageRoman[list2_ptr] == "VI") {
       list2Stage[list2_ptr] = 5
    }
    if (list2StageRoman[list2_ptr] == "VII") {
       list2Stage[list2_ptr] = 6
    }

#wavelength in nm starts at position 23 and is in %8.3f format 
#not expecting anything greater than 9999.999 nm
     list2Mass[list2_ptr] = getMass(list2Element[list2_ptr])
     species = glue(list2Element[list2_ptr], "I")
     list2ChiI1[list2_ptr] = getIonE(species)
     species = glue(list2Element[list2_ptr], "II")
     list2ChiI2[list2_ptr] = getIonE(species)
     species = glue(list2Element[list2_ptr], "III")
     list2ChiI3[list2_ptr] = getIonE(species)
     species = glue(list2Element[list2_ptr], "IV")
     list2ChiI4[list2_ptr] = getIonE(species)
     species = glue(list2Element[list2_ptr], "V")
     list2ChiI5[list2_ptr] = getIonE(species)
     species = glue(list2Element[list2_ptr], "VI")
     list2ChiI6[list2_ptr] = getIonE(species)

     list2LogGammaCol[list2_ptr] = logGammaCol

     list2_ptr = list2_ptr + 1
}

numLines2 <- list2_ptr

##END FILE I/O SECTION##

#initialize the accumulator 
gaussLineCntr <- 0 

#array of pointers to lines which make the cut in the array of
#pointers to lines that make the cut in the array of pointers to
#lines that make the cut in the master line list
gaussLine_ptr <- vector("numeric", numLines21)
isFirstLine <- TRUE #initialization
firstLine <- 0 #default initialization
#This holds 2-element temperature-dependent base 10 log partition fn
thisUwV <- vector("numeric", numAtmPrtTmps)
#Below created a loop to initialize each value to zero for the five 
#temperature lburns

for (iLine in 1:len(numLines2)) {
    if ((list2Element[iLine] == "H") || (list2Element[iLine] == "He"))
    {
	zScaleList = 1.0
	if (list2Lam0[iLine] <= 657.0) {
	  list2GwL[iLine] <- 8.0 #fix for Balmer lines
	} else {
	  list2GwL[iLine] <- 18.0 #fix for Paschen lines
	}
    } else {
      zScaleList <- zScale
    }

    list2Lam0[iLine] <- list2Lam0[iLine] * nm2cm #nm to cm
    iAbnd <- 0 #initialization
    logNums_ptr <- 0
    for (jj in 1:len(nelemAbnd)) {
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 0) {
	      species <- glue(cname[jj], "I")
	      logNums_ptr <- 0 
	      }
	   }
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 1) {
	      species <- glue(cname[jj], "II")
	      logNums_ptr <- 1 
	      }
	   }
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 2) {
	      species <- glue(cname[jj], "III")
	      logNums_ptr <- 4 
	      }
	   }
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 3) {
	      species <- glue(cname[jj], "IV")
	      logNums_ptr <- 5 
	      }
	   }
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 4) {
	      species <- glue(cname[jj], "V")
	      logNums_ptr <- 6 
	      }
	   }
    	if (list2Element[iLine] == cname[jj]) {
	   if (list2Stage[iLine] == 5) {
	      species <- glue(cname[jj], "VI")
	      logNums_ptr <- 7 
	      }
	   }
	   thisUwV <- getPartFn2(species) #base e log_e U
	   break
	}
    	iAbnd <- iAbnd + 1    	
    }
    list2LogNums <- matrix(0.0, numDeps, numStages+2)
    list2LogNums[1] <- masterStagePops[1,iAbnd]
    list2LogNums[2] <- masterStagePops[2,iAbnd]
    list2LogNums[5] <- masterStagePops[3,iAbnd]
    list2LogNums[6] <- masterStagePops[4,iAbnd]
    list2LogNums[7] <- masterStagePops[5,iAbnd]
    list2LogNums[8] <- masterStagePops[6,iAbnd]

    numHelp <- levelpops(lsit2Lam0[iLine], list2LogNums[logNums_ptr],
    list2ChiL[iLine], thisUwV, list2GwL[iLine], numDeps, temp)

    list2LogNums[3] <- numHelp
    list2LogNums[4] <- numHelp/2.0

    #linePoints: row 0 in cm (will need to be in nm for
    #Plack.planck). Row 1 in Doppler widths for now initial strength
    #check with delta fn profiles at line center for triage
    listNumPointsDelta <- 1
    listLinePointsDelta <- lineGridDelta(list2Lam0[iLine], list2Mass[iLine], xiT, numDeps, teff)
    listLineProfileDelta <- 
}