inst/doc/skyscapeR.R

In skyscapeR: Data Analysis and Visualization for Skyscape Archaeology

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## -----------------------------------------------------------------------------
# install.packages(skyscapeR)

## -----------------------------------------------------------------------------
# install.packages('devtools', repos='https://cloud.r-project.org')
# devtools::install_github('f-silva-archaeo/skyscapeR')

## -----------------------------------------------------------------------------
# install.packages("swephRdata", repos = "https://rstub.github.io/drat/", type = "source")

## -----------------------------------------------------------------------------
library(skyscapeR)

## -----------------------------------------------------------------------------
BC.AD(1)
BC.AD(0)
BC.AD(-1)
BC.AD(-99)
BC.AD(-100)

## -----------------------------------------------------------------------------
skyscapeR.vars() # shows current values
skyscapeR.vars(timezone='CET') # to set timezone to Central European Time
skyscapeR.vars(calendar='J') # to set calendar to Julian
skyscapeR.vars(refraction=F) # to switch refraction calculations off

## -----------------------------------------------------------------------------
skyscapeR.vars(timezone='GMT') # to set timezone to Central European Time
skyscapeR.vars(calendar='G') # to set calendar to Julian
skyscapeR.vars(refraction=T) # to switch refraction calculations off


## -----------------------------------------------------------------------------
az2dec(az=92, loc=c(35,-8,100), alt=2)

## -----------------------------------------------------------------------------
# geocentric declination for june solstice at 2501 BC
jS(-2500)

# topocentric declination for december solstice at 2501 BC from Stonehenge
dS(-2500, loc=c(51.17889,1.826111,10))

# geocentric declination for northern minor lunar extreme at 3001 BC
nmnLX(-3000)

# topocentric declination for northern minor lunar extreme at 3001 BC from Stonehenge
sMjLX(-3000, loc=c(51.17889,1.826111,10))

# zenith sun for Teotihuacan
zenith(loc=c(19.6925,98.84389,200))

## -----------------------------------------------------------------------------
mag.dec(loc=c(35,-8), date="2020/06/07")

## -----------------------------------------------------------------------------
loc <- c(35,-8,100)
mag.az <- c(89.5, 105, 109.5)
data <- reduct.compass(loc, mag.az, "2016/04/02", alt=c(1,2,0))
data

## -----------------------------------------------------------------------------
loc <- c(52,-3,100)
sunAz(loc, '2017-10-04 12:32:14', 'Europe/London')

## -----------------------------------------------------------------------------
lat <- ten(35,50,37.8) # latitude
lon <- ten(14,34,6.4) # longitude
elev <- 100
az <- c(ten(298,24,10), ten(302,20,40))
alt <- c(2, 5)
az.sun <- ten(327,29,50)
date <- "2016/02/20"
time <- "11:07:17"
data <- reduct.theodolite(c(lat,lon,elev), az, date , time, tz= "Europe/Malta", az.sun, alt=alt)
data

## ---- fig.width = 6, fig.height=4---------------------------------------------
az <- c(0,90,180,270,360)
alt <- c(0,5,5,0,0)
loc <- c(51.17889,1.826111,10) 
hor <- createHor(az, alt, alt.unc=0.5, loc=loc, name='Stonehenge Test')
plot(hor)

## ---- fig.width = 6, fig.height=4---------------------------------------------
hor <- createHWT(lat=ten(51,30,45), lon=ten(0,5,26.1), name='London Mithraeum')
plot(hor)

## ---- fig.width = 5, fig.height=5---------------------------------------------
az <- c(120, 100, 93, 97, 88, 115, 112, 67)
plotAzimuth(az)

## ---- fig.width = 5, fig.height=5---------------------------------------------
tt <- sky.objects('solar extremes', epoch=-2000, loc=c(35,-8), col='red')
plotAzimuth(az, obj=tt)

## ---- fig.width = 5, fig.height=5---------------------------------------------
tt <- sky.objects(c('solar extremes', 'lunar extremes'), epoch=-2000, loc=c(35,-8), col=c('red','blue'))
plotAzimuth(obj=tt)

## ---- fig.width = 5, fig.height=5---------------------------------------------
plotAzimuth(az, obj=tt, lwd=2, col='black')

## ---- fig.width = 6, fig.height=4---------------------------------------------
az <- c(120, 100, 93, 97, 88, 115, 112, 67)
hor <- createHWT(51.17889, 1.826111, name='Stonehenge')
dec <- az2dec(az, loc=hor)
kde <- density(dec)
plot(kde)

## ---- fig.width = 6, fig.height=4---------------------------------------------
kde <- density(dec, bw=2) # forces uncertainty to be a given value, in this case 2 degrees
plot(kde)

kde <- density(dec, bw=2, kernel='epanechnikov') # changes kernel to epanechnikov
plot(kde)

## -----------------------------------------------------------------------------
bernoulli.trial(n=8, p=12/180, r=4)

## ---- fig.width = 6, fig.height=4---------------------------------------------
az <- c(120, 100, 93, 97, 88, 115, 112, 67) # same azimuths as before
unc <- c(2, 3, 2, 2.5, 1.5, 3, 4, 3.5)
az.prob <- az.pdf(az=az, unc=unc)
plot(az.prob)

## ---- fig.width = 6-----------------------------------------------------------
dec.prob <- coordtrans(az.prob, hor) # same horizon as before
plot(dec.prob)

## ---- fig.width = 6, fig.height=4---------------------------------------------
ss <- spd(dec.prob)
plot(ss)

## ---- fig.width = 6, fig.height=4---------------------------------------------
# on azimuth data
st.az <- randomTest(az.prob, nsims=10, ncores=1) 
plot(st.az)

# on declination data
st.dec <- randomTest(dec.prob, nsims=10, ncores=1)
plot(st.dec)

## ---- fig.width = 6, fig.height=4---------------------------------------------
plot(st.dec, show.local=T)
st.dec$metadata$local.pval

## -----------------------------------------------------------------------------
findTargets(c(-25,-18), c(-2499,-1749))

## -----------------------------------------------------------------------------
findTargets(c(-25,-18), c(-2499,-1749), max.mag=3)

## ---- results='hide'----------------------------------------------------------
sp <- star.phases('Sirius', -2999, loc=c(30.0, 31.2, 25), alt.hor=2)

## -----------------------------------------------------------------------------
sp$metadata$type
sp$metadata$events
sp$metadata$seasons

## ---- fig.width = 6-----------------------------------------------------------
plot(sp)

## ---- fig.width = 6-----------------------------------------------------------
sky.sketch(time='2019/01/07 09:30', loc=c(35,-8,100))