163 lines
4.9 KiB
Matlab
163 lines
4.9 KiB
Matlab
% [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
|