renamed and stuff

particle-track-and-trace/src/CMakeLists.txt

@@ -0,0 +1,92 @@
+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+
+project(ParticleTrackTrace)
+
+set(CMAKE_CXX_STANDARD 23)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+find_package(VTK COMPONENTS 
+  CommonColor
+  CommonColor
+  CommonCore
+  CommonDataModel
+  FiltersGeneral
+  FiltersGeometry
+  FiltersProgrammable
+  FiltersSources
+  ImagingSources
+  InteractionStyle
+  IOImage
+  RenderingContextOpenGL2
+  RenderingCore
+  RenderingCore
+  RenderingFreeType
+  RenderingGL2PSOpenGL2
+  RenderingOpenGL2)
+
+
+if (NOT VTK_FOUND)
+  message(FATAL_ERROR "VtkBase: Unable to find the VTK build folder.")
+endif()
+
+# netcdf setup
+find_package(netCDF REQUIRED)
+
+add_executable(ParticleTrackTrace MACOSX_BUNDLE main.cpp
+        layers/BackgroundImage.cpp
+        layers/BackgroundImage.h
+        layers/EGlyphLayer.cpp
+        layers/EGlyphLayer.h
+        layers/Layer.cpp
+        layers/Layer.h
+        layers/LGlyphLayer.cpp
+        layers/LGlyphLayer.h
+        Program.cpp
+        Program.h
+  commands/TimerCallbackCommand.h
+  commands/TimerCallbackCommand.cpp
+        commands/SpawnPointCallback.h
+        commands/SpawnPointCallback.cpp
+        CartographicTransformation.cpp
+        advection/AdvectionKernel.h
+        advection/EulerAdvectionKernel.cpp
+        advection/EulerAdvectionKernel.h
+        advection/interpolate.cpp
+        advection/interpolate.h
+        advection/readdata.cpp
+        advection/readdata.h
+        advection/RK4AdvectionKernel.cpp
+        advection/RK4AdvectionKernel.h
+        advection/UVGrid.cpp
+        advection/UVGrid.h
+        advection/Vel.cpp
+        advection/Vel.h
+)
+
+execute_process(
+    COMMAND nc-config --includedir
+    OUTPUT_VARIABLE NETCDF_INCLUDE_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+execute_process(
+    COMMAND ncxx4-config --libdir
+    OUTPUT_VARIABLE NETCDFCXX_LIB_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE
+)
+
+target_include_directories(ParticleTrackTrace PUBLIC ${netCDF_INCLUDE_DIR})
+
+find_library(NETCDF_LIB NAMES netcdf-cxx4 netcdf_c++4 PATHS ${NETCDFCXX_LIB_DIR} NO_DEFAULT_PATH)
+
+# Prevent a "command line is too long" failure in Windows.
+set(CMAKE_NINJA_FORCE_RESPONSE_FILE "ON" CACHE BOOL "Force Ninja to use response files.")
+target_link_libraries(ParticleTrackTrace ${NETCDF_LIB} ${VTK_LIBRARIES})
+
+# vtk_module_autoinit is needed
+vtk_module_autoinit(
+  TARGETS ParticleTrackTrace
+  MODULES ${VTK_LIBRARIES}
+)

particle-track-and-trace/src/CartographicTransformation.cpp

@@ -0,0 +1,46 @@
+#include "CartographicTransformation.h"
+#include <vtkMatrix4x4.h>
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+
+vtkSmartPointer<vtkCamera> createNormalisedCamera() {
+    vtkSmartPointer<vtkCamera> camera = vtkSmartPointer<vtkCamera>::New();
+    camera->ParallelProjectionOn();  // Enable parallel projection
+
+    camera->SetPosition(0, 0, 1000);  // Place the camera above the center
+    camera->SetFocalPoint(0, 0, 0);  // Look at the center
+    camera->SetViewUp(0, 1, 0); // Set the up vector to be along the Y-axis
+    camera->SetParallelScale(1); // x,y in [-1, 1]
+
+    return camera;
+}
+
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid) {
+    const double XMin = uvGrid->lons.front();
+    const double XMax = uvGrid->lons.back();
+    const double YMin = uvGrid->lats.front();
+    const double YMax = uvGrid->lats.back();
+
+    double eyeTransform[] = {
+            2/(XMax-XMin), 0, 0, -(XMax+XMin)/(XMax-XMin),
+            0, 2/(YMax-YMin), 0, -(YMax+YMin)/(YMax-YMin),
+            0, 0, 1, 0,
+            0, 0, 0, 1
+    };
+
+    auto matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+    matrix->DeepCopy(eyeTransform);
+    return matrix;
+}
+
+// Assumes Normalised camera is used
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid) {
+    vtkNew<vtkTransform> transform;
+
+    transform->SetMatrix(getCartographicTransformMatrix(uvGrid));
+
+    vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
+    transformFilter->SetTransform(transform);
+
+    return transformFilter;
+}

particle-track-and-trace/src/CartographicTransformation.h

@@ -0,0 +1,31 @@
+#include <vtkCamera.h>
+#include <vtkTransformFilter.h>
+#include "advection/UVGrid.h"
+
+#ifndef NORMALISEDCARTOGRAPHICCAMERA_H
+#define NORMALISEDCARTOGRAPHICCAMERA_H
+
+
+/**
+ * Constructs a orthographically projected camera that looks at the square x,y in [-1, 1] with z = 0 and w = 1.
+ * The space [-1,1] x [-1,1] x {0} will be referred to as the normalised space.
+ * @return pointer to camera
+ */
+vtkSmartPointer<vtkCamera> createNormalisedCamera();
+
+/**
+ * Constructs a 4x4 projection matrix that maps homogenious (longitude, latitude, 0, 1) points
+ * to the normalised space.
+ * TODO: This transformation has room for improvement see:
+ * https://github.com/MakeNEnjoy/interactive-track-and-trace/issues/12
+ * @return pointer to 4x4 matrix
+ */
+vtkSmartPointer<vtkMatrix4x4> getCartographicTransformMatrix(const std::shared_ptr<UVGrid> uvGrid);
+
+/**
+ * Convenience function that converts the 4x4 projection matrix into a vtkTransformFilter
+ * @return pointer to transform filter
+ */
+vtkSmartPointer<vtkTransformFilter> createCartographicTransformFilter(const std::shared_ptr<UVGrid> uvGrid);
+
+#endif //NORMALISEDCARTOGRAPHICCAMERA_H

particle-track-and-trace/src/Program.cpp

@@ -0,0 +1,88 @@
+#include <vtkRenderWindow.h>
+#include <vtkPointData.h>
+#include <vtkDoubleArray.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkRegularPolygonSource.h>
+#include <vtkGlyph2D.h>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <netcdf>
+#include <vtkArrowSource.h>
+#include <vtkNew.h>
+#include <vtkCallbackCommand.h>
+#include <vtkInteractorStyleUser.h>
+
+#include "Program.h"
+#include "commands/TimerCallbackCommand.h"
+#include "commands/SpawnPointCallback.h"
+
+void Program::setWinProperties() {
+  this->win->SetWindowName("Simulation");
+  this->win->SetSize(661, 661);
+  this->win->SetDesiredUpdateRate(60);
+
+  this->interact->SetRenderWindow(this->win);
+  this->interact->Initialize();
+
+  vtkNew<vtkInteractorStyleUser> style;
+  interact->SetInteractorStyle(style);
+}
+
+void Program::setupTimer(int dt) {
+  auto callback = vtkSmartPointer<TimerCallbackCommand>::New(this);
+  callback->SetClientData(this);
+  callback->setDt(dt);
+  this->interact->AddObserver(vtkCommand::TimerEvent, callback);
+  this->interact->AddObserver(vtkCommand::KeyPressEvent, callback);
+  this->interact->CreateRepeatingTimer(17); // 60 fps == 1000 / 60 == 16.7 ms per frame
+}
+
+Program::Program(int timerDT) {
+  this->win = vtkSmartPointer<vtkRenderWindow>::New();
+  this->interact = vtkSmartPointer<vtkRenderWindowInteractor>::New();
+
+  this->win->SetNumberOfLayers(0);
+  setWinProperties();
+  setupTimer(timerDT);
+}
+
+
+void Program::addLayer(Layer *layer) {
+  this->layers.push_back(layer);
+  this->win->AddRenderer(layer->getLayer());
+  this->win->SetNumberOfLayers(this->win->GetNumberOfLayers() + 1);
+}
+
+void Program::removeLayer(Layer *layer) {
+  this->win->RemoveRenderer(layer->getLayer());
+
+  auto it = std::find(this->layers.begin(), this->layers.end(), layer);
+  if (it != this->layers.end()) {
+    this->layers.erase(it);
+    this->win->SetNumberOfLayers(this->win->GetNumberOfLayers() - 1);
+  }
+}
+
+void Program::updateData(int t) {
+  win->Render();
+  for (Layer *l: layers) {
+    l->updateData(t);
+  }
+}
+
+void Program::setupInteractions() {
+  for (Layer *l: layers) {
+    l->addObservers(interact);
+  }
+}
+
+void Program::render() {
+  setupInteractions();
+  win->Render();
+  interact->Start();
+}

particle-track-and-trace/src/Program.h

@@ -0,0 +1,68 @@
+#ifndef PROGRAM_H
+#define PROGRAM_H
+
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+
+#include "layers/Layer.h"
+#include "commands/SpawnPointCallback.h"
+
+/** This class manages the upper levels of the vtk pipeline; it has attributes for the vtkrenderWindow and a vector of Layers to represent a variable number of vtkRenderers.
+  * It can also set up a vtkTimer by connecting an instance of TimerCallbackCommand with its contained vtkRenderWindowInteractor.
+  */ 
+class Program {
+private:
+  /** This attribute models a variable number of vtkRenderers, managed through the abstract Layer class.
+    */ 
+  std::vector<Layer *> layers;
+
+  /** The window this program's layers render to.
+    */ 
+  vtkSmartPointer<vtkRenderWindow> win;
+
+  /** The interactor through which the layers can interact with the window.
+    */ 
+  vtkSmartPointer<vtkRenderWindowInteractor> interact;
+
+  /** This function sets some default properties on the vtkRenderWindow. Extracted to its' own function to keep the constructor from becoming cluttered.
+    */ 
+  void setWinProperties();
+
+  /** This function sets up and connects a TimerCallbackCommand with the program.
+    */ 
+  void setupTimer(int dt);
+
+  void setupInteractions();
+
+public:
+  /** Constructor.
+    */ 
+  Program(int timerDT);
+
+  /** This function adds a new layer (and thus vtkRenderer) to the program.
+    * The layer is expected to set its own position in the vtkRenderWindow layer system.
+    * @param layer : pointer to the layer to add.
+    */ 
+  void addLayer(Layer *layer);
+
+  /** This function removes a given layer from the vtkRenderWindow and layers vector.
+    * If the given layer is not actually in the program, nothing happens.
+    * @param layer : the layer to removeLayer
+    */
+  void removeLayer(Layer *layer);
+
+  /** This function updates the data for the associated layers to the given timestamp.
+    * Also updates the renderWindow.
+    * @param t : the timestamp to update the data to.
+    */
+  void updateData(int t);
+
+  /**
+   * This function renders the vtkRenderWindow for the first time.
+   * Only call this function once!
+   */
+  void render();
+};
+
+#endif

particle-track-and-trace/src/advection/AdvectionKernel.h

@@ -0,0 +1,31 @@
+#ifndef ADVECTIONKERNEL_H
+#define ADVECTIONKERNEL_H
+
+#include <tuple>
+
+#include "Vel.h"
+
+/*
+ * Implement this class for every integration method.
+ */
+class AdvectionKernel {
+public:
+    /**
+     * This function must take a time, latitude and longitude of a particle and must output
+     * a new latitude and longitude after being advected once for AdvectionKernel::DT time as defined above.
+     * @param time Time since the beginning of the data
+     * @param latitude Latitude of particle
+     * @param longitude Longitude of particle
+     * @return A pair of latitude and longitude of particle.
+     */
+    virtual std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const = 0;
+
+    // Taken from Parcels https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/tools/converters.py#L155
+    const static double metreToDegrees(double metre) {
+        return metre / 1000. / 1.852 / 60.;
+    }
+
+    virtual ~AdvectionKernel() = default; // Apparently I need this, idk why
+};
+
+#endif //ADVECTIONKERNEL_H

particle-track-and-trace/src/advection/EulerAdvectionKernel.cpp

@@ -0,0 +1,13 @@
+
+#include "EulerAdvectionKernel.h"
+#include "interpolate.h"
+
+using namespace std;
+
+EulerAdvectionKernel::EulerAdvectionKernel(std::shared_ptr<UVGrid> grid) : grid(grid) {}
+
+std::pair<double, double> EulerAdvectionKernel::advect(int time, double latitude, double longitude, int dt) const {
+  auto [u, v] = bilinearinterpolate(*grid, time, latitude, longitude);
+
+  return {latitude + metreToDegrees(v * dt), longitude + metreToDegrees(u * dt)};
+}

particle-track-and-trace/src/advection/EulerAdvectionKernel.h

@@ -0,0 +1,25 @@
+#ifndef EULERADVECTIONKERNEL_H
+#define EULERADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+
+/**
+ * Implementation of AdvectionKernel.
+ * The basic equation is:
+ * new_latitude = latitude + v*DT
+ * new_longitude = longitude + u*DT
+ *
+ * Uses bilinear interpolation as implemented in interpolate.h
+ */
+class EulerAdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit EulerAdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+
+};
+
+
+#endif //EULERADVECTIONKERNEL_H

particle-track-and-trace/src/advection/RK4AdvectionKernel.cpp

@@ -0,0 +1,35 @@
+#include "RK4AdvectionKernel.h"
+#include "interpolate.h"
+
+using namespace std;
+
+RK4AdvectionKernel::RK4AdvectionKernel(std::shared_ptr<UVGrid> grid): grid(grid) { }
+
+std::pair<double, double> RK4AdvectionKernel::advect(int time, double latitude, double longitude, int dt) const {
+    auto [u1, v1] = bilinearinterpolate(*grid, time, latitude, longitude);
+//    lon1, lat1 = (particle.lon + u1*.5*particle.dt, particle.lat + v1*.5*particle.dt);
+    double lon1 = longitude + metreToDegrees(u1 * 0.5*dt);
+    double lat1 = latitude + metreToDegrees(v1 * 0.5*dt);
+
+//    (u2, v2) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat1, lon1, particle]
+    auto [u2, v2] = bilinearinterpolate(*grid, time + 0.5 * dt, lat1, lon1);
+
+//    lon2, lat2 = (particle.lon + u2*.5*particle.dt, particle.lat + v2*.5*particle.dt)
+    double lon2 = longitude + metreToDegrees(u2 * 0.5 * dt);
+    double lat2 = latitude + metreToDegrees(v2 * 0.5 * dt);
+
+//    (u3, v3) = fieldset.UV[time + .5 * particle.dt, particle.depth, lat2, lon2, particle]
+    auto [u3, v3] = bilinearinterpolate(*grid, time + 0.5 * dt, lat2, lon2);
+
+//    lon3, lat3 = (particle.lon + u3*particle.dt, particle.lat + v3*particle.dt)
+    double lon3 = longitude + metreToDegrees(u3 * dt);
+    double lat3 = latitude + metreToDegrees(v3 * dt);
+
+//    (u4, v4) = fieldset.UV[time + particle.dt, particle.depth, lat3, lon3, particle]
+    auto [u4, v4] = bilinearinterpolate(*grid, time + dt, lat3, lon3);
+
+    double lonFinal = longitude + metreToDegrees((u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt);
+    double latFinal = latitude + metreToDegrees((v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt);
+
+    return {latFinal, lonFinal};
+}

particle-track-and-trace/src/advection/RK4AdvectionKernel.h

@@ -0,0 +1,22 @@
+#ifndef RK4ADVECTIONKERNEL_H
+#define RK4ADVECTIONKERNEL_H
+
+#include "AdvectionKernel.h"
+#include "UVGrid.h"
+
+/**
+ * Implementation of Advection kernel using RK4 integration
+ * See https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods for more details.
+ * Uses bilinear interpolation as implemented in interpolate.h
+ */
+class RK4AdvectionKernel: public AdvectionKernel {
+private:
+    std::shared_ptr<UVGrid> grid;
+public:
+    explicit RK4AdvectionKernel(std::shared_ptr<UVGrid> grid);
+    std::pair<double, double> advect(int time, double latitude, double longitude, int dt) const override;
+
+};
+
+
+#endif //RK4ADVECTIONKERNEL_H

particle-track-and-trace/src/advection/UVGrid.cpp

@@ -0,0 +1,66 @@
+#include <ranges>
+
+#include "UVGrid.h"
+#include "readdata.h"
+
+#define sizeError2 "The sizes of the hydrodynamic data files are different"
+#define sizeError "The sizes of the hydrodynamicU or -V files does not correspond with the sizes of the grid file"
+
+using namespace std;
+
+UVGrid::UVGrid() {
+  auto us = readHydrodynamicU();
+  auto vs = readHydrodynamicV();
+  if (us.size() != vs.size()) {
+    throw domain_error(sizeError2);
+  }
+
+  tie(times, lats, lons) = readGrid();
+
+  timeSize = times.size();
+  latSize = lats.size();
+  lonSize = lons.size();
+
+  size_t gridSize = timeSize * latSize * lonSize;
+  if (gridSize != us.size()) {
+    throw domain_error(sizeError);
+  }
+
+  uvData.reserve(gridSize);
+
+  for (auto vel: views::zip(us, vs)) {
+    uvData.push_back(Vel(vel));
+  }
+}
+
+const Vel &UVGrid::operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const {
+  if (timeIndex < 0 or timeIndex >= timeSize
+      or latIndex < 0 or latIndex >= latSize
+      or lonIndex < 0 or lonIndex >= lonSize) {
+    throw std::out_of_range("Index out of bounds");
+  }
+  size_t index = timeIndex * (latSize * lonSize) + latIndex * lonSize + lonIndex;
+  return uvData[index];
+}
+
+double UVGrid::lonStep() const {
+  return lons[1] - lons[0];
+}
+
+double UVGrid::latStep() const {
+  return lats[1] - lats[0];
+}
+
+int UVGrid::timeStep() const {
+  return times[1] - times[0];
+}
+
+void UVGrid::streamSlice(ostream &os, size_t t) {
+  for (int x = 0; x < latSize; x++) {
+    for (int y = 0; y < lonSize; y++) {
+      auto vel = (*this)[t, x, y];
+      os << vel << " ";
+    }
+    os << endl;
+  }
+}

particle-track-and-trace/src/advection/UVGrid.h

@@ -0,0 +1,65 @@
+#ifndef UVGRID_H
+#define UVGRID_H
+
+#include <vector>
+#include "Vel.h"
+
+class UVGrid {
+private:
+    /**
+     * 1D data vector of all the us and vs
+     */
+    std::vector<Vel> uvData;
+public:
+    UVGrid();
+
+    /**
+     * The matrix has shape (timeSize, latSize, lonSize)
+     */
+    size_t timeSize;
+    size_t latSize;
+    size_t lonSize;
+
+    /**
+     * Assuming grid is a regular grid, gives the longitudinal spacing of grid.
+     * @return longitudinal spacing
+     */
+    double lonStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the latitudinal spacing of grid.
+     * @return latitudinal spacing
+     */
+    double latStep() const;
+
+    /**
+     * Assuming grid is a regular grid, gives the time spacing of grid.
+     * @return time spacing
+     */
+    int timeStep() const;
+
+    /**
+     * times, lats, lons are vector of length timeSize, latSize, lonSize respectively.
+     * The maintain the following invariant:
+     * grid[timeIndex,latIndex,lonIndex] gives the u,v at the point with latitude at lats[latIndex],
+     * with longitude at lons[lonIndex], and with time at times[timeIndex].
+     */
+    std::vector<int> times;
+    std::vector<double> lats;
+    std::vector<double> lons;
+
+    /**
+     * The 3D index into the data. The array is sized by [8761][67][116]
+     * @return Velocity at that index
+     */
+    const Vel& operator[](size_t timeIndex, size_t latIndex, size_t lonIndex) const;
+
+    /**
+     * Streams a slice at timeIndex t of the matrix to the outstream given by os
+     * @param os outstream
+     * @param t index with which to slice matrix
+     */
+    void streamSlice(std::ostream &os, size_t t);
+};
+
+#endif //UVGRID_H

particle-track-and-trace/src/advection/Vel.cpp

@@ -0,0 +1,40 @@
+#include "Vel.h"
+#include <stdexcept>
+#include <iomanip>
+
+using namespace std;
+
+Vel::Vel(double u, double v) : u(u), v(v) {}
+
+Vel::Vel(const std::pair<double, double>& p) : u(p.first), v(p.second) {}
+
+Vel& Vel::operator=(const std::pair<double, double>& p) {
+    u = p.first;
+    v = p.second;
+    return *this;
+}
+
+Vel Vel::operator+(const Vel& other) const {
+    return Vel(u + other.u, v + other.v);
+}
+
+Vel& Vel::operator+=(const Vel& other) {
+    u += other.u;
+    v += other.v;
+    return *this;
+}
+
+template<typename Scalar>
+Vel Vel::operator/(Scalar scalar) const {
+    if (scalar == 0) throw std::runtime_error("Division by zero");
+    return Vel(u / scalar, v / scalar);
+}
+
+std::ostream& operator<<(ostream& os, const Vel& vel) {
+    os << "(";
+    os << fixed << setprecision(2) << setw(5) << vel.u;
+    os << ", ";
+    os << fixed << setprecision(2) << setw(5) << vel.v;
+    os << ")";
+    return os;
+}

particle-track-and-trace/src/advection/Vel.h

@@ -0,0 +1,44 @@
+#ifndef VEL_H
+#define VEL_H
+
+#include <utility>
+#include <stdexcept>
+#include <iostream>
+#include <format>
+
+class Vel {
+public:
+    double u; // Eastward Current Velocity in the Water Column
+    double v; // Northward Current Velocity in the Water Column
+
+    Vel(double u, double v);
+    Vel(const std::pair<double, double>& p);  // Conversion constructor
+    Vel& operator=(const std::pair<double, double>& p);
+
+    // Operator + to add two Vel objects
+    Vel operator+(const Vel& other) const;
+
+    // Operator += to add another Vel object to this object
+    Vel& operator+=(const Vel& other);
+
+    // Operator * to multiply Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator*(Scalar scalar) const {
+        return Vel(u * scalar, v * scalar);
+    }
+
+    // Operator / to divide Vel by a scalar, defined as a member template
+    template<typename Scalar>
+    Vel operator/(Scalar scalar) const;
+
+    // Friend declaration for the stream insertion operator
+    friend std::ostream& operator<<(std::ostream& os, const Vel& vel);
+};
+
+// Non-member function for scalar multiplication on the left
+template<typename Scalar>
+Vel operator*(Scalar scalar, const Vel& p) {
+    return Vel(p.u * scalar, p.v * scalar);
+}
+
+#endif //VEL_H

particle-track-and-trace/src/advection/interpolate.cpp

@@ -0,0 +1,47 @@
+#include "interpolate.h"
+
+using namespace std;
+
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon) {
+    double latStep = uvGrid.latStep();
+    double lonStep = uvGrid.lonStep();
+    int timeStep = uvGrid.timeStep();
+
+    int latIndex = (lat - uvGrid.lats[0]) / latStep;
+    int lonIndex = (lon - uvGrid.lons[0]) / lonStep;
+    int timeIndex = (time - uvGrid.times[0]) / timeStep;
+
+    double timeRatio = (static_cast<double>(time) - uvGrid.times[timeIndex]) / timeStep;
+    double latRatio = (lat - uvGrid.lats[latIndex]) / latStep;
+    double lonRatio = (lon - uvGrid.lons[lonIndex]) / lonStep;
+
+    Vel point = {0, 0};
+    for (int timeOffset = 0; timeOffset <= 1; timeOffset++) {
+        for (int latOffset = 0; latOffset <= 1; latOffset++) {
+            for (int lonOffset = 0; lonOffset <= 1; lonOffset++) {
+                auto vertex = uvGrid[
+                        timeIndex + 1 < uvGrid.timeSize ? timeIndex + timeOffset : timeIndex,
+                                latIndex + 1 < uvGrid.latSize ? latIndex + latOffset : latIndex,
+                                lonIndex + 1 < uvGrid.lonSize ? lonIndex + lonOffset : lonIndex
+                ];
+
+                double timeRation = (1 - timeOffset) * (1 - timeRatio) + timeOffset * timeRatio;
+                double latRation = (1 - latOffset) * (1 - latRatio) + latOffset * latRatio;
+                double lonRation = (1 - lonOffset) * (1 - lonRatio) + lonOffset * lonRatio;
+                point += timeRation * latRation * lonRation * vertex;
+            }
+        }
+    }
+
+    return point;
+}
+
+vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, vector<tuple<int, double, double>> points) {
+    vector<Vel> result;
+    result.reserve(points.size());
+    for (auto [time, lat, lon]: points) {
+        result.push_back(bilinearinterpolate(uvGrid, time, lat, lon));
+    }
+
+    return result;
+}

particle-track-and-trace/src/advection/interpolate.h

@@ -0,0 +1,28 @@
+#ifndef INTERPOLATE_H
+#define INTERPOLATE_H
+
+#include <vector>
+
+#include "UVGrid.h"
+
+/**
+ * Bilinearly interpolate the point (time, lat, lon) to produce the interpolated velocity.
+ * Since it is in 3D, this means that it interpolates against 8 points (excluding edges).
+ * As described in https://numerical.recipes/book.html Chapter 3.6
+ * @param uvGrid velocity grid
+ * @param time time of point
+ * @param lat latitude of point
+ * @param lon longitude of point
+ * @return interpolated velocity
+ */
+Vel bilinearinterpolate(const UVGrid &uvGrid, int time, double lat, double lon);
+
+/**
+ * Helper function for bilnearly interpolating a vector of points
+ * @param uvGrid velocity grid
+ * @param points vector of points
+ * @return interpolated velocities
+ */
+std::vector<Vel> bilinearinterpolation(const UVGrid &uvGrid, std::vector<std::tuple<int, double, double>> points);
+
+#endif //INTERPOLATE_H

particle-track-and-trace/src/advection/readdata.cpp

@@ -0,0 +1,50 @@
+#include <stdexcept>
+
+#include <netcdf>
+
+#include "readdata.h"
+
+using namespace std;
+using namespace netCDF;
+
+template <typename T>
+vector<T> getVarVector(const NcVar &var) {
+    int length = 1;
+    for (NcDim dim : var.getDims()) {
+        length *= dim.getSize();
+    }
+
+    vector<T> vec(length);
+
+    var.getVar(vec.data());
+
+    return vec;
+}
+
+vector<double> readHydrodynamicU() {
+    // Vs and Us flipped cause the files are named incorrectly
+    netCDF::NcFile data("../../../../data/hydrodynamic_V.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("vo")->second);
+}
+
+vector<double> readHydrodynamicV() {
+    // Vs and Us flipped cause the files are named incorrectly
+    netCDF::NcFile data("../../../../data/hydrodynamic_U.h5", netCDF::NcFile::read);
+
+    multimap< string, NcVar > vars = data.getVars();
+
+    return getVarVector<double>(vars.find("uo")->second);
+}
+
+tuple<vector<int>, vector<double>, vector<double>> readGrid() {
+    netCDF::NcFile data("../../../../data/grid.h5", netCDF::NcFile::read);
+    multimap< string, NcVar > vars = data.getVars();
+    vector<int> time = getVarVector<int>(vars.find("times")->second);
+    vector<double> longitude = getVarVector<double>(vars.find("longitude")->second);
+    vector<double> latitude = getVarVector<double>(vars.find("latitude")->second);
+
+    return {time, latitude, longitude};
+}

particle-track-and-trace/src/advection/readdata.h

@@ -0,0 +1,22 @@
+#ifndef READDATA_H
+#define READDATA_H
+
+/**
+ * reads the file hydrodynamic_U.h5
+ * @return the data vector of us
+ */
+std::vector<double> readHydrodynamicU();
+
+/**
+ * reads the file hydrodynamic_V.h5
+ * @return the data vector of vs
+ */
+std::vector<double> readHydrodynamicV();
+
+/**
+ * Reads the file grid.h5
+ * @return a tuple of (times, latitude, longitude)
+ */
+std::tuple<std::vector<int>, std::vector<double>, std::vector<double>> readGrid();
+
+#endif //READDATA_H

particle-track-and-trace/src/commands/SpawnPointCallback.cpp

@@ -0,0 +1,82 @@
+#include "SpawnPointCallback.h"
+
+#include <vtkVertex.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkSmartPointer.h>
+#include <vtkCommand.h>
+#include <vtkRenderWindow.h>
+
+#include "../CartographicTransformation.h"
+
+void convertDisplayToWorld(vtkRenderer *renderer, int x, int y, double *worldPos) {
+  double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
+  renderer->SetDisplayPoint(displayPos);
+  renderer->DisplayToWorld();
+  renderer->GetWorldPoint(worldPos);
+}
+
+void SpawnPointCallback::Execute(vtkObject *caller, unsigned long evId, void *callData) {
+  // Note the use of reinterpret_cast to cast the caller to the expected type.
+  auto interactor = reinterpret_cast<vtkRenderWindowInteractor *>(caller);
+
+  if (evId == vtkCommand::LeftButtonPressEvent) {
+    dragging = true;
+  }
+  if (evId == vtkCommand::LeftButtonReleaseEvent) {
+    dragging = false;
+  }
+  if (!dragging) {
+    return;
+  }
+
+  int x, y;
+  interactor->GetEventPosition(x, y);
+
+  double worldPos[4] = {2, 0, 0, 0};
+  double displayPos[3] = {static_cast<double>(x), static_cast<double>(y), 0.0};
+  ren->SetDisplayPoint(displayPos);
+  ren->DisplayToWorld();
+  ren->GetWorldPoint(worldPos);
+  inverseCartographicProjection->MultiplyPoint(worldPos, worldPos);
+  cout << "clicked on lon = " << worldPos[0] << " and lat = " << worldPos[1] << endl;
+
+  vtkIdType id = points->InsertNextPoint(worldPos[0], worldPos[1], 0);
+  data->SetPoints(points);
+
+  vtkSmartPointer<vtkVertex> vertex = vtkSmartPointer<vtkVertex>::New();
+  vertex->GetPointIds()->SetId(0, id);
+
+  vtkSmartPointer<vtkCellArray> vertices = vtkSmartPointer<vtkCellArray>::New();
+  vertices->InsertNextCell(vertex);
+  data->SetVerts(vertices);
+  ren->GetRenderWindow()->Render();
+}
+
+
+SpawnPointCallback::SpawnPointCallback() : data(nullptr),
+                                           points(nullptr),
+                                           inverseCartographicProjection(nullptr),
+                                           uvGrid(nullptr) { }
+
+SpawnPointCallback *SpawnPointCallback::New() {
+    return new SpawnPointCallback;
+}
+
+void SpawnPointCallback::setData(const vtkSmartPointer<vtkPolyData> &data) {
+  this->data = data;
+}
+
+void SpawnPointCallback::setPoints(const vtkSmartPointer<vtkPoints> &points) {
+  this->points = points;
+}
+
+void SpawnPointCallback::setRen(const vtkSmartPointer<vtkRenderer> &ren) {
+  this->ren = ren;
+}
+
+void SpawnPointCallback::setUVGrid(const std::shared_ptr<UVGrid> &uvGrid) {
+  this->uvGrid = uvGrid;
+  inverseCartographicProjection = getCartographicTransformMatrix(uvGrid);
+  inverseCartographicProjection->Invert();
+}

particle-track-and-trace/src/commands/SpawnPointCallback.h

@@ -0,0 +1,40 @@
+#ifndef SPAWNPOINTCALLBACK_H
+#define SPAWNPOINTCALLBACK_H
+
+
+#include <vtkCallbackCommand.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkMatrix4x4.h>
+#include "../advection/UVGrid.h"
+
+class SpawnPointCallback : public vtkCallbackCommand {
+
+public:
+  static SpawnPointCallback *New();
+
+  SpawnPointCallback();
+
+  void setPoints(const vtkSmartPointer<vtkPoints> &points);
+
+  void setData(const vtkSmartPointer<vtkPolyData> &data);
+
+  void setRen(const vtkSmartPointer<vtkRenderer> &ren);
+
+  void setUVGrid(const std::shared_ptr<UVGrid> &uvGrid);
+
+private:
+  vtkSmartPointer<vtkPolyData> data;
+  vtkSmartPointer<vtkPoints> points;
+  vtkSmartPointer<vtkRenderer> ren;
+  std::shared_ptr<UVGrid> uvGrid;
+  vtkSmartPointer<vtkMatrix4x4> inverseCartographicProjection;
+
+  void Execute(vtkObject *caller, unsigned long evId, void *callData) override;
+
+  bool dragging = false;
+};
+
+
+#endif //SPAWNPOINTCALLBACK_H

particle-track-and-trace/src/commands/TimerCallbackCommand.cpp

@@ -0,0 +1,44 @@
+#include "TimerCallbackCommand.h"
+#include "../Program.h"
+
+
+TimerCallbackCommand::TimerCallbackCommand() : dt(3600), maxTime(3600*24*365), time(0) {}
+
+TimerCallbackCommand* TimerCallbackCommand::New(Program *program) {
+  TimerCallbackCommand *cb = new TimerCallbackCommand();
+  cb->setProgram(program);
+  cb->setPaused(false);
+  return cb;
+}
+
+void TimerCallbackCommand::Execute(vtkObject* caller, unsigned long eventId, void* vtkNotUsed(callData)) {
+  auto intr = reinterpret_cast<vtkRenderWindowInteractor *>(caller);
+  
+  if (eventId == vtkCommand::KeyPressEvent and not strcmp("space", intr->GetKeySym())) {
+    this->paused = ! this->paused;
+  } else if (eventId == vtkCommand::TimerEvent and not this->paused) {
+  this->time += this->dt;
+
+  if (this->time >= this->maxTime) {
+    return;
+    // TODO: how do we deal with reaching the end of the simulated dataset? Do we just stop simulating, loop back around? What about the location of the particles in this case? Just some ideas to consider, but we should iron this out pretty soon.
+  }
+
+  this->program->updateData(this->time);
+  }
+
+}
+
+
+void TimerCallbackCommand::setProgram(Program *program) {
+  this->program = program;
+}
+
+
+void TimerCallbackCommand::setPaused(const bool val) {
+  this->paused = val;
+}
+
+void TimerCallbackCommand::setDt(int dt) {
+  this->dt = dt;
+}

particle-track-and-trace/src/commands/TimerCallbackCommand.h

@@ -0,0 +1,26 @@
+#ifndef TIMERCALLBACKCOMMAND_H
+#define TIMERCALLBACKCOMMAND_H
+
+#include <vtkCallbackCommand.h>
+#include "../Program.h"
+
+class TimerCallbackCommand : public vtkCallbackCommand {
+public:
+  TimerCallbackCommand();
+  static TimerCallbackCommand* New(Program *program);
+  void Execute(vtkObject* caller, unsigned long eventId, void* vtkNotUsed(callData)) override;
+
+  void setProgram(Program *program);
+  void setPaused(const bool val);
+
+  void setDt(int dt);
+
+private:
+  int time;
+  int dt;
+  int maxTime;
+  bool paused;
+  Program *program;
+};
+
+#endif

particle-track-and-trace/src/layers/BackgroundImage.cpp

@@ -0,0 +1,62 @@
+#include "BackgroundImage.h"
+#include <vtkCamera.h>
+#include <vtkImageActor.h>
+#include <vtkImageData.h>
+#include <vtkImageReader2.h>
+
+using std::string;
+
+BackgroundImage::BackgroundImage(string imagePath) : imagePath(imagePath) {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(0);
+  this->ren->InteractiveOff();
+  updateImage();
+}
+
+
+void BackgroundImage::updateImage() {
+  vtkSmartPointer<vtkImageData> imageData;
+
+  vtkSmartPointer<vtkImageReader2> imageReader;
+  
+  imageReader.TakeReference(this->readerFactory->CreateImageReader2(this->imagePath.c_str()));
+  imageReader->SetFileName(this->imagePath.c_str());
+  imageReader->Update();
+  imageData = imageReader->GetOutput();
+
+  vtkNew<vtkImageActor> imageActor;
+  imageActor->SetInputData(imageData);
+
+  this->ren->AddActor(imageActor);
+
+
+  // camera stuff
+  // essentially sets the camera to the middle of the background, and points it at the background
+  // TODO: extract this to its own function, separate from the background class.
+  double origin[3], spacing[3];
+  int extent[6];
+  imageData->GetOrigin(origin);
+  imageData->GetSpacing(spacing);
+  imageData->GetExtent(extent);
+
+  vtkCamera *camera = this->ren->GetActiveCamera();
+  camera->ParallelProjectionOn();
+
+  double xc = origin[0] + 0.5 * (extent[0] + extent[1]) * spacing[0];
+  double yc = origin[1] + 0.5 * (extent[2] + extent[3]) * spacing[1];
+  double yd = (extent[3] - extent[2] + 1) * spacing[1];
+  double d = camera->GetDistance();
+  camera->SetParallelScale(0.5 * yd);
+  camera->SetFocalPoint(xc, yc, 0.0);
+  camera->SetPosition(xc, yc, d);
+}
+
+
+string BackgroundImage::getImagePath() {
+  return this->imagePath;
+}
+
+void BackgroundImage::setImagePath(string imagePath) {
+  this->imagePath = imagePath;
+  updateImage();
+}

particle-track-and-trace/src/layers/BackgroundImage.h

@@ -0,0 +1,37 @@
+#ifndef BACKGROUND_H
+#define BACKGROUND_H
+
+#include "Layer.h"
+#include <vtkImageReader2Factory.h>
+
+/** Implements the Layer class for the case of a background image.
+  * Specifically, reads a backgroundImage given by the imagePath attribute and puts it on layer 0.
+  */
+class BackgroundImage : public Layer {
+private:
+  std::string imagePath;
+  vtkSmartPointer<vtkImageReader2Factory> readerFactory;
+
+  /** This private function updates the background image using the imagePath attribute.
+   */
+  void updateImage();
+
+
+public:
+  /** Constructor.
+    * @param imagePath : String to the path of the image to use as background.
+    */
+  BackgroundImage(std::string imagePath);
+
+  /** Getter.
+    * @return the imagePath attribute.
+    */
+  std::string getImagePath();
+
+  /** Setter. Can be used to change the background image 
+    * @param imagePath : String to the path of the new image to use.
+    */
+  void setImagePath(std::string imagePath);
+};
+
+#endif

particle-track-and-trace/src/layers/EGlyphLayer.cpp

@@ -0,0 +1,109 @@
+#include "EGlyphLayer.h"
+#include <vtkPointData.h>
+#include <vtkDoubleArray.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkRegularPolygonSource.h>
+#include <vtkGlyph2D.h>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkProperty.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <vtkArrowSource.h>
+#include "../CartographicTransformation.h"
+#include "../advection/readdata.h"
+#include "../advection/interpolate.h"
+
+using namespace std;
+
+EGlyphLayer::EGlyphLayer(std::shared_ptr<UVGrid> uvGrid) {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(1);
+  this->ren->InteractiveOff();
+
+  this->uvGrid = uvGrid;
+
+  this->data = vtkSmartPointer<vtkPolyData>::New();
+  this->direction = vtkSmartPointer<vtkDoubleArray>::New();
+  this->direction->SetName("direction");
+  readCoordinates();
+}
+
+
+void EGlyphLayer::readCoordinates() {
+  vtkNew<vtkPoints> points;
+  this->numLats = uvGrid->lats.size();
+  this->numLons = uvGrid->lons.size();
+
+  this->direction->SetNumberOfComponents(3);
+  this->direction->SetNumberOfTuples(numLats * numLons);
+  points->Allocate(numLats * numLons);
+
+  auto camera = createNormalisedCamera();
+  ren->SetActiveCamera(camera);
+
+  int i = 0;
+  int latIndex = 0;
+  for (double lat: uvGrid->lats) {
+    int lonIndex = 0;
+    for (double lon: uvGrid->lons) {
+      auto [u, v] = (*uvGrid)[0, latIndex, lonIndex];
+      direction->SetTuple3(i, 5*u, 5*v, 0);
+      points->InsertPoint(i++, lon, lat, 0);
+      lonIndex++;
+    }
+    latIndex++;
+  }
+  this->data->SetPoints(points);
+  this->data->GetPointData()->AddArray(this->direction);
+  this->data->GetPointData()->SetActiveVectors("direction");
+
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
+  transformFilter->SetInputData(data);
+
+  vtkNew<vtkGlyphSource2D> arrowSource;
+  arrowSource->SetGlyphTypeToArrow();
+  arrowSource->SetScale(0.2); //TODO: set this properly
+  arrowSource->Update();
+
+  vtkNew<vtkGlyph2D> glyph2D;
+  glyph2D->SetSourceConnection(arrowSource->GetOutputPort());
+  glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+  glyph2D->OrientOn();
+  glyph2D->ClampingOn();
+  glyph2D->SetScaleModeToScaleByVector();
+  glyph2D->SetVectorModeToUseVector();
+  glyph2D->Update();
+
+//  vtkNew<vtkCoordinate> coordinate;
+//  coordinate->SetCoordinateSystemToWorld();
+
+  vtkNew<vtkPolyDataMapper>(mapper);
+  // mapper->SetTransformCoordinate(coordinate);
+  mapper->SetInputConnection(glyph2D->GetOutputPort());
+  mapper->Update();
+
+  vtkNew<vtkActor> actor;
+  actor->SetMapper(mapper);
+
+  actor->GetProperty()->SetColor(0, 0, 0);
+  actor->GetProperty()->SetOpacity(0.2);
+
+  this->ren->AddActor(actor);
+}
+
+
+void EGlyphLayer::updateData(int t) {
+  int i = 0;
+  for (int lat = 0; lat < uvGrid->latSize; lat++) {
+    for (int lon = 0; lon < uvGrid->lonSize; lon++) {
+      auto [u, v] = (*uvGrid)[t/3600, lat, lon];
+      // TODO: 5*v scaling stuff should really be a filter transform
+      this->direction->SetTuple3(i, 5*u, 5*v, 0);
+      i++;
+    }
+  }
+  this->direction->Modified();
+}

particle-track-and-trace/src/layers/EGlyphLayer.h

@@ -0,0 +1,38 @@
+#ifndef EGLYPHLAYER_H
+#define EGLYPHLAYER_H
+
+#include "Layer.h"
+#include <vtkPolyData.h>
+
+#include "../advection/UVGrid.h"
+
+/** Implements the Layer class for the case of a Eulerian visualization.
+  * Specifically, this class models the eulerian flow-fields of the simulation using the 'glyph' mark and 'direction' and 'form' channels to denote direction and strength of velocities.
+  */
+class EGlyphLayer : public Layer {
+private:
+  vtkSmartPointer<vtkPolyData> data;
+  vtkSmartPointer<vtkDoubleArray> direction;
+  std::shared_ptr<UVGrid> uvGrid;
+  int numLats;
+  int numLons;
+
+  /** This private function sets up the initial coordinates for the glyphs in the dataset.
+    * It also reads some initial data to actually display.
+    */
+  void readCoordinates();
+
+public:
+  /** Constructor.
+    */
+  EGlyphLayer(std::shared_ptr<UVGrid> uvGrid);
+
+  /** updates the glyphs to reflect the given timestamp in the dataset.
+    * @param t : the time at which to fetch the data.
+    */
+  void updateData(int t);
+
+};
+
+
+#endif

particle-track-and-trace/src/layers/LGlyphLayer.cpp

@@ -0,0 +1,101 @@
+#include "LGlyphLayer.h"
+#include "../commands/SpawnPointCallback.h"
+#include <vtkActor2D.h>
+#include <vtkGlyph2D.h>
+#include <vtkGlyphSource2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkProperty2D.h>
+#include <vtkVertexGlyphFilter.h>
+#include <vtkInteractorStyle.h>
+#include <vtkInteractorStyleUser.h>
+#include <vtkTransform.h>
+#include <vtkTransformFilter.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkRenderWindow.h>
+#include <vtkCamera.h>
+
+#include "../CartographicTransformation.h"
+
+vtkSmartPointer<SpawnPointCallback> LGlyphLayer::createSpawnPointCallback() {
+  auto newPointCallBack = vtkSmartPointer<SpawnPointCallback>::New();
+  newPointCallBack->setData(data);
+  newPointCallBack->setPoints(points);
+  newPointCallBack->setRen(ren);
+  newPointCallBack->setUVGrid(uvGrid);
+  return newPointCallBack;
+}
+
+// Further notes; current thinking is to allow tracking a particle's age by using a scalar array in the VtkPolyData. This would be incremented for every tick/updateData function call.
+// Another challenge is the concept of beaching; dead particles must not be included in the advect function call (wasted computations), but they should not be outright deleted from the vtkPoints either (we still want to display them). Working Solution: have another array of ints in the vtkPolyData, which tracks for how many calls of UpdateData a given particle has not had its position changed. If this int reaches some treshold (5? 10? 3? needs some testing), exclude the particle from the advect call.
+// TODO: modelling all this in vtkClasses is workable, but ideally i would want to work with a native C++ class. See if this is doable and feasible.
+
+LGlyphLayer::LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel) {
+  this->ren = vtkSmartPointer<vtkRenderer>::New();
+  this->ren->SetLayer(2);
+
+  this->points = vtkSmartPointer<vtkPoints>::New();
+  this->data = vtkSmartPointer<vtkPolyData>::New();
+  this->data->SetPoints(this->points);
+
+  advector = std::move(advectionKernel);
+  this->uvGrid = uvGrid;
+
+  auto camera = createNormalisedCamera();
+  ren->SetActiveCamera(camera);
+
+  vtkSmartPointer<vtkTransformFilter> transformFilter = createCartographicTransformFilter(uvGrid);
+  transformFilter->SetInputData(data);
+
+  vtkNew<vtkGlyphSource2D> circleSource;
+  circleSource->SetGlyphTypeToCircle();
+  circleSource->SetScale(0.05);
+  circleSource->Update();
+
+  vtkNew<vtkGlyph2D> glyph2D;
+  glyph2D->SetSourceConnection(circleSource->GetOutputPort());
+  glyph2D->SetInputConnection(transformFilter->GetOutputPort());
+  glyph2D->SetColorModeToColorByScalar();
+  glyph2D->Update();
+
+  vtkNew<vtkPolyDataMapper> mapper;
+  mapper->SetInputConnection(glyph2D->GetOutputPort());
+  mapper->Update();
+
+  vtkNew<vtkActor> actor;
+  actor->SetMapper(mapper);
+
+  this->ren->AddActor(actor);
+}
+
+// creates a few points so we can test the updateData function
+void LGlyphLayer::spoofPoints() {
+  this->points->InsertNextPoint(-4.125, 61.375, 0);
+  this->points->InsertNextPoint(6.532949683882039, 53.24308582564463, 0); // Coordinates of Zernike
+  this->points->InsertNextPoint(5.315307819255385, 60.40001057122271, 0); // Coordinates of Bergen
+  this->points->InsertNextPoint(6.646210231365825, 46.52346296009023, 0); // Coordinates of Lausanne
+  this->points->InsertNextPoint(-6.553894313570932, 62.39522131195857,0); // Coordinates of the top of the Faroe islands
+
+  this->points->Modified();
+}
+
+void LGlyphLayer::updateData(int t) {
+  const int SUPERSAMPLINGRATE = 4;
+  double point[3];
+  for (vtkIdType n = 0; n < this->points->GetNumberOfPoints(); n++) {
+    this->points->GetPoint(n, point);
+    for (int i = 0; i < SUPERSAMPLINGRATE; i++) {
+      std::tie(point[1], point[0]) = advector->advect(t, point[1], point[0], (t-lastT)/SUPERSAMPLINGRATE);
+    }
+    this->points->SetPoint(n, point[0], point[1], 0);
+  }
+  lastT = t;
+  this->points->Modified();
+}
+
+void LGlyphLayer::addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) {
+  auto newPointCallBack = createSpawnPointCallback();
+  interactor->AddObserver(vtkCommand::LeftButtonPressEvent, newPointCallBack);
+  interactor->AddObserver(vtkCommand::LeftButtonReleaseEvent, newPointCallBack);
+  interactor->AddObserver(vtkCommand::MouseMoveEvent, newPointCallBack);
+}

particle-track-and-trace/src/layers/LGlyphLayer.h

@@ -0,0 +1,40 @@
+#ifndef LGLYPHLAYER_H
+#define LGLYPHLAYER_H
+
+#include "Layer.h"
+#include "../advection/AdvectionKernel.h"
+#include "../commands/SpawnPointCallback.h"
+#include <vtkPolyData.h>
+#include <vtkInteractorStyle.h>
+
+/** Implements the Layer class for the case of a Lagrangian visualization.
+  * Specifically, this class models the Lagrangian particles in the simulation using the 'glyph' mark and 'transparency' channel to denote age.
+  */
+class LGlyphLayer : public Layer {
+private:
+  vtkSmartPointer<vtkPoints> points;
+  vtkSmartPointer<vtkPolyData> data;
+  std::unique_ptr<AdvectionKernel> advector;
+  std::shared_ptr<UVGrid> uvGrid;
+  int lastT = 1000;
+
+public:
+  /** Constructor.
+    */
+  LGlyphLayer(std::shared_ptr<UVGrid> uvGrid, std::unique_ptr<AdvectionKernel> advectionKernel);
+
+  /** This function spoofs a few points in the dataset. Mostly used for testing.
+    */
+  void spoofPoints();
+
+  /** updates the glyphs to reflect the given timestamp in the dataset.
+    * @param t : the time at which to fetch the data.
+    */
+  void updateData(int t) override;
+
+  vtkSmartPointer<SpawnPointCallback> createSpawnPointCallback();
+
+  void addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) override;
+};
+
+#endif

particle-track-and-trace/src/layers/Layer.cpp

@@ -0,0 +1,17 @@
+#include "Layer.h"
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+
+using std::string;
+
+vtkSmartPointer<vtkRenderer> Layer::getLayer() {
+  return this->ren;
+}
+
+void Layer::updateData(int t) {
+  // By default, do nothing
+}
+
+void Layer::addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor) {
+    // By default, do nothing
+}

particle-track-and-trace/src/layers/Layer.h

@@ -0,0 +1,32 @@
+#ifndef LAYER_H
+#define LAYER_H
+
+#include <vtkInteractorStyle.h>
+#include <vtkRenderer.h>
+
+/** This class represents one abstract layer to be rendered to VTK.
+  * It exists to manage multiple different layers under the Program class.
+  */
+class Layer {
+protected:
+  vtkSmartPointer<vtkRenderer> ren;
+
+public:
+  /** gets the vtkRenderer to assign it to the vtkRenderWindow of the program class.
+   * @return pointer to the vtkRenderer of this class.
+   */
+  virtual vtkSmartPointer<vtkRenderer> getLayer();
+
+
+  /** updates the data in the layer to reflect the given timestamp.
+    * @param t : the timestamp which the data should reflect.
+    */
+  virtual void updateData(int t);
+
+  /** Adds observers to the renderWindowinteractor within which this layer is active.
+    * @param interactor : pointer to the interactor that observers can be added to.
+    */
+  virtual void addObservers(vtkSmartPointer<vtkRenderWindowInteractor> interactor);
+};
+
+#endif

particle-track-and-trace/src/main.cpp

@@ -0,0 +1,38 @@
+#include <netcdf>
+#include <vtkActor2D.h>
+#include <vtkNamedColors.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataMapper2D.h>
+#include <vtkProperty2D.h>
+#include <vtkRenderer.h>
+#include <vtkVertexGlyphFilter.h>
+#include <memory>
+
+#include "layers/BackgroundImage.h"
+#include "layers/EGlyphLayer.h"
+#include "layers/LGlyphLayer.h"
+#include "Program.h"
+#include "advection/UVGrid.h"
+#include "advection/RK4AdvectionKernel.h"
+
+using namespace std;
+
+#define DT 60 * 60 // 60 sec/min * 60 mins
+
+int main() {
+  shared_ptr<UVGrid> uvGrid = std::make_shared<UVGrid>();
+  auto kernelRK4 = make_unique<RK4AdvectionKernel>(uvGrid);
+
+  auto l = new LGlyphLayer(uvGrid, std::move(kernelRK4));
+
+  Program *program = new Program(DT);
+  program->addLayer(new BackgroundImage("../../../../data/map_661-661.png"));
+  program->addLayer(new EGlyphLayer(uvGrid));
+  program->addLayer(l);
+
+  program->render();
+
+  return EXIT_SUCCESS;
+}
+

particle-track-and-trace/src/script.py

@@ -0,0 +1,226 @@
+#!/usr/bin/env python
+
+import collections
+import json
+import os
+import re
+from pathlib import Path
+
+
+def get_program_parameters(argv):
+    import argparse
+    description = 'Generate a find_package(VTK COMPONENTS ...) command for CMake.'
+    epilogue = '''
+Uses modules.json and your source files to generate a
+  find_package(VTK COMPONENTS ...) command listing all the vtk modules
+  needed by the C++ source and header files in your code.
+
+Paths for more than one source path can be specified.
+
+Note than if there are spaces in the paths, enclose the path in quotes.
+
+If it is unable to find modules for your headers then
+  a list of these, along with the files they are in, is produced
+  so you can manually add the corresponding modules or rebuild VTK
+  to include the missing modules.
+
+You will need to manually add any third-party modules
+   (if used) to the find_package command.
+    '''
+    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
+                                     formatter_class=argparse.RawTextHelpFormatter)
+    parser.add_argument('json', default=['modules.json'], help='The path to the VTK JSON file (modules.json).')
+    parser.add_argument('sources', nargs='+', help='The path to the source files.')
+    parser.add_argument('-f', '--file', help='The file name to write the output too.')
+    args = parser.parse_args()
+    return args.json, args.sources, args.file
+
+
+class Patterns:
+    header_pattern = re.compile(r'^#include *[<\"](\S+)[>\"]')
+    vtk_include_pattern = re.compile(r'^(vtk\S+)')
+    vtk_qt_include_pattern = re.compile(r'^(QVTK\S+)')
+
+
+def get_headers_modules(json_data):
+    """
+    From the parsed JSON data file make a dictionary whose key is the
+     header filename and value is the module.
+    :param json_data: The parsed JSON file modules.json.
+    :return:
+    """
+
+    # The headers should be unique to a module, however we will not assume this.
+    res = collections.defaultdict(set)
+    for k, v in json_data['modules'].items():
+        if 'headers' in v:
+            for k1 in v['headers']:
+                res[k1].add(k)
+    return res
+
+
+def get_vtk_components(jpath, paths):
+    """
+    Get the VTK components
+    :param jpath: The path to the JSON file.
+    :param paths: The C++ file paths.
+    :return:
+    """
+
+    with open(jpath) as data_file:
+        json_data = json.load(data_file)
+    vtk_headers_modules = get_headers_modules(json_data)
+
+    modules = set()
+    inc_no_mod = set()
+    inc_no_mod_headers = collections.defaultdict(set)
+    mod_implements = collections.defaultdict(set)
+    headers = collections.defaultdict(set)
+
+    for path in paths:
+        if path.is_file():
+            content = path.read_text().split('\n')
+            for line in content:
+                m = Patterns.header_pattern.match(line.strip())
+                if m:
+                    # We have a header name, split it from its path (if the path exists).
+                    header_parts = os.path.split(m.group(1))
+                    m = Patterns.vtk_include_pattern.match(header_parts[1])
+                    if m:
+                        headers[m.group(1)].add(path)
+                        continue
+                    m = Patterns.vtk_qt_include_pattern.match(header_parts[1])
+                    if m:
+                        headers[m.group(1)].add(path)
+    for incl in headers:
+        if incl in vtk_headers_modules:
+            m = vtk_headers_modules[incl]
+            for v in m:
+                modules.add(v)
+        else:
+            inc_no_mod.add(incl)
+            inc_no_mod_headers[incl] = headers[incl]
+
+    if headers:
+        for m in modules:
+            if not json_data['modules'][m]['implementable']:
+                continue
+            for i in json_data['modules']:
+                if i in modules:
+                    continue
+                if m in json_data['modules'][i]['implements']:
+                    # Suggest module i since it implements m
+                    mod_implements[i].add(m)
+
+    return modules, mod_implements, inc_no_mod, inc_no_mod_headers
+
+
+def disp_components(modules, module_implements):
+    """
+    For the found modules display them in a form that the user can
+     copy/paste into their CMakeLists.txt file.
+    :param modules: The modules.
+    :param module_implements: Modules implementing other modules.
+    :return:
+    """
+    res = ['find_package(VTK\n COMPONENTS']
+    for m in sorted(modules):
+        res.append('    {:s}'.format(m.split('::')[1]))
+    if module_implements:
+        keys = sorted(module_implements)
+        max_width = len(max(keys, key=len).split('::')[1])
+        comments = [
+            '    #',
+            '    # These modules are suggested since they implement an existing module.',
+            '    # You may need to uncomment one or more of these.',
+            '    # If vtkRenderWindow is used and you want to use OpenGL,',
+            '    #   you also need the RenderingOpenGL2 module.',
+            '    # If vtkRenderWindowInteractor is used,',
+            '    #    uncomment RenderingUI and possibly InteractionStyle.',
+            '    # If text rendering is used, uncomment RenderingFreeType',
+            '    #'
+        ]
+        res.extend(comments)
+        for key in keys:
+            res.append(
+                f'    # {key.split("::")[1]:<{max_width}} # implements {", ".join(sorted(module_implements[key]))}')
+    res.append(')\n')
+
+    return res
+
+
+def disp_missing_components(inc_no_mod, inc_no_mod_headers):
+    """
+    Display the headers along with the missing VTK modules.
+
+    :param inc_no_mod: Missing modules.
+    :param inc_no_mod_headers: Headers with missing modules.
+    :return:
+    """
+    if inc_no_mod:
+        res = [''
+               '*' * 64,
+               'You will need to manually add the modules that',
+               '  use these headers to the find_package command.',
+               'These could be external modules not in the modules.json file.',
+               'Or you may need to rebuild VTK to include the missing modules.',
+               '',
+               'Here are the vtk headers and corresponding files:']
+        sinmd = sorted(inc_no_mod)
+        for i in sinmd:
+            sh = sorted(list(inc_no_mod_headers[i]))
+            res.append(f'in {i}:')
+            for j in sh:
+                res.append(f'   {j}')
+        res.append('*' * 64)
+
+        return res
+    else:
+        return None
+
+
+def main(json_path, src_paths, ofn):
+    jpath = Path(json_path)
+    if jpath.is_dir():
+        jpath = jpath / 'modules.json'
+    if not jpath.is_file():
+        print(f'Non existent JSON file: {jpath}')
+        return
+
+    paths = list()
+    valid_ext = ['.h', '.hxx', '.cxx', '.cpp', '.txx']
+    path_list = list()
+    for fn in src_paths:
+        path = Path(fn)
+        if path.is_file() and path.suffix in valid_ext:
+            paths.append(path)
+        elif path.is_dir():
+            for e in valid_ext:
+                path_list += list(Path(fn).rglob(f'*{e}'))
+            program_path = Path(__file__)
+            for path in path_list:
+                if path.resolve() != program_path.resolve():
+                    paths.append(path)
+        else:
+            print(f'Non existent path: {path}')
+
+    modules, mod_implements, inc_no_mod, inc_no_mod_headers = get_vtk_components(jpath, paths)
+
+    res = '\n'.join(disp_components(modules, mod_implements))
+    if inc_no_mod:
+        res += '\n'.join(disp_missing_components(inc_no_mod, inc_no_mod_headers))
+
+    if ofn:
+        path = Path(ofn)
+        if path.suffix == '':
+            path = Path(ofn).with_suffix('.txt')
+        path.write_text(res)
+    else:
+        print(res)
+
+
+if __name__ == '__main__':
+    import sys
+
+    json_paths, src_paths, ofn = get_program_parameters(sys.argv)
+    main(json_paths, src_paths, ofn)