Weather routing: Difference between revisions
Content added Content deleted
(→{{header|Wren}}: Now deals with Windows line separators when reading the file.) |
|||
Line 915: | Line 915: | ||
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]] |
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. |
which has duration 4.0 hours, 43.697879668707344 minutes. |
||
</pre> |
|||
=={{header|Phix}}== |
|||
{{trans|Go}} |
|||
<lang Phix>-- demo/rosetta/Weather_Routing.exw |
|||
function to_numbers(sequence s) |
|||
for i=1 to length(s) do |
|||
s[i] = to_number(s[i]) |
|||
end for |
|||
return s |
|||
end function |
|||
function getpolardata(string s) |
|||
-- |
|||
-- 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 |
|||
-- |
|||
sequence lines = split_any(s,"\r\n",no_empty:=true), |
|||
winds = to_numbers(split(lines[1],";")[2..$]), |
|||
degrees = {}, speeds = {} |
|||
for i=2 to length(lines) do |
|||
sequence l = to_numbers(split(lines[i],";")) |
|||
if length(l)!=length(winds)+1 then ?9/0 end if |
|||
degrees = append(degrees,l[1]) |
|||
speeds = append(speeds,l[2..$]) |
|||
end for |
|||
return {winds, degrees, speeds} |
|||
end function |
|||
-- |
|||
-- 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] |
|||
-- |
|||
constant {winds, degrees, speeds} = getpolardata(get_text("polar.csv")) |
|||
-- (note the distributed version uses a literal string constant instead) |
|||
constant R = 6372800 -- Earth's approximate radius in meters |
|||
constant timeinterval = 10 -- the minutes duration for each TimeSlice |
|||
function deg2rad(atom deg) return remainder(deg*PI/180+2*PI,2*PI) end function |
|||
function rad2deg(atom rad) return remainder (rad*(180/PI)+360,360) end function |
|||
function sind(atom deg) return sin(deg2rad(deg)) end function |
|||
function cosd(atom deg) return cos(deg2rad(deg)) end function |
|||
function asind(atom deg) return rad2deg(arcsin(deg)) end function |
|||
function atand(atom x,y) return rad2deg(atan2(x,y)) end function |
|||
function haversine(atom 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. |
|||
-- |
|||
atom dlat = lat2 - lat1, |
|||
dlon = lon2 - lon1, |
|||
a = power(sind(dlat/2),2) + cosd(lat1)*cosd(lat2)*power(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 = remainder(theta+360, 360) |
|||
return {R*c*0.5399565, theta} |
|||
end function |
|||
function findfirst(atom v, sequence s) |
|||
-- Returns the index of the first element of s >= v |
|||
for i=1 to length(s) do |
|||
if s[i]>=v then return i end if |
|||
end for |
|||
return -1 |
|||
end function |
|||
function findlast(atom v, sequence s) |
|||
-- Returns the index of the last element of s <= v |
|||
for i=length(s) to 1 by -1 do |
|||
if s[i]<=v then return i end if |
|||
end for |
|||
return -1 |
|||
end function |
|||
function boatspeed(atom pointofsail, windspeed) |
|||
-- |
|||
-- 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) |
|||
-- |
|||
integer udeg = findlast(pointofsail, degrees), |
|||
odeg = findfirst(pointofsail, degrees), |
|||
uvel = findlast(windspeed, winds), |
|||
ovel = findfirst(windspeed, winds) |
|||
if find(-1,{udeg, odeg, uvel, ovel}) then return -1 end if |
|||
atom wu = winds[uvel], |
|||
wh = winds[ovel], |
|||
du = degrees[udeg], |
|||
dh = degrees[odeg], |
|||
f |
|||
if odeg==udeg then |
|||
f = iff(uvel==ovel?1:(windspeed-wu)/(wh-wu)) |
|||
elsif uvel==ovel then |
|||
f = (pointofsail-du)/(dh-du) |
|||
else |
|||
f = ((pointofsail-du)/(dh-du)+(windspeed-wu)/(wh-wu))/2 |
|||
end if |
|||
atom su = speeds[udeg,uvel], |
|||
sh = speeds[odeg,ovel], |
|||
res = su + f*(sh-su) |
|||
return res |
|||
end function |
|||
function sailingspeed(atom azimuth, pointofsail, 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 (pointofsail is the ship direction from wind), |
|||
-- and velocity of wind (ws) is in knots. |
|||
-- |
|||
return boatspeed(pointofsail, ws) * cosd(abs(pointofsail-azimuth)) |
|||
end function |
|||
struct SurfaceParameters |
|||
-- |
|||
-- wind and surface current, direction and velocity, for a given position |
|||
-- directions are in degrees from north, and velocities are in knots. |
|||
-- |
|||
public atom winddirection, |
|||
windvelocity, |
|||
currentdirection, |
|||
currentvelocity |
|||
end struct |
|||
function bestvectorspeed(atom dirtravel, SurfaceParameters p) |
|||
-- |
|||
-- Calculate the net direction and velocity of a sailing ship. |
|||
-- |
|||
atom wd = p.winddirection, |
|||
wv = p.windvelocity, |
|||
cd = p.currentdirection, |
|||
cv = p.currentvelocity, |
|||
azimuth = remainder(dirtravel-wd,360) |
|||
if azimuth<0 then azimuth += 360 end if |
|||
if azimuth>180 then azimuth = 360-azimuth end if |
|||
atom vmg = boatspeed(azimuth, wv), |
|||
other = -1 |
|||
integer idx = -1 |
|||
for i=1 to length(degrees) do |
|||
atom ss = sailingspeed(azimuth, degrees[i], wv) |
|||
if ss>other then {other,idx} = {ss,i} end if |
|||
end for |
|||
if other>vmg then |
|||
azimuth = degrees[idx] |
|||
vmg = other |
|||
end if |
|||
atom dirchosen = deg2rad(wd + azimuth), |
|||
dircurrent = deg2rad(cd), |
|||
wx = vmg * sin(dirchosen), |
|||
wy = vmg * cos(dirchosen), |
|||
curx = cv * sin(dircurrent), |
|||
cury = cv * cos(dircurrent) |
|||
return sqrt(power(wx+curx,2) + power(wy+cury,2)) |
|||
end function |
|||
function sailsegmenttime(SurfaceParameters p, atom lat1, lon1, lat2, lon2) |
|||
-- |
|||
-- Calculate the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2). |
|||
-- Uses the data in SurfaceParameters p for wind and current velocity and direction. |
|||
-- |
|||
atom {distance, direction} = haversine(lat1, lon1, lat2, lon2), |
|||
velocity = bestvectorspeed(direction, p) |
|||
-- minutes/s * m / (knots * (m/s / knot)) = minutes |
|||
atom res = (1/60) * distance / (velocity * 1.94384) |
|||
return res |
|||
end function |
|||
-- |
|||
-- 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. |
|||
-- |
|||
constant chart = split(""" |
|||
...S..... |
|||
x......x. |
|||
....x..x. |
|||
x...xx.x. |
|||
x...xx.x. |
|||
..xxxx.xx |
|||
x..xxx... |
|||
.......xx |
|||
x..F..x.x""",'\n') |
|||
function minimum_time_route(sequence timeframe, start, finish) |
|||
-- |
|||
-- 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. |
|||
-- |
|||
atom t0 = time(), |
|||
mintime = 1000.0 |
|||
integer xmax = length(chart[1]), |
|||
ymax = length(chart), |
|||
{py,px} = start |
|||
sequence todo = {start}, |
|||
costs = repeat(repeat(-1,xmax),ymax), -- (lowest durations) |
|||
paths = repeat(repeat(0,xmax),ymax), -- (single backlinks) |
|||
minpath = {} |
|||
costs[py,px] = 0 |
|||
while length(todo) do |
|||
{py,px} = todo[1] |
|||
todo = todo[2..$] |
|||
atom duration = costs[py,px] |
|||
integer sdx = remainder(floor(round(duration)/timeinterval),length(timeframe))+1 |
|||
sequence s = timeframe[sdx] |
|||
for nx=px-1 to px+1 do |
|||
for ny=py-1 to py+1 do |
|||
if (nx!=px or ny!=py) |
|||
and nx>=1 and nx<=xmax |
|||
and ny>=1 and ny<=ymax |
|||
and chart[ny,nx]!='x' then |
|||
sequence gp1 = s[py,px], -- {pt,SurfaceParameters} |
|||
gp2 = s[ny,nx] -- "" |
|||
atom {lat1, lon1} = gp1[1], |
|||
{lat2, lon2} = gp2[1] |
|||
SurfaceParameters sp = gp1[2] |
|||
atom nt = duration + sailsegmenttime(sp, lat1, lon1, lat2, lon2) |
|||
if costs[ny,nx]=-1 or nt<costs[ny,nx] then |
|||
costs[ny,nx] = nt |
|||
paths[ny,nx] = {py,px} |
|||
if not find({ny,nx},todo) then |
|||
todo = append(todo,{ny,nx}) |
|||
end if |
|||
elsif nt==costs[ny,nx] then |
|||
-- (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.) |
|||
?9/0 |
|||
end if |
|||
end if |
|||
end for |
|||
end for |
|||
s = {} -- (simplify debugging) |
|||
end while |
|||
timeframe = {} -- (simplify debugging) |
|||
{py,px} = finish |
|||
mintime = costs[py,px] |
|||
minpath = {finish} |
|||
while true do |
|||
object pyx = paths[py,px] |
|||
if pyx=0 then exit end if |
|||
minpath = prepend(minpath,pyx) |
|||
paths[py,px] = 0 |
|||
{py,px} = pyx |
|||
end while |
|||
if minpath[1]!=start then ?9/0 end if |
|||
return {minpath,elapsed(mintime*60),elapsed(time()-t0)} |
|||
end function |
|||
function surfacebylongitude(atom lon) |
|||
-- Create regional wind patterns on the map. |
|||
sequence p = iff(lon < -155.03 ? { -5.0, 8, 150, 0.5} : |
|||
iff(lon < -155.99 ? {-90.0, 20, 150, 0.4} : |
|||
{180.0, 25, 150, 0.3})) |
|||
SurfaceParameters res = new(p) |
|||
return res |
|||
end function |
|||
procedure mutatetimeslices(sequence slices) |
|||
-- Vary wind speeds over time. |
|||
for i=1 to length(slices) do |
|||
sequence s = slices[i] |
|||
for j=1 to length(s) do |
|||
sequence sj = s[j] |
|||
for k=1 to length(sj) do |
|||
object sjk = sj[k] |
|||
SurfaceParameters p = sjk[2] |
|||
p.windvelocity = p.windvelocity * (1+0.002*i) |
|||
end for |
|||
end for |
|||
end for |
|||
end procedure |
|||
sequence slices = repeat(null,200) |
|||
for s=1 to length(slices) do |
|||
sequence gpoints = repeat(null,9) |
|||
for i=1 to 9 do |
|||
atom lat = 19.78 - 2/60 + i/60 |
|||
gpoints[i] = repeat(null,9) |
|||
for j=1 to 9 do |
|||
atom lon = -155.0 - 6/60 + j/60 |
|||
gpoints[i][j] = {{lat,lon}, surfacebylongitude(lon)} |
|||
end for |
|||
end for |
|||
slices[s] = gpoints |
|||
end for |
|||
mutatetimeslices(slices) |
|||
constant fmt = """ |
|||
The route taking the least time found was: |
|||
%v |
|||
which has duration %s [route found in %s] |
|||
""" |
|||
printf(1,fmt,minimum_time_route(slices,{1,4},{9,4}))</lang> |
|||
{{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> |
</pre> |
||