lazySailor8

function [thrustR, thrustC] = solver8(chart, aIndex, bIndex, maxThrottle)

  windR = chart(:,:,1);
  windC = chart(:,:,2);
  
  [nR,nC] = size(windR);
  [AR,AC] = ind2sub([nR,nC],aIndex);
  [BR,BC] = ind2sub([nR,nC],bIndex);

  [X,Y] = meshgrid(-maxThrottle:maxThrottle);
  pick = abs(X(:))+abs(Y(:)) < maxThrottle;
  tryThrustR = X(pick);
  tryThrustC = Y(pick);
  
  [thrustR, thrustC, score, turnAround] = quickSolver([], [], false);
  
  for j=1:120

      if isempty(thrustR)
        break
      end
      
      index = randi(length(thrustR));
      
      thrustRT = thrustR(1:index);
      thrustCT = thrustC(1:index);
      
      thrustRT(end) = thrustRT(end) + randi([-1 1]);
      thrustCT(end) = thrustCT(end) + randi([-1 1]);
      
      [thrustRT, thrustCT, scoreT, turnAround] ...
        = quickSolver(thrustRT, thrustCT, turnAround*(turnAround<index));
         
      if scoreT<score
        score = scoreT;
        thrustR = thrustRT;
        thrustC = thrustCT;
      end
    
  end
      
  function [thrR, thrC, score, turnAround] = quickSolver(thrR, thrC, turnAround)

    posR = AR;
    posC = AC;
    velR = 0;
    velC = 0;
    dB   = inf;
    thrS = 0;
    
    for ii=1:length(thrR)
      velRi = velR + thrR(ii) + windR(posR,posC);
      velCi = velC + thrC(ii) + windC(posR,posC);
      posRi = posR + velRi;
      posCi = posC + velCi;
      thrS  = thrS + abs(thrR(ii)) + abs(thrC(ii));
      if posRi>nR || posRi<1 || posCi>nC || posCi<1
        break
      end
      posR = posRi;
      posC = posCi;
      velR = velRi;
      velC = velCi;
      dBTemp = (posR-BR)^2 + (posC-BC)^2;
      if dBTemp < dB
        dB = dBTemp;
      end
    end
    
    for ii=1:1000

      tryVelR = velR + tryThrustR + windR(posR,posC);
      tryVelC = velC + tryThrustC + windC(posR,posC);
      tryPosR = posR + tryVelR;
      tryPosC = posC + tryVelC;

      inside = find(tryPosR<size(windR,1) ...
                  & tryPosR>1 ...
                  & tryPosC<size(windR,2) ...
                  & tryPosC>1);
      
      if turnAround
        distance    = (posR-AR).^2 + (posC-AC).^2;
        tryDistance = (tryPosR(inside)-AR).^2 + (tryPosC(inside)-AC).^2;
      else
        distance    = (posR-BR).^2 + (posC-BC).^2;
        tryDistance = (tryPosR(inside)-BR).^2 + (tryPosC(inside)-BC).^2;
      end

      improveDistance = inside(tryDistance < distance);
      if isempty(improveDistance)
        if ~turnAround
          turnAround = length(thrR);
          continue
        else
          break
        end
      end
      [dummy,I] = min(abs(tryThrustR(improveDistance))+abs(tryThrustC(improveDistance)));
      I = improveDistance(I);

      posR = tryPosR(I);
      posC = tryPosC(I);

      velR = tryVelR(I);
      velC = tryVelC(I);

      dBTemp = (posR-BR)^2 + (posC-BC)^2;
      if dBTemp < dB
        dB = dBTemp;
      end
      
      thrR(end+1) = tryThrustR(I);
      thrC(end+1) = tryThrustC(I);
      
      thrS  = thrS + abs(thrR(end)) + abs(thrC(end));

    end
    
    dA = (posR-AR)^2 + (posC-AC)^2;
    
    score = dB + dA + thrS;

  end
  
  
end