%% PART 1 LOAD THE DATA in ZMAP FORMAT
a=load('tutorial_100day_tohokuzmap.mat');

mtx_data=a.mtx_tok100dayzmap;
mtx_times= [mtx_data(:,3:5) mtx_data(:,8:10)]; 

%% PART 2 ADD A COLUMN GIVING THE ELAPSED TIME FROM THE MAINSHOCK
% datenum of the mainshock
var_otime = datenum(2011, 3, 11, 5, 46, 18.12);
% datenum of the other events
vec_evtimes=datenum(mtx_times);
% elapsed time in seconds between the mainshock and aftershocks
vec_elaptime=etime( datevec(vec_evtimes) , datevec(var_otime*ones(size(mtx_data,1),1)) );
% add the column back onto the matrix
mtx_data100=[mtx_data vec_elaptime];

% year column dates need to be in dec.year format  year.DayMonthHourMinSec
vec_decdates=decyear(mtx_times);
% set year column to decimal year format
mtx_data100(:,3)=vec_decdates;

% find the first day events
%%% create a logical array of the 1st day of events 86400 seconds = 1 day 
var_tmp=mtx_data100(:,11) <= 86400;
%%% filter the 1st day using the logical array
mtx_data1=mtx_data100(var_tmp,:);

%% PART 2 MAGNITUDE HISTOGRAM
% grab the unrounded magnitude values and store
a=mtx_data100(:,6);
c=mtx_data1(:,6);
unrnd100=a;
unrnd1=c;

% round the values
anew=a*10;
anew=ceil(anew)/10;
mtx_tok100dayzmap(:,6)=anew;

cnew=c*10;
cnew=ceil(cnew)/10;
mtx_tok1dayzmap(:,6)=cnew;

% create histogram vectors
x=0:.1:9;
[n1day,bin1day]= hist(mtx_tok1dayzmap(:,6),x');
[n100day,bin100day]= hist(mtx_tok100dayzmap(:,6),x');
xp=x';

[n1dayNORND,bin1dayNORND]= hist(c,x');
[n100dayNORND,bin100dayNORND]= hist(a,x');

% plot the histogram
%%% WARNING USES PANEL function,  either use your own SUBPLOT code or
%%% panel from http://www.mathworks.com/matlabcentral/fileexchange/20003
%%%%%%%%%%%%%%%%%%%%%%%%%%%% BEGIN PANEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%
phist = panel();
phist.pack(2, 2);
phist.marginleft=25;
phist.de.margin = [5 10 5 10];
phist.fontsize=14;

phist(1, 1).select();
h=bar(xp,n100day,'histc');
phist(1,1).hold('on');
h1=bar(xp,n1day,'histc');
set(h1,'FaceColor','r','EdgeColor','k');
xlim([0 9]);
set(gca,'XMinorTick','on');
set(gca,'TickDir','out');
ylabel('Number (100 day)');
axis('square');

phist(1, 2).select();
h=bar(xp,n100dayNORND,'histc');
phist(1,2).hold('on');
h1=bar(xp,n1dayNORND,'histc');
set(h1,'FaceColor','r','EdgeColor','k');
xlim([0 9]);
set(gca,'XMinorTick','on');
set(gca,'TickDir','out');
ylabel('Number (100 day)');
axis('square');


phist(2, 1).select();
h=bar(xp,n1day,'histc');
set(h,'FaceColor','r','EdgeColor','k');
xlim([0 9]);
set(gca,'XMinorTick','on');
set(gca,'TickDir','out');
ylabel('Number (1 day)');
axis('square');
xlabel('Ceil Applied');

phist(2, 2).select();
h=bar(xp,n1dayNORND,'histc');
set(h,'FaceColor','r','EdgeColor','k');
xlim([0 9]);
set(gca,'XMinorTick','on');
set(gca,'TickDir','out');
ylabel('Number (1 day)');
xlabel('No Ceil Applied');
axis('square');
%%%%%%%%%%%%%%%%%%%%%%%%%%%% END PANEL %%%%%%%%%%%%%%%%%%%%%%%%%%%%
% close the plot
close all;

%% PART 3 ESTIMATE THE Mc PART I  (BOOTSTRAPPING)
% clear stuff we do not need
a100new=anew;
c1new=cnew;
clearvars -except mtx_data1 mtx_data100 c1new a100new;
% setup Mc calculation parameters which will not change
%%% fbinning is the binning used for the magnitudes, typically 0.1
%%% nSample is the number of bootstraps to use
%%% mCatalog is the catalog in ZMAP format
fBinning=0.1;
nSample=1;
%%% change this
mCatalog=mtx_data1(:,1:10);
mCatalog(:,6)=c1new;
mCatalog=mCatalog(2:end,:);

% MAXIMUM CURVATURE AVERAGE
    fMcCorr = .2; 
    nMinNum = 20;
    nMethod =1;
    [fMcMAX, fStd_McMAX, fBvalueMAX, fStd_BMAX, fAvalueMAX, fStd_AMAX, vMcMAX, mBvalueMAX] = calc_McBboot(mCatalog, fBinning, nSample, nMethod, nMinNum, fMcCorr);
    fMcMAX= ceil(fMcMAX*10)/10;
    t1=[fMcMAX fStd_McMAX fBvalueMAX fStd_BMAX fAvalueMAX fStd_AMAX];
    display(['MAXC AVERAGE VALUE IS ' num2str(fMcMAX)])

% BEST COMBINATION AVERAGE
    fMcCorr = 0; 
    nMinNum = 20;
    nMethod =5;
    [fMcMAX, fStd_McMAX, fBvalueMAX, fStd_BMAX, fAvalueMAX, fStd_AMAX, vMcMAX, mBvalueMAX] = calc_McBboot(mCatalog, fBinning, nSample, nMethod, nMinNum, fMcCorr);
    fMcMAX= ceil(fMcMAX*10)/10;
    t2=[fMcMAX fStd_McMAX fBvalueMAX fStd_BMAX fAvalueMAX fStd_AMAX];
     display(['BEST COMBINATION AVERAGE VALUE IS ' num2str(fMcMAX)]) 
     
 % EMR COMBINATION AVERAGE
    fMcCorr = 0; 
    nMinNum = 20;
    nMethod =6;
    [fMcMAX, fStd_McMAX, fBvalueMAX, fStd_BMAX, fAvalueMAX, fStd_AMAX, vMcMAX, mBvalueMAX] = calc_McBboot(mCatalog, fBinning, nSample, nMethod, nMinNum, fMcCorr);
    fMcMAX= ceil(fMcMAX*10)/10;
    t3=[fMcMAX fStd_McMAX fBvalueMAX fStd_BMAX fAvalueMAX fStd_AMAX];
%     [hPlot] = plot_FMDboth(mCatalog, fMcMAX, fAvalueMAX, fBvalueMAX, fBinning, nMinNum);
%     xlabel('magnitude','FontSize',14,'FontWeight','bold');
%     ylabel('frequency','FontSize',14,'FontWeight','bold');
    display(['EMR AVERAGE VALUE IS ' num2str(fMcMAX)]) 
   
   var_mcaverages=[t1;t2;t3];
  
%% PART 3 ESTIMATE THE Mc PART II  (NO BOOTSTRAPPING)
% Calculate Mc using Maximum Curvature
nMethod =1; % Maximum curvature
fBinning = 0.1; 
fMcCorrection = .2; % Used to add a correction because maximum curvature often underestimates Mc
[fMc] = calc_Mc(mCatalog, nMethod, fBinning, fMcCorrection);
% Calculate a and b-value according to Shi & Bolt; a-value representative
% for the Catalog length
% Select catalog above Mc
vSel = mCatalog(:,6) >= fMc;
mCat = mCatalog(vSel,:);
[fMeanMag, fBvalue, fStdDev, fAvalue] =  calc_bmemag(mCat, fBinning);
t4=[fMc,fMeanMag, fBvalue, fStdDev, fAvalue];
display(['MAXC VALUE IS ' num2str(fMc)]) 

% Calculate Mc using Best Combination
nMethod =5; % BEST COMBINATION
fBinning = 0.1; 
fMcCorrection = .2; % Used to add a correction because maximum curvature often underestimates Mc
[fMc] = calc_Mc(mCatalog, nMethod, fBinning, fMcCorrection);
% Calculate a and b-value according to Shi & Bolt; a-value representative
% for the Catalog length
% Select catalog above Mc
vSel = mCatalog(:,6) >= fMc;
mCat = mCatalog(vSel,:);
[fMeanMag, fBvalue, fStdDev, fAvalue] =  calc_bmemag(mCat, fBinning);
t5=[fMc,fMeanMag, fBvalue, fStdDev, fAvalue];
display(['BEST COMBINATION VALUE IS ' num2str(fMc)]) 

% Calculate Mc using EMR
nMethod =6; % EMR
fBinning = 0.1; 
fMcCorrection = 0; % Used to add a correction because maximum curvature often underestimates Mc
[fMc] = calc_Mc(mCatalog, nMethod, fBinning, fMcCorrection);
% Calculate a and b-value according to Shi & Bolt; a-value representative
% for the Catalog length
% Select catalog above Mc
vSel = mCatalog(:,6) >= fMc;
mCat = mCatalog(vSel,:);
[fMeanMag, fBvalue, fStdDev, fAvalue] =  calc_bmemag(mCat, fBinning);
t6=[fMc,fMeanMag, fBvalue, fStdDev, fAvalue];   
display(['EMR VALUE IS ' num2str(fMc)]) 

var_mcs=[t4;t5;t6;];
   
%% PART 4 ESTIMATE THE b-value
clearvars -except mtx_data1 mtx_data100 c1new a100new mCatalog var_mcs var_mcaverages;
% Estimate the b-value
var_magnitudes=mCatalog(:,6);
bin=.1;
var_tmp=6;
mtx_lsq=[];
mtx_mlm=[];
mtx_zm=[];

for i=1:.1:var_tmp
var_tmplogi= var_magnitudes >= i;
var_tempmag=var_magnitudes(var_tmplogi);
   
[blsq, alsq, std_errlsq] = bval_lsqreg(var_tempmag, bin);
[bmlm, amlm, std_errmlm] = bval_maxlkh2(var_tempmag, bin);

[fMeanMag, fBValue, fStdDev, fAValue] =  calc_bmemag(var_tempmag, bin);

var_tmprow=[blsq alsq std_errlsq i size(var_tempmag,1) ];
mtx_lsq=[mtx_lsq;var_tmprow];

var_tmprow=[bmlm amlm std_errmlm i size(var_tempmag,1) ];
mtx_mlm=[mtx_mlm;var_tmprow];  

var_tmprow=[fBValue fAValue fStdDev i size(var_tempmag,1) ];
mtx_zm=[mtx_mlm;var_tmprow];  

end

figureBvalues=figure;
h1=plot([1:.1:var_tmp],mtx_lsq(:,1),'k');
hold on
h2=plot([1:.1:var_tmp],mtx_lsq(:,1),'.k');
h3=plot([1:.1:var_tmp],mtx_mlm(:,1),'b');
h4=plot([1:.1:var_tmp],mtx_mlm(:,1),'.b');

h5=plot([1:.1:var_tmp],mtx_mlm(:,1),'g');
h6=plot([1:.1:var_tmp],mtx_mlm(:,1),'.g');
hold off
hline(1,'r');
legend([h1,h3,h5],'Least Squares','Maximum Like','ZMAP')
ylabel('b-value');


figureNumvalues=figure;
[ax1,h1,i1]=plotyy(1:.1:var_tmp,mtx_lsq(:,1),1:.1:var_tmp,mtx_lsq(:,5));
hold on
[ax2,h2,i2]=plotyy(1:.1:var_tmp,mtx_mlm(:,1),1:.1:var_tmp,mtx_mlm(:,5));
hold off

set(h1,'Color','k')
set(h1,'Marker','.')

set(h2,'Color','b')
set(h2,'Marker','.')

set(i2,'Marker','.')
h8=hline(1,'r');
ylabel('b-value');
xlabel('Mc value');
hh=legend([h1,h2,i2,h8],'Least Squares','Maximum Like','# Events','b=1');

%% PART 5 ESTIMATE THE MOL parameters
%%% The modified Omori law   n(t) = K/(t+c)^p
%%% all parameters are real and positive.  
%%% It is usual to normalize the units in days
%%% K is the productivity constant
%%% c is the time adjustment of catalog incompleteness
%%% t is the time in days
%%% p is the decay constant usually between 0.7 and 1.5

%%% Script parameters
%%% windowInDays - the number of days of data to use. script will start at
%%% first data and use data up to this value
%%% bootStraps - number of bootstraps to use
%%% model - number of model to use. 1 unless you are familiar
%%% evt_time -  vector of event times cut at Mc normalized to days.  hence .9999 is
%%% just under 1 day
%%% evt_mag - vector of event magnitudes cut at Mc 
%%% blankfigurehandle - pass in a blank figure handle
close all;

Mc=4;

catalogtouse=mtx_data100;

% filter by Mc
var_ind=catalogtouse(:,6) >= Mc;
mtx_tempcatalog=catalogtouse(var_ind,:);
% convert times to days since ms
mtx_tempcatalog(:,11)=mtx_tempcatalog(:,11)/86400;

windowInDays=100;
bootstraps=300;
model=1;
evt_time=mtx_tempcatalog(1:end,11);
evt_mag=mtx_tempcatalog(1:end,6);
blankfigurehandle=figure;

[p,k,c,timevec,winsize,time_as,cummodel,cumobs,ksstuff,fAIC,fRMS,nMod,mRes]= f_calc_omori_auto(windowInDays,bootstraps,model,evt_time,evt_mag,blankfigurehandle);



