---
title: "Introduction to the *earthtide* package"
author: "Jonathan Kennel, Beth Parker"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{introduction}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


## Code background 

**earthtide** is a port of the 'Fortran ETERNA 3.4' (Wenzel, 1996) predict and part of the analyze codes with the Kudryavtsev 2004 update.  The original 'Fortran' code was rewritten in **R**, and **C++** using the great **Rcpp**, **RcppArmadillo**, and **RcppParallel**. The package is useful for generating synthetic earth tides using highly accurate tidal catalogs for prediction and regression. Attempts were made to ensure that results were consistent with the 'ETERNA 3.4', however, there is always the possibility that a bug was introduced in during the conversion and update.  For the most feature rich and up-to-date version of 'ETERNA' please see http://ggp.bkg.bund.de/eterna/ maintained by Klaus Schueller.


Hartmann, T., Wenzel, H.-G., 1995. The HW95 tidal potential catalogue. Geophys. Res. Lett. 22, 3553–3556. \url(https://doi.org/10.1029/95GL03324)

Kudryavtsev, S.M., 2004. Improved harmonic development of the Earth tide-generating potential. J. Geod. 77, 829–838. \url(https://doi.org/10.1007/s00190-003-0361-2)

Wenzel, H.G. 1996: The nanogal software: Earth tide data processing package ETERNA 3.30. Bull. Inf. Marges Terrestres. 124, 9425-9439. \url(http://www.eas.slu.edu/GGP/ETERNA34/MANUAL/ETERNA33.HTM)


## Setup

The first step to generate earthtides is to create the object with the prediction times, station, catalog and frequency details. The following are the most commonly used parameters:

- *utc* is a vector of POSIXct times for earth tide prediction.  This should typically be a regular time series (equally spaced) as the method is most efficient for this case, however, irregular time series are also accepted.
- *latitude* numeric of length one.  The latitude of the monitoring location.
- *longitude* numeric of length one.  The longitude of the monitoring location.
- *cutoff* numeric of length one. Cutoff determines how many waves to include in the analysis.  A lower number yields more precise results but at the cost of runtime. A value of 1e-10 includes all waves in the catalog and may be quite slow for long series. The default is 1e-6.
- *wave_groups* is a `data.frame` of of the frequency ranges in cycles per day. It requires a column named start and an column named end which give the ranges for each wave group.

See `?Earthtide` for the complete list of available arguments.

```{r newearthtide}
library(earthtide)

# One month of hourly data
tms <- as.POSIXct("2015-01-01", tz = "UTC") + 0:(24*31) * 3600

et <- Earthtide$new(
  utc = tms,
  latitude = 52.3868,
  longitude = 9.7144,
  wave_groups = data.frame(start = 0.0, end = 6.0))

```

## LOD (length of day) tide and pole tide

The Length of Day (LOD) and Pole tides are calculated in this initial step.

```{r lodpole, echo = FALSE, fig.width = 6.5, fig.height = 5}
tide <-et$lod_tide()$
  pole_tide()$
  tide()


layout(matrix(1:2, ncol=1, nrow = 2))
par(mai = c(0.4, 0.9, 0.1, 0.1))

# Plot the results

plot(lod_tide~datetime, tide,
     xlab = '',
     ylab = expression('LOD tide nm/s' ^ 2),
     type='l',
     lwd = 2,
     col = '#5696BC',
     las = 1)

plot(pole_tide~datetime, tide,
     xlab = '',
     ylab = expression('Pole tide nm/s' ^ 2),
     type='l',
     lwd = 2,
     col = '#5696BC',
     las = 1)

```

## Predict

The _**predict**_ method of the earthtide class computes the tidal response. There are two options for this function:

- *method* determines the which response is desired. Currently valid options are "gravity", "tidal_potential", "tidal_tilt", "vertical_displacement", "vertical_strain", "areal_strain", "volume_strain", and "ocean_tides". 
- *astro_update* determines how often phases are updated. For long time series (greater than 100,000 times) a value greater than 1 can speed up the results. In general, values greater than 1 will only be necessary for predicting at the minute or second scale or if you are predicting for many years.  There may be a slight loss of precision when the cutoff is small.  For small cutoff values use a small value for *astro_update* and we suggest testing to achieve the desired precision.
  

```{r predict}

et$predict(method = 'gravity', astro_update = 24)

```


```{r predictplot, fig.width = 6.5, fig.height = 3, fig.ext='png', echo = FALSE}
# Plot the results
par(mai = c(0.6, 0.9, 0.1, 0.1))
plot(gravity~datetime, et$tide(),
     ylab = expression('Gravity nm/s' ^ 2),
     xlab = '',
     type='l',
     lwd = 2,
     col = '#5696BC',
     xaxs = 'i',
     las = 1)
```


## Analyze

To generate normalized curves for use in regression equations we use the _**analyze**_ method instead of the _**predict**_ method.  All other arguments remain the same.  This will generate normalized sin and cos terms for each wavegroup.  In the example we plot the first 5 of a 22 wave groups set. These may be used in your preferred analyses.

```{r wg, fig.width = 6.5, fig.height = 8, fig.ext='png', echo = TRUE}
wave_groups <- na.omit(eterna_wavegroups[eterna_wavegroups$time == '1 month', 
                                        c('start', 'end')])

et <- Earthtide$new(utc = tms, 
                    latitude = 49.00937,
                    longitude = 8.40444,
                    elevation = 120,
                    cutoff = 1e-10,
                    catalog = 'hw95s',
                    wave_groups = wave_groups)

print(wave_groups[1:5,], row.names = FALSE)
```

```{r analyze}

et$analyze(method = 'gravity', astro_update = 1)

```

```{r analyzeplot, fig.width = 6.5, fig.height = 8, fig.ext='png', echo = FALSE}

tides <- et$tide()

layout(matrix(1:5, ncol=1, nrow = 5))
par(mai = c(0.3, 0.9, 0.1, 0.1))

for (i in seq(2, 11, 2)) {

  plot(tides[,1], tides[,i],
       ylab = expression('Gravity nm/s' ^ 2),
       xlab = '',
       type = 'l',
       lwd = 2,
       col = '#AAB6A2',
       las = 1)
  points(tides[,1], tides[,i+1],
         type = 'l',
         lwd = 2,
         col  = '#5696BC')

}


```


## Get the tidal values

The _**tide**_ method of the earthtide class returns a `data.frame` of tidal values. The following example in the chaining section shows how to calculate the gravity, LOD (length of day), and pole tides and return a `data.frame` of the values. 

## Chaining

It is also possible to chain operations together for more compact code.  

```{r r6chain, echo = TRUE}


tides <- Earthtide$
  new(utc = as.POSIXct("2015-01-01", tz = "UTC") + 0:(24*31) * 3600,
      latitude = 52.3868,
      longitude = 9.7144,
      wave_groups = data.frame(start = 0.0, end = 6.0))$
  predict(method = "gravity", astro_update = 1)$
  lod_tide()$
  pole_tide()$
  tide()
  

print(tides[1:5,], row.names = FALSE)

```


## TODO: 

This package is still in development.  The following changes are planned:

- Speed enhancements
- Print methods and getter/setters for love numbers, catalog, station variables
- Optional wave group names
- Ocean loading


