Ship

function [rT, cT] = ship(chart, a, b, mT)
%

% results: 170577
% time: 0.89
% complexity: 9
% node count: 519?

global verbose
verbose = false;
if verbose
    fprintf('\n')
end

% winds
rW = chart(:,:,1);
cW = chart(:,:,2);

% dimensions
[m, n] = size(rW);

% start
[ai, aj] = ind2sub([m n], a);
% island
[bi, bj] = ind2sub([m n], b);

% current pos
pi = ai;
pj = aj;

% current vel
vi = 0;
vj = 0;

% 1st target: try to reach b
[k, pi, pj, vi, vj, rT1, cT1] = reachTarget(0, pi, pj, vi, vj, ...
    bi, bj, rW, cW, mT);
if k < 1000
    % 2st target: try to return to a
    [k, pi, pj, vi, vj, rT2, cT2] = reachTarget(k, pi, pj, vi, vj, ...
        ai, aj, rW, cW, mT);
else
    rT2 = [];
    cT2 = [];
end

% thrust
rT = [rT1; rT2];
cT = [cT1; cT2];

end

function [k, pi, pj, vi, vj, rT, cT] = reachTarget(k0, pi, pj, vi, vj, ...
    ti, tj, rW, cW, mT)

global verbose

N = 40;
rT = zeros(N,1);
cT = zeros(N,1);

mT = fix(mT/2);

% dimensions
[m, n] = size(rW);

for k = 1:N % N turns one away, N the other
    % wind
    wi = rW(pi, pj);
    wj = cW(pi, pj);

    % desired direction
    di = sat(ti - pi, 2);
    dj = sat(tj - pj, 2);

    % thrust
    r = sat(di - (vi + wi), mT);
    c = sat(dj - (vj + wj), mT);

    % update position
    [pi, vi, r, qi] = keepInMap(pi, vi, r, wi, m, mT);
    [pj, vj, c, qj] = keepInMap(pj, vj, c, wj, n, mT);

    % check bounds
    if qi ~= 0 || qj ~= 0
        if verbose
            fprintf(['next move would place ship out of bounds: ' ...
                'turn %d, stopped at (%d,%d)\n'], ...
                k, pi, pj)
        end
        rT = rT(1:k-1);
        cT = cT(1:k-1);
        k = 1000;
        return;
    end

    % save thrust
    rT(k) = r;
    cT(k) = c;
    
    % reached target?
    if pi == ti && pj == tj
        rT = rT(1:k);
        cT = cT(1:k);
        if verbose
            fprintf('reached target (%d,%d)\n', ti, tj)
        end
        return;
    end
    
end

if verbose
    fprintf('out of turns (%d): stopped at (%d,%d)\n', ...
        k, pi, pj)
end

end

function u = sat(u, mT)

u = sign(u)*min(mT, abs(u));

end

function [p, v, t, q] = keepInMap(p0, v0, t, w, d, m)

v = v0 + t + w;
p = p0 + v;
q = 1*(p < 1) - 1*(p > d); % = 1 if p < 1, -1 if p > n
while (t < m && q == 1) || (t > 0 && q == -1)
    t = t + q;
    v = v0 + t + w;
    p = p0 + v;
    q = 1*(p < 1) - 1*(p > d);
end

end