clc;
clear;

% a proposed model used for cs18.6:
% The number of simulations are defined based on the rupture magnitude
model_Mw = [5.5 6.0 8 8.5]; 
model_NumSim = [10 10 50 50];

% % -------------------------------------------------------------------------
% % read mag for all faults

% read all fault mag and length
fid=fopen('Fault_LeonardMag_Length_Width_Rake.txt','r'); % read the first line
read_format='%s %s %s %s %s';
C = textscan(fgetl(fid),read_format,1);
LastLine=false;
Num=0;

while LastLine==false
    
    Num=Num+1;

    fault.name(Num)=C{1};
    fault.mag(Num)=str2num(cell2mat(C{2}));
    fault.length(Num)=str2num(cell2mat(C{3}));
    fault.width(Num)=str2num(cell2mat(C{4}));
    fault.rake(Num)=str2num(cell2mat(C{5}));

    str=fgetl(fid);
    if length(str)>1
    C = textscan(str,read_format,1);
    LastLine=false;
    else
    LastLine=true;
    end
end
fclose(fid);

%% read source names to create srf

fid=fopen('CS_list_v18p6.txt','r'); % read the first line
format='%s %s %s';
C = textscan(fgetl(fid),format,1);
LastLine=false;
Num=0;

while LastLine==false
    
    Num=Num+1;

    srf.name(Num)=C{1};
    
    str=fgetl(fid);
    if length(str)>1
    C = textscan(str,format,1);
    LastLine=false;
    else
    LastLine=true;
    end
end
fclose(fid);

%% -------------------------------------------------------------------------
% % create comparison data
k=0;
for i=1:1:length(srf.name)
    
    for j=1:1:length(fault.name)
        
        if strcmp(srf.name(i),fault.name(j))==1
            
            k=k+1;
            srf.name(i);
            
            data.name(k) = fault.name(i);         
            data.mag(k) = fault.mag(j);  

%             if data.mag(k) > 7.8
%                 data.name(k)
%             end 
        end
    end
end


%% calc number of simulations for given fault name

for k=1:1:length(data.name)
    
    fault_Mw = data.mag(k);
    
    % using multi-linear model
    if (fault_Mw >= model_Mw(1)) && (fault_Mw < model_Mw(2))
        
        func = @(x) ((model_NumSim(2)-model_NumSim(1)) / (model_Mw(2)-model_Mw(1)))*(x - model_Mw(1)) + model_NumSim(1);
        num_sim = func(fault_Mw);
        
    elseif (fault_Mw >= model_Mw(2)) && (fault_Mw < model_Mw(3))
        
        func = @(x) ((model_NumSim(3)-model_NumSim(2)) / (model_Mw(3)-model_Mw(2)))*(x - model_Mw(2)) + model_NumSim(2);
        num_sim = func(fault_Mw);
        
    elseif (fault_Mw >= model_Mw(3)) && (fault_Mw < model_Mw(4))

        func = @(x) ((model_NumSim(4)-model_NumSim(3)) / (model_Mw(4)-model_Mw(3)))*(x - model_Mw(3)) + model_NumSim(3);
        num_sim = func(fault_Mw);
    end
        
    data.num_sim(k) = ceil(num_sim);
    
end

disp(['total number of realizations = ',num2str(sum(data.num_sim))])

% % write file used for SRF generation

outFile = fopen(['fault_selection_c18p6.txt'],'w');

for k=1:1:length(data.name)
    
    fprintf(outFile,'%s %d%s \n',cell2mat(data.name(k)),data.num_sim(k),'r');
end
fclose(outFile);

%% Mw - NumSim
A(:,1) = data.mag;
A(:,2) = data.num_sim;
B=sortrows(A,1);

figure('Position',[600 80 700 600]);
h1=plot(B(:,1),B(:,2),'Color','b','LineWidth',1.5,'MarkerSize',8); hold on;
set(gca,'FontSize',5+18,'FontUnits','pixels','XMinorTick','on','YMinorTick','on','XScale','linear','YScale','linear'); grid on;  hold on;
set(gca,'GridLineStyle','-', 'MinorGridLineStyle', '-'); hold on;
xlabel('Rupture magnitude, M_w','FontSize',5+ 20);
ylabel('Number of realizations','FontSize',5+ 20);
legend ([h1],{'v18.6'},'Location','NorthWest');
set(legend,'FontSize',5+ 14);
% xlim([0 ceil(max(T))]);
ylim([1 100]);
% annotation('textbox',[0.68 0.62 0.30 0.04],'String',{'Wairau'},'HorizontalAlignment','right','Margin',0,'FontSize',5+14,'LineStyle','none','BackgroundColor',[1 1 1]);
saveas(gcf,['fig_NumSim_vs_Mw'],'fig'); %close;


            