Weather routing: Difference between revisions
→{{header|Phix}}
(→{{header|Wren}}: Now deals with Windows line separators when reading the file.) |
|||
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]]
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>
| |||