---
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)
library(knitr)
library(earthtide)
```


## 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 **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 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.

## Capabilities

The following tidal componenents are implemented in **earthtide**.

```{r capabilities, echo = FALSE, results = 'asis'}

tidal_component <- c(
  "tidal_potential", "gravity",  "tidal_tilt", "vertical_displacement", 
  "horizontal_displacement", "n_s_displacement", "e_w_displacement",
  "vertical_strain", "areal_strain", "volume_strain", "horizontal_strain",
  "ocean_tides")

output_units <- c('$meters^2/second^2$','$nanometers/second^2$', '$milliarcsec$', '$millimeter$',
                  '$millimeter$', '$millimeter$', '$millimeter$', 
                  '$nanostrain$', '$nanostrain$', '$nanostrain$', '$nanostrain$',
                  '$millimeter$')

status <- c("tested", "tested", "tested", "tested", 
            "preliminary", "preliminary","preliminary", 
            "tested","tested","tested","tested",
            "preliminary")
dat <- data.frame(tidal_component, status, output_units)
names(dat) <- c('Tidal component', 'Status', 'Output units')
kable(dat)
```

## Example

The primary inputs are the date-time in UTC, the component name from the previous table, and the latitude and longitude.  For most cases these are the minimum requirements necessary. For the full list of options see the documentation for _calc_earthtide_.

```{r standardmethod, echo = TRUE}
tms <- as.POSIXct("2015-01-01", tz = "UTC") + 0:(24*31) * 3600

gravity_tide <- calc_earthtide(utc = tms, 
                               method = 'gravity',
                               latitude = 52.3868,
                               longitude = 9.7144)


```


There are two main methods of generating Earth tides: predict and analyze. Predict returns the combined tidal signal, and analyze returns a set of sin and cos curves for each wave group that is specified. This option is set using the **do_predict** parameter which defaults to TRUE.

### Predict

```{r predict, echo = TRUE}

gravity_tide <- calc_earthtide(utc = tms,
                               do_predict = TRUE,
                               method = 'gravity',
                               latitude = 52.3868,
                               longitude = 9.7144)

```


```{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, gravity_tide,
     ylab = expression('Gravity nm/s' ^ 2),
     xlab = '',
     type='l',
     lwd = 2,
     col = '#5696BC',
     xaxs = 'i',
     las = 1)
```

### Analyze

In analyze mode, results are separated by wave group into sin and cos curves.  The resulting sin and cos curves can be used in further analysis such as least squares models.  The first five constituents are plotted in the following example.  The **wave_groups** parameter is specified using a *data.frame* having the start and end frequencies for each component.

```{r analyze, echo = TRUE}
wg <- eterna_wavegroups
wg <- na.omit(wg[wg$time=='1 month',])

tides <- calc_earthtide(utc = tms,
                        do_predict = FALSE,
                        method = 'gravity',
                        latitude = 52.3868,
                        longitude = 9.7144,
                        wave_groups = wg)

```


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

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')

}


```


## Wave grouping

The choice of wave groups is important.  Wave groups are specified using a *data.frame* of start and end values for each group.  Example groupings are provided in the dataset **eterna_wavegroups**.  The choice of the appropriate wave groups is dependent purpose of the study and the duration of the dataset to be analyzed. For example, if the goal is to generate tidal harmonics to analyze one (1) month of data you would select the "1 month" dataset.

```{r analyze1month, echo = TRUE}

tms <- as.POSIXct("2015-01-01", tz = "UTC") + 0:(24*31) * 3600

wg <- eterna_wavegroups
wg <- na.omit(wg[wg$time=='1 month',])

head(wg)
```

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

The Length of Day (LOD) and Pole tides can also be calculated. These results differ from ETERNA in that we interpolate using splines.

```{r lodpolecalc, echo = TRUE}
tide <- calc_earthtide(utc = tms,
                       method = c('lod_tide', 'pole_tide'),
                       latitude = 52.3868,
                       longitude = 9.7144)
```

```{r lodpoleplot, echo = FALSE, fig.width = 6.5, fig.height = 5}

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)

```

## Speed for large datasets

There are two primary ways to increase the speed of generating the tidal datasets: increasing the **astro_update** parameter; and increasing the **cutoff** parameter. In general, the speed should be as good as or better than the fortran version of ETERNA, given the parallel computation.  The **cutoff** parameter determines the number of waves used in the analysis.  A larger cutoff value means fewer waves will be used leading to a faster but less accurate result.  The **astro_update** parameter determines how often the astronomical arguments are updated.  In ETERNA this value ranges from daily to yearly. The current default in the earthtide package for **astro_update** is 1, which means that the parameters are updated for every time provided, and therefore if you are predicting in hourly increments it will be updated 24 times a day.  Increasing this value typically speeds up computation, and in general values larger than 1 can be safely used unless the time interval between predictions is large. 


## References

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)


## 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
- Updated examples


