{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Detailed Look\n",
    "\n",
    "\n",
    "A more detailed look at spatial interpolation, integration through time, and I/O.\n",
    "For additional documentation e.g. see\n",
    "[1](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),\n",
    "[2](https://JuliaClimate.github.io/MeshArrays.jl/dev/),\n",
    "[3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),\n",
    "[4](https://en.wikipedia.org/wiki/Displacement_(vector)).\n",
    "Here we illustrate a few things in more detail:\n",
    "\n",
    "- reading velocities from file.\n",
    "  - gridded velocity output (U*data, V*data)\n",
    "  - pre-computed trajectory output (`float_traj*data`)\n",
    "- interpolating `U,V` from gridded output to individual locations\n",
    "  - compared with `u,v` from `float_traj*data`\n",
    "- computing trajectories (location v time) using `OrdinaryDiffEq.jl`\n",
    "  - compared with `x(t),y(t)` from `float_traj*data`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 1. Import Software"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using IndividualDisplacements, MITgcm\n",
    "import IndividualDisplacements.OrdinaryDiffEq as OrdinaryDiffEq\n",
    "import IndividualDisplacements.DataFrames as DataFrames\n",
    "p=dirname(pathof(IndividualDisplacements))\n",
    "include(joinpath(p,\"../examples/more/example123.jl\"))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 2. Read Trajectory Output\n",
    "\n",
    "from `MITgcm/pkg/flt`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "IndividualDisplacements.flt_example_download()\n",
    "dirIn=IndividualDisplacements.flt_example_path\n",
    "prec=Float32\n",
    "df=read_flt(dirIn,prec);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 3. Read Gridded Variables\n",
    "\n",
    "using `MeshArrays.jl` and e.g. a NamedTyple"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "𝑃,Γ=example2_setup();"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 4. Visualize Velocity Fields\n",
    "\n",
    "```\n",
    "plt=heatmap(Γ.mskW[1,1].*𝑃.u0,title=\"U at the start\")\n",
    "plt=heatmap(Γ.mskW[1,1].*𝑃.u1-𝑃.u0,title=\"U end - U start\")\n",
    "```"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 5. Visualize Trajectories\n",
    "\n",
    "(select one trajectory)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmp=df[df.ID .== 200, :]\n",
    "tmp[1:4,:]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Super-impose trajectory over velocity field (first for u ...)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "x=Γ.XG.f[1][:,1]\n",
    "y=Γ.YC.f[1][1,:]\n",
    "z=transpose(Γ.mskW[1].*𝑃.u0);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Super-impose trajectory over velocity field (... then for v)"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "x=Γ.XC.f[1][:,1]\n",
    "y=Γ.YG.f[1][1,:]\n",
    "z=transpose(Γ.mskW[1].*𝑃.v0);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 6. Interpolate Velocities"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dx=Γ.dx\n",
    "uInit=[tmp[1,:lon];tmp[1,:lat]]./dx\n",
    "nSteps=Int32(tmp[end,:time]/3600)-2\n",
    "du=fill(0.0,2);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Visualize and compare with actual grid point values -- jumps on the tangential component are expected with linear scheme:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmpu=fill(0.0,100)\n",
    "tmpv=fill(0.0,100)\n",
    "tmpx=fill(0.0,100)\n",
    "for i=1:100\n",
    "    tmpx[i]=500.0 *i./dx\n",
    "    dxdt!(du,[tmpx[i];0.499./dx],𝑃,0.0)\n",
    "    tmpu[i]=du[1]\n",
    "    tmpv[i]=du[2]\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "And similarly in the other direction"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tmpu=fill(0.0,100)\n",
    "tmpv=fill(0.0,100)\n",
    "tmpy=fill(0.0,100)\n",
    "for i=1:100\n",
    "    tmpy[i]=500.0 *i./dx\n",
    "    dxdt!(du,[0.499./dx;tmpy[i]],𝑃,0.0)\n",
    "    tmpu[i]=du[1]\n",
    "    tmpv[i]=du[2]\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compare recomputed velocities with those from `pkg/flt`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nSteps=2998\n",
    "tmpu=fill(0.0,nSteps); tmpv=fill(0.0,nSteps);\n",
    "tmpx=fill(0.0,nSteps); tmpy=fill(0.0,nSteps);\n",
    "refu=fill(0.0,nSteps); refv=fill(0.0,nSteps);\n",
    "for i=1:nSteps\n",
    "    dxy_dt_replay(du,[tmp[i,:lon],tmp[i,:lat]],tmp,tmp[i,:time])\n",
    "    refu[i]=du[1]./dx\n",
    "    refv[i]=du[2]./dx\n",
    "    dxdt!(du,[tmp[i,:lon],tmp[i,:lat]]./dx,𝑃,Float64(tmp[i,:time]))\n",
    "    tmpu[i]=du[1]\n",
    "    tmpv[i]=du[2]\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## 6. Compute Trajectories\n",
    "\n",
    "Solve through time using `OrdinaryDiffEq.jl` with\n",
    "\n",
    "- `dxdt!` is the function computing `d(position)/dt`\n",
    "- `uInit` is the initial condition `u @ tspan[1]`\n",
    "- `tspan` is the time interval\n",
    "- `𝑃` are parameters for `dxdt!`\n",
    "- `Tsit5` is the time-stepping scheme\n",
    "- `reltol` and `abstol` are tolerance parameters"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tspan = (0.0,nSteps*3600.0)\n",
    "#prob = OrdinaryDiffEq.ODEProblem(dxy_dt_replay,uInit,tspan,tmp)\n",
    "prob = OrdinaryDiffEq.ODEProblem(dxdt!,uInit,tspan,𝑃)\n",
    "sol = OrdinaryDiffEq.solve(prob,OrdinaryDiffEq.Tsit5(),reltol=1e-8,abstol=1e-8)\n",
    "sol[1:4]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Compare recomputed trajectories with originals from `MITgcm/pkg/flt`"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])\n",
    "maxLon=80*5.e3\n",
    "maxLat=42*5.e3\n",
    "#show(size(ref))\n",
    "for i=1:nSteps-1\n",
    "    ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing\n",
    "    ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing\n",
    "    ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing\n",
    "    ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing\n",
    "end\n",
    "ref=ref./dx;"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.2"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}