# User:Andrew Nikitin/drnav

## Dead Reckoning

Positions are represented as 2-element array of *lat lon*.
Both quantities are in radians.
*lat* must be between -π/2 and π/2. *lon* can be anything,
but canonical representations is -π<*lon*≤π.

Vectors are represented as 2-element array *az* (azimuth, also referred to as
*TC*, *true course*) and distance. Both quantities, including distance, are
in angular units, radians. We provide verbs that convert vector representations
from more natural units.

Usually vectors represent offset between 2 points on a sphere.
Sometimes vectors represent speeds. In this case *dist* component is
interpreted as *radians per hour*.

Depending on context, vectors can be interpreted as

- familiar 2D vectors on a plane
- arcs of great circle (in this case
*az*component refers to an angle that this arc makes with meridian at the origin) - part of loxodrome between 2 points on a sphere (in this case,
*dist*is measured along loxodrome); such "vector" looks like a straight line segment on a Mercator chart

[{{#file: "drnav.ijs"}} Download script: drnav.ijs ]

require '~nsg/geo/angular.ijs'

This library (angular.ijs) provides parsing and formatting verbs for various representations of angular values. [{{#file: "drnav.ijs"}} Download script: drnav.ijs ]

coinsert 'geo' cocurrent 'geo' «loxodrome» «vector» «chart» «misc»

### True mercator, loxodrome

[{{#file: "loxodrome"}} Download script: loxodrome ]

NB.*mlvecfpos v returns vector from x to y mlvecfpos=:(4 : 0)&posfpos 'lat lon'=.x 'lat2 lon2'=.y lon2=.lon ((= <./)@:(|@-~) # ]) lon2+/_1 0 1*2p1 tc=.12 o. j.~/(lon2-lon) , lat2 -&mercatorstretch lat if. 0.002<|c=.2 o. tc do. tc, |(lat2-lat)%c else. tc, |(2 o. -:lat2+lat)*(lon2-lon)%1 o. tc end. )

Loxodrome length is equal to . The formula does not work very well when (TC is horizontal) and may give roundoff errors when TC is close to horizontal. In this limiting case we use midlatitude formula. [{{#file: "loxodrome"}} Download script: loxodrome ]

NB.*mlposfvec v returns new position after moving along vector NB. x is old position NB. y is vector NB. result is new position NB. Movement is along loxodrome, line of equal azimuths mlposfvec=:4 : 0 'lat lon'=.posfpos x 'tc dist'=.y lat2=.lat+dist*c=.2 o. tc if. 0.002<|c do. lon2=.lon+(3 o. tc)*lat2 -&mercatorstretch lat else. lon2=.lon+dist * (1 o. tc)%2 o.-:lat+lat2 end. posfpos lat2,lon2 )

#### Shorthand

[{{#file: "loxodrome"}} Download script: loxodrome ]

NB.*pfv v = mlposfvec, new position from old and direction vector NB.*vfp v = mlvecfpos, direction vector from old position to new pfv=:mlposfvec vfp=:mlvecfpos

In order to simplify typing we give shorter (and, hopefully, still mnemonic) aliases to these verbs.

### Vector operations

[{{#file: "vector"}} Download script: vector ]

NB.*getgist v extract distance component from direction NB. vector and return it as nm getdist=:mfr@(1&{)

1 nautical mile = 1 minute of arc by definition. [{{#file: "vector"}} Download script: vector ]

NB.*newdist v x is direction vector, y is distance in nm NB. returns new direction vector with same TC and new NB. distance component newdist=:1}~ rfm

[{{#file: "vector"}} Download script: vector ]

NB.*rectify v rectify vector into (horizontal,vertical) rectify=:([: +. *.^:_1)&.|."1

`rectify` interprets vector as regular planar 2D vector and converts it into
"horizontal" and "vertical" components.
Horizontal means along latitude, vertical means along meridian.
The unit is still "radian", or rather *180·60/π* nautical miles per unit.

This simplification is useful for manipulation on small scale, such as current consideration, radar ranges, lines of position, etc.

*Note:* rectify&:_1 does not work at the origin
[{{#file: "vector"}} Download script: vector ]

NB.*scalevec v x is direction vector, y is scaling factor"1 scalevec=:* 1&,

If vector represents speed, then scaling it by number of hours will give actual offset. [{{#file: "vector"}} Download script: vector ]

NB.*Vec a applies operations on complex numbers to vectors NB. in (az,dist) format. Most useful are +Vec and -Vec Vec=:&.(*.^:_1@|.) ("1)

Vec converts argument vector array(s) to complex number(s) and performs operation on the number(s).

Vec=:&.rectify

Alternative for `scalevec`:

scalevec=:3 : 'x&*Vec y'

### Grid paper charting

Grid paper can be used as convenient charting medium. While Mercator grid paper is available commercially, regular (1/4", 1/5", 5 mm or any other kind) may be used with only little bit more effort.

Grid imposes rectangular coordinates on a piece of paper, what is needed is to find way to convert them to and from geographic coordinates. [{{#file: "chart"}} Download script: chart ]

MXY=:1 0 0

The conversion depends on 3 parameters: *scale* and the anchor for *lower
left corner*. These 3 parameters will be referred to as *M*, *X* and *Y*
and will be stored as 3 element vector in global `MXY`.
[{{#file: "chart"}} Download script: chart ]

islon=:, ,&1 0 islat=:, mercatorstretch , 0 1"_ islatlon=:((islon&{. ,: islat&{:)|.)"1

These parameters can be calculated by a library in several ways.
One of horizontal lines of a grid must be identified with certain
geographic latitude by means of verb `islat`.
One of the vertical lines must be identified with certain geographic longitude
(verb `islon`).
Verb `x_y islatlon lat_lon` combines these 2 verbs by assigning given point with
(x,y) coordinates (lat,lon) point.

In addition to these, one more line (either horizontal
or vertical, choice depends on convenience) should be identified with
appropriate geographic line (verb `islon` or `islat`, depending on kind
of a line). These 3 identification, together with the assumption that
geographic lines are spaced like on Mercator chart, are enough to
determine M, X and Y.
[{{#file: "chart"}} Download script: chart ]

setchart=: 3 : 0 MXY=:({."1 %. }."1) y if. 0>{.MXY do. i=.(<a:;1) y=. ((i{y) * 1 _1{~2{"1 y) i} y MXY=:({."1 %. }."1) y end. )

The verbs `islat`, `islon`, `islatlon` generate individual linear
equations for parameters `M`, `X` and `Y`. When there are enough of them
(3), these equations can be solved and `MXY` can be determined from
them. `setchart` does just that. If `M` turns out negative, it is
interpreted as Y axis direction reversal. This allows to use coordinate
system used in the bitmaps so that latitude increases towards top of the
image (y coordinate decreases) while east longitude increases left
to right.
[{{#file: "chart"}} Download script: chart ]

xyfpos=:3 : 0 MXY xyfpos y : 'lat lon'=.y 'M X Y'=.x (X+(|M)*lon),(Y+M*mercatorstretch lat) )

[{{#file: "chart"}} Download script: chart ]

posfxy=:3 : 0 MXY posfxy y : 'xx yy'=.y 'M X Y'=.x (mercatorstretch inv M%~yy-Y),(|M)%~xx-X )

({."1 %. }."1) (0 islat rfd 43),(35 islat rfd 45),(0 islon rfd _70),i.0 0

### Drift calculator

[{{#file: "misc"}} Download script: misc ]

NB.*driftsc v calculates course to steer and speed made good NB. for x=current speed vector and y=(desired_tc,boat_speed) array driftsc=:4 : 0 delta=.1 o. x -&{. y sa=.y (%~ *&delta)&{: x if. 1>:|sa do. nc=.2p1 | ({.y) - _1 o. sa nc, {: x +Vec nc, {:y else. _ 0 end. )

*boat_speed* is in units of radian/hour. Note that `y` and *result*
superficially resemble vectors in having direction as their 0
component and distance/speed as 1 component. However, in both `y`
and *result*, their components are not related to each other.

**WARNING**: The formula assumes that boat speed relative to the water does
not depend on course steered. This is not the case for sailboats.

### LOP intersection

[{{#file: "misc"}} Download script: misc ]

NB.*mpp v least squares D OFF from LOPS defined by orthogonal vectors NB. point of maximum probability=minimum sum of square distances from NB. point to lines mpp_rect=:({: %. [: |: 1 2 o./ {.)@|: mpp_off=:rectify^:_1@mpp_rect

Takes n×2 array of (Az,intercept) pairs. Each pair defines line, which is an approximation to a circle of equal altitudes. Returns pair (Az,dist), which is an offset to an MPP -- maximum probability position, given the observations.

(Used to be implemented as

doff=:((] %. %"_1) +/&.:*:"1)&.:rectify doff=: (%.~ +/&:*:"1)&.:rectify

which all have deficiencies or outright wrong)

### Utilities

Validators, converters, etc. [{{#file: "misc"}} Download script: misc ]

NB.*posfpos v Check the validity of position vector NB. and coerse lon into -π..π posfpos=:3 : 0 assert. 2=#y assert. 1r2p1>:|0{y NB. ({.y),(] - 2p1 * 1p1 < ])@(2p1&|) {:y ({.y),2p1 ([ toneg |) {:y )

[{{#file: "misc"}} Download script: misc ]

NB.*mercatorstretch v latitude to y mercator coordinate mercatorstretch=:^.@(3&o.)@(1r4p1 + -:)

### Mercator midlatitude

`mc...` procedures are deprecated, their algorithms are included as
limiting cases into `ml...` family.

The original reason to include them was to allow to reproduce/check calculations made using midlatitude method. [{{#file: "midlatitude"}} Download script: midlatitude ]

NB.*mcvecfpos ? returns vector from x to y NB. DEPRECATED by mlvecfpos mcvecfpos=:4 : 0 'lat lon'=.x 'lat2 lon2'=.y |.*.(2 o. -:lat+lat2)*(lat2 -&mercatorstretch lat)j.(lon2-lon) )

[{{#file: "midlatitude"}} Download script: midlatitude ]

NB.*mcposfvec ? returns new pos from x by y NB. DEPRECATED by mlposfvec mcposfvec=: 4 : 0 'lat lon'=.x 'tc dist'=.y lat2=.lat+dist * 2 o. tc lat2,lon+dist * (1 o. tc)%2 o.-:lat+lat2 )

Contributed by Andrew Nikitin