%% Matlab pseudocode for Retention Clock Matrix (RCM) Analysis
%   
% For the accompanying manuscript, please see:
% Defne, Z, NK Ganju, A Aretxabaleta (2016), Estimating time-dependent connectivity in marine systems.
% Geophysical Research Letters
%
% October 2015
% zdefne@usgs.gov
% Zafer Defne
%
--------------------------------------------------------------------------
%   START of copying permission statement
%--------------------------------------------------------------------------
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   at your option) any later version.
% 
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%   <http://www.gnu.org/licenses/>
%--------------------------------------------------------------------------
%   END of copying permission statement
%--------------------------------------------------------------------------

clear
%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Release times shifted to t=0
% read LTRANS output from 24 individual runs
nrelease=24;
% for each release (release label icd)
 
for icd =1:nrelease
    % calculate shift. 2 represents 1 hr shift with 0.5 hr output
    shft=2*icd;
    shft_last=2*nrelease;
    % load lon lat values in LTRANS output 
    LON{icd+1}=lon( 1+shft : end+1+shft - shft_last , : );
    LAT{icd+1}=lat( 1+shft : end+1+shft - shft_last , : );
end
%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create a stack of connectivity matrices
 
% read the initial release coordinates
x0=LON(1,:)
y0=LAT(1,:)
 
% Set number of days to use from LTRANS output
dayn=28;
% number of LTRANS output time steps per day
outPerDay=48;
% total number of time steps
ti_tot=dayn*outPerDay
% LTRANS output time step (1 for 30 mins)
ti_step=1; 
 
% load ESA boundaries
% polyListESA is a list of ESA polygons with x, y vertices for each polygon 
% numPolys = total number of Polys 
polyListESA{numESAS}(x, y)
 
% mark the particles in the source polygons initially
for i=1:numPolys
    m{1,i}= inpolygon(x0, y0, polyListESA{i}(x), polyListESA{i}(y)); 
    % number of particles in source
    S(i)=sum(m{1,i});
end
 
% mark the particles in the destination polygons at each time step  
k=0;
% for each time step
for ti=1 : ti_step : ti_tot
    x2=LON(ti,:);
    y2=LAT(ti,:);
    % mark the particles in the destination polygons 
    for i=1:numPolys
        m{2,i}= inpolygon(x2, y2, polyListESA{i}(x), polyListESA{i}(y)); 
    end
    for i=1:numPolys %Source on y axis
        for j=1:numPolys %Destination on x axis
            % number of particles in destiantion
            N(i,j)=sum(m{1,i}.*m{2,j}); 
            N(N==0)=nan;
            % connectivity at each time step calculated with eqn: (N_di / N_so)
            NbyS(i,j)=N(i,j)/S(i); % 
        end
    end
    
    k=k+1;
    % create a stack of connectivity matrices. Pn(time, source, destination)  
    Pn(k,:,:)=NbyS(:,:);
end
 
%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot Retention Clock Matrix (RCM)
% this part of the code is presented mostly as is for ease of reproducibility
 
% maximum time scale in days
T=28
% time granularity (size of slices in each RCM)
dt=48
 
%calculate number of timesteps
hrsPerDay=24;
numtstep= T*hrsPerDay / dt;
 
% calculate averaging window length
wlen=dt*(outPerDay/24)
 
% calculate average P in each slice
Pavg=[];
for k=1:numtstep
    ind=(k-1)*wlen +1 : k*wlen-1;
    Pavg(k, :, :)=nanmean(Pn(ind,:,:),1);
end
 
% set minimum treshold for P (optional)
Pmin=0.001;
Pavg(Pavg<Pmin)=nan;
 
t=[0:numtstep]/numtstep*2*pi;
figure
hp=gcf
set(gcf, 'position', [500 200 800 800])
for di=1:numPolys
    % RCM orientation:
    % Source Polys on y-axis, destination Polys on x-axis
    % top left corner:(ESA1,ESA1), bottom right corner:(numPolys,numPolys)
    x=(di-1)/numPolys;
    for si=1:numPolys
        y=(numPolys-si)/numPolys; 
        subplot('position', [x y 1/numPolys 1/numPolys])
        hold on
        P=squeeze(Pavg(:, si, di));
        P(isnan(P))=0;
        if max(P)>0 % print Retention Clock if P is larger than a treshold value 
            for i=1:numtstep
                if(P(i)~=0)
                    patch([sin(t(i:i+1)) 0], [cos(t(i:i+1)) 0],  P(i));
                else
                    patch([sin(t(i:i+1)) 0], [cos(t(i:i+1)) 0], [1 1 1]);
                end
            end
            if di==si
                text(0,0, sprintf('%d', si), 'background', 'w', ...
                    'HorizontalAlignment', 'center', 'fontName', 'Arial')
            end
        end 
        set(gca, 'xColor', [1 1 1 ])
        set(gca, 'yColor', [1 1 1 ])
        set(gca, 'xtick', [])
        set(gca, 'ytick', [])
        set(gca, 'xtickLabel', [])
        set(gca, 'ytickLabel', [])
        caxis([0 cmax])
        set(gca, 'box', 'on')
        if si==di
            set(gca,'color', [.5 .5 .5])
        else
            set(gca,'color', [.8 .8 .8])
        end
    end
end
% colormap(flipud(copper))
set(gcf, 'paperPositionMode','auto')
set(gcf, 'renderer', 'painters')
 
% set colormap
colormap(w2c2k(32))
 
hc = colorbar
set(hc,'position', [ 0.94         0.33        0.015         0.42])
 
figname=fullfile('connectivity',dirname,sprintf('RCM_%dd_%dhr', T, dt));
% save as Matlab figure
saveas(gcf, figname)
% print eps
print(gcf, '-r600', '-depsc', figname)