clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Program to reduce OpenSees free field data generated using freefield.tcl
% Created by Lindsay Baynes (lbaynes@u.washington.edu)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load input acceleration and select depth for profile plots
  gConvert = 9.8/1;                             % Factor needed to scale max acc and convert acc unit to m/s^2
  baseAcc = load('quak.dat')*gConvert;          % Make sure earthquake motion is same as that used for analysis
  depMax = 20;                                  % Choose max depth (m) for acceleration colormap and other depth plots
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analysis options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Stress/strain analysis (stress vs strain and q v p)
analOpt1 = 1;               % 0 = off, 1 = on 
elements1 = [185 164 159];  % Specify elements to analyze (increase with elevation)

% Displacement analysis (time histories)
analOpt2 = 1;           % 0 = off, 1 = on 
nodes2 = [221];         % Specify nodes to analyze

% Acceleration analysis (time histories, response spectra)
analOpt3 = 1;           % 0 = off, 1 = on.  Base acceleration always plotted when analysis is performed.
nodes3 = [221 1];       % Specify nodes to analyze
                        
% Pressure analysis (time histories)
analOpt4 = 1;               % 0 = off, 1 = on 
elements4 = [185 164 159];  % Specify elements to analyze (increase with elevation)
                          
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Data import
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Import parameters from OpenSees
param = load('outputParam.out');
numLayer = length(param(1,:))-2;         % number of soil layers 
drnLayer = param(1,1);                   % number of layers above GWT
layerDepth = param(2,:);                 % depths of layers top to bottom
deltaT = param(3,1);

% Build arrays of OpenSees output
stressArray = load([num2str(numLayer),'.stress']);      % Initialize array of element stresses.  Column 1 = time;
[numSteps,sDim] = size(stressArray);                    % remaining columns = sigma_x, _y, _z, _xy, and nu_r for all elements

strainArray = load([num2str(numLayer),'.strain']);      % Initialize array of element strains.  Column 1 = time;
eDim = size(strainArray,2);                             % remaining columns = epsilon_xx, _yy, and gamma_xy for all elements

if drnLayer > 0                                         % Initialize array of element pressures. Column 1 = time;
    pressArray = load([num2str(numLayer),'.pressure']); % remaining columns = excess pwp and excess pwp ratio for all elements
    pDim = size(pressArray,2);                          % below the GWT
end

numEle(1:numLayer) = zeros(1,numLayer);              % Initialize array containing number of elements in each layer
numEle(numLayer) = (sDim-1)/5;   

PDMYArray = load([num2str(numLayer),'.pdmy']);      % Initialize PDMY parameter array. Row 1 = relative density; row 2 = 
                                                    % saturated mass density; row 3 = ref shear mod; row 4 = ref bulk mod;
                                                    % row 5 = friction angle; row 6 = peak shear strain; row 7 = ref pressure; 
                                                    % row 8 = press dep coe; row 9 = PT angle; row 10 = contrac; row 11 = 
                                                    % dilat 1; row 12 = dilat 2; row 13 = liquefac 1; row 14 = liquefac 2;
                                                    % row 15 = liquefac 3; row 16 = void ratio

% Create time vector
time = stressArray(:,1);

for ii = numLayer-1:-1:1
%   Stress array
    tempStress = load([num2str(ii),'.stress']);         
    b = size(tempStress,2);
    tempStress = tempStress(:,2:b);
    sDim_new = b-1;
    stressArray(:,sDim+1:sDim+sDim_new) = tempStress;
    sDim = sDim+sDim_new;
%   Strain array
    tempStrain = load([num2str(ii),'.strain']);         
    b = size(tempStrain,2);
    tempStrain = tempStrain(:,2:b);
    eDim_new = b-1;
    strainArray(:,eDim+1:eDim+eDim_new) = tempStrain;
    eDim = eDim+eDim_new;
%   Pressure array
    if ii > drnLayer
        tempPress = load([num2str(ii),'.pressure']);
        b = size(tempPress,2);
        tempPress = tempPress(:,2:b);
        pDim_new = b-1;
        pressArray(:,pDim+1:pDim+pDim_new) = tempPress;
        pDim = pDim+pDim_new;
    end
%   Element per layer array    
    numEle(ii) = sDim_new/5;
%   PDMY parameter array 
    PDMYArray(:,ii) = load([num2str(ii),'.pdmy']);    
end

% Load displacement and acceleration arrays
dispArray = load('disp.out');
accArray = load('accel.out');
% Eliminate unnecessary data (data duplicated for coupled nodes)
[a,b] = size(dispArray);
dispArray = dispArray(:,2:b);
accArray = accArray(:,2:b);
tempD = zeros(a,(b-1)/2);
dispHor = zeros(a,(b-1)/4);
tempA = tempD;
accHor = dispHor;
j = 1;
k = 1;
for ii = 1:2:((b-1)/2)
    tempD(:,ii:ii+1) = dispArray(:,j:j+1);
    dispHor(:,k) = dispArray(:,j);
    tempA(:,ii:ii+1) = accArray(:,j:j+1);
    accHor(:,k) = accArray(:,j)+baseAcc;        % Add input motion for total acceleration
    j = j+4;
    k = k+1;
end
dispArray = tempD;
accArray = tempA;

% Construct array identifying elements by layer and element number
ii = 1;
numDrnEl = 0;
while ii <= sum(numEle)
    for j = numLayer:-1:1
        layerThick(j) = layerDepth(j+1)-layerDepth(j);
        for k = 1:numEle(j)
            eleDepth(ii) = layerDepth(j+1)-layerThick(j)/numEle(j)*(k-1);
            eleID(ii,:) = [j,ii,eleDepth(ii)];      % Characterize element by layer, element #, depth
            if eleDepth(ii) == depMax
                dMaxID = ii;
            end
            ii = ii+1;
        end
        if j <= drnLayer
            numDrnEl = numDrnEl+numEle(j);
        end
    end
end



    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stress/strain analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Determine elements for analysis
    
    figure(n); clf;
    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperPosition', [1 1 5.8 3.8]);
    for j = 1:sum(numEle)
    %   Compose array of shear strains for selected elements
        exy(:,j) = strainArray(:,(j-1)*3+4);
    %   Compose arrays of normal and shear stresses for selected elements       
        sxx(:,j) = stressArray(:,(j-1)*5+2);
        syy(:,j) = stressArray(:,(j-1)*5+3);
        szz(:,j) = stressArray(:,(j-1)*5+4);
        sxy(:,j) = stressArray(:,(j-1)*5+5);
    %   Confinement stress array        
        po(:,j) = (sxx(:,j)+syy(:,j)+szz(:,j))/3;
    %   Deviatoric stress array        
        for ii = 1:length(sxx)
            qo(ii,j) = ((sxx(ii,j)-syy(ii,j))^2+(syy(ii,j)-szz(ii,j))^2+(sxx(ii,j)-szz(ii,j))^2) + 6.0*sxy(ii,j)^2;
            qo(ii,j) = sign(sxy(ii,j))*1/sqrt(2)*qo(ii,j)^0.5;
        end
    end
    if analOpt1 > 0
    for j = 1:length(elements1)
    %   Plot relationships
        % Stress-strain
        subplot(length(elements1),2,j*2-1)
        plot(2*exy(:,elements1(j))*100,sxy(:,elements1(j)));
        if j == length(elements1)
            xLabel('\gamma (%)');
        end
        yLabel('\tau (kPa)');
        xmin = -5; xmax = 5; ymin = -150; ymax = 150;
        axis([xmin xmax ymin ymax])
        t1 = xmin+3/4*(xmax-xmin); t2 = ymin+4/5*(ymax-ymin);
        text(t1,t2,[num2str((eleID(elements1(j),3)),'%-1.1f'),' m'])
                
        % Deviatoric vs. confinement stress
        subplot(length(elements1),2,j*2), plot(-po(:,elements1(j)),qo(:,elements1(j))); 
        if j == length(elements1)
            xLabel('p (kPa)');
        end
        yLabel('q (kPa)');
        xmin = 0; xmax = 150; ymin = -150; ymax = 150;
        axis([xmin xmax ymin ymax])
        t1 = xmin+3/4*(xmax-xmin); t2 = ymin+4/5*(ymax-ymin);
        text(t1,t2,[num2str((eleID(elements1(j),3)),'%-1.1f'),' m'])
    end
    n = n+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Displacement analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if analOpt2 > 0
%   Determine nodes for analysis
    x = zeros(length(dispArray),length(nodes2));
    for j = 1:length(nodes2)
    %   Compose arrays of horizontal displacements for selected nodes
        x(:,j) = dispArray(:,(nodes2(j)*2));
    %   lateral displacement at element top
        figure(n); clf;
        plot(time,x(:,j)*100,'r');
        title (['Lateral displacement (',num2str((eleID(nodes2(j)-1,3)-layerThick(eleID((nodes2(j)-1),1))/numEle(eleID((nodes2(j)-1),1))),'%-1.1f'),' m)']);
        xLabel('Time (s)'); yLabel('Disp. (cm)'); 
        n = n+1;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Acceleration analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if analOpt3 > 0
    figure(n); clf;
    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperPosition', [1 1 5.8 3.8]);
    xdd = zeros(length(accArray),length(nodes3));
    for j = 1:length(nodes3)
    %   Compose arrays of horizontal accelerations for selected nodes
        xdd(:,j) = accArray(:,(nodes3(j)*2-1))+baseAcc;
    %   Compute repsonse spectra for selected nodes
        N=length(xdd(:,j));
        for jj=2:N+1                 % add initial value of zero and rename
            aa(jj)=xdd(jj-1,j);
        end
        aa(1)=0.;
        dw=2.*3.14159/(N*deltaT);
        nn=100;                     % number of periods at which spectral values are to be computed       
        delta_ee = 4/nn;
        ee(1) = -3;
        p(1) = 0.001;               % set initial period just greater than zero
        for jj=2:nn
            ee(jj) = ee(jj-1)+delta_ee;
            p(jj) = 10^ee(jj);
        end

        k=1000.;                 % choose arbitrary value of k
        xe=0.05;                 % specify damping ratio
        afft=fft(aa);             % compute fft of input motion
        t=0.:deltaT:N*deltaT;    % set up time vector
        w=0.:dw:N*dw;            % set up vector of circular frequencies
        for kk=1:nn  
            m(kk)=((p(kk)/(2*3.14159))^2)*k;            % compute mass values required to produce desired periods	
            c(kk)=xe*2*(k*m(kk))^0.5;                   % compute dashpot coefficients required to produce desired damping ratios  
            qfft=-m(kk)*afft;	                        % define fft(q)	
            for jj=2:N/2+1
                h(jj)=1./(-m(kk)*w(jj)*w(jj)+i*c(kk)*w(jj)+k);      % compute transfer function
                h(N+3-jj)=conj(h(jj));					            % mirror image of transfer function
            end
            h(1)=1./k;
            for jj=1:N+1
                u(jj)=h(jj)*qfft(jj);	                    % compute displacement in frequency domain using transfer function
            end
            utime=real(ifft(u));		                    % compute displacement in time domain (neglect small imaginary part)
            umax(kk)=max(abs(utime));			            % spectral displacement
            vmax(kk)=2.*3.14159*umax(kk)/p(kk);		        % spectral velocity
            amax(kk,j)=2.*3.14159*vmax(kk)/p(kk);		    % spectral acceleration
        end
    end
    for j = 1:length(nodes3)
    %   plot time histories and response spectra
        % time history
        subplot(length(nodes3),2,j*2-1)
        plot(time,xdd(:,j));
        yLabel(['Acc (m/s^{2})']);
        if j == length(nodes3)
            xLabel('Time (s)'); 
        end
        xmin = 0; xmax = max(time); ymin = -9; ymax = 5;
        axis([xmin xmax ymin ymax])
        t1 = xmin+3/4*(xmax-xmin); t2 = ymin+4/5*(ymax-ymin);
        if nodes3(j) > 1
            text(t1,t2,[num2str((eleID(nodes3(j)-1,3)-layerThick(eleID((nodes3(j)-1),1))/numEle(eleID((nodes3(j)-1),1))),'%-1.1f'),' m'])
        else
            text(t1,t2,[num2str((eleID(1,3)),'%-1.1f'),' m'])
        end
        % spectral acceleration
        subplot(length(nodes3),2,j*2)
        plot(p,amax(:,j))
        if j == length(nodes3)
            xlabel('Period (s)'); 
        end
        yLabel('S_a (m/s^{2})');
        xmin = 0; xmax = max(p); ymin = 0; ymax = 25;
        axis([xmin xmax ymin ymax])
        t1 = xmin+3/4*(xmax-xmin); t2 = ymin+4/5*(ymax-ymin);
        if nodes3(j) > 1
            text(t1,t2,[num2str((eleID(nodes3(j)-1,3)-layerThick(eleID((nodes3(j)-1),1))/numEle(eleID((nodes3(j)-1),1))),'%-1.1f'),' m'])
        else
            text(t1,t2,[num2str((eleID(1,3)),'%-1.1f'),' m'])
        end
    end
    n = n+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pressure analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for j = 1:sum(numEle)-numDrnEl
%   Compose array of excess pwp and excess pwp ratio for selected elements
    epp(:,j) = pressArray(:,(j-1)*2+2);
    eppr(:,j) = pressArray(:,(j-1)*2+3);
end
if analOpt4 > 0
    figure(n); clf;
    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'PaperUnits', 'inches');
    set(gcf, 'PaperPosition', [1 1 2.9 5.4]);
    for j = 1:length(elements4)
        if eleID(elements4(j),1) > drnLayer
            subplot(length(elements4),1,j),plot(time,eppr(:,elements4(j)));
            if j == length(elements4)
                xLabel('Time (s)'); 
            end
            yLabel('r_u');
            xmin = 0; xmax = max(time); ymin = 0; ymax = 1.5;
            axis([xmin xmax ymin ymax])
            t1 = xmin+3/4*(xmax-xmin); t2 = ymin+4/5*(ymax-ymin);
            text(t1,t2,[num2str((eleID(elements4(j),3)),'%-1.1f'),' m'])
        end
    end
    n = n+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vs depth analysis - to be performed regardless of node/element selections at beginning of script
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Maximum excess pore pressure ratio, strain, and horizonal displacement with depth ---------------------------------------

figure(n);clf;
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [1 1 5.8 3.8]);
% Excess pore pressure ratio
subplot(1,3,1)
plot([max(eppr(:,dMaxID:(sum(numEle)-numDrnEl))) zeros(1,numDrnEl)],eleID(dMaxID:sum(numEle),3),'LineWidth',2);
axis ij
xmin = 0; xmax = 1.5; ymin = 0; ymax = depMax;
axis([xmin xmax ymin ymax])
xLabel('r_u'); yLabel('Depth (m)');

% Plot maximum shear strain vs depth
subplot(1,3,2)
plot(max(exy(:,dMaxID:size(exy,2)))*100,eleID(dMaxID:length(eleID),3),'LineWidth',2);
axis ij
xmin = 0; xmax = 6; ymin = 0; ymax = depMax;
axis([xmin xmax ymin ymax])
xLabel('\gamma_m_a_x (%)');

% Plot maximum horizontal displacement versus depth
dMax = max(abs(dispHor));
subplot(1,3,3)
plot(dMax(dMaxID:length(dMax))*100,[eleID(dMaxID:length(eleID),3)' 0],'LineWidth',2);
axis ij
xmin = 0; xmax = 30; ymin = 0; ymax = depMax;
axis([xmin xmax ymin ymax])
xLabel('d_m_a_x (cm)');
n = n+1;

% Vertical deformation with time -----------------------------------------------------------------
figure(n);clf;
set(gcf, 'PaperPositionMode', 'manual');
set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPosition', [1 1 5.8 2.5]);
plot(time,dispArray(:,size(dispArray,2)))
xmin = 0; xmax = max(time); ymin = -0.1; ymax = 0;
axis([xmin xmax ymin ymax])
xLabel('Time (s)');ylabel('Vert disp (m)');
n = n+1;