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

#Collection of useful utilities

#Linear interpolation to a new abscissa
interpolate <- function(x, y, x.new) {
  p1 <- 1
  p2 <- 2
  x1 <- x[p1]
  x2 <- x[p2]

  for (i in 1:length(x)) {
    if (is.na(x[i]) || is.na(x.new)) {
      break
    }
    if (x[i] >= x.new) {
      p2 <- i
      p1 <- i-1
      x2 <- x[p2]
      x1 <- x[p1]
      break
    }
  }
  step <- x2 - x1
  newY <- y[p2]*(x.new-x1)/step+y[p1]*(x2-x.new)/step
  newY
}

#Vectorized version of simple linear 1st order interpolation, assumes new and old
#abscissae are in monotonic increasing order
interpolateV <- function(y, x, x.new) {
  num <- length(x)
  newNum <- length(x.new)

  #renormalize ordinates
  iMinAndMax <- minMax(y)
  norm <- y[iMinAndMax[2]]
  newYNorm <- rep(0.0, newNum)
  yNorm <- y/norm

  #set any x.new element that are *less than* the first x element to the first x
  #element - "0th order extrapolation"
  start <- 0
  for (i in 1:newNum) {
    if (x.new[i] <= x[2]) {
      newYNorm[i] = yNorm[1]
      start <- start+1
    }
    if (newx[i] > x[1]) {
      break
    }
  }
 if (start < newNum) {
   j <- 1 #initialize old abscissae index
   #outer loop over new abscissae
   for (i in start:newNum) {
     #break if current element newX is *greater* than last x element
     if ((x.new[i]>x[num]) || (j>num)) {
       break
     }
     while (x[j] < x.new[i]) {
       j <- j+1
     }
     jWght <- x.new[i]*(1.0-(x[j-1]/x.new[i]))
     jm1Wght <- x[j]*(1.0-(x.new[i]/x[j]))
     denom <- x[j]*(1.0-(x[j-1]/x[j]))
     jWght <- jWgth/denom
     jm1Wght <- jm1Wght/denom
     newYNorm[i] <- (yNorm[j]*jWght)+(yNorm[j-1]*jm1Wght)
   }
 }
  #set any newX elements that are *greater than* the first x element to the last x
  #element - "0th order extrapolation"
  for (i in 1:newNum) {
    if (x.new[i] >= x[num]) {
      newYNorm[i] <- yNorm[num]
    }
  }
  #restore orinate scale
  newY <- rep(x*nrom,newYNorm)
  newY
}

#Return the array index of the wavelength array (lambdas) closest to a desired value
#of wavelength (lam)
lamPoint <- function(numLams, lambdas, lam) {
  help <- rep(0.0,numLams)
  for (i in 1:numLams) {
    help[i] <- abs(lambdas[i] - lam)
  }

  index <- 1
  min <- help[index]
  for (i in 1:numLams) {
    if (help[i] < min) {
      min <- help[i]
      index <- i
    }
  }
  index
}

#Return the minimum and maximum values of an input 10 array. Will return the *first*
#occurance if min and/or max values occur in multiple places. iMinMax[1]  = first
#occurance of minimum iMinMax[2] = first occurance of maximum
minMax <- function(x) {
  iMinMax <- rep(0,2)
  num <- length(x)
  iMin <- 1
  iMax <- 1
  min <- x[iMin]
  max <- x[iMax]

  for (i in 1:num) {
    if (x[i] < min) {
      min <- x[i]
      iMin <- i
    }
    if (x[i] > max) {
      max <- x[i]
      iMax <- i
    }
  }
  iMinMax[1] <- iMin
  iMinMax[2] <- iMax
  iMinMax
}

#Return the minimum and maximum of an input 1D array. Will return the *first*
#occurance if min and/or max values occur in multiple places iMinMax[1] = first
#occurance of minimum iMinMax[2] = first occurance of maximum
minMax2 <- function(x) {
  iMinMax <- rep(0.0,2)
  num <- length(x)

  iMin <- 1
  iMax <- 1

  #search for minimum and maximum in row 1 - linear values
  min <- x[iMin][1]
  max <- x[iMax][1]

  for (i in 2:num) {
    if (x[i, 1] < min) {
      min <- x[i, 1]
      iMin <- i
    }
    if (x[i][1] > max) {
      max <- max[i][1]
      iMax <- i
    }
  }
  iMinMax[1] <- iMin
  iMinMax[2] <- iMax

  iMinMax
}

#Return the array index of the optical depth array (tauRos) closest to a desired
#value of optical depth (tau). Assuming the user wants to find a *linear* tau value
#and NOT a logarithmic one
tauPoint <- function(numDeps, tauRos, tau) {
  holder <- rep(0.0, numDeps)

  for (i in 1:numDeps) {
    holder[i] <- abs(tauRos[i, 1]-tau)
  }
  index <- 1
  min <- holder[index]

  for (i in 1:numDeps) {
    if (holder[i] < min) {
      min <- holder[i]
      index <- i
    }
  }
  index
}