R/lineGrid.R

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

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

useful <- readRDS("./data/useful.rds")
c <- useful[,"c"]
logC <- useful[,"logC"]
logK <- useful[,"logK"]
amu <- useful[,"amu"]
ln10 <- log(10)
ln2 <- log(2.0)
logE <- log10(exp(1))

lineGridDelta <- function(lam0In, massIn, xiTIn, numDeps, teff) {
  #put input parameters into linear cgs units
  logTeff <- log(teff)
  xiT <- xiTIn*1E5
  lam0 <- lam0In
  logLam0 <- log(lam0)
  logMass <- log(massIn * amu)

  #compute depth-independent Doppler width, Delta_lambda_D
  logHelp <- ln2+logK+logTeff-logMass #M-B dist, square of v_mode
  help <- exp(logHelp)+xiT*xiT #quadratic sum of thermal v and turbulent v
  logHelp <- 0.5*log(help)
  logDopp <- logHelp+logLam0-logC
  doppler <- exp(logDopp)

  numPoints <- 1
  linePoints <- matrix(0.0, numPoints, 2)
  v <- vector("numeric", numPoints)
  il <- 0
  ii <- il

  #in core, space v poits linearly
  v[il] <- ii
  linePoints[il,1] <- doppler*v[il]
  linePoints[il,2] <- v[il]

  return linePoints
}

lineGridGauss <- function(lam0In, massIn, xiTIn, numDeps, teff, numCore) {
  #put input parameters into linear cgs units
  logTeff <- log(teff)
  xiT <- xiTIn*1E5
  lam0 <- lam0In
  logLam0 <- log(lam0)
  logMass <- log(massIn * amu)

  #compute depth-independent Doppler width, Delta_lambda_D
  logHelp <- ln2+logK+logTeff-logMass #M-B dist, square of v_mode
  help <- exp(logHelp)+xiT*xiT #quadratic sum of thermal v and turbulent v
  logHelp <- 0.5*log(help)
  logDopp <- logHelp+logLam0-logC
  doppler <- exp(logDopp)

  #set up a half-profile Delta_lambda grid in Doppler width units from line
  #centre to wing
  numPoints <- numCores

  #a 2D 2 X numPoints array of Delta Lambdas
  linePoints <- matrix(0.0, numPoints, 2)

  #line profile points in Doppler widths - needed for Voigt function, H(a,v)
  v <- vector("numeric", numPoints)
  maxCoreV <- 3.5 #core half-width - in Doppler widths

  minWingDeltaLogV <- log(maxCoreV + 1.5)
  maxWingDeltaLogV <- 9.0 + minWingDeltaLogV

  for (il in 1:numPoints) {
      ii <- il

      #in core, space v points linearily
      v[il] <- ii*maxCoreV/(numCore-1)
      linePoints[il,1] <- doppler*v[il]
      linePoints[il,2] <- v[il]
  }
  #add the negative DeltaLambda half of the line
  numPoints2 <- (2*numPoints)-1

  #return a 2D 2 X (2XnumPoints-1) array of Delta Lambdas
  #row 1: delta lambdas in cm - will need to be in nm for Planck and Rad Trans
  #row 2: delta lambdas in Doppler widths
  linePoints2 <- matrix(0.0, numPoints2, 2)

  #wavelentsh are depth-independent - just put them in the 0th depth slot
  for (il2 in 1:numPoints2) {
    if (il2 < numPoints - 1) {
      il <- (numPoints - 1) - il2
      linePoints2[il2,1] <- -1.0*linePoints[il,1]
      linePoints2[il2,2] <- -1.0*linePoints[il,2]
    } else {
      il <- il2-(numPoints-1)
      linePoints2[il2,1] <- linePoints[il,1]
      linePoints2[il2,2] <- linePoints[il,2]
    }
  }
  return linePoints2
}

lineGridVoigt <- function(lam0In, massIn, xiTIn, numDeps, teff, numCore, numWing, species) {
  c <- useful[,"c"]
  logC <- useful[,"logC"]
  logK <- useful[,"logK"]
  amu <- useful[,"amu"]

  # put input parameters into linear cgs units
  logTeff <- log(teff)
  xiT <- xiTIn * 1.0E5
  lam0 <- lam0In
  logLam0 <- log(lam0)
  logMass <- log(massIn * amu)

  logHelp <- ln2 + logK + logTeff - logMass
  help <- exp(logHelp) + xiT * xiT # quadratic sum of thermal v and turbulent v
  logHelp <- 0.5 * log(help)
  logDopp <- logHelp + logLam0 - logC

  doppler <- exp(logDopp)

  # set up a half-profile Delta_lambda grid in Doppler width units
  numPoints <- numCore + numWing

  # a 2D 2 x numpoints array of Delta Lambdas
  linePoints = matrix(0, numPoints, 2)

  # line profile points in Doppler widths - needed for Voigt function, H(a,v)
  v = vector("numeric", numpoints)

  macCoreV <- 3.5 # core half-widths - in Doppler widths
  minWingDeltaLogV <- log(maxCoreV + 1.5)
  maxWingDeltaLogV <- 9.0 + minWingDeltaLogV

  if (species == "HI" && teff >= 7000) {
     maxCorev <- 3.5
     minWingDeltaLogV <- log(maxCoreV + 1.5)
     maxWingDeltaLogV <- 12.0 + minWingDeltaLogV
  }

  for (il in 1:numPoints) {
      if (il < numCore) {
      	 # in core, space v points linearly
	 v[il] <- il * maxCoreV / (numCore - 1)
	 linePoints[il][1] <- doppler * v[il]
	  linePoints[il][2] <- v[il]
      } else {
      	 # space v points logarithmically in wing
	 jj <- ii - numCore
	 logV <- (jj * (maxWingDeltaLogV - minWingDeltaLogV) / (numPoints - 1)) + minwingDeltaLogV
	 v[il] <- exp(logV)
	 linePoints[il][1] <- doppler * v[il]
	 linePoints[il][2] <- v[il]
      }
  }

  # add the negative DeltaLambda half of the line
  numPoints2 <- (2 * numPoints) - 1

  # return a 2D 2 x (2 x numPoints - 1) array of Delta Lambdas
  # row 0: delta lambdas in cm
  # row 1: delta lambdas in Doppler widths
  linePoints2 <- matrix(0, numPoints2, 2)

  # wavelength are depth-independent - just put them in the 0th depth slot
  for (il2 in 1:numPoints2) {
      if (il2 < numPoints - 1) {
      	 il <- (numPoints - 1) - il2
	 linePoints2[il2][1] <- -1.0 * linePoints[il][1]
	 linePoints2[il2][2] <- -1.0 * linePoints[il][2]
      } else {
      	il <- il2 - (numPoints - 1)
	linePoints2[il2][1] <- linePoints[il][1]
	linePoints2[il2][2] <- linePoints[il][2]
      }
  }

linePoints2
 
}