The Liz grid.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with a great arc passing from $B$ to $C$ and know that $ABCD$ will form a box with $AB$ and $CD$ also great arcs and $ABC$ and $BCD$ being right angles.\n",
    "\n",
    "- What curve is $AD$ is perpendicular to $AB$ and $CD$?\n",
    "- What are the scale factors and grid metrics for a grid laid down within $ABCD$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following dequence of rotations that map the three known sides of the box (great arcs, $AB$, $BC$, and $CD$) to a regular box in a familiar orthogonal coordinate system.\n",
    "\n",
    "1. What rotation $R_1$ brings $B$ to the Greenwich Meridian (i.e. longitude of 0$^\\circ$E)? (This is for convenience)\n",
    "\n",
    "  - A rotation about the North pole of $-\\lambda_B$.\n",
    "\n",
    "2. What rotation $R_2$ then brings $R_1 B$ on to the Equator?\n",
    "\n",
    "  - A rotation of $-\\phi_B$ about the pole passing through the Equator at longitudes $\\pm 90^\\circ$.\n",
    "\n",
    "3. What rotation about $R_3 B$ then brings $R_2 R_1 C$ on to the Equator so that $R_3 R_2 R_1 AB$ is in the Northern hemisphere?\n",
    "\n",
    "  - A rotation about the x-axis by some angle that has to be determined.\n",
    "\n",
    "$B'C'$ is now a great arc on the Equator. Label the resulting box $A'B'C'D' = R_3 R_2 R_1 ABCD$.\n",
    "\n",
    "- If $A'B'$ is perpendicular to $B'C'$ and is a great arc then it falls on a Meridian.\n",
    "- If $C'D'$ is perpendicular to $B'C'$ and is a great arc then it falls on a Meridian.\n",
    "\n",
    "Requiring $AD$ to be perpendicular to both $AB$ and $CD$ is equivalent to requiring $A'D'$ to be perpendicular to both $A'B'$ and $C'D'$ since rotations do not change angles between lines. The curve that is perpendicular to a Meridian is a line of constant latitude. Thus $A'D'$ is a line of constant latitude and $AD$ is an small arc obtained by transforming $R_1^{-1}R_2^{-1}R_3^{-1}$ the line of constant latitude, $A'D'$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid generation\n",
    "\n",
    "A grid that fits $A'B'C'D'$ has nodes given by\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda'(i,j) &= i \\Delta \\lambda \\;\\; \\mbox{for} \\;\\; i = 0, 1, \\ldots, N \\\\\n",
    "\\phi'(i,j) &= j \\Delta \\phi \\;\\; \\mbox{for} \\;\\; j = 0, 1, \\ldots, M\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working in 3D coordinates (x,y,z)\n",
    "\n",
    "Use the right-hand coordiante system $x,y,z$ where the x-axis ($y=z=0$) passes through the Greenwich Meridian on the equator, y-axis ($x=0,z=0$) is at right angles and on the EQuator and the z-axis ($x=0,y=0$) passes through the North pole.\n",
    "\n",
    "The position in space of a point with geographical coordinates $\\lambda,\\phi$ (corresponding to longitude and latitude measured in radians) is\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\n",
    "x \\\\ y \\\\ z\n",
    "\\end{pmatrix}\n",
    "&=\n",
    "\\begin{pmatrix}\n",
    "R_e \\cos{\\phi} \\cos{\\lambda} \\\\\n",
    "R_e \\cos{\\phi} \\sin{\\lambda} \\\\\n",
    "R_e \\sin{\\phi}\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "where $R_e$ is the (spherical) radius of the Earth.\n",
    "\n",
    "The geographic coordiantes of a point $\\mathbf{r}=(x,y,z)$ projected to the Earth's surface from the origin is\n",
    "\n",
    "\\begin{align}\n",
    "\\phi &= \\sin^{-1} \\frac{z}{|\\mathbf{r}|} \\\\\n",
    "\\lambda &= \\tan^{-1} \\frac{y}{x}\n",
    "\\end{align}\n",
    "\n",
    "Note that for points on the Earth's surface, $|r|=Re$. The above rotations can now be described in terms of rotations matrices that operate on the position vectors of points we are concerned with.\n",
    "\n",
    "A rotation about the North pole of angle $\\alpha$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\n",
    "\\cos \\alpha & -\\sin \\alpha & 0 \\\\\n",
    "\\sin \\alpha & \\cos \\alpha & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "so $R_1$ is\n",
    "\n",
    "\\begin{align}\n",
    "R_1 =\n",
    "\\begin{pmatrix}\n",
    "\\cos \\lambda_B & \\sin \\lambda_B & 0 \\\\\n",
    "-\\sin \\lambda_B & \\cos \\lambda_B & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "A rotation about the y-axis of angle $\\beta$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\n",
    "\\cos \\beta & 0 & \\sin \\beta \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "-\\sin \\beta & 0 & \\cos \\beta\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "so $R_2$ is\n",
    "\n",
    "\\begin{align}\n",
    "R_2 &=\n",
    "\\begin{pmatrix}\n",
    "\\cos \\phi_B & 0 & -\\sin \\phi_B \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "\\sin \\phi_B & 0 & \\cos \\phi_B\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "A rotation about the x-axis of angle $\\gamma$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\begin{pmatrix}\n",
    "1 & 0 & 0 \\\\\n",
    "0 & \\cos \\gamma & -\\sin \\gamma \\\\\n",
    "0 & \\sin \\gamma & \\cos \\gamma\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Now $\\gamma$ is determined by the position of $C$ after $R_2 R_1$ which is\n",
    "\n",
    "\\begin{align}\n",
    "R_2 R_1 C &=\n",
    "\\begin{pmatrix}\n",
    "\\cos \\phi_B & 0 & -\\sin \\phi_B \\\\\n",
    "0 & 1 & 0 \\\\\n",
    "\\sin \\phi_B & 0 & \\cos \\phi_B\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "\\cos \\lambda_B & \\sin \\lambda_B & 0 \\\\\n",
    "-\\sin \\lambda_B & \\cos \\lambda_B & 0 \\\\\n",
    "0 & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "R_e \\cos{\\phi_C} \\cos{\\lambda_C} \\\\\n",
    "R_e \\cos{\\phi_C} \\sin{\\lambda_C} \\\\\n",
    "R_e \\sin{\\phi_C}\n",
    "\\end{pmatrix}\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Haversine formula should be used (for accuracy) to calculate the length $B'C'=BC$:\n",
    "\n",
    "\\begin{align}\n",
    "B'C'/R_e &= BC/R_e \\\\\n",
    "&= 2 \\sin^{-1} \\left( \\sqrt{ \\sin^2 \\left( \\frac{\\phi_B - \\phi_C}{2} \\right) +\n",
    "\\cos \\phi_B \\cos \\phi_C \\sin^2 \\left( \\frac{\\lambda_B - \\lambda_C}{2} \\right)}\\right)\n",
    "\\end{align}\n",
    "\n",
    "That angle can be found by dropping a perpendicular (great arc) from $C'$ to the Equator. Label the intersection point $P$. From the spherical law of cosines we then know angle $\\gamma_{C'}$ between $B'P$ (along the Equator) and $B'C'$ is\n",
    "\n",
    "\\begin{align}\n",
    "\\cos \\gamma_{C'}\n",
    "&= \\frac{\\cos PC'/R_e - \\cos B'P/R_e \\cos B'C'/R_e}{ \\sin B'P/R_e \\sin B'C'/R_e } \\\\\n",
    "&= \\frac{\\cos \\phi_{C'} - \\cos \\lambda_{C'} \\cos B'C'/R_e}{ \\sin \\lambda_{C'} \\sin B'C'/R_e }\n",
    "\\end{align}"
   ]
  }
 ],
 "metadata": {
  "gist": {
   "data": {
    "description": "The Liz grid.ipynb",
    "public": false
   },
   "id": ""
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}