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 $\Delta \phi /\cos {\rm {TC}}$. The formula does not work very well when $\cos {\rm {TC}}=0$ (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 + -:)

$y=\ln \tan(\phi /2+45^{\circ })$

• http://www.erikdeman.de/html/sail020m.htm

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