%   [lamb_all,lamb,lmab_err,unit,eps,b]=UnlimitGR(t,M,iop,Mmin)
%
% ESTIMATES THE MEAN ACTIVITY RATE WITHIN THE WHOLE SAMPLE AND WITHIN THE 
% COMPLETE PART OF THE SAMPLE, THE ROUND-OFF ERROR OF MAGNITUDE AND THE 
% GUTENBERG-RICHTER B-VALUE USING THE UNLIMITED G-R LED MAGNITUDE 
% DISTRIBUTION MODEL 
%
% !! THIS FUNCTION MUST BE EXECUTED AT START-UP OF THE UNBOUNDED 
%       GUTENBERG-RICHETR HAZARD ESTIMATION MODE !!
%
% AUTHOR: S. Lasocki ver 2 01/2015 within IS-EPOS project.
%
% DESCRIPTION: The assumption on the unlimited Gutenberg-Richter relation 
% leads to the exponential distribution model of magnitude distribution 
% from and above the catalog completness level Mmin. The shape parameter of 
% this distribution and consequently the G-R b-value is estimated by 
% maximum likelihood method (Aki-Utsu procedure).  
% The mean activity rate, lamb, is the number of events >=Mmin into the
% length of the period in which they occurred. Upon the value of the input 
% parameter, iop, the used unit of time can be either day ot month or year.    
% The round-off interval length - eps if the least non-zero difference 
% between sample data or 0.1 is the least difference is greater than 0.1. 
% 
% INPUT:
%   t - vector of earthquake occurrence times
%   M - vector of magnitudes from a user selected catalog
%   iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month', 
%       iop=2 - 'year'
%   Mmin - catalog completeness level. Must be determined externally.
%   can take any value from [min(M), max(M)].
%
% OUTPUT:
%   lamb_all - mean activity rate for all events
%   lamb - mean activity rate for events >= Mmin
%   lamb_err - error paramter on the number of events >=Mmin. lamb_err=0
%       for 7 or more events >=Mmin and the parameter estimation is
%       continued, lamb_err=1 otherwise, all output paramters except 
%       lamb_all and lamb are set to zero and the function execution is 
%       terminated.  
%   unit - string with name of time unit used ('year' or 'month' or 'day').
%   eps - length of the round-off interval of magnitudes.
%   b - Gutenberg-Richter b-value
%
% LICENSE
%     This file is a part of the IS-EPOS e-PLATFORM.
%
%     This 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.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%


function [lamb_all,lamb,lamb_err,unit,eps,b]=UnlimitGR(t,M,iop,Mmin)
if isempty(t) || numel(t)<3 || isempty(M(M>=Mmin)) %K03OCT
    t=[1 2];M=[1 2];   end   %K30SEP


lamb_err=0;
n=length(M);
t1=t(1);
for i=1:n
    if M(i)>=Mmin; break; end
    t1=t(i+1);
end
    t2=t(n);
for i=n:1
    if M(i)>=Mmin; break; end
    t2=t(i-1);
end
nn=0;
for i=1:n
    if M(i)>=Mmin
        nn=nn+1;
    end
end

% SL 03MAR2015 ----------------------------------
[NM,unit]=time_diff(t(1),t(n),iop);    
lamb_all=n/NM;
[NM,unit]=time_diff(t1,t2,iop);        
lamb=nn/NM;
% SL 03MAR2015 ----------------------------------

if nn<7
    eps=0;b=0;
    lamb_err=1;
    return;
end

eps=magn_accur(M);
xx=M(M>=Mmin); %K21OCT2014
% x=sort(M,'descend');
% for i=1:n
%      if x(i)<Mmin; break; end
%      xx(i)=x(i);                                                      %
% end
clear x;
beta=1/(mean(xx)-Mmin+eps/2);
b=beta/log(10);
clear xx
end

function [NM,unit]=time_diff(t1,t2,iop)              % SL 03MAR2015

% TIME DIFFERENCE BETWEEEN t1,t2 EXPRESSED IN DAY, MONTH OR YEAR UNIT
%
% t1 - start time (in MATLAB numerical format)
% t2 - end time (in MATLAB numerical format)  t2>=t1
%   iop - determines the used unit of time. iop=0 - 'day', iop=1 - 'month', 
%       iop=2 - 'year'
%
%   NM - number of time units from t1 to t2
%   unit - string with name of time unit used ('year' or 'month' or 'day').

if iop==0
    NM=(t2-t1);
    unit='day';
elseif iop==1
    V1=datevec(t1);
    V2=datevec(t2);
    NM=V2(3)/eomday(V2(1),V2(2))+V2(2)+12-V1(2)-V1(3)/eomday(V1(1),V1(2))...
        +(V2(1)-V1(1)-1)*12;
    unit='month';
else
    V1=datevec(t1);
    V2=datevec(t2);
    NM2=V2(3);
    if V2(2)>1
        for k=1:V2(2)-1
            NM2=NM2+eomday(V2(1),k);
        end
    end
    day2=365; if eomday(V2(1),2)==29; day2=366; end;
    NM2=NM2/day2;
    NM1=V1(3);
    if V1(2)>1
        for k=1:V1(2)-1
            NM1=NM1+eomday(V1(1),k);
        end
    end
    day1=365; if eomday(V1(1),2)==29; day1=366; end;
    NM1=(day1-NM1)/day1;
    NM=NM2+NM1+V2(1)-V1(1)-1;
    unit='year';
end
end

function [eps]=magn_accur(M)
x=sort(M);
d=x(2:length(x))-x(1:length(x)-1);
eps=min(d(d>0));
if eps>0.1; eps=0.1;end
end
