% Script to perform particle tracking using deep-sea ocean current data
% generated by "generateVectorFields.m"

% Each day the particles in the simulation experience a new vector field
% drawn from a pool of size nVectorFields

% Jon Yearsley & Julia Sigwart (Aug 2010)
% jon.yearsley@ucd.ie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

nReps = 500; % Number of replicate simulations to run
nDays = 500; % Number of days of simulation
nVectorFields = 200; % Number of different vector fields to use in each simulation

inputPrefix = 'Shallow_800m_1400m_Data'; % Prefix for inout vector field data
vectorFieldDir = '.'; % Directory where vector field files are kept
outputPrefix = 'shallowSims'; % Prefix for results file
outputDir = 'particleTracks'; % Directory where simulation results will be saved

% Define a box of starting trajectories
chitonFile = 'deepseaChitons.xls'; % Excel containgin the observed distribution of chitons
longlatSD = 0.1; % Std. Dev. of initial particle positions about observed chiton locations
island = 1; % Index defining island where simulate particles will start
% island = 1  start at Vanuatu
% island = 2  start at New Caledonia
% island = 3  start at Solomon Islands
% island = 4  start at a random location anywhere in the simulation arena

numParticles = 20; % Number of particles to throw into the box
% (all these particles will experience the same set of vector fields, whilst
% particles in different replicate simulations will experience different
% sets of vector fields)

displayOn = false;  % If true plot figures
options = odeset('AbsTol', 1e-2,'RelTol',1e-6,'NormControl','on','MaxOrder',5); % ODE solver soptions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% No more parameters to change below this line

%% Read chiton data and coastline data 

load('coastline.mat') % Load coastline data

% Read in observed distribution of deep-sea chitons 
[obsNUM,obsTXT] = xlsread(chitonFile);
popIslandID = obsNUM(:,1);
popLongLat = obsNUM(:,2:3);
popAbundance = obsNUM(:,4);
popAbundance(popIslandID~=island) = 0; % Set islands not of interest to zero abundance
popIslandName = obsTXT(2:end,1);
islandNames = {'Vanuatu', 'NewCaledonia', 'Solomon','RandomStart'};

for rep=1:nReps
    disp(['Simulation Rep number ' num2str(rep)])
    
    % Load randomly chosen vector fields
    inputFile = dir([vectorFieldDir '\' inputPrefix '*.mat']);
    nFiles = length(inputFile);
    if nFiles==0
        error('No vector field files found')
    elseif nFiles<nVectorFields
        warning('Fewer vector field files than required for simulation');
        useFiles = randi(nFiles,nVectorFields,1);        
    else
        tmp = randperm(nFiles);
        useFiles = tmp(1:nVectorFields);
    end
    odefun = cell(nVectorFields,1);
    for f=1:nVectorFields
        inFile = [vectorFieldDir '\' inputFile(useFiles(f)).name];
        load(inFile);
                
        % Define an ODE function for each vector field
        odefun{f} =  @(t,y) dlonglatdt(y, longGrid, latGrid, v1Grid, v2Grid);
    end
    
    longVec = longGrid(1,:);
    latVec = latGrid(:,1)';
    gridStep = longVec(2)-longVec(1);
    
    % Clear some space in memory
    clear 'timestep' 'v1Data' 'v2Data' 'longlatEnd' 'longlatStart' 'xEnd' 'xStart' 
    clear 'distStep' 'depth' 'year' 'month' 'velocity' 'timeStep' 'xData' 'dlatdt' 'dlongdt' 'dt'
    
    %% Simulate the tracks of some particles
    tic
    
    latTrack = nan(nDays+1,numParticles);
    longTrack = nan(nDays+1,numParticles);
    startEpicentre = zeros(numParticles,2);
    
    % Find limits of each island
    ind = find(isnan(longlatCoast(:,1)));
    ind = ind(ind(1:(length(ind)-1))+1 ~= ind(2:end)); 
       
    for p = 1:numParticles % Simulating particles 1 by 1 appears as quick as simulating them all together
        % Set the initial positions
        overLand = true;
        y0 = [nan, nan];       
        while overLand
            if island==4 % If a random starting location for a particle
                minLong = min(longlatCoast(:,1))+2;
                maxLong = max(longlatCoast(:,1))-2;
                minLat = min(longlatCoast(:,2))+2;
                maxLat = max(longlatCoast(:,2))-2;
                y0(1) = minLong + (maxLong-minLong)*rand;
                y0(2) = minLat + (maxLat-minLat)*rand;
            else % Else a partcile starts off near the locality of an observation
                % Pick locality proportional to relative abundance
                tmp = cumsum(popAbundance/sum(popAbundance))-rand;
                tmpInd = find(abs(tmp-min(tmp(tmp>0)))<1e-6,1);
                
                startEpicentre(p,:) = popLongLat(tmpInd,:);

                % Define starting position of a particle
                y0(1) = popLongLat(tmpInd,1) + longlatSD*randn;
                y0(2) = popLongLat(tmpInd,2) + longlatSD*randn;
            end
            
            % If particle is over land, try again
            overLand = false;
            for i=1:(length(ind)-1),
                polyx = longlatCoast((ind(i)+1):(ind(i+1)-1),1);
                polyy = longlatCoast((ind(i)+1):(ind(i+1)-1),2);
                
                if inpolygon(y0(1),y0(2),polyx,polyy);
                    overLand = true;
                end
            end
        end
        
        track = nan(nDays+1,2);
        track(1,:) = y0;
        tspan = [0,1];
        for t=1:nDays
            % Simulate particles for 1 day using a randomly selected vector field
            [T,Y] = ode15s(odefun{randi(nVectorFields)},tspan,y0,options); % ode15s was the quickest routine for the accuracy required, but other ode solvers can be used here
            track(t+1,:) = Y(end,:); % Save particle position at the end of one day
            y0 = Y(end,:);           % Reset starting point and pick a new vector field
        end
        % Save the particle's track  before simulating next particle
        longTrack(:,p) = track(:,1);
        latTrack(:,p) = track(:,2);
    end
    tspan = (0:nDays);
    toc
    %% Plot the track
    if displayOn
        figure(2)
        map = colormap(jet);
        col = map(round(linspace(1,64,numParticles)),:);
        v1Grid = squeeze(v1GridMulti(randi(nVectorFields),:,:));
        v2Grid = squeeze(v2GridMulti(randi(nVectorFields),:,:));
        
        for i=1:numParticles
            plot(longTrack(:,i),latTrack(:,i),'-','Color',col(i,:))
            hold on
            plot(longTrack(1:20:end,i),latTrack(1:20:end,i),'x','Color',col(i,:),'MarkerSize',10)
            plot(longTrack(1,i), latTrack(1,i),'s','MarkerSize',10,'MarkerFaceColor',col(i,:),'MarkerEdgeColor','k')
            plot(longTrack(end,i), latTrack(end,i),'ok','MarkerSize',10,'MarkerFaceColor',col(i,:))
        end
        ind = find(isnan(longlatCoast(:,1)));
        for i=3:length(ind)-1,
            fill(longlatCoast((ind(i)+1):(ind(i+1)-1),1),longlatCoast((ind(i)+1):(ind(i+1)-1),2),[0,0.4,0])
            hold on
        end
        
        dlon = 1;
        dlat = 1;
        axis([min(longTrack(:))-dlon, max(longTrack(:))+dlon, min(latTrack(:))-dlat, max(latTrack(:))+dlat])
        hold off
        xlabel('Longitude','FontSize',25)
        ylabel('Latitude','FontSize',25)
        title(['Simulated tracks of ' num2str(numParticles) ' particles for ' num2str(nDays) ' days (x''s mark 20 day intervals)'])
    end
    %% Save the particle Tracks in a unique filename
    
    if ~exist(outputDir,'dir') % Check output folder exists, and create it if it doesn't
        mkdir(outputDir)
    end
    x = 1;
    filename = [outputDir '\' outputPrefix, '_' islandNames{island} '_' num2str(x) '.mat'];
    while exist(filename,'file')==2
        x = x+1;
    filename = [outputDir '\' outputPrefix, '_' islandNames{island} '_' num2str(x) '.mat'];
    end
    
    save(filename,'longTrack','latTrack','tspan','options','inputFile','useFiles','nVectorFields','nDays','popLongLat','popAbundance','longlatSD','startEpicentre')
    
end

