%% Introduction

% This code gets the 1st cycle backbone curve of the digitized PZ hysteretic loops
% in terms of Vpz-gamma and calculates the elastic stiffness (Ke) and the 1st  
% yield according to 3 methods. (a) at gamma_y, (b) at 20% drop in the
% stiffness (for this to be accurately calculated, there is need for
% interpolation and (c) according to the Japanese provisions, that define yield 
% by 0.35% plastic deformations. That way, a parallel line to the initial stiffness,
% is brought and the intersection with the envelope defined that yield.
clear; clc; close all;

%% Input Parameters

stiffness_drop = 0.7; % i.e. 70% drop in the elastic stiffness means that we search for 0.3Ke
points = 10000; % points for interpolation curve
gamma_y = 0.0021; % According to Krawinkler (1978)
multiplier = 100; % Present gamma values in different scale (either 100, or 1)
error = 0.001; % For the Vy_0_35 calculation.

% First row: panel zone shear distortion, gamma
% Second row: panel zone shear force, V_pz
% This example uses the 1st positive cycle envelope of A2 specimen, Krawinkler (1971)
first_cycle_envelope = [
0	0
0.001856964	1643.415689
0.002559598	1863.400467
0.00336261	2044.564401
0.005570891	2277.489459
0.009184442	2523.354799
0.013500627	2730.399295
]; 

%% Data processing parameters

shear_yield_offset = 0.0035; % According to 0.2% yield strain

data_length = length(first_cycle_envelope);
data(:,1) = first_cycle_envelope(1:data_length,1);
data(:,2) = first_cycle_envelope(1:data_length,2);

%% Vy,drop calculation

figure('position',[550 100 520 300],'color','white');
hold on; grid on; box on;

plot(data(:,1).*multiplier,data(:,2),'LineStyle','-','color','black','LineWidth', 1.2, 'HandleVisibility', 'On');
Ke = (data(2,2)-data(1,2))/(data(2,1)-data(1,1));
data_length = length(data(:,1));
xplastic = zeros(data_length - 1,1);
yplastic = zeros(data_length - 1,1);

for i = 1:data_length-1
    xplastic(i) = data(i+1,1);
    yplastic(i) = data(i+1,2);
end

step = (xplastic(end)-xplastic(1))/points;
xxplastic = xplastic(1):step:xplastic(end);
yyplastic = pchip(xplastic,yplastic,xxplastic); % 3rd order polynomial fit
xxelastic = data(1,1):step:data(2,1);
yyelastic = interp1q([data(1,1) data(2,1)]', [data(1,2) data(2,2)]', xxelastic');

plot(xxelastic.*multiplier,yyelastic,'LineStyle','--','color','red','LineWidth', 1.5, 'HandleVisibility', 'On');
plot(xxplastic.*multiplier,yyplastic,'LineStyle','--','color','red','LineWidth', 1.5, 'HandleVisibility', 'Off');
for i = 1:length(xxplastic)-1  
    Klocal(i) = (yyplastic(i+1)-yyplastic(i))/(xxplastic(i+1)-xxplastic(i));
    if Klocal(i)<(1-stiffness_drop)*Ke
        V_drop = yyplastic(i);
        x_drop = xxplastic(i);
        scatter(x_drop.*multiplier,V_drop,30,'o', 'MarkerEdgeColor','black',...
                          'MarkerFaceColor','white',...
                          'LineWidth',0.1, 'HandleVisibility', 'Off');
        break
    end
end

%% Vy (@ gamma_y) & Vp (@ 4gamma_y) calculation

max_y = 10000; % a relatively big value to account for the error
max_4y = 10000;
gamma_4y = 4*gamma_y;
if xxplastic(1) < gamma_y
    for i = 1:length(xxplastic)
        if abs(gamma_y-xxplastic(i)) < max_y
            max_y = abs(gamma_y-xxplastic(i));
            V_gy = yyplastic(i);
            x_gy = xxplastic(i);
        end
    end
else
    for i = 1:length(xxelastic)
        if abs(gamma_y-xxelastic(i)) < max_y
            max_y = abs(gamma_y-xxelastic(i));
            V_gy = yyelastic(i);
            x_gy = xxelastic(i);
        end
    end
end
for i = 1:length(xxplastic)
    if abs(gamma_4y-xxplastic(i)) < max_4y
        max_4y = abs(gamma_4y-xxplastic(i));
        V_4gy = yyplastic(i);
        x_4gy = xxplastic(i);
    end
end
plot([x_gy x_gy].*multiplier,[0 V_gy],'LineStyle','--','color',[0.6 0.6 0.6],'LineWidth', 1.0, 'HandleVisibility', 'Off'); hold on;
plot([x_4gy x_4gy].*multiplier,[0 V_4gy],'LineStyle','--','color',[0.6 0.6 0.6],'LineWidth', 1.0, 'HandleVisibility', 'Off'); hold on;
scatter(x_gy.*multiplier,V_gy,30,'o', 'MarkerEdgeColor','black',...
                      'MarkerFaceColor','white',...
                      'LineWidth',0.1, 'HandleVisibility', 'Off');
scatter(x_4gy.*multiplier,V_4gy,30,'o', 'MarkerEdgeColor','black',...
                      'MarkerFaceColor','white',...
                      'LineWidth',0.1, 'HandleVisibility', 'Off');
                  
%% Vy_0_35 according to 0.2% yield strain (0.35% shear yield strain)

flag = 0;
step2 = 1;
japanx(1) = shear_yield_offset;
japany(1) = 0;
while flag == 0
    step2 = step2+1;
    japanx(step2) = japanx(step2-1)+0.000001;
    japany(step2) = japany(step2-1)+Ke*(japanx(step2)-japanx(step2-1));
    japanvec = sqrt((japany(step2))^2+(japanx(step2))^2);
    for i = 1:length(xxplastic)
        pointvec = sqrt((yyplastic(i))^2+(xxplastic(i))^2);
        if abs((japanvec-pointvec)/pointvec) < error && abs((japanx(step2)-xxplastic(i))/xxplastic(i)) < error
            V_0_35 = yyplastic(i);
            x_0_35 = xxplastic(i);
            plot(japanx.*multiplier,japany,'LineStyle','--','color',[0.6 0.6 0.6],'LineWidth', 1.0, 'HandleVisibility', 'Off'); hold on;
            scatter(x_0_35.*multiplier,V_0_35,30,'o', 'MarkerEdgeColor','black',...
                      'MarkerFaceColor','white',...
                      'LineWidth',0.1, 'HandleVisibility', 'Off');
            flag = 1;
            break
        end
    end
end  
if multiplier == 100;
   xlabel('Panel zone shear distortion, \gamma [% rad]');
else
   xlabel('Panel zone shear distortion, \gamma [rad]');
end 
ylabel('Panel zone shear force, \itV_{pz}\rm [kN]');
set(gca,'FontSize',12);set(gca,'FontName','Times New Roman');
set(gca,'GridLineStyle','--')
set(gca,'color','white'); % transparent background
lgd = legend('First cycle envelope', 'Fitted curve');
lgd.Location = 'bestoutside';
ax = gca;
xlim([0 0.01].*multiplier);
ylim([0 3000]);
ax.YTick = [0:500:3000];
ax.XTick = [0:0.0025:0.01].*multiplier;

%% Tags in the figure

htext = text(shear_yield_offset.*multiplier+0.05,5,'0.35%');          %// Add text label
set(htext,'VerticalAlignment','bottom', 'HorizontalAlignment','left', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_gy.*multiplier,5,'\gamma_{\ity}');          %// Add text label
set(htext,'VerticalAlignment','bottom', 'HorizontalAlignment','right', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_4gy.*multiplier,5,'4\gamma_{\ity}');          %// Add text label
set(htext,'VerticalAlignment','bottom', 'HorizontalAlignment','right', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_drop.*multiplier,V_drop*1.2,5,'\itV_{y,m,drop}');          %// Add text label
set(htext,'VerticalAlignment','top', 'HorizontalAlignment','center', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_4gy.*multiplier,V_4gy*1.2,5,'\itV_{p,m}');          %// Add text label
set(htext,'VerticalAlignment','top', 'HorizontalAlignment','center', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_gy.*multiplier,V_gy,5,'\itV_{y,m}');          %// Add text label
set(htext,'VerticalAlignment','top', 'HorizontalAlignment','center', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
htext = text(x_0_35.*multiplier,V_0_35*1.2,5,'\itV_{0.35%,m}');          %// Add text label
set(htext,'VerticalAlignment','top', 'HorizontalAlignment','center', 'Rotation', 0, 'FontSize', 12,'FontName','Times New Roman', 'FontWeight', 'normal')
