Weather routing: Difference between revisions
m
→{{header|Wren}}: Changed to Wren S/H
No edit summary |
m (→{{header|Wren}}: Changed to Wren S/H) |
||
| (28 intermediate revisions by 5 users not shown) | |||
Line 10:
Given the above information and a specific path, progress and arrival time are determined. The weather routing problem, conversely, is to determine the path which results in the earliest arrival time.
=={{header|Go}}==
{{trans|Julia}}
This runs in only 37 seconds which is surprisingly quick compared to Julia. However, I've just noticed that I'm using an out of date version of Julia (1.0.4) so hopefully the latest version will be able to close the gap.
<syntaxhighlight lang="go">package main
import (
"fmt"
"io/ioutil"
"log"
"math"
"strconv"
"strings"
)
type matrixF = [][]float64
type pred = func(float64) bool
/*
Structure that represents a polar CSV file's data.
Note 0 degrees is directly into the wind, 180 degrees is directly downwind.
*/
type SailingPolar struct {
winds []float64 // vector of windspeeds
degrees []float64 // vector of angles in degrees of direction relative to the wind
speeds matrixF // matrix of sailing speeds indexed by wind velocity and angle of boat to wind
}
/*
Structure that represents wind and surface current direction and velocity for a given position.
Angles in degrees, velocities in knots.
*/
type SurfaceParameters struct{ windDeg, windKts, currentDeg, currentKts float64 }
// Checks for fatal errors.
func check(err error) {
if err != nil {
log.Fatal(err)
}
}
// Reads a sailing polar CSV file and returns a SailingPolar struct containing the file data.
// A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
// The first line of file contains labels for the wind velocities that make up columns, and
// the first entry of each row makes up a column of angle of sailing direction from wind in degrees.
func getPolarData(fileName string) *SailingPolar {
content, err := ioutil.ReadFile(fileName)
check(err)
lines := strings.Split(string(content), "\n")
line0 := strings.TrimSpace(lines[0])
header := strings.Split(line0, ";")
var winds, degrees []float64
var speeds matrixF
for _, col := range header[1:] {
t, err := strconv.ParseFloat(col, 64)
check(err)
winds = append(winds, t)
}
for _, line := range lines[1:] {
line = strings.TrimSpace(line)
if line == "" {
break // ignore final blank line if there is one
}
cols := strings.Split(line, ";")
f, err := strconv.ParseFloat(cols[0], 64)
check(err)
degrees = append(degrees, f)
var temp []float64
for _, col := range cols[1:] {
t, err := strconv.ParseFloat(col, 64)
check(err)
temp = append(temp, t)
}
speeds = append(speeds, temp)
}
return &SailingPolar{winds, degrees, speeds}
}
const R = 6372800.0 // Earth's approximate radius in meters
/* various helper methods which work with degrees rather than radians. */
// Converts degrees to radians.
func deg2Rad(deg float64) float64 { return math.Mod(deg*math.Pi/180+2*math.Pi, 2*math.Pi) }
// Converts radians to degrees.
func rad2Deg(rad float64) float64 { return math.Mod(rad*180/math.Pi+360, 360) }
// Trig functions.
func sind(d float64) float64 { return math.Sin(deg2Rad(d)) }
func cosd(d float64) float64 { return math.Cos(deg2Rad(d)) }
func asind(d float64) float64 { return rad2Deg(math.Asin(d)) }
func atand(x, y float64) float64 { return rad2Deg(math.Atan2(x, y)) }
// Calculates the haversine function for two points on the Earth's surface.
// Given two latitude, longitude pairs in degrees for a point on the Earth,
// get distance in meters and the initial direction of travel in degrees for
// movement from point 1 to point 2.
func haversine(lat1, lon1, lat2, lon2 float64) (float64, float64) {
dlat := lat2 - lat1
dlon := lon2 - lon1
a := math.Pow(sind(dlat/2), 2) + cosd(lat1)*cosd(lat2)*math.Pow(sind(dlon/2), 2)
c := 2 * asind(math.Sqrt(a))
theta := atand(sind(dlon)*cosd(lat2), cosd(lat1)*sind(lat2)-sind(lat1)*cosd(lat2)*cosd(dlon))
theta = math.Mod(theta+360, 360)
return R * c * 0.5399565, theta
}
// Returns the index of the first element of 'a' for which 'p' returns true or -1 otherwise.
func findFirst(a []float64, p pred) int {
for i, e := range a {
if p(e) {
return i
}
}
return -1
}
// Returns the index of the last element of 'a' for which 'p' returns true or -1 otherwise.
func findLast(a []float64, p pred) int {
for i := len(a) - 1; i >= 0; i-- {
if p(a[i]) {
return i
}
}
return -1
}
// Calculate the expected sailing speed in a specified direction in knots,
// given sailing polar data, a desired point of sail in degrees, and wind speed in knots.
func boatSpeed(sp *SailingPolar, pointOfSail, windSpeed float64) float64 {
winds := sp.winds
degrees := sp.degrees
speeds := sp.speeds
udeg := findLast(degrees, func(t float64) bool { return t <= pointOfSail })
odeg := findFirst(degrees, func(t float64) bool { return t >= pointOfSail })
uvel := findLast(winds, func(t float64) bool { return t <= windSpeed })
ovel := findFirst(winds, func(t float64) bool { return t >= windSpeed })
if udeg == -1 || odeg == -1 || uvel == -1 || ovel == -1 {
return -1
}
var frac float64
switch {
case odeg == udeg && uvel == ovel:
frac = 1
case odeg == udeg:
frac = (windSpeed - winds[uvel]) / (winds[ovel] - winds[uvel])
case uvel == ovel:
frac = (pointOfSail - degrees[udeg]) / (degrees[odeg] - degrees[udeg])
default:
frac = ((pointOfSail-degrees[udeg])/(degrees[odeg]-degrees[udeg]) +
(windSpeed-winds[uvel])/(winds[ovel]-winds[uvel])) / 2
}
return speeds[udeg][uvel] + frac*(speeds[odeg][ovel]-speeds[udeg][uvel])
}
// Calculates the expected net boat speed in a desired direction versus the wind ('azimuth').
// This is generally different from the actual boat speed in its actual direction.
// Directions are in degrees ('pointos' is point of sail the ship direction from the wind),
// and velocity of wind ('ws') is in knots.
func sailingSpeed(sp *SailingPolar, azimuth, pointos, ws float64) float64 {
return boatSpeed(sp, pointos, ws) * cosd(math.Abs(pointos-azimuth))
}
// Calculates the net direction and velocity of a sailing ship.
// Arguments are sailing polar data, direction of travel in degrees from north, wind direction in
// degrees from north, wind velocity in knots, surface current direction in degrees, and
// current velocity in knots.
func bestVectorSpeed(sp *SailingPolar, dirTravel, dirWind, windSpeed, dirCur, velCur float64) (float64, float64) {
azimuth := math.Mod(dirTravel-dirWind, 360)
if azimuth < 0 {
azimuth += 360
}
if azimuth > 180 {
azimuth = 360 - azimuth
}
vmg := boatSpeed(sp, azimuth, windSpeed)
other := -1.0
idx := -1
for i, d := range sp.degrees {
ss := sailingSpeed(sp, azimuth, d, windSpeed)
if ss > other {
other = ss
idx = i
}
}
if other > vmg {
azimuth = sp.degrees[idx]
vmg = other
}
dirChosen := deg2Rad(dirWind + azimuth)
wx := vmg * math.Sin(dirChosen)
wy := vmg * math.Cos(dirChosen)
curX := velCur * math.Sin(deg2Rad(dirCur))
curY := velCur * math.Cos(deg2Rad(dirCur))
return rad2Deg(math.Atan2(wy+curY, wx+curX)), math.Sqrt(math.Pow(wx+curX, 2) + math.Pow(wy+curY, 2))
}
// Calculates the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2).
// Uses the data in SurfaceParameters for wind and current velocity and direction.
func sailSegmentTime(sp *SailingPolar, p SurfaceParameters, lat1, lon1, lat2, lon2 float64) float64 {
distance, dir := haversine(lat1, lon1, lat2, lon2)
_, vel := bestVectorSpeed(sp, dir, p.windDeg, p.windKts, p.currentDeg, p.currentKts)
// minutes/s * m / (knots * (m/s / knot)) = minutes
return (1.0 / 60.0) * distance / (vel * 1.94384)
}
/* Structure that represents a point in 2-D space. */
type Point2 struct{ x, y int }
func (p Point2) add(p2 Point2) Point2 { return Point2{p.x + p2.x, p.y + p2.y} }
func (p Point2) equals(p2 Point2) bool { return p.x == p2.x && p.y == p2.y }
func (p Point2) String() string { return fmt.Sprintf("[%d, %d]", p.x, p.y) }
/*
Structure that consists of a tuple of latitude and longitude in degrees.
NB: This uses latitude (often considered to be y) first then longitude (often considered to be x).
This latitude, then longitude ordering is as per ISO 6709 (en.wikipedia.org/wiki/ISO_6709).
*/
type Position struct{ lat, lon float64 }
/* Structure that represents a Position with the SurfaceParameters of wind and current at the Position. */
type GridPoint struct {
pt Position
sp SurfaceParameters
}
type MatrixG = [][]*GridPoint
/*
Type alias for a matrix of GridPoints, each Position point with their SurfaceParameters.
A Vector of TimeSlice can give the surface characteristics for an ocean region over time.
*/
type TimeSlice = MatrixG
/* Structure that represents a routing problem. */
type RoutingProblem struct {
timeInterval float64 // the minutes duration for each TimeSlice
timeFrame []TimeSlice // a vector of sequential timeslices for the ocean region
obstacleIndices []Point2 // the Cartesian indices in each TimeSlice for
// obstacles, such as land or shoals, where the ship may not go
startIndex int // the TimeSlice position for time of starting
start Point2 // starting location on grid of GridPoints
finish Point2 // destination / finish location on grid of GridPoints
allowRepeatVisits bool // whether the vessel may overlap its prior path, usually false
}
/* Structure that represents a timed path. */
type TimedPath struct {
duration float64 // minutes total to travel the path
path []Point2 // vector of Cartesian indices of points in grid for path to travel
}
func (t TimedPath) String() string { return fmt.Sprintf("%g %v", t.duration, t.path) }
func (t TimedPath) equals(t2 TimedPath) bool { return t.String() == t2.String() }
func findMin(a []float64) (float64, int) {
min := a[0]
idx := 0
for i, e := range a[1:] {
if e < min {
min = e
idx = i + 1
}
}
return min, idx
}
var ntuples = [][2]int{{-1, -1}, {-1, 0}, {-1, 1}, {0, -1}, {0, 1}, {1, -1}, {1, 0}, {1, 1}}
var neighbors = make([]Point2, len(ntuples))
func init() {
for i := 0; i < len(ntuples); i++ {
neighbors[i] = Point2{ntuples[i][0], ntuples[i][1]}
}
}
func contains(points []Point2, point Point2) bool {
for _, p := range points {
if p.equals(point) {
return true
}
}
return false
}
// Returns a list of points surrounding 'p' which are not otherwise excluded.
func surround(p Point2, mat TimeSlice, excluded []Point2) []Point2 {
xmax := len(mat)
ymax := len(mat[0])
var res []Point2
for _, x := range neighbors {
q := x.add(p)
if (0 <= q.x && q.x < xmax) && (0 <= q.y && q.y < ymax) && !contains(excluded, q) {
res = append(res, q)
}
}
return res
}
// Get the route (as a TimedPath) that minimizes time from start to finish for a given
// RoutingProblem (sea parameters) given sailing polar data (ship parameters).
func minimumTimeRoute(rp *RoutingProblem, sp *SailingPolar, verbose bool) *TimedPath {
timedPaths := []*TimedPath{&TimedPath{0, []Point2{rp.start}}}
completed := false
minPath := &TimedPath{1000, []Point2{}}
for i := 0; i < 1000; i++ {
var newPaths []*TimedPath
if verbose {
fmt.Printf("Checking %d paths of length %d\n", len(timedPaths), len(timedPaths[0].path))
}
for _, tpath := range timedPaths {
le := len(tpath.path)
if tpath.path[le-1] == rp.finish {
completed = true
newPaths = append(newPaths, tpath)
} else {
p1 := tpath.path[le-1]
num := int(math.Round(tpath.duration))
den := int(math.Round(rp.timeInterval))
slice := rp.timeFrame[(num/den)%len(rp.timeFrame)]
for _, p2 := range surround(p1, slice, rp.obstacleIndices) {
if !rp.allowRepeatVisits && contains(tpath.path, p2) {
continue
}
gp1 := slice[p1.x][p1.y]
gp2 := slice[p2.x][p2.y]
lat1 := gp1.pt.lat
lon1 := gp1.pt.lon
lat2 := gp2.pt.lat
lon2 := gp2.pt.lon
t := sailSegmentTime(sp, gp1.sp, lat1, lon1, lat2, lon2)
path := make([]Point2, len(tpath.path))
copy(path, tpath.path)
path = append(path, p2)
newPaths = append(newPaths, &TimedPath{tpath.duration + t, path})
}
}
}
set := make(map[string]*TimedPath)
for _, np := range newPaths {
set[np.String()] = np
}
timedPaths = timedPaths[:0]
for k := range set {
timedPaths = append(timedPaths, set[k])
}
if completed {
var durations []float64
for _, x := range timedPaths {
durations = append(durations, x.duration)
}
minDur, _ := findMin(durations)
var finished []*TimedPath
for _, x := range timedPaths {
le := len(x.path)
if x.path[le-1] == rp.finish {
finished = append(finished, x)
}
}
durations = durations[:0]
for _, x := range finished {
durations = append(durations, x.duration)
}
minFinDur, idx := findMin(durations)
if verbose {
fmt.Printf("Current finished minimum: %v, others %v\n", finished[idx], minDur)
}
if minDur == minFinDur {
minPath = finished[idx]
break
}
}
}
return minPath
}
/*
The data is selected so the best time path is slightly longer than the
shortest length path. The forbidden regions are x, representing land or reef.
The allowed sailing points are . and start and finish are S and F.
x . . F . . x . x
. . . . . . . x x
x . . x x x . . .
. . x x x x . x x
x . . . x x . x .
x . . . x x . x .
. . . . x . . x .
x . . . . . . x .
. . . S . . . . .
*/
// These need to be changed to 0-based for Go.
var ftuples = [][2]int{
{1, 8}, {2, 1}, {2, 8}, {3, 5}, {3, 8}, {4, 1}, {4, 5}, {4, 6}, {4, 8}, {5, 1},
{5, 5}, {5, 6}, {5, 8}, {6, 3}, {6, 4}, {6, 5}, {6, 6}, {6, 8}, {6, 9}, {7, 1},
{7, 4}, {7, 5}, {7, 6}, {8, 8}, {8, 9}, {9, 1}, {9, 7}, {9, 9},
}
var forbidden = make([]Point2, len(ftuples))
func init() {
for i := 0; i < len(ftuples); i++ {
forbidden[i] = Point2{ftuples[i][0] - 1, ftuples[i][1] - 1}
}
}
// Create regional wind patterns on the map.
func surfaceByLongitude(lon float64) SurfaceParameters {
switch {
case lon < -155.03:
return SurfaceParameters{-5, 8, 150, 0.5}
case lon < -155.99:
return SurfaceParameters{-90, 20, 150, 0.4}
default:
return SurfaceParameters{180, 25, 150, 0.3}
}
}
// Vary wind speeds over time.
func mutateTimeSlices(slices []TimeSlice) {
i := 1
for _, slice := range slices {
for j := 0; j < len(slice); j++ {
for k := 0; k < len(slice[j]); k++ {
x := slice[j][k]
x.sp = SurfaceParameters{x.sp.windDeg, x.sp.windKts * (1 + 0.002*float64(i)),
x.sp.currentDeg, x.sp.currentKts}
}
}
i++
}
}
func main() {
startPos := Point2{0, 3} // 0-based
endPos := Point2{8, 3} // ditto
slices := make([]MatrixG, 200)
for s := 0; s < 200; s++ {
gpoints := make([][]*GridPoint, 9)
for i := 0; i < 9; i++ {
gpoints[i] = make([]*GridPoint, 9)
for j := 0; j < 9; j++ {
pt := Position{19.78 - 1.0/60.0 + float64(i)/60, -155.0 - 5.0/60.0 + float64(j)/60}
gpoints[i][j] = &GridPoint{pt, surfaceByLongitude(pt.lon)}
}
}
slices[s] = gpoints
}
mutateTimeSlices(slices)
routeProb := &RoutingProblem{10, slices, forbidden, 0, startPos, endPos, false}
fileName := "polar.csv"
sp := getPolarData(fileName)
tp := minimumTimeRoute(routeProb, sp, false)
fmt.Println("The route taking the least time found was:\n", tp.path, "\nwhich has duration",
int(tp.duration/60), "hours,", int(math.Round(math.Mod(tp.duration, 60))), "minutes.")
}</syntaxhighlight>
{{out}}
<pre>
The route taking the least time found was:
[[0, 3] [0, 4] [1, 5] [2, 6] [3, 6] [4, 6] [5, 6] [6, 6] [7, 5] [7, 4] [8, 3]]
which has duration 4 hours, 44 minutes.
</pre>
=={{header|Julia}}==
Brute force optimization search, practical for shorter path lengths, but would require a better algorithm for paths over twice this size.
<syntaxhighlight lang="julia">module SailingPolars
using DelimitedFiles
export SailingPolar, SurfaceParameters, getpolardata, deg2rad, rad2deg, cartesian2polar
export polar2cartesian, haversine, inverse_haversine, boatspeed, bestvectorspeed
export sailingspeed, sailsegmenttime
"""
Structure to represent a polar CSV file's data.
Contains a matrix, speeds, of sailing speeds indexed by wind velocity and angle of boat to wind
winds is a list of wind speeds
degrees is a list of angles in degrees of direction relative to the wind
Note 0.0 degrees is directly into the wind, 180 degrees is directly downwind.
"""
struct SailingPolar
winds::Vector{Float32}
degrees::Vector{Float32}
speeds::Matrix{Float32} # speeds[wind direction degrees, windspeed knots]
end
"""
struct SurfaceParameters
Structure that represents wind and surface current direction and velocity for a given position
Angles in degrees, velocities in knots
"""
struct SurfaceParameters
winddeg::Float32
windkts::Float32
currentdeg::Float32
currentkts::Float32
end
"""
function getpolardata(filename)
Read a sailing polar CSV file and return a SailingPolar containing the file data.
A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
The first line of file contains labels for the wind velocities that make up columns, and
the first entry of each row makes up a column of angle of sailing direction from wind in degrees
"""
function getpolardata(filename)
datacells, headercells = readdlm(filename, ';', header=true)
winds = map(x -> parse(Float32, x), headercells[2:end])
degrees = datacells[:, 1]
speeds = datacells[:, 2:end]
return SailingPolar(winds, degrees, speeds)
end
const R = 6372800 # Earth's approximate radius in meters
"""
deg2rad(deg)
Convert degrees to radians
"""
deg2rad(deg) = (deg * π / 180.0 + 2π) % 2π
"""
rad2deg(rad)
Convert radians to degrees
"""
rad2deg(rad) = (rad * (180.0 / π) + 360.0) % 360.0
"""
cartesian2polard(x, y)
Convert x, y coordinates to polar coordinates with angle in degrees
"""
cartesian2polard(x, y) = sqrt(x * x + y * y), atand(x, y)
"""
polard2cartesian(r, deg)
Convert polar coordinates in degrees to cartesian x, y coordinates
"""
polard2cartesian(r, deg) = r .* sincosd(deg)
"""
function haversine(lat1, lon1, lat2, lon2)
Calculate the haversine function for two points on the Earth's surface.
Given two latitude, longitude pairs in degrees for a point on the Earth,
get distance in meters and the initial direction of travel in degrees for
movement from point 1 to point 2.
"""
function haversine(lat1, lon1, lat2, lon2)
dlat = lat2 - lat1
dlon = lon2 - lon1
a = sind(dlat / 2)^2 + cosd(lat1) * cosd(lat2) * sind(dlon / 2)^2
c = 2.0 * asind(sqrt(a))
theta = atand(sind(dlon) * cosd(lat2),
cosd(lat1) * sind(lat2) - sind(lat1) * cosd(lat2) * cosd(dlon))
theta = (theta + 360) % 360
return R * c * 0.5399565, theta
end
"""
function inverse_haversine(lat1, lon1, distance, direction)
Calculate an inverse haversine function.
Takes the point of origin in degrees (latitude, longitude), distance in meters, and
initial direction in degrees, and returns the latitude and longitude of the endpoint
in degrees after traveling the specified distance.
"""
function inverse_haversine(lat1, lon1, distance, direction)
lat2 = asind(sind(lat1) * cos(distance / R) + cosd(lat1) * sin(distance / R) * cosd(direction))
lon2 = lon1 + atand(sind(direction) * sin(distance / R) * cosd(lat1),
cos(distance / R) - sind(lat1) * sind(lat2))
return lat2, lon2
end
"""
function boatspeed(sp::SailingPolar, pointofsail, windspeed)
Calculate the expected sailing speed in a specified direction in knots,
given sailing polar data, a desired point of sail in degrees, and wind speed in knots
"""
function boatspeed(sp::SailingPolar, pointofsail, windspeed)
winds, degrees, speeds = sp.winds, sp.degrees, sp.speeds
udeg = findlast(t -> t <= pointofsail, degrees)
odeg = findfirst(t -> t >= pointofsail, degrees)
uvel = findlast(t -> t <= windspeed, winds)
ovel = findfirst(t -> t >= windspeed, winds)
if any(t -> t == nothing, [udeg, odeg, uvel, ovel])
return -1.0
end
frac = (odeg == udeg && uvel == ovel) ? 1.0 :
(odeg == udeg) ? (windspeed - winds[uvel]) / (winds[ovel] - winds[uvel]) :
(uvel == ovel) ? (pointofsail - degrees[udeg]) / (degrees[odeg] - degrees[udeg]) :
((pointofsail - degrees[udeg]) / (degrees[odeg] - degrees[udeg]) +
(windspeed - winds[uvel]) / (winds[ovel] - winds[uvel])) / 2
return speeds[udeg, uvel] + frac * (speeds[odeg, ovel] - speeds[udeg, uvel])
end
"""
sailingspeed(sp::SailingPolar, azimuth, pointos, ws)
Calculate the expected net boat speed in a desired direction versus the wind (azimuth).
This is generally different from the actual boat speed in its actual direction.
Directions are in degrees (pointos is point of sail, the ship direction from wind),
and velocity of wind (ws) is in knots.
"""
sailingspeed(sp, azimuth, pointos, ws) = boatspeed(sp, pointos, ws) * cosd(abs(pointos - azimuth))
"""
function bestvectorspeed(sp::SailingPolar, dirtravel, dirwind, windspeed, dircur, velcur)
Calculate the net direction and velocity of a sailing ship.
Arguments are sailing polar data, direction of travel in degrees from north, wind direction in
degrees from north, wind velocity in knots, surface current direction in degrees, and
current velocity in knots.
"""
function bestvectorspeed(sp::SailingPolar, dirtravel, dirwind, windspeed, dircur, velcur)
azimuth = (dirtravel - dirwind) % 360.0
azimuth = azimuth < 0 ? azimuth + 360.0 : azimuth
azimuth = azimuth > 180.0 ? 360.0 - azimuth : azimuth
VMG = boatspeed(sp, azimuth, windspeed)
other, idx = findmax([sailingspeed(sp, azimuth, x, windspeed) for x in sp.degrees])
if other > VMG
azimuth = sp.degrees[idx]
VMG = other
end
dirchosen = deg2rad(dirwind + azimuth)
wx, wy = VMG * sin(dirchosen), VMG * cos(dirchosen)
curx, cury = velcur * sin(deg2rad(dircur)), velcur * cos(deg2rad(dircur))
return rad2deg(atan(wy + cury, wx + curx)), sqrt((wx + curx)^2 + (wy + cury)^2)
end
"""
function sailsegmenttime(sp::SailingPolar, p::SurfaceParameters, lat1, lon1, lat2, lon2)
Calculate the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2).
Uses the data in SurfaceParameters for wind and current velocity and direction.
"""
function sailsegmenttime(sp::SailingPolar, p::SurfaceParameters, lat1, lon1, lat2, lon2)
distance, dir = haversine(lat1, lon1, lat2, lon2)
dir2, vel = bestvectorspeed(sp, dir, p.winddeg, p.windkts, p.currentdeg, p.currentkts)
# minutes/s * m / (knots * (m/s / knot)) = minutes
return (1 / 60) * distance / (vel * 1.94384)
end
end # module
module SailingNavigation
export Position, lat, lon, GridPoint, TimeSlice, TimedPath, closestpoint, surround
export RoutingProblem, minimumtimeroute
using GeometryTypes
using ..SailingPolars
# NB: This uses latitude (often considered to be y) first then longitude (often considered to be x).
# This latitude, then longitude ordering is as per ISO 6709 (en.wikipedia.org/wiki/ISO_6709)
# Position is a Float32 2-tuple of latitude and longitude in degrees
Position = Point2f0
# latitude from Position
lat(p::Position) = p[1]
# longitude from Position
lon(p::Position) = p[2]
# A GridPoint is a Position with the SurfaceParameters of wind and current at the Position
mutable struct GridPoint
pt::Position
sp::SurfaceParameters
end
"""
TimeSlice
A TimeSlice is a matrix of GridPoints, each Position point with their SurfaceParameters
A Vector of TimeSlice can give the surface characteristics for an ocean region over time.
"""
TimeSlice = Matrix{GridPoint}
"""
mutable struct RoutingProblem
timeinterval: the minutes duration for each TimeSlice
timeframe: a vector of sequential timeslices for the ocean region
obstacleindices: the Cartesian indices in each timeslice for
obstacles, such as land or shoals, where the ship may not go
startindex: the timeslice position for time of starting
start: starting location on grid of GridPoints
finish: destination / finish location on grid of GridPoints
allowrepeatvisits: whether the vessel may overlap its prior path, usually false
"""
mutable struct RoutingProblem
timeinterval::Float64 # minutes between timeframe slices
timeframe::Vector{TimeSlice}
obstacleindices::Vector{Point2{Int}}
startindex::Int
start::Point2{Int}
finish::Point2{Int}
allowrepeatvisits::Bool
end
"""
mutable struct TimedPath
duration: minutes total to travel the path
path: vector of Cartesian indices of points in grid for path to travel
"""
mutable struct TimedPath
duration::Float64
path::Vector{Point2{Int}}
end
"""
closestpoint(p, mat)
Get the closest GridPoint in matrix mat to a given position p.
p: Cartesian indices of a Position (latitude, longitude in degrees) in grid of GridPoints
mat: matrix of Gridpoints
"""
closestpoint(p, mat) = findmin(gp -> haversine(p[1], p[2], gp.pt[1], gp.pt[2])[1], mat)[2]
function surround(p, mat, excluded)
neighbors = Point2{Int}[(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)]
(xmax, ymax) = size(mat)
return filter(q -> 1 <= q[1] <= xmax && 1 <= q[2] <= ymax && !(q in excluded),
[x + p for x in neighbors])
end
"""
function minimumtimeroute(rp::RoutingProblem, sp::SailingPolar, verbose=false)
Get the route (as a TimedPath) that minimizes time from start to finish for a given
RoutingProblem (sea parameters) given sailing polar data (ship parameters).
"""
function minimumtimeroute(rp::RoutingProblem, sp::SailingPolar, verbose=false)
timedpaths = [TimedPath(0.0, [rp.start])]
completed, mintime, minpath = false, 1000.0, TimedPath(1000.0, [])
for i in 1:1000
newpaths = TimedPath[]
verbose && println("Checking ", length(timedpaths), " paths of length ",
length(timedpaths[1].path) - 1)
for tpath in timedpaths
if tpath.path[end] == rp.finish
completed = true
push!(newpaths, tpath)
else
p1 = tpath.path[end]
slice = rp.timeframe[div(Int(round(tpath.duration)),
Int(round(rp.timeinterval))) %
length(rp.timeframe) + 1]
for p2 in surround([p1[1], p1[2]], slice, rp.obstacleindices)
!rp.allowrepeatvisits && p2 in tpath.path && continue
gp1, gp2 = slice[p1[1], p1[2]], slice[p2[1], p2[2]]
lat1, lon1, lat2, lon2 = gp1.pt[1], gp1.pt[2], gp2.pt[1], gp2.pt[2]
t = sailsegmenttime(sp, gp1.sp, lat1, lon1, lat2, lon2)
path = deepcopy(tpath.path)
push!(path, p2)
push!(newpaths, TimedPath(tpath.duration + t, path))
end
end
end
timedpaths = unique(newpaths)
if completed
mindur = minimum(map(x -> x.duration, timedpaths))
finished = filter(x -> x.path[end] == rp.finish, timedpaths)
minfindur, idx = findmin(map(x -> x.duration, finished))
verbose && println("Current finished minimum: ", finished[idx], ", others $mindur")
if mindur == minfindur
minpath = finished[idx]
break
end
end
end
return minpath
end
end # module
using GeometryTypes
using .SailingNavigation, .SailingPolars
#=
The data is selected so the best time path is slightly longer than the
shortest length path. The forbidden regions are x, representing land or reef.
The allowed sailing points are . and start and finish are S and F.
x . . F . . x . x
. . . . . . . x x
x . . x x x . . .
. . x x x x . x x
x . . . x x . x .
x . . . x x . x .
. . . . x . . x .
x . . . . . . x .
. . . S . . . . .
=#
const forbidden = Point2{Int}.([
[1, 8], [2, 1], [2, 8], [3, 5], [3, 8], [4, 1], [4, 5], [4, 6], [4, 8], [5, 1],
[5, 5], [5, 6], [5, 8], [6, 3], [6, 4], [6, 5], [6, 6], [6, 8], [6, 9], [7, 1],
[7, 4], [7, 5], [7, 6], [8, 8], [8, 9], [9, 1], [9, 7], [9, 9],
])
# Create regional wind patterns on the map.
function surfacebylongitude(lon)
return lon < -155.03 ? SurfaceParameters(-5.0, 8, 150, 0.5) :
lon < -155.99 ? SurfaceParameters(-90.0, 20, 150, 0.4) :
SurfaceParameters(180.0, 25, 150, 0.3)
end
# Vary wind speeds over time.
function mutatetimeslices!(slices)
for (i, slice) in enumerate(slices), x in slice
x.sp = SurfaceParameters(x.sp.winddeg, x.sp.windkts * (1 + 0.002 * i),
x.sp.currentdeg, x.sp.currentkts)
end
end
const startpos = Point2{Int}(1, 4)
const endpos = Point2{Int}(9, 4)
const pmat = [Position(19.78 - 1/60 + i/60, -155.0 - 5/60 + j/60) for i in 0:8, j in 0:8]
const gpoints = map(pt -> GridPoint(pt, surfacebylongitude(lon(pt))), pmat)
const slices = [deepcopy(gpoints) for _ in 1:200]
mutatetimeslices!(slices)
const routeprob = RoutingProblem(10.0, slices, forbidden, 1, startpos, endpos, false)
const filename = "polar.csv"
const sp = getpolardata(filename)
const tp = minimumtimeroute(routeprob, sp)
println("The route taking the least time found was:\n ", tp.path,
"\nwhich has duration $(div(tp.duration, 60)) hours, $(rem(tp.duration, 60)) minutes.")
</syntaxhighlight>
The polar CSV file used for this solution, named polar.csv, is as follows. Note that this is a very detailed
polar, chosen to stress the testing of the code. Most polar files are far smaller,
with fewer choices for angle and wind speed.
<pre style="height:80ex;overflow:scroll;">
TWA\TWS;0;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;35;40;60;70
40;0;0.53;0.54;0.49;0.4;0.31;0.21;0.16;0.11;0.08;0.05;0.03;0.02;0.01;0;0;0;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.05;-0.1;-0.11
41;0;0.61;0.62;0.56;0.47;0.36;0.25;0.19;0.14;0.1;0.07;0.04;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;0;0;-0.04;-0.09;-0.1
44;0;0.89;0.91;0.82;0.69;0.56;0.42;0.33;0.24;0.18;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;-0.02;-0.06;-0.06
45;0;0.99;1.02;0.92;0.78;0.64;0.49;0.39;0.29;0.22;0.15;0.1;0.07;0.04;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.05;-0.05
46;0;1.11;1.14;1.02;0.87;0.73;0.57;0.45;0.35;0.26;0.18;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.04;-0.05
47;0;1.23;1.25;1.14;0.97;0.82;0.66;0.53;0.41;0.31;0.22;0.15;0.1;0.07;0.04;0.02;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.01;-0.03;-0.04
48;0;1.37;1.37;1.26;1.08;0.93;0.76;0.61;0.48;0.36;0.26;0.19;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.02;-0.03
49;0;1.5;1.5;1.39;1.2;1.05;0.87;0.71;0.56;0.42;0.31;0.22;0.15;0.1;0.07;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.02;-0.02
50;0;1.65;1.64;1.52;1.33;1.18;1;0.81;0.65;0.49;0.37;0.26;0.19;0.13;0.08;0.05;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;-0.01;-0.02
51;0;1.79;1.77;1.67;1.46;1.32;1.13;0.92;0.74;0.57;0.43;0.31;0.22;0.15;0.1;0.07;0.05;0.03;0.02;0.02;0.01;0.01;0;0;0;0;0;0;0;0;-0.01;-0.02
53;0;2.1;2.07;1.99;1.76;1.62;1.4;1.14;0.95;0.74;0.57;0.43;0.31;0.22;0.16;0.1;0.08;0.06;0.04;0.03;0.02;0.01;0.01;0.01;0;0;0;0;0;0;-0.01;-0.01
54;0;2.26;2.22;2.16;1.92;1.78;1.55;1.28;1.06;0.84;0.65;0.5;0.37;0.27;0.19;0.13;0.1;0.07;0.06;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;-0.01
55;0;2.43;2.39;2.34;2.09;1.95;1.7;1.42;1.18;0.95;0.74;0.57;0.43;0.32;0.23;0.16;0.12;0.09;0.07;0.05;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;-0.01
60;0;3.29;3.33;3.33;3.08;2.93;2.64;2.29;1.98;1.66;1.36;1.1;0.88;0.68;0.53;0.39;0.32;0.26;0.21;0.17;0.13;0.1;0.07;0.05;0.04;0.03;0.02;0.01;0;0;0;0
70;0;5.2;5.53;5.74;5.59;5.5;5.22;4.84;4.46;3.94;3.51;3.08;2.65;2.26;1.9;1.55;1.38;1.22;1.06;0.92;0.78;0.66;0.55;0.46;0.37;0.3;0.24;0.18;0.03;0;0;0
80;0;6.8;7.43;7.97;8.02;8.23;8.34;8.2;7.9;7.37;6.91;6.43;5.9;5.32;4.72;4.12;3.83;3.55;3.25;2.96;2.67;2.4;2.13;1.88;1.65;1.43;1.22;1.04;0.37;0.09;0.01;0
90;0;7.59;8.5;9.4;9.73;10.4;11.16;11.53;11.56;11.3;11.05;10.77;10.44;9.83;9.07;8.34;8;7.65;7.27;6.88;6.46;6.04;5.61;5.15;4.74;4.33;3.88;3.51;1.72;0.67;0.12;0.03
100;0;7.34;8.25;9.16;9.86;10.5;11.95;12.79;13.5;14.02;14.4;14.37;14.5;14.4;13.92;13.52;13.19;12.79;12.51;12.1;11.66;11.22;10.77;10.26;9.72;9.2;8.58;8.01;4.87;2.51;0.7;0.23
110;0;7.09;7.97;8.84;9.74;10.09;11.85;12.75;13.84;14.99;16.02;16.33;17.1;17.83;17.99;18.32;18.14;17.81;17.84;17.6;17.3;17.05;16.83;16.53;16.03;15.59;15.03;14.37;10.26;6.41;2.32;0.86
120;0;6.59;7.42;8.3;9.1;9.56;10.83;11.6;13.1;13.87;14.66;15.75;16.67;17.63;18.43;19.62;20.17;20.6;21.12;21.55;21.75;21.91;22.07;21.9;21.58;21.29;20.92;20.29;16.47;12.03;5.49;2.26
129;0;6.14;6.93;7.83;8.52;9.09;9.89;10.57;12.42;12.87;13.43;15.23;16.16;17.08;18.07;19.48;20.35;21.22;21.93;22.85;23.44;23.98;24.55;24.59;24.55;24.51;24.46;24;21.56;17.75;9.64;4.25
130;0;6.07;6.87;7.76;8.44;9.02;9.8;10.48;12.29;12.73;13.27;15.08;16.03;16.97;17.96;19.36;20.25;21.15;21.88;22.82;23.44;24.03;24.6;24.66;24.68;24.67;24.64;24.24;22;18.33;10.11;4.5
135;0;5.72;6.57;7.36;8.02;8.65;9.38;10.11;11.52;11.97;12.55;13.85;15.31;16.31;17.33;18.54;19.48;20.35;21.28;22.3;23.08;24.09;24.63;24.69;24.78;24.79;24.91;24.82;23.74;20.98;12.39;5.78
136;0;5.66;6.5;7.28;7.93;8.57;9.3;10.04;11.34;11.82;12.42;13.62;15.06;16.17;17.2;18.35;19.29;20.15;21.12;22.15;22.96;24.07;24.6;24.67;24.76;24.75;24.85;24.81;23.98;21.45;12.8;6.03
139;0;5.42;6.31;6.92;7.67;8.34;9.08;9.86;10.86;11.32;12.03;12.99;14.3;15.73;16.76;17.76;18.71;19.53;20.6;21.66;22.54;23.92;24.44;24.53;24.64;24.58;24.65;24.67;24.47;22.68;13.79;6.73
140;0;5.35;6.22;6.79;7.59;8.26;9;9.8;10.72;11.16;11.89;12.79;14.06;15.5;16.62;17.57;18.51;19.32;20.43;21.49;22.4;23.84;24.36;24.46;24.58;24.51;24.57;24.61;24.56;23.02;14.08;6.96
141;0;5.29;6.12;6.67;7.48;8.18;8.93;9.74;10.57;11.02;11.77;12.62;13.82;15.26;16.47;17.38;18.32;19.04;20.28;21.31;22.07;23.53;24;24.21;24.29;24.43;24.48;24.55;24.61;23.33;14.31;7.18
142;0;5.23;6.02;6.57;7.39;8.1;8.86;9.67;10.43;10.88;11.64;12.45;13.59;15.03;16.24;17.14;18.06;18.77;19.98;21.01;21.75;23.18;23.65;23.86;23.95;24.34;24.39;24.48;24.61;23.61;14.54;7.4
143;0;5.16;5.93;6.45;7.3;8;8.78;9.54;10.27;10.75;11.5;12.28;13.36;14.8;16.01;16.9;17.81;18.5;19.69;20.72;21.43;22.84;23.31;23.52;23.61;24.05;24.27;24.41;24.57;23.85;14.8;7.6
144;0;5.09;5.83;6.33;7.23;7.92;8.66;9.41;10.13;10.62;11.39;12.13;13.13;14.57;15.78;16.65;17.56;18.24;19.41;20.43;21.12;22.5;22.97;23.19;23.28;23.73;24.08;24.33;24.49;24.04;15;7.8
145;0;5.02;5.73;6.23;7.15;7.85;8.55;9.28;9.98;10.51;11.27;11.98;12.92;14.35;15.56;16.42;17.31;17.97;19.13;20.14;20.82;22.17;22.64;22.87;22.96;23.42;23.81;24.23;24.41;24.19;15.14;7.98
146;0;4.96;5.64;6.12;7.07;7.77;8.43;9.15;9.84;10.38;11.16;11.83;12.71;14.12;15.35;16.19;17.07;17.72;18.86;19.86;20.51;21.84;22.31;22.56;22.65;23.12;23.48;23.94;24.33;24.3;15.3;8.16
148;0;4.82;5.45;5.91;6.9;7.59;8.21;8.89;9.55;10.14;10.89;11.55;12.29;13.7;14.92;15.74;16.6;17.23;18.32;19.3;19.91;21.2;21.67;21.95;22.05;22.53;22.87;23.38;24.13;24.39;15.52;8.46
149;0;4.76;5.35;5.81;6.78;7.49;8.09;8.78;9.42;10.01;10.76;11.41;12.1;13.48;14.71;15.52;16.36;16.98;18.06;19.03;19.63;20.89;21.37;21.67;21.77;22.26;22.61;23.12;23.98;24.37;15.57;8.58
150;0;4.69;5.26;5.7;6.67;7.37;7.96;8.64;9.26;9.86;10.6;11.24;11.89;13.26;14.48;15.29;16.11;16.73;17.79;18.74;19.33;20.55;21.04;21.37;21.48;21.98;22.34;22.83;23.69;24.11;15.6;8.67
155;0;4.33;4.74;5.16;6.16;6.79;7.33;7.96;8.51;9.15;9.81;10.4;10.85;12.14;13.37;14.1;14.87;15.48;16.42;17.3;17.8;18.88;19.39;19.86;20;20.54;20.97;21.45;22.25;22.77;15.38;8.89
160;0;4.09;4.41;4.83;5.77;6.39;6.94;7.55;8.04;8.67;9.28;9.83;10.24;11.46;12.69;13.39;14.11;14.73;15.6;16.41;16.87;17.85;18.4;18.97;19.15;19.72;20.2;20.65;21.35;21.84;14.95;8.74
162;0;4;4.29;4.69;5.62;6.23;6.77;7.38;7.86;8.48;9.07;9.6;10;11.18;12.42;13.1;13.81;14.43;15.27;16.06;16.5;17.43;18;18.62;18.81;19.39;19.89;20.33;20.99;21.48;14.76;8.61
168;0;3.74;3.93;4.35;5.15;5.75;6.31;6.93;7.34;7.92;8.45;8.95;9.35;10.44;11.68;12.32;12.99;13.63;14.39;15.11;15.5;16.31;16.94;17.7;17.92;18.53;19.08;19.49;19.99;20.46;14.34;8.3
170;0;3.69;3.85;4.27;5.04;5.65;6.22;6.82;7.23;7.8;8.31;8.8;9.22;10.27;11.51;12.15;12.81;13.45;14.19;14.9;15.28;16.06;16.7;17.51;17.73;18.34;18.91;19.31;19.77;20.22;14.24;8.24
174;0;3.57;3.69;4.11;4.83;5.43;6.01;6.62;7;7.55;8.03;8.5;8.93;9.95;11.19;11.81;12.45;13.11;13.81;14.48;14.84;15.57;16.24;17.11;17.35;17.98;18.57;18.95;19.33;19.77;14.03;8.13
180;0;3.51;3.6;4.03;4.71;5.31;5.91;6.51;6.88;7.41;7.88;8.33;8.79;9.78;11.02;11.63;12.26;12.93;13.61;14.26;14.61;15.31;15.99;16.9;17.15;17.79;18.39;18.77;19.09;19.52;13.87;8.07
</pre>
{{out}}
<pre>
The route taking the least time found was:
Point{2,Int64}[[1, 4], [1, 5], [2, 6], [3, 7], [4, 7], [5, 7], [6, 7], [7, 7], [8, 6], [8, 5], [9, 4]]
which has duration 4.0 hours, 43.697879668707344 minutes.
</pre>
=={{header|Nim}}==
{{trans|Go}}
The Go version runs in about 44 seconds on my computer. This Nim version, compiled in release mode (which includes runtime checks), runs in 22 seconds (18 seconds if link time optimization is activated). When compiled without checks (“danger” mode), it runs in 15.5 seconds.
<syntaxhighlight lang="nim">import hashes, math, parsecsv, sequtils, sets, strutils, sugar
type
MatrixF = seq[seq[float]]
Pred = float -> bool
# Structure that represents a polar CSV file's data.
# Note 0 degrees is directly into the wind, 180 degrees is directly downwind.
SailingPolar = object
winds: seq[float] # Vector of windspeeds.
degrees: seq[float] # Vector of angles in degrees of direction relative to the wind.
speeds: MatrixF # Matrix of sailing speeds indexed by wind velocity and angle of boat to wind.
# Structure that represents wind and surface current direction and velocity for a given position.
# Angles in degrees, velocities in knots.
SurfaceParameters = tuple[windDeg, windKts, currentDeg, currentKts: float]
proc getPolarData(filename: string): SailingPolar =
## Reads a sailing polar CSV file and returns a SailingPolar struct containing the file data.
## A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
## The first line of file contains labels for the wind velocities that make up columns, and
## the first entry of each row makes up a column of angle of sailing direction from wind in degrees.
var parser: CsvParser
parser.open(filename, separator = ';')
parser.readHeaderRow()
for col in 1..parser.headers.high:
result.winds.add parser.headers[col].parseFloat()
while parser.readRow():
if parser.row.len == 0: break # Ignore final blank line if there is one.
result.degrees.add parser.row[0].parseFloat()
result.speeds.add @[]
for col in 1..parser.row.high:
result.speeds[^1].add parser.row[col].parseFloat()
const R = 6372800.0 # Earth's approximate radius in meters.
template sind(d: float): float = sin(degToRad(d))
template cosd(d: float): float = cos(degToRad(d))
template asind(x: float): float = radToDeg(arcsin(x))
template atand(x, y: float): float = radToDeg(arctan2(x, y))
func haversine(lat1, long1, lat2, long2: float): (float, float) =
## Calculates the Haversine function for two points on the Earth's surface.
## Given two latitude, longitude pairs in degrees for a point on the Earth,
## get distance in meters and the initial direction of travel in degrees for
## movement from point 1 to point 2.
let dlat = lat2 - lat1
let dlong = long2 - long1
let a = sind(dlat/2)^2 + cosd(lat1) * cosd(lat2) * sind(dlong/2)^2
let c = 2 * asind(sqrt(a))
var theta = atand(sind(dlong) * cosd(lat2),
cosd(lat1) * sind(lat2) - sind(lat1) * cosd(lat2) * cosd(dlong))
theta = (theta + 360) mod 360
result = (R * c * 0.5399565, theta)
func findFirst(a: seq[float]; p: Pred): int =
## Returns the index of the first element of 'a' for which 'p' returns true or -1 otherwise.
for i in 0..a.high:
if p(a[i]): return i
result = -1
func findLast(a: seq[float]; p: Pred): int =
## Returns the index of the last element of 'a' for which 'p' returns true or -1 otherwise.
for i in countdown(a.high, 0):
if p(a[i]): return i
result = -1
func boatSpeed(sp: SailingPolar; pointOfSail, windSpeed: float): float =
## Calculate the expected sailing speed in a specified direction in knots,
## given sailing polar data, a desired point of sail in degrees, and wind speed in knots.
let
udeg = sp.degrees.findLast(t => t <= pointOfSail)
odeg = sp.degrees.findFirst(t => t >= pointOfSail)
uvel = sp.winds.findLast(t => t <= windSpeed)
ovel = sp.winds.findFirst(t => t >= windSpeed)
if udeg == -1 or odeg == -1 or uvel == -1 or ovel == -1: return -1
let frac = if odeg == udeg and uvel == ovel:
1.0
elif odeg == udeg:
(windSpeed - sp.winds[uvel]) / (sp.winds[ovel] - sp.winds[uvel])
elif uvel == ovel:
(pointOfSail - sp.degrees[udeg]) / (sp.degrees[odeg] - sp.degrees[udeg])
else:
((pointOfSail - sp.degrees[udeg]) / (sp.degrees[odeg] - sp.degrees[udeg]) +
(windSpeed - sp.winds[uvel]) / (sp.winds[ovel] - sp.winds[uvel])) / 2
result = sp.speeds[udeg][uvel] + frac * (sp.speeds[odeg][ovel] - sp.speeds[udeg][uvel])
func sailingSpeed(sp: SailingPolar; azimuth, pointos, ws: float): float =
## Calculates the expected net boat speed in a desired direction versus the wind ('azimuth').
## This is generally different from the actual boat speed in its actual direction.
## Directions are in degrees ('pointos' is point of sail the ship direction from the wind),
## and velocity of wind ('ws') is in knots.
sp.boatSpeed(pointos, ws) * cosd(abs(pointos - azimuth))
func bestVectorSpeed(sp: SailingPolar;
dirTravel, dirWind, windSpeed, dirCur, velCur: float): (float, float) =
## Calculates the net direction and velocity of a sailing ship.
## Arguments are sailing polar data, direction of travel in degrees from north, wind direction in
## degrees from north, wind velocity in knots, surface current direction in degrees, and
## current velocity in knots.
var azimuth = (dirTravel - dirWind) mod 360
if azimuth < 0: azimuth += 360
if azimuth > 180: azimuth = 360 - azimuth
var vmg = sp.boatSpeed(azimuth, windSpeed)
var other = -1.0
var idx = -1
for i, d in sp.degrees:
let ss = sp.sailingSpeed(azimuth, d, windSpeed)
if ss > other:
other = ss
idx = i
if other > vmg:
azimuth = sp.degrees[idx]
vmg = other
let
dirChosen = dirWind + azimuth
wx = vmg * sind(dirChosen)
wy = vmg * cosd(dirChosen)
curX = velCur * sind(dirCur)
curY = velCur * cosd(dirCur)
result = (atand(wy + curY, wx + curX), sqrt((wx + curX)^2 + (wy + curY)^2))
func sailSegmentTime(sp: SailingPolar; p: SurfaceParameters;
lat1, long1, lat2, long2: float): float =
## Calculates the trip time in minutes from (lat1, long1) to the destination (lat2, long2).
## Uses the data in SurfaceParameters for wind and current velocity and direction.
let (distance, dir) = haversine(lat1, long1, lat2, long2)
let (_, vel) = sp.bestVectorSpeed(dir, p.windDeg, p.windKts, p.currentDeg, p.currentKts)
## minutes/s * m / (knots * (m/s / knot)) = minutes
result = (1 / 60) * distance / (vel * 1.94384)
# Structure that represents a point in 2-D space.
type Point2 = tuple[x, y: int]
func `+`(p1, p2: Point2): Point2 = (p1.x + p2.x, p1.y + p2.y)
func `$`(p: Point2): string = "($1, $2)".format(p.x, p.y)
type
# Tuple that consists of a tuple of latitude and longitude in degrees.
# NB: This uses latitude (often considered to be y) first then longitude (often considered to be x).
# This latitude, then longitude ordering is as per ISO 6709 (en.wikipedia.org/wiki/ISO_6709).
Position = tuple[lat, long: float]
# Tuple that represents a Position with the SurfaceParameters of wind and current at the Position.
GridPoint = tuple[pt: Position; sp: SurfaceParameters]
MatrixG = seq[seq[GridPoint]]
# Type alias for a matrix of GridPoints, each Position point with their SurfaceParameters.
# A vector of TimeSlice can give the surface characteristics for an ocean region over time.
TimeSlice = MatrixG
# Structure that represents a routing problem.
RoutingProblem = object
timeInterval: float # The minutes duration for each TimeSlice.
timeFrame: seq[TimeSlice] # A vector of sequential timeslices for the ocean region.
obstacleIndices: seq[Point2] # The cartesian indices in each TimeSlice for obstacles
# such as land or shoals, where the ship may not go.
startIndex: int # The TimeSlice position for time of starting.
start: Point2 # Starting location on grid of GridPoints.
finish: Point2 # Destination / finish location on grid of GridPoints.
allowRepeatVisits: bool # Whether the vessel may overlap its prior path, usually false.
# Structure that represents a timed path.
TimedPath = object
duration: float # Minutes total to travel the path.
path: seq[Point2] # Vector of cartesian indices of points in grid for path to travel.
func hash(t: TimedPath): Hash =
## Hash function to allow building a set of TimedPath values.
result = t.duration.hash !& t.path.hash
result = !$result
const Neighbors: seq[Point2] = @[(-1, -1), (-1, 0), (-1, 1), (0, -1),
( 0, 1), ( 1, -1), ( 1, 0), (1, 1)]
func surround(p: Point2; mat: TimeSlice; excluded: openArray[Point2]): seq[Point2] =
## Returns a list of points surrounding 'p' which are not otherwise excluded.
let xmax = mat.len
let ymax = mat[0].len
for x in Neighbors:
let q = p + x
if q.x >= 0 and q.x < xmax and q.y >= 0 and q.y < ymax and q notin excluded:
result.add q
proc minimumTimeRoute(rp: RoutingProblem; sp: SailingPolar; verbose: bool): TimedPath =
## Get the route (as a TimedPath) that minimizes time from start to finish for a given
## RoutingProblem (sea parameters) given sailing polar data (ship parameters).
var timedPaths = @[TimedPath(duration: 0, path: @[rp.start])]
var completed = false
result = TimedPath(duration: 1000)
for _ in 1..1000:
var newPaths: seq[TimedPath]
if verbose:
echo "Checking $1 paths of length $2".format(timedPaths.len, timedPaths[0].path.len)
for tpath in timedPaths:
if tpath.path[^1] == rp.finish:
completed = true
newPaths.add tpath
else:
let p1 = tpath.path[^1]
let num = tpath.duration.toInt
let den = rp.timeInterval.toInt
let slice = rp.timeFrame[num div den mod rp.timeFrame.len]
for p2 in p1.surround(slice, rp.obstacleIndices):
if not rp.allowRepeatVisits and p2 in tpath.path:
continue
let gp1 = slice[p1.x][p1.y]
let gp2 = slice[p2.x][p2.y]
let (lat1, long1) = gp1.pt
let (lat2, long2) = gp2.pt
let t = sp.sailSegmentTime(gp1.sp, lat1, long1, lat2, long2)
let path = tpath.path & p2
newPaths.add TimedPath(duration: tpath.duration + t, path: path)
timedPaths = newPaths.toHashSet().toSeq()
if completed:
var durations = collect(newSeq, for t in timedPaths: t.duration)
let minDur = min(durations)
let finished = collect(newSeq):
for t in timedPaths:
if t.path[^1] == rp.finish: t
durations = collect(newSeq, for f in finished: f.duration)
let idx = minIndex(durations)
let minFinDur = durations[idx]
if verbose:
echo "Current finished minimum: $1, others $2".format(finished[idx], minDur)
if minDur == minFinDur:
result = finished[idx]
break
#[ The data is selected so the best time path is slightly longer than the
shortest length path. The forbidden regions are x, representing land or reef.
The allowed sailing points are . and start and finish are S and F.
x . . F . . x . x
. . . . . . . x x
x . . x x x . . .
. . x x x x . x x
x . . . x x . x .
x . . . x x . x .
. . . . x . . x .
x . . . . . . x .
. . . S . . . . .
]#
const Forbidden: seq[Point2] = @[(0, 7), (1, 0), (1, 7), (2, 4), (2, 7), (3, 0), (3, 4),
(3, 5), (3, 7), (4, 0), (4, 4), (4, 5), (4, 7), (5, 2),
(5, 3), (5, 4), (5, 5), (5, 7), (5, 8), (6, 0), (6, 3),
(6, 4), (6, 5), (7, 7), (7, 8), (8, 0), (8, 6), (8, 8)]
func surfaceByLongitude(long: float): SurfaceParameters =
## Create regional wind patterns on the map.
if long < -155.03:
(-5.0, 8.0, 150.0, 0.5)
elif long < -155.99:
(-90.0, 20.0, 150.0, 0.4)
else:
(180.0, 25.0, 150.0, 0.3)
func mutateTimeSlices(slices: var seq[TimeSlice]) =
var i = 1
for slice in slices.mitems:
for j in 0..slice.high:
for x in slice[j].mitems:
x.sp = (x.sp.windDeg, x.sp.windKts * (1 + 0.002 * float64(i)),
x.sp.currentDeg, x.sp.currentKts)
inc i
let startPos: Point2 = (0, 3)
let endPos: Point2 = (8, 3)
var slices = newSeq[MatrixG](200)
for s in 0..slices.high:
var gpoints = newSeq[seq[GridPoint]](9)
for i in 0..<9:
gpoints[i].setLen(9)
for j in 0..<9:
let pt: Position = (19.78 - 1/60 + i/60, -155.0 - 5/60 + j/60)
gpoints[i][j] = (pt, surfaceByLongitude(pt.long))
slices[s] = move(gpoints)
slices.mutateTimeSlices()
let routeProb = RoutingProblem(timeInterval: 10, timeFrame: slices,
obstacleIndices: Forbidden, startIndex: 0,
start: startPos, finish: endPos, allowRepeatVisits: false)
let fileName = "polar.csv"
let sp = getPolarData(fileName)
let tp = routeProb.minimumTimeRoute(sp, false)
echo "The route taking the least time found was:"
echo tp.path
echo "which has duration ", int(tp.duration / 60), " hours ", toInt(tp.duration mod 60), " minutes."</syntaxhighlight>
{{out}}
<pre>The route taking the least time found was:
@[(0, 3), (0, 4), (1, 5), (2, 6), (3, 6), (4, 6), (5, 6), (6, 6), (7, 5), (7, 4), (8, 3)]
which has duration 4 hours 44 minutes.</pre>
=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #000080;font-style:italic;">--
-- demo\rosetta\Weather_Routing.exw
-- ================================
--</span>
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">polar_csv</span> <span style="color: #0000FF;">=</span> <span style="font-size: 75%; color: #008000;">"""
TWA\TWS;0;4;5;6;7;8;9;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;35;40;60;70
40;0;0.53;0.54;0.49;0.4;0.31;0.21;0.16;0.11;0.08;0.05;0.03;0.02;0.01;0;0;0;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.05;-0.1;-0.11
41;0;0.61;0.62;0.56;0.47;0.36;0.25;0.19;0.14;0.1;0.07;0.04;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;0;0;-0.04;-0.09;-0.1
44;0;0.89;0.91;0.82;0.69;0.56;0.42;0.33;0.24;0.18;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;-0.02;-0.06;-0.06
45;0;0.99;1.02;0.92;0.78;0.64;0.49;0.39;0.29;0.22;0.15;0.1;0.07;0.04;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.05;-0.05
46;0;1.11;1.14;1.02;0.87;0.73;0.57;0.45;0.35;0.26;0.18;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;0;-0.01;-0.04;-0.05
47;0;1.23;1.25;1.14;0.97;0.82;0.66;0.53;0.41;0.31;0.22;0.15;0.1;0.07;0.04;0.02;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.01;-0.03;-0.04
48;0;1.37;1.37;1.26;1.08;0.93;0.76;0.61;0.48;0.36;0.26;0.19;0.13;0.08;0.05;0.03;0.02;0.01;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.02;-0.03
49;0;1.5;1.5;1.39;1.2;1.05;0.87;0.71;0.56;0.42;0.31;0.22;0.15;0.1;0.07;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;0;-0.02;-0.02
50;0;1.65;1.64;1.52;1.33;1.18;1;0.81;0.65;0.49;0.37;0.26;0.19;0.13;0.08;0.05;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;0;0;-0.01;-0.02
51;0;1.79;1.77;1.67;1.46;1.32;1.13;0.92;0.74;0.57;0.43;0.31;0.22;0.15;0.1;0.07;0.05;0.03;0.02;0.02;0.01;0.01;0;0;0;0;0;0;0;0;-0.01;-0.02
53;0;2.1;2.07;1.99;1.76;1.62;1.4;1.14;0.95;0.74;0.57;0.43;0.31;0.22;0.16;0.1;0.08;0.06;0.04;0.03;0.02;0.01;0.01;0.01;0;0;0;0;0;0;-0.01;-0.01
54;0;2.26;2.22;2.16;1.92;1.78;1.55;1.28;1.06;0.84;0.65;0.5;0.37;0.27;0.19;0.13;0.1;0.07;0.06;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;0;-0.01
55;0;2.43;2.39;2.34;2.09;1.95;1.7;1.42;1.18;0.95;0.74;0.57;0.43;0.32;0.23;0.16;0.12;0.09;0.07;0.05;0.04;0.03;0.02;0.01;0.01;0;0;0;0;0;0;-0.01
60;0;3.29;3.33;3.33;3.08;2.93;2.64;2.29;1.98;1.66;1.36;1.1;0.88;0.68;0.53;0.39;0.32;0.26;0.21;0.17;0.13;0.1;0.07;0.05;0.04;0.03;0.02;0.01;0;0;0;0
70;0;5.2;5.53;5.74;5.59;5.5;5.22;4.84;4.46;3.94;3.51;3.08;2.65;2.26;1.9;1.55;1.38;1.22;1.06;0.92;0.78;0.66;0.55;0.46;0.37;0.3;0.24;0.18;0.03;0;0;0
80;0;6.8;7.43;7.97;8.02;8.23;8.34;8.2;7.9;7.37;6.91;6.43;5.9;5.32;4.72;4.12;3.83;3.55;3.25;2.96;2.67;2.4;2.13;1.88;1.65;1.43;1.22;1.04;0.37;0.09;0.01;0
90;0;7.59;8.5;9.4;9.73;10.4;11.16;11.53;11.56;11.3;11.05;10.77;10.44;9.83;9.07;8.34;8;7.65;7.27;6.88;6.46;6.04;5.61;5.15;4.74;4.33;3.88;3.51;1.72;0.67;0.12;0.03
100;0;7.34;8.25;9.16;9.86;10.5;11.95;12.79;13.5;14.02;14.4;14.37;14.5;14.4;13.92;13.52;13.19;12.79;12.51;12.1;11.66;11.22;10.77;10.26;9.72;9.2;8.58;8.01;4.87;2.51;0.7;0.23
110;0;7.09;7.97;8.84;9.74;10.09;11.85;12.75;13.84;14.99;16.02;16.33;17.1;17.83;17.99;18.32;18.14;17.81;17.84;17.6;17.3;17.05;16.83;16.53;16.03;15.59;15.03;14.37;10.26;6.41;2.32;0.86
120;0;6.59;7.42;8.3;9.1;9.56;10.83;11.6;13.1;13.87;14.66;15.75;16.67;17.63;18.43;19.62;20.17;20.6;21.12;21.55;21.75;21.91;22.07;21.9;21.58;21.29;20.92;20.29;16.47;12.03;5.49;2.26
129;0;6.14;6.93;7.83;8.52;9.09;9.89;10.57;12.42;12.87;13.43;15.23;16.16;17.08;18.07;19.48;20.35;21.22;21.93;22.85;23.44;23.98;24.55;24.59;24.55;24.51;24.46;24;21.56;17.75;9.64;4.25
130;0;6.07;6.87;7.76;8.44;9.02;9.8;10.48;12.29;12.73;13.27;15.08;16.03;16.97;17.96;19.36;20.25;21.15;21.88;22.82;23.44;24.03;24.6;24.66;24.68;24.67;24.64;24.24;22;18.33;10.11;4.5
135;0;5.72;6.57;7.36;8.02;8.65;9.38;10.11;11.52;11.97;12.55;13.85;15.31;16.31;17.33;18.54;19.48;20.35;21.28;22.3;23.08;24.09;24.63;24.69;24.78;24.79;24.91;24.82;23.74;20.98;12.39;5.78
136;0;5.66;6.5;7.28;7.93;8.57;9.3;10.04;11.34;11.82;12.42;13.62;15.06;16.17;17.2;18.35;19.29;20.15;21.12;22.15;22.96;24.07;24.6;24.67;24.76;24.75;24.85;24.81;23.98;21.45;12.8;6.03
139;0;5.42;6.31;6.92;7.67;8.34;9.08;9.86;10.86;11.32;12.03;12.99;14.3;15.73;16.76;17.76;18.71;19.53;20.6;21.66;22.54;23.92;24.44;24.53;24.64;24.58;24.65;24.67;24.47;22.68;13.79;6.73
140;0;5.35;6.22;6.79;7.59;8.26;9;9.8;10.72;11.16;11.89;12.79;14.06;15.5;16.62;17.57;18.51;19.32;20.43;21.49;22.4;23.84;24.36;24.46;24.58;24.51;24.57;24.61;24.56;23.02;14.08;6.96
141;0;5.29;6.12;6.67;7.48;8.18;8.93;9.74;10.57;11.02;11.77;12.62;13.82;15.26;16.47;17.38;18.32;19.04;20.28;21.31;22.07;23.53;24;24.21;24.29;24.43;24.48;24.55;24.61;23.33;14.31;7.18
142;0;5.23;6.02;6.57;7.39;8.1;8.86;9.67;10.43;10.88;11.64;12.45;13.59;15.03;16.24;17.14;18.06;18.77;19.98;21.01;21.75;23.18;23.65;23.86;23.95;24.34;24.39;24.48;24.61;23.61;14.54;7.4
143;0;5.16;5.93;6.45;7.3;8;8.78;9.54;10.27;10.75;11.5;12.28;13.36;14.8;16.01;16.9;17.81;18.5;19.69;20.72;21.43;22.84;23.31;23.52;23.61;24.05;24.27;24.41;24.57;23.85;14.8;7.6
144;0;5.09;5.83;6.33;7.23;7.92;8.66;9.41;10.13;10.62;11.39;12.13;13.13;14.57;15.78;16.65;17.56;18.24;19.41;20.43;21.12;22.5;22.97;23.19;23.28;23.73;24.08;24.33;24.49;24.04;15;7.8
145;0;5.02;5.73;6.23;7.15;7.85;8.55;9.28;9.98;10.51;11.27;11.98;12.92;14.35;15.56;16.42;17.31;17.97;19.13;20.14;20.82;22.17;22.64;22.87;22.96;23.42;23.81;24.23;24.41;24.19;15.14;7.98
146;0;4.96;5.64;6.12;7.07;7.77;8.43;9.15;9.84;10.38;11.16;11.83;12.71;14.12;15.35;16.19;17.07;17.72;18.86;19.86;20.51;21.84;22.31;22.56;22.65;23.12;23.48;23.94;24.33;24.3;15.3;8.16
148;0;4.82;5.45;5.91;6.9;7.59;8.21;8.89;9.55;10.14;10.89;11.55;12.29;13.7;14.92;15.74;16.6;17.23;18.32;19.3;19.91;21.2;21.67;21.95;22.05;22.53;22.87;23.38;24.13;24.39;15.52;8.46
149;0;4.76;5.35;5.81;6.78;7.49;8.09;8.78;9.42;10.01;10.76;11.41;12.1;13.48;14.71;15.52;16.36;16.98;18.06;19.03;19.63;20.89;21.37;21.67;21.77;22.26;22.61;23.12;23.98;24.37;15.57;8.58
150;0;4.69;5.26;5.7;6.67;7.37;7.96;8.64;9.26;9.86;10.6;11.24;11.89;13.26;14.48;15.29;16.11;16.73;17.79;18.74;19.33;20.55;21.04;21.37;21.48;21.98;22.34;22.83;23.69;24.11;15.6;8.67
155;0;4.33;4.74;5.16;6.16;6.79;7.33;7.96;8.51;9.15;9.81;10.4;10.85;12.14;13.37;14.1;14.87;15.48;16.42;17.3;17.8;18.88;19.39;19.86;20;20.54;20.97;21.45;22.25;22.77;15.38;8.89
160;0;4.09;4.41;4.83;5.77;6.39;6.94;7.55;8.04;8.67;9.28;9.83;10.24;11.46;12.69;13.39;14.11;14.73;15.6;16.41;16.87;17.85;18.4;18.97;19.15;19.72;20.2;20.65;21.35;21.84;14.95;8.74
162;0;4;4.29;4.69;5.62;6.23;6.77;7.38;7.86;8.48;9.07;9.6;10;11.18;12.42;13.1;13.81;14.43;15.27;16.06;16.5;17.43;18;18.62;18.81;19.39;19.89;20.33;20.99;21.48;14.76;8.61
168;0;3.74;3.93;4.35;5.15;5.75;6.31;6.93;7.34;7.92;8.45;8.95;9.35;10.44;11.68;12.32;12.99;13.63;14.39;15.11;15.5;16.31;16.94;17.7;17.92;18.53;19.08;19.49;19.99;20.46;14.34;8.3
170;0;3.69;3.85;4.27;5.04;5.65;6.22;6.82;7.23;7.8;8.31;8.8;9.22;10.27;11.51;12.15;12.81;13.45;14.19;14.9;15.28;16.06;16.7;17.51;17.73;18.34;18.91;19.31;19.77;20.22;14.24;8.24
174;0;3.57;3.69;4.11;4.83;5.43;6.01;6.62;7;7.55;8.03;8.5;8.93;9.95;11.19;11.81;12.45;13.11;13.81;14.48;14.84;15.57;16.24;17.11;17.35;17.98;18.57;18.95;19.33;19.77;14.03;8.13
180;0;3.51;3.6;4.03;4.71;5.31;5.91;6.51;6.88;7.41;7.88;8.33;8.79;9.78;11.02;11.63;12.26;12.93;13.61;14.26;14.61;15.31;15.99;16.9;17.15;17.79;18.39;18.77;19.09;19.52;13.87;8.07
"""</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">to_numbers</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">to_number</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">s</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">getpolardata</span><span style="color: #0000FF;">(</span><span style="color: #004080;">string</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
-- The first line of the file contains labels for the wind velocities that make up columns, and
-- the first entry of each row makes up a column of angle of sailing direction from wind in degrees
--</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">lines</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split_any</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\r\n"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">winds</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">to_numbers</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span><span style="color: #008000;">";"</span><span style="color: #0000FF;">)[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$]),</span>
<span style="color: #000000;">degrees</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{},</span> <span style="color: #000000;">speeds</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">l</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">to_numbers</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lines</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span><span style="color: #008000;">";"</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">l</span><span style="color: #0000FF;">)!=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">winds</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">degrees</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">degrees</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">])</span>
<span style="color: #000000;">speeds</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">speeds</span><span style="color: #0000FF;">,</span><span style="color: #000000;">l</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$])</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">winds</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">speeds</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000080;font-style:italic;">--
-- winds is a list of wind speeds
-- degrees is a list of angles in degrees of direction relative to the wind
-- (note 0 degrees is directly into the wind, 180 degrees is directly downwind)
-- each speeds[i] is an array of length(winds) for each degrees[i]
--</span>
<span style="color: #008080;">constant</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">winds</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">speeds</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">getpolardata</span><span style="color: #0000FF;">(</span><span style="color: #000000;">polar_csv</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--constant {winds, degrees, speeds} = getpolardata(get_text("polar.csv")) -- alt</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">R</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">6372800</span> <span style="color: #000080;font-style:italic;">-- Earth's approximate radius in meters</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">timeinterval</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">10</span> <span style="color: #000080;font-style:italic;">-- the minutes duration for each TimeSlice</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">deg2rad</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">deg</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">/</span><span style="color: #000000;">180</span><span style="color: #0000FF;">+</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">*</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">rad2deg</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">rad</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">remainder</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">rad</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">180</span><span style="color: #0000FF;">/</span><span style="color: #004600;">PI</span><span style="color: #0000FF;">)+</span><span style="color: #000000;">360</span><span style="color: #0000FF;">,</span><span style="color: #000000;">360</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">deg</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg2rad</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">deg</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg2rad</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">asind</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">deg</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">rad2deg</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">arcsin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">deg</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">atand</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #000000;">rad2deg</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">atan2</span><span style="color: #0000FF;">(</span><span style="color: #000000;">x</span><span style="color: #0000FF;">,</span><span style="color: #000000;">y</span><span style="color: #0000FF;">))</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">haversine</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lat2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon2</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Calculate the haversine function for two points on the Earth's surface.
--
-- Given two latitude, longitude pairs in degrees for a point on the Earth,
-- get distance in meters and the initial direction of travel in degrees for
-- movement from point 1 to point 2.
--</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">dlat</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lat2</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">dlon</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">lon2</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">a</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlat</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat2</span><span style="color: #0000FF;">)*</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlon</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span><span style="color: #000000;">2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">c</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2.0</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">asind</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)),</span>
<span style="color: #000000;">theta</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">atand</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlon</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat2</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">sind</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat2</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dlon</span><span style="color: #0000FF;">))</span>
<span style="color: #000000;">theta</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">theta</span><span style="color: #0000FF;">+</span><span style="color: #000000;">360</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">360</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">R</span><span style="color: #0000FF;">*</span><span style="color: #000000;">c</span><span style="color: #0000FF;">*</span><span style="color: #000000;">0.5399565</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">theta</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">find_range</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- Returns the indexes of s of the first >=v, and the last <=v</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]>=</span><span style="color: #000000;">v</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">to</span> <span style="color: #000000;">1</span> <span style="color: #008080;">by</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]<=</span><span style="color: #000000;">v</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span><span style="color: #000000;">j</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">exit</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">boatspeed</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">windspeed</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Calculate the expected sailing speed in a specified direction in knots,
-- given a desired point of sail in degrees, and wind speed in knots (for
-- the previously loaded sailing polar data)
--</span>
<span style="color: #004080;">integer</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">ld</span><span style="color: #0000FF;">,</span><span style="color: #000000;">ud</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">find_range</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">),</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">lv</span><span style="color: #0000FF;">,</span><span style="color: #000000;">uv</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">find_range</span><span style="color: #0000FF;">(</span><span style="color: #000000;">windspeed</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">winds</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">ld</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ud</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lv</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">uv</span><span style="color: #0000FF;">})</span> <span style="color: #008080;">then</span> <span style="color: #008080;">return</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">wu</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">winds</span><span style="color: #0000FF;">[</span><span style="color: #000000;">uv</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">wl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">winds</span><span style="color: #0000FF;">[</span><span style="color: #000000;">lv</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">du</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ud</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">dl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ld</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">f</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">ld</span><span style="color: #0000FF;">==</span><span style="color: #000000;">ud</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">uv</span><span style="color: #0000FF;">==</span><span style="color: #000000;">lv</span> <span style="color: #0000FF;">?</span> <span style="color: #000000;">1</span> <span style="color: #0000FF;">:</span>
<span style="color: #0000FF;">(</span><span style="color: #000000;">wu</span><span style="color: #0000FF;">-</span><span style="color: #000000;">windspeed</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">wu</span><span style="color: #0000FF;">-</span><span style="color: #000000;">wl</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">elsif</span> <span style="color: #000000;">uv</span><span style="color: #0000FF;">==</span><span style="color: #000000;">lv</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">du</span><span style="color: #0000FF;">-</span><span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">du</span><span style="color: #0000FF;">-</span><span style="color: #000000;">dl</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">else</span>
<span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">((</span><span style="color: #000000;">du</span><span style="color: #0000FF;">-</span><span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">du</span><span style="color: #0000FF;">-</span><span style="color: #000000;">dl</span><span style="color: #0000FF;">)+</span>
<span style="color: #0000FF;">(</span><span style="color: #000000;">wu</span><span style="color: #0000FF;">-</span><span style="color: #000000;">windspeed</span><span style="color: #0000FF;">)/(</span><span style="color: #000000;">wu</span><span style="color: #0000FF;">-</span><span style="color: #000000;">wl</span><span style="color: #0000FF;">))/</span><span style="color: #000000;">2</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">su</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">speeds</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ud</span><span style="color: #0000FF;">,</span><span style="color: #000000;">uv</span><span style="color: #0000FF;">],</span>
<span style="color: #000000;">sl</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">speeds</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ld</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lv</span><span style="color: #0000FF;">],</span>
<span style="color: #000080;font-style:italic;">-- res = su + f*(sl-su) -- (original)</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">su</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">su</span><span style="color: #0000FF;">-</span><span style="color: #000000;">sl</span><span style="color: #0000FF;">)</span> <span style="color: #000080;font-style:italic;">-- (equivalent)
-- res = sl + (1-f)*(su-sl) -- (also equivalent)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">sailingspeed</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">azimuth</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ws</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Calculate the expected net boat speed in a desired direction versus the wind (azimuth).
-- This is generally different from the actual boat speed in its actual direction.
-- Directions are in degrees (pointofsail is the ship direction from wind),
-- and velocity of wind (ws) is in knots.
--</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">boatspeed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ws</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">cosd</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">pointofsail</span><span style="color: #0000FF;">-</span><span style="color: #000000;">azimuth</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">bestvectorspeed</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">dirtravel</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">sequence</span> <span style="color: #000000;">surfaceparameters</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Calculate the net direction and velocity of a sailing ship.
--</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">winddirection</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">windvelocity</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">currentdirection</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">currentvelocity</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">surfaceparameters</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">azimuth</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dirtravel</span><span style="color: #0000FF;">-</span><span style="color: #000000;">winddirection</span><span style="color: #0000FF;">,</span><span style="color: #000000;">360</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">azimuth</span><span style="color: #0000FF;"><</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #000000;">azimuth</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">360</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">azimuth</span><span style="color: #0000FF;">></span><span style="color: #000000;">180</span> <span style="color: #008080;">then</span> <span style="color: #000000;">azimuth</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">360</span><span style="color: #0000FF;">-</span><span style="color: #000000;">azimuth</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">vmg</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">boatspeed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">azimuth</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">windvelocity</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">other</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">idx</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">degrees</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">ss</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sailingspeed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">azimuth</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">windvelocity</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">ss</span><span style="color: #0000FF;">></span><span style="color: #000000;">other</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">other</span><span style="color: #0000FF;">,</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">ss</span><span style="color: #0000FF;">,</span><span style="color: #000000;">i</span><span style="color: #0000FF;">}</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">other</span><span style="color: #0000FF;">></span><span style="color: #000000;">vmg</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">azimuth</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">degrees</span><span style="color: #0000FF;">[</span><span style="color: #000000;">idx</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">vmg</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">other</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">dirchosen</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">deg2rad</span><span style="color: #0000FF;">(</span><span style="color: #000000;">winddirection</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">azimuth</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">dircurrent</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">deg2rad</span><span style="color: #0000FF;">(</span><span style="color: #000000;">currentdirection</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">wx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">vmg</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dirchosen</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">wy</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">vmg</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dirchosen</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">curx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">currentvelocity</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dircurrent</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">cury</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">currentvelocity</span> <span style="color: #0000FF;">*</span> <span style="color: #7060A8;">cos</span><span style="color: #0000FF;">(</span><span style="color: #000000;">dircurrent</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">wx</span><span style="color: #0000FF;">+</span><span style="color: #000000;">curx</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span> <span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">wy</span><span style="color: #0000FF;">+</span><span style="color: #000000;">cury</span><span style="color: #0000FF;">,</span><span style="color: #000000;">2</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">sailsegmenttime</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">surfaceparameters</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lat2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon2</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Calculate the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2).
-- Uses the data in surfaceparameters for wind and current velocity and direction.
--</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">distance</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">direction</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">haversine</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lat2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon2</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">velocity</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">bestvectorspeed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">direction</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">surfaceparameters</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- minutes/s * m / (knots * (m/s / knot)) = minutes</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">60</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">distance</span> <span style="color: #0000FF;">/</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">velocity</span> <span style="color: #0000FF;">*</span> <span style="color: #000000;">1.94384</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #000080;font-style:italic;">--
-- The data is selected so the best time path is slightly longer than the
-- shortest length path. The forbidden regions are x, representing land or reef.
-- The allowed sailing points are . and start and finish are S and F.
-- </span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">chart</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">split</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"""
...S.....
x......x.
....x..x.
x...xx.x.
x...xx.x.
..xxxx.xx
x..xxx...
.......xx
x..F..x.x"""</span><span style="color: #0000FF;">,</span><span style="color: #008000;">'\n'</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">minimum_time_route</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">timeframe</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">start</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">finish</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">--
-- Get the fastest route from start to finish for some detailed sea/ship parameters.
-- timeframe is a massive 200 * 9x9 * {pt,surfaceparameters}
-- note that polar data (ie winds, degrees, speeds) is static here, for simplicity.
--</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">(),</span>
<span style="color: #000000;">mintime</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1000.0</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">xmax</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">chart</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]),</span>
<span style="color: #000000;">ymax</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">chart</span><span style="color: #0000FF;">),</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">start</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">todo</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">start</span><span style="color: #0000FF;">},</span>
<span style="color: #000000;">costs</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xmax</span><span style="color: #0000FF;">),</span><span style="color: #000000;">ymax</span><span style="color: #0000FF;">),</span> <span style="color: #000080;font-style:italic;">-- (lowest durations)</span>
<span style="color: #000000;">paths</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">xmax</span><span style="color: #0000FF;">),</span><span style="color: #000000;">ymax</span><span style="color: #0000FF;">),</span> <span style="color: #000080;font-style:italic;">-- (single backlinks)</span>
<span style="color: #000000;">minpath</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span>
<span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">while</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">todo</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">todo</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">todo</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">todo</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">..$]</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">duration</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">sdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">round</span><span style="color: #0000FF;">(</span><span style="color: #000000;">duration</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">timeinterval</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">timeframe</span><span style="color: #0000FF;">))+</span><span style="color: #000000;">1</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">timeframe</span><span style="color: #0000FF;">[</span><span style="color: #000000;">sdx</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">nx</span><span style="color: #0000FF;">=</span><span style="color: #000000;">px</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">px</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">ny</span><span style="color: #0000FF;">=</span><span style="color: #000000;">py</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">py</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">px</span> <span style="color: #008080;">or</span> <span style="color: #000000;">ny</span><span style="color: #0000FF;">!=</span><span style="color: #000000;">py</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">and</span> <span style="color: #000000;">nx</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">nx</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">xmax</span>
<span style="color: #008080;">and</span> <span style="color: #000000;">ny</span><span style="color: #0000FF;">>=</span><span style="color: #000000;">1</span> <span style="color: #008080;">and</span> <span style="color: #000000;">ny</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">ymax</span>
<span style="color: #008080;">and</span> <span style="color: #000000;">chart</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]!=</span><span style="color: #008000;">'x'</span> <span style="color: #008080;">then</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">gp1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">],</span> <span style="color: #000080;font-style:italic;">-- {pt,surfaceparameters}</span>
<span style="color: #000000;">gp2</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]</span> <span style="color: #000080;font-style:italic;">-- ""</span>
<span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gp1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">],</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">lat2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon2</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gp2</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">surfaceparameters</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gp1</span><span style="color: #0000FF;">[</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">nt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">duration</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">sailsegmenttime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">surfaceparameters</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lat1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon1</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lat2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lon2</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]=-</span><span style="color: #000000;">1</span> <span style="color: #008080;">or</span> <span style="color: #000000;">nt</span><span style="color: #0000FF;"><</span><span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #000080;font-style:italic;">-- a larger (than 9x9) simulation might benefit from not
-- putting any already-too-long routes back on the todo
-- list and/or processing todo lowest duration first.</span>
<span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">nt</span>
<span style="color: #000000;">paths</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">if</span> <span style="color: #008080;">not</span> <span style="color: #7060A8;">find</span><span style="color: #0000FF;">({</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">},</span><span style="color: #000000;">todo</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">todo</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">todo</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">elsif</span> <span style="color: #000000;">nt</span><span style="color: #0000FF;">==</span><span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">ny</span><span style="color: #0000FF;">,</span><span style="color: #000000;">nx</span><span style="color: #0000FF;">]</span> <span style="color: #008080;">then</span>
<span style="color: #000080;font-style:italic;">-- (Should multiple same-time routes exist, we could store
-- multiple back-links and whip up a simple [recursive]
-- routine to rebuild them all. Or just ignore them.)</span>
<span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #000080;font-style:italic;">-- (simplify debugging)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #000000;">timeframe</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{}</span> <span style="color: #000080;font-style:italic;">-- (simplify debugging)</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">finish</span>
<span style="color: #000000;">mintime</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">costs</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">minpath</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">finish</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">while</span> <span style="color: #004600;">true</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">object</span> <span style="color: #000000;">pyx</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">paths</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">]</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">pyx</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">minpath</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">prepend</span><span style="color: #0000FF;">(</span><span style="color: #000000;">minpath</span><span style="color: #0000FF;">,</span><span style="color: #000000;">pyx</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">paths</span><span style="color: #0000FF;">[</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span> <span style="color: #000080;font-style:italic;">-- (be safe, why not)</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">py</span><span style="color: #0000FF;">,</span><span style="color: #000000;">px</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">pyx</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">minpath</span><span style="color: #0000FF;">[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]!=</span><span style="color: #000000;">start</span> <span style="color: #008080;">then</span> <span style="color: #0000FF;">?</span><span style="color: #000000;">9</span><span style="color: #0000FF;">/</span><span style="color: #000000;">0</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">minpath</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #000000;">mintime</span><span style="color: #0000FF;">*</span><span style="color: #000000;">60</span><span style="color: #0000FF;">),</span><span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">)}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">surfacebylongitude</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">lon</span><span style="color: #0000FF;">)</span>
<span style="color: #000080;font-style:italic;">-- Create regional wind patterns on the map.</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">surfaceparameters</span> <span style="color: #0000FF;">=</span> <span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lon</span> <span style="color: #0000FF;"><</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">155.03</span> <span style="color: #0000FF;">?</span> <span style="color: #0000FF;">{</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">5.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">150</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.5</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:</span>
<span style="color: #008080;">iff</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lon</span> <span style="color: #0000FF;"><</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">155.99</span> <span style="color: #0000FF;">?</span> <span style="color: #0000FF;">{-</span><span style="color: #000000;">90.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">20</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">150</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.4</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">180.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">25</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">150</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">0.3</span><span style="color: #0000FF;">}))</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">surfaceparameters</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">slices</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">null</span><span style="color: #0000FF;">,</span><span style="color: #000000;">200</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">mutatetimeslices</span><span style="color: #0000FF;">()</span>
<span style="color: #000080;font-style:italic;">-- Vary wind speeds over time.</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">slices</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">deep_copy</span><span style="color: #0000FF;">(</span><span style="color: #000000;">slices</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">sj</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">sj</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">windvelocity</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sj</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span>
<span style="color: #000000;">windvelocity</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">1</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">0.002</span><span style="color: #0000FF;">*</span><span style="color: #000000;">i</span>
<span style="color: #000000;">sj</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">][</span><span style="color: #000000;">2</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">windvelocity</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">sj</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">slices</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">s</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">slices</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">gpoints</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">null</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">lat</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">19.78</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">/</span><span style="color: #000000;">60</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">/</span><span style="color: #000000;">60</span>
<span style="color: #000000;">gpoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #004600;">null</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">9</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">lon</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">155.0</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">/</span><span style="color: #000000;">60</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">/</span><span style="color: #000000;">60</span>
<span style="color: #000000;">gpoints</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">][</span><span style="color: #000000;">j</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{{</span><span style="color: #000000;">lat</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lon</span><span style="color: #0000FF;">},</span> <span style="color: #000000;">surfacebylongitude</span><span style="color: #0000FF;">(</span><span style="color: #000000;">lon</span><span style="color: #0000FF;">)}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">slices</span><span style="color: #0000FF;">[</span><span style="color: #000000;">s</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">gpoints</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">mutatetimeslices</span><span style="color: #0000FF;">()</span>
<span style="color: #008080;">constant</span> <span style="color: #000000;">fmt</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">"""
The route taking the least time found was:
%v
which has duration %s [route found in %s]
"""</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">fmt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">minimum_time_route</span><span style="color: #0000FF;">(</span><span style="color: #000000;">slices</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">},{</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">4</span><span style="color: #0000FF;">}))</span>
<span style="color: #0000FF;">{}</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">wait_key</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
The route taking the least time found was:
{{1,4},{1,5},{2,6},{3,7},{4,7},{5,7},{6,7},{7,7},{8,6},{8,5},{9,4}}
which has duration 4 hours, 43 minutes and 41s [route found in 0.0s]
</pre>
=={{header|Wren}}==
{{trans|Julia}}
A reasonably faithful translation though I haven't bothered to split the code up into modules (which would mean separate files in Wren) and have dispensed altogether with four functions which aren't actually called.
As Wren uses 0-based indexing the points in the minimum path have coordinates one less than those in the Julia results.
As you'd expect, this takes many times longer than Julia to run (about 24.5 minutes versus 3 minutes 20 seconds) but gets there in the end :)
<syntaxhighlight lang="wren">import "io" for File
/*
Class that represents a polar CSV file's data.
Contains a matrix, 'speeds', of sailing speeds indexed by wind velocity and angle of boat to wind.
'winds' is a list of wind speeds.
'degrees' is a list of angles in degrees of direction relative to the wind.
Note 0 degrees is directly into the wind, 180 degrees is directly downwind.
*/
class SailingPolar {
construct new(winds, degrees, speeds) {
_winds = winds
_degrees = degrees
_speeds = speeds // speeds[wind direction degrees, windspeed knots]
}
winds { _winds }
degrees { _degrees }
speeds {_speeds }
}
/*
Class that represents wind and surface current direction and velocity for a given position.
Angles in degrees, velocities in knots.
*/
class SurfaceParameters {
construct new(windDeg, windKts, currentDeg, currentKts) {
_windDeg = windDeg
_windKts = windKts
_currentDeg = currentDeg
_currentKts = currentKts
}
windDeg { _windDeg }
windKts { _windKts }
currentDeg { _currentDeg }
currentKts { _currentKts }
}
// Reads a sailing polar CSV file and returns a SailingPolar object containing the file data.
// A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
// The first line of file contains labels for the wind velocities that make up columns, and
// the first entry of each row makes up a column of angle of sailing direction from wind in degrees.
var getPolarData = Fn.new { |fileName|
var lines = File.read(fileName).split("\n")
var header = lines[0].trim().split(";")
var winds = header[1..-1].map { |x| Num.fromString(x) }.toList
var degrees = []
var speeds = []
for (line in lines[1..-1]) {
line = line.trim()
if (line == "") break // ignore final blank line if there is one
var cols = line.split(";")
degrees.add(Num.fromString(cols[0]))
speeds.add(cols[1..-1].map{ |x| Num.fromString(x) }.toList)
}
return SailingPolar.new(winds, degrees, speeds)
}
var R = 6372800 // Earth's approximate radius in meters
/* Class containing various helper methods which work with degrees rather than radians. */
class D {
// Converts degrees to radians.
static deg2Rad(deg) { (deg*Num.pi/180 + 2*Num.pi) % (2*Num.pi) }
// Converts radians to degrees.
static rad2Deg(rad) { (rad*180/Num.pi + 360) % 360 }
// Trig functions.
static sin(d) { deg2Rad(d).sin }
static cos(d) { deg2Rad(d).cos }
static asin(d) { rad2Deg(d.asin) }
static atan(x, y) { rad2Deg(x.atan(y)) }
}
// Calculates the haversine function for two points on the Earth's surface.
// Given two latitude, longitude pairs in degrees for a point on the Earth,
// get distance in meters and the initial direction of travel in degrees for
// movement from point 1 to point 2.
var haversine = Fn.new { |lat1, lon1, lat2, lon2|
var dlat = lat2 - lat1
var dlon = lon2 - lon1
var a = D.sin(dlat/2).pow(2) + D.cos(lat1) * D.cos(lat2) * (D.sin(dlon/2).pow(2))
var c = 2 * D.asin(a.sqrt)
var theta = D.atan(D.sin(dlon) * D.cos(lat2),
D.cos(lat1)*D.sin(lat2) - D.sin(lat1) * D.cos(lat2) * D.cos(dlon))
theta = (theta + 360) % 360
return [R * c * 0.5399565, theta]
}
// Returns the index of the first element of 'a' for which 'pred' returns true or -1 otherwise.
var findFirst = Fn.new { |a, pred|
for (i in 0...a.count) if (pred.call(a[i])) return i
return -1
}
// Returns the index of the last element of 'a' for which 'pred' returns true or -1 otherwise.
var findLast = Fn.new { |a, pred|
for (i in a.count-1..0) if (pred.call(a[i])) return i
return -1
}
// Calculate the expected sailing speed in a specified direction in knots,
// given sailing polar data, a desired point of sail in degrees, and wind speed in knots.
var boatSpeed = Fn.new { |sp, pointOfSail, windSpeed|
var winds = sp.winds
var degrees = sp.degrees
var speeds = sp.speeds
var udeg = findLast.call(degrees) { |t| t <= pointOfSail }
var odeg = findFirst.call(degrees) { |t| t >= pointOfSail }
var uvel = findLast.call(winds) { |t| t <= windSpeed }
var ovel = findFirst.call(winds) { |t| t >= windSpeed }
if ([udeg, odeg, uvel, ovel].any { |t| t == -1 }) return -1
var frac = (odeg == udeg && uvel == ovel) ? 1 :
(odeg == udeg) ? (windSpeed - winds[uvel])/(winds[ovel] - winds[uvel]) :
(uvel == ovel) ? (pointOfSail - degrees[udeg])/(degrees[odeg] - degrees[udeg]) :
((pointOfSail - degrees[udeg])/(degrees[odeg] - degrees[udeg]) +
(windSpeed - winds[uvel])/(winds[ovel] - winds[uvel]))/2
return speeds[udeg][uvel] + frac * (speeds[odeg][ovel] - speeds[udeg][uvel])
}
// Calculates the expected net boat speed in a desired direction versus the wind ('azimuth').
// This is generally different from the actual boat speed in its actual direction.
// Directions are in degrees ('pointos' is point of sail the ship direction from the wind),
// and velocity of wind ('ws') is in knots.
var sailingSpeed = Fn.new { |sp, azimuth, pointos, ws|
return boatSpeed.call(sp, pointos, ws) * D.cos((pointos - azimuth).abs)
}
// Calculates the net direction and velocity of a sailing ship.
// Arguments are sailing polar data, direction of travel in degrees from north, wind direction in
// degrees from north, wind velocity in knots, surface current direction in degrees, and
// current velocity in knots.
var bestVectorSpeed = Fn.new { |sp, dirTravel, dirWind, windSpeed, dirCur, velCur|
var azimuth = (dirTravel - dirWind) % 360
azimuth = (azimuth < 0) ? azimuth + 360 : azimuth
azimuth = (azimuth > 180) ? 360 - azimuth : azimuth
var VMG = boatSpeed.call(sp, azimuth, windSpeed)
var other = -1
var idx = -1
for (i in 0...sp.degrees.count) {
var ss = sailingSpeed.call(sp, azimuth, sp.degrees[i], windSpeed)
if (ss > other) {
other = ss
idx = i
}
}
if (other > VMG) {
azimuth = sp.degrees[idx]
VMG = other
}
var dirChosen = D.deg2Rad(dirWind + azimuth)
var wx = VMG * (dirChosen.sin)
var wy = VMG * (dirChosen.cos)
var curX = velCur * (D.deg2Rad(dirCur).sin)
var curY = velCur * (D.deg2Rad(dirCur).cos)
return [D.rad2Deg((wy + curY).atan(wx + curX)), ((wx + curX).pow(2) + (wy + curY).pow(2)).sqrt]
}
// Calculates the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2).
// Uses the data in SurfaceParameters for wind and current velocity and direction.
var sailSegmentTime = Fn.new { |sp, p, lat1, lon1, lat2, lon2|
var h = haversine.call(lat1, lon1, lat2, lon2)
var distance = h[0]
var dir = h[1]
var vel = bestVectorSpeed.call(sp, dir, p.windDeg, p.windKts, p.currentDeg, p.currentKts)[1]
// minutes/s * m / (knots * (m/s / knot)) = minutes
return (1 / 60) * distance / (vel * 1.94384)
}
/* Class that represents a point in 2-D space. Need value type semantics for comparisons etc. */
class Point2 {
construct new(x, y) {
_x = x
_y = y
}
x { _x }
y { _y }
+ (other) { Point2.new(x + other.x, y + other.y) }
== (other) { x == other.x && y == other.y }
!= (other) { !(this == other) }
toString { "[%(_x), %(_y)]" }
}
/*
Class that consists of a tuple of latitude and longitude in degrees.
NB: This uses latitude (often considered to be y) first then longitude (often considered to be x).
This latitude, then longitude ordering is as per ISO 6709 (en.wikipedia.org/wiki/ISO_6709).
*/
class Position {
construct new(lat, lon) {
_lat = lat
_lon = lon
}
lat { _lat }
lon { _lon }
}
/* Class that represents a Position with the SurfaceParameters of wind and current at the Position. */
class GridPoint {
construct new(pt, sp) {
_pt = pt
_sp = sp
}
pt { _pt }
pt=(value) { _pt = value }
sp { _sp }
sp=(value) { _sp = value }
}
/*
Class that consists of a matrix of GridPoints, each Position point with their SurfaceParameters.
A Vector of TimeSlice can give the surface characteristics for an ocean region over time.
*/
class TimeSlice {
construct new(gridPoints) {
_gridPoints = gridpoints
}
gridPoints { _gridPoints }
}
/*
Class that represents a routing problem and requiring the following parameters:
* timeinterval: the minutes duration for each TimeSlice
* timeframe: a vector of sequential timeslices for the ocean region
* obstacleindices: the Cartesian indices in each timeslice for
obstacles, such as land or shoals, where the ship may not go
* startindex: the timeslice position for time of starting
* start: starting location on grid of GridPoints
* finish: destination / finish location on grid of GridPoints
* allowrepeatvisits: whether the vessel may overlap its prior path, usually false.
*/
class RoutingProblem {
construct new(timeInterval, timeFrame, obstacleIndices, startIndex, start, finish, allowRepeatVisits) {
_timeInterval = timeInterval // minutes between timeFrame slices
_timeFrame = timeFrame
_obstacleIndices = obstacleIndices
_startIndex = startIndex
_start = start
_finish = finish
_allowRepeatVisits = allowRepeatVisits
}
timeInterval { _timeInterval }
timeFrame { _timeFrame }
obstacleIndices { _obstacleIndices }
startIndex { _startIndex }
start { _start }
finish { _finish }
allowRepeatVisits { _allowRepeatVisits }
}
/*
Class that represents a timed path and requires the following parameters:
* duration: minutes total to travel the path
* path: vector of Cartesian indices of points in grid for path to travel.
*/
class TimedPath {
construct new(duration, path) {
_duration = duration
_path = path
}
duration { _duration }
path { _path }
toString { "%(_duration) %(_path)" }
== (other) { this.toString == other.toString }
!= (other) { this.toString != other.toString }
}
var findMin = Fn.new { |a|
var min = a[0]
var idx = 0
for (i in 1...a.count) {
if (a[i] < min) {
min = a[i]
idx = i
}
}
return [min, idx]
}
var ntuples = [ [-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1] ]
var neighbors = List.filled(ntuples.count, null)
(0...ntuples.count).each { |i| neighbors[i] = Point2.new(ntuples[i][0], ntuples[i][1]) }
// Returns a list of points surrounding 'p' which are not otherwise excluded.
var surround = Fn.new { |p, mat, excluded|
var xmax = mat.count
var ymax = mat[0].count
return neighbors.map { |x| x + p }.where { |q|
return (0 <= q.x && q.x < xmax) && (0 <= q.y && q.y < ymax) && !excluded.contains(q)
}.toList
}
// Get the route (as a TimedPath) that minimizes time from start to finish for a given
// RoutingProblem (sea parameters) given sailing polar data (ship parameters).
var minimumTimeRoute = Fn.new { |rp, sp, verbose|
var timedPaths = [TimedPath.new(0, [rp.start])]
var completed = false
var minPath = TimedPath.new(1000, [])
for (i in 0...1000) {
var newPaths = []
verbose && System.print("Checking %(timedPaths.count) paths of length %(timedPaths[0].path.count)")
for (tpath in timedPaths) {
if (tpath.path[-1] == rp.finish) {
completed = true
newPaths.add(tpath)
} else {
var p1 = tpath.path[-1]
var num = tpath.duration.round
var den = rp.timeInterval.round
var slice = rp.timeFrame[(num/den).truncate % rp.timeFrame.count]
for (p2 in surround.call(p1, slice, rp.obstacleIndices)) {
if (rp.allowRepeatVisits || !tpath.path.contains(p2)) {
var gp1 = slice[p1.x][p1.y]
var gp2 = slice[p2.x][p2.y]
var lat1 = gp1.pt.lat
var lon1 = gp1.pt.lon
var lat2 = gp2.pt.lat
var lon2 = gp2.pt.lon
var t = sailSegmentTime.call(sp, gp1.sp, lat1, lon1, lat2, lon2)
var path = tpath.path.toList
path.add(p2)
newPaths.add(TimedPath.new(tpath.duration + t, path))
}
}
}
}
var set = {}
for (np in newPaths) set[np.toString] = np
timedPaths = set.values.toList
if (completed) {
var minDur = findMin.call(timedPaths.map { |x| x.duration }.toList)[0]
var finished = timedPaths.where { |x| x.path[-1] == rp.finish }.toList
var mi = findMin.call(finished.map { |x| x.duration }.toList)
var minFinDur = mi[0]
var idx = mi[1]
if (verbose) {
System.print("Current finished minimum: %(finished[idx]), others %(minDur)")
}
if (minDur == minFinDur) {
minPath = finished[idx]
break
}
}
}
return minPath
}
/*
The data is selected so the best time path is slightly longer than the
shortest length path. The forbidden regions are x, representing land or reef.
The allowed sailing points are . and start and finish are S and F.
x . . F . . x . x
. . . . . . . x x
x . . x x x . . .
. . x x x x . x x
x . . . x x . x .
x . . . x x . x .
. . . . x . . x .
x . . . . . . x .
. . . S . . . . .
*/
// These need to be changed to 0-based for Wren.
var ftuples = [
[1, 8], [2, 1], [2, 8], [3, 5], [3, 8], [4, 1], [4, 5], [4, 6], [4, 8], [5, 1],
[5, 5], [5, 6], [5, 8], [6, 3], [6, 4], [6, 5], [6, 6], [6, 8], [6, 9], [7, 1],
[7, 4], [7, 5], [7, 6], [8, 8], [8, 9], [9, 1], [9, 7], [9, 9]
]
var forbidden = List.filled(ftuples.count, null)
(0...ftuples.count).each { |i| forbidden[i] = Point2.new(ftuples[i][0]-1, ftuples[i][1]-1) }
// Create regional wind patterns on the map.
var surfaceByLongitude = Fn.new { |lon|
return (lon < -155.03) ? SurfaceParameters.new(-5, 8, 150, 0.5) :
(lon < -155.99) ? SurfaceParameters.new(-90, 20, 150, 0.4) :
SurfaceParameters.new(180, 25, 150, 0.3)
}
// Vary wind speeds over time.
var mutateTimeSlices = Fn.new { |slices|
var i = 1
for (slice in slices) {
for (j in 0...slice.count) {
for (k in 0...slice[j].count) {
var x = slice[j][k]
x.sp = SurfaceParameters.new(x.sp.windDeg, x.sp.windKts * (1 + 0.002 * i),
x.sp.currentDeg, x.sp.currentKts)
}
}
i = i + 1
}
}
var startPos = Point2.new(0, 3) // 0-based
var endPos = Point2.new(8, 3) // ditto
var slices = List.filled(200, null)
for (s in 0...200) {
var gpoints = List.filled(9, null)
for (i in 0..8) {
gpoints[i] = List.filled(9, null)
for (j in 0..8) {
var pt = Position.new(19.78 - 1/60 + i/60, -155.0 - 5/60 + j/60)
gpoints[i][j] = GridPoint.new(pt, surfaceByLongitude.call(pt.lon))
}
}
slices[s] = gpoints
}
mutateTimeSlices.call(slices)
var routeProb = RoutingProblem.new(10, slices, forbidden, 0, startPos, endPos, false)
var fileName = "polar.csv"
var sp = getPolarData.call(fileName)
var tp = minimumTimeRoute.call(routeProb, sp, false)
System.print("The route taking the least time found was:\n %(tp.path) \nwhich has duration " +
"%((tp.duration/60).truncate) hours, %((tp.duration%60).round) minutes.")</syntaxhighlight>
{{out}}
<pre>
The route taking the least time found was:
[[0, 3], [0, 4], [1, 5], [2, 6], [3, 6], [4, 6], [5, 6], [6, 6], [7, 5], [7, 4], [8, 3]]
which has duration 4 hours, 44 minutes.
</pre>
| |||