function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_ver1(mCatalog,nMethod)
% function [mCatDecluster, mCatAfter, vCluster, vCl, vMainCluster] = calc_decluster_ver1(mCatalog,nMethod)
% ----------------------------------------------------------------------------------------------------------
%
% Function to decluster earthquake catalog using the windowing technique in space and time by
% Knopoff & Gardner, GJR astr. Soc, 28, 311-313, 1972
% Gardner & Knopoff, BSSA, 64,5, 1363-1367, 1974
% using different windows
%
% Incoming variables
% mCatalog : Incoming earthquake catalog (ZMAP format)
% nMethod  : Window length for declustering (see calc_windows.m)
%            1: Gardener & Knopoff, 1974
%            2: Gruenthal pers. communication
%            3: Urhammer, 1986
%
% Outgoing variables:
% mCatDecluster : Declustered earthquake catalog
% mCatAfter     : Catalog of aftershocks (foreshocks)
% vCluster      : Vector indicating only aftershocks/foreshocls in cluster using a cluster number
% vCl           : Vector indicating all events in clusters using a cluster number
% vMainCluster  : Vector indicating mainshocks of clusters using a cluster number
%
% J. Woessner, woessner@seismo.ifg.ethz.ch
% last update: 31.07.02

%% Added:
% 31.07.02 Correction for problem of mainshocks with a cluster number as aftershocks belong to two sequences
% 31.07.02 Corrected fMaxClusterMag(nMagCount) to fMaxClusterMag, since counting variable not needed

%%% Remember: Improve zero length cluster problem, no waitbars yet, time difference mainshock-original event not implemented
%% This is the version before calc_decluster.m

%% Initialize Vectors and Matrices
mCatDecluster = [];
mCatAfter = [];
vCluster = zeros(length(mCatalog),1); % Initialize all events as mainshock
vCl = zeros(length(mCatalog),1); % Initialize all events as mainshock
vSel = zeros(length(mCatalog),1); % Initialize all events as mainshock
vMainCluster = zeros(length(mCatalog),1); % Initialize

[nXSize, nYsize] = size(mCatalog);
if nXSize == 0
    disp('Load new catalog');
    return
end

vDecDate = mCatalog(:,3);
nCount = 0;    % Variable of cluster number

fMagThreshold = min(mCatalog(:,6)); % Set Threshold to minimum magnitude of catalog
for nEvent=1:length(mCatalog(:,6))
    %nEvent
    if vCluster(nEvent) == 0
        fMagnitude(nEvent) = mCatalog(nEvent, 6);
        if fMagnitude(nEvent) >= fMagThreshold
            %% Define first aftershock zone and determine magnitude of strongest aftershock
            fMag = fMagnitude(nEvent);
            [fSpace, fTime] = calc_windows(fMagnitude(nEvent), nMethod);
            fSpaceDeg = km2deg(fSpace);
            %% This first if is for events with no location given
            if isnan(mCatalog(nEvent,1))
                %vSel = (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime  & vCluster(nEvent) == 0);
                vSel = (mCatalog(:,3) == mCatalog(nEvent,3));
            else
                vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) & (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                    & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vCluster(nEvent) == 0);
            end
            mTmp = mCatalog(vSel,:);
            %nMagCount = 1; % Variable to determine max. magnitude in a cluster
            if length(mTmp(:,1)) == 1  % Only one event thus no cluster
                fMaxClusterMag = fMag;
            else
                fMaxClusterMag = max(mTmp(:,6));
                % Search for event with bigger magnitude in cluster and resize windows if necessary
                while fMaxClusterMag-fMag > 0
                    [fSpace, fTime] = calc_windows(fMaxClusterMag, nMethod);
                    fSpaceDeg = km2deg(fSpace);
                    vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) &...
                        (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                        & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vSel == 1);
                    mTmp = mCatalog(vSel,:);
                    if isempty(mTmp) % no events in aftershock zone
                        break;
                    end
                    fMag = fMaxClusterMag;
                    %nMagCount = nMagCount+1;
                    fMaxClusterMag = max(mTmp(:,6));
                    if fMaxClusterMag - fMag == 0 % no bigger event in aftershock zone
                        break;
                    end
                end  % End of while
                nCount = nCount + 1; % Set cluster number
            end % End of if
            % Define final aftershock zone if different from original one
            if fMaxClusterMag > fMagnitude(nEvent)
                [fSpace, fTime] = calc_windows(fMaxClusterMag, nMethod);
                fSpaceDeg = km2deg(fSpace);
                % Select events with the Event having the number nEvent, since it is a foreshock (time condition)
                vSel = (abs(mCatalog(:,1) - mCatalog(nEvent,1)) <= fSpaceDeg) & (abs(mCatalog(:,2) - mCatalog(nEvent,2)) <= fSpaceDeg)...
                    & (vDecDate(:,1)-vDecDate(nEvent,1) >= 0) & (vDecDate(:,1)-vDecDate(nEvent,1) <= fTime & vCluster(nEvent) == 0);
            end
            [vIndice]=find(vSel); % Vector of indices with Clusters
            vTmpCluster(vIndice,:) = nCount;
            length(vTmpCluster(vIndice,:));
            nI=1; % Variable counting the length of the cluster
            % Select the right numbers for the cluster using the indice vector vIndice
            % First: Insert cluster number after check for length
            % Second: Check if it's a mainshock
            % Third: Keep the former cluster indice;
            while nI <= length(vIndice)
                if (~isempty(vTmpCluster(vIndice(nI))) & length(vTmpCluster(vIndice,:)) > 1 & vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = vTmpCluster(vIndice(nI));
                    %vEventnr(vIndice,:) = nEvent;
                elseif  (~isempty(vTmpCluster(vIndice(nI)))  &&  length(vTmpCluster(vIndice,:)) == 1  &&  vCluster(vIndice(nI)) == 0)
                    vCluster(vIndice(nI)) = 0;
                else
                    vCluster(vIndice(nI)) = vCluster(vIndice(nI));
                end
                nI=nI+1;
            end %End of while nI
            %%% Check if the Cluster is not just one event which can happen in case of keeping the former
            %%% cluster number in preceeding while-Loop
            vSelSingle = (vCluster == nCount);
            [vIndiceSingle] = find(vSelSingle);
            %vTmpSingle(vIndiceSingle,:);
            if length(vIndiceSingle) == 1
                nCount
                vIndiceSingle
                vCluster(vIndiceSingle)=0; % Set the event as mainsock
                nCount = nCount-1; % Correct the cluster number down by one
            end
        end % End of if checking magnitude threshold fMagThreshold
    end % End of if chekcing if vCluster == 0
end % End of FOR over mCatalog
%nCount
%% vCL Cluster vector with mainshocks in it; vCluster is now modified to get rid of mainshocks
vCl = vCluster;

%% Matrix with cluster indice, magnitude and time
mTmpCat = [vCluster mCatalog(:,6) mCatalog(:,3)];
%% Delete largest event from cluster series and add to mainshock catalog
for nCevent = 1:nCount
    %nCevent
    vSel4 = (mTmpCat(:,1) == nCevent); % Select cluster
    mTmpCat2 = mCatalog(vSel4,:);
    fTmpMaxMag = max(mTmpCat2(:,6)); % Select max magnitude of cluster
    vSelMag = (mTmpCat2(:,6) == fTmpMaxMag);
    [nMag] = find(vSelMag);
    if length(nMag) == 1
        vSel5 = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag); % Select the event
        [vIndiceMag] = find(vSel5); % Find indice
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent; % Set mainshock vector to cluster number
    elseif isempty(nMag)
        disp('Nothing in ')
        nCevent
    else
        vSel = (mTmpCat(:,1) == nCevent & mTmpCat(:,2) == fTmpMaxMag);
        mTmpCat3 = mCatalog(vSel,:);
        [vIndiceMag] = min(find(vSel)); % Find minimum indice of event with max magnitude in cluster
        vCluster(vIndiceMag) = 0;  % Set cluster value to zero, so it is a mainshock
        vMainCluster(vIndiceMag) = nCevent;  % Set mainshock vector to cluster number
    end
end % End of For nCevent

%% Create a catalog of aftershocks (mCatAfter) and of declustered catalog (mCatDecluster)
vSel = (vCluster(:,1) > 0);
mCatDecluster=mCatalog(~vSel,:);
mCatAfter = mCatalog(vSel,:);

