function [SDOFIDAresult]=SDOFIDAanalysis(OSpath,SDOFdata,ACCdata,IDAdata)

% function [SDOFIDAresult]=SDOFIDAanalysis(OSpath,SDOFdata,ACCdata,IDAdata)
% perform IDA analysis for SDOF system. The result IDA curves and
% summarized IDA curves
% matjaž dolšek 26.10.2007
% 
% INPUT DATA
% 
%   - OSPath=is a path and exe file of OpenSees
%            for example: 'D:\opensees162\opensees.exe'
%   - SDOFdata is a structure containing following field
%     m    = mass of SDOF system
%     F    = forces of the F-U envelope of SDOF system
%     U    = displalement of the F-U envelope of SDOF system
%     hysteretic = parameters of hysteretic behaviour
%            Hysteretic material : [pinchX pinchY d1 d2 beta]
%            Takeda material     : [beta]
%     dynamic=[tInteg alfaM betaK betaKcomm betaKinit]
%             !!! Newmark integration scheme is used (gama=0.5, beta 0.25)
% 
%  - ACCdata is a structure containing following field
%    name= {'acc file name'...}
%    path= {'acc path'...}, e.g. 'D:/SDOF/acc/' (!!! tcl notation)
%    dt  = [dt1...] time step in acc file
% 
%  - IDAdata is a strucutre
%    dPGA      = initial PGA value for which IDA is performed (e.g. 0.05). 
%    maxPoints = maximum number of additional seismic responses (between the points which were 
%                calculated for the determiantion of capacity
%                point) for construction of IDA curves (e.g. 10)
%    dPGAc     = tolerance (in g) for the estimation of capacity point in terms
%                of PGA (e.g 0.02)
%
% OUTPUT DATA
%
%  - SDOFIDAresult is a structure containing following field:
%    maxD    - max displacement of IDA analysis
%    IM      - pga (g) of IDA analysis
%    maxDall - max displacement of all analyses
%    Mall    - pga of all analyses in (g)
%    collapse- collapse indicator (1: collapse, 2: non-convergency, 0: non-collapse 
%    tA      - time until the analysis is performed
%
%
%
% Referenced manuscripts:
%
% Dolsek M. (2010). Development of computing environment for the seismic
% performance assessment of reinforced concrete frames by using simplified
% nonlinear models. Bulletin of Earthquake Engineering (Published online
% 16 May 2010), DOI: 10.1007/s10518-010-9184-8. 
%
% Contact: mdolsek@ikpir.fgg.uni-lj.si
%
% Copyright (C) 2010 Matjaz Dolsek
%
% This program 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, version 2.
% 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.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% accelerograms

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input data for SDOF model and dynamic analysis
m=SDOFdata.m;
F=SDOFdata.F;
U=SDOFdata.D;
Cdrift=SDOFdata.Collapse;
%Cdrift=max(max(abs(SDOFdata.envelope.UNC)))*1.01;
hysteretic=SDOFdata.hysteretic;
tInteg=SDOFdata.dynamic(1);
alfaM=SDOFdata.dynamic(2);
betaK=SDOFdata.dynamic(3);
betaKcomm=SDOFdata.dynamic(4);
betaKinit=SDOFdata.dynamic(5);


    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% input file for OpenSees SDOF model
file = fopen('SDOF.tcl', 'w+');

fprintf(file,'# Model for SDOF system \n');
fprintf(file,'model BasicBuilder -ndm 3 -ndf 6 \n');
% definition of nodes and mass
fprintf(file,'#Nodes, restraints, material end element \n');
fprintf(file,'node 1  0.0 0.0 0.0 \n');
fprintf(file,'node 2  0.0 0.0 0.0 -mass %15.3e 0.0 0.0 0.0 0.0 0.0 \n',m);
fprintf(file,'fix  1  1 1 1 1 1 1\n'); %fixed node
fprintf(file,'fix  2  0 1 1 1 1 1\n'); 
%                                               Fp1	Up1	Fp2	Up2	Fp3 Up3 Fm1 Um1 Fm2	Um2 Fm3 Um3 pincX pinchY d1 d2 beta
fprintf(file,'uniaxialMaterial Hysteretic	1   %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e %15.3e \n',...
                                                F(1),U(1),F(2),U(2),F(3),U(3),-F(1),-U(1),-F(2),-U(2),-F(3),-U(3),hysteretic(1),hysteretic(2),hysteretic(3),hysteretic(4),hysteretic(5));


% definition of section aggregator and zero length element
%                                       tag mat1 code1 
fprintf(file,'section Aggregator        1     1   P  \n');
%                                      eleTag iNode jNode secTag 
fprintf(file,'element zeroLengthSection 1 	  1      2      1   \n');

% definition of recorders
fprintf(file,'# definition of recorders \n');
fprintf(file,'recorder Node -file TDsdof.out -time -node 2 -dof 1 2 disp \n');
fprintf(file,'recorder EnvelopeElement -time -file EFsdof.out -ele  1 section force \n');
fprintf(file,'recorder EnvelopeElement -time -file EDsdof.out -ele  1 section deformation \n');
fprintf(file,'recorder Element -time -file Fsdof.out -ele  1 section force \n');
fprintf(file,'recorder Element -time -file Dsdof.out -ele  1 section deformation \n');

% parameters
fprintf(file,'test NormDispIncr 1.0e-6 12 0 \n');
fprintf(file,'algorithm Newton \n');
fprintf(file,'system SparseGeneral -piv \n');
fprintf(file,'numberer RCM \n');
fprintf(file,'constraints Transformation \n');
fprintf(file,'integrator Newmark   0.5  0.25  \n');
fprintf(file,'rayleigh %15.3e %15.3e %15.3e %15.3e \n',alfaM,betaK,betaKcomm,betaKinit);
fprintf(file,'analysis Transient \n');
fprintf(file,'set ok 0 \n');
fprintf(file,'set currentTime 0.0 \n');
fprintf(file,'set napaka 1e-4; \n');
fprintf(file,'set napaka2 1e-4; \n');
fprintf(file,'set dt %15.3e; \n',tInteg);
fprintf(file,'set dt2 %15.3e; \n',tInteg);
fprintf(file,'set dt3 %15.3e; \n',tInteg/2);
fprintf(file,'set dt4 %15.3e; \n',tInteg/2);
fclose(file);



%%%%%%%%%%%%%%%%%%%%%%%%%%%
% IDA ANALYSIS
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% na=max(size(ACCdata.name));
% a=1;

% variables
maxPoints=IDAdata.maxPoints; % minmum number of points for IDA analysis
dPGAg=IDAdata.dPGA;
dPGAc=IDAdata.dPGAc;
accName=ACCdata.name;
acc=load(accName);
PGA=max(abs(acc(:,1)));
g=9.81;
dPGA=dPGAg*g; % puts dPGA in m/s2
dt=ACCdata.dt;
tFinal=dt*max(size(acc));
% accPath=char(ACCdata.path);
%T=2*pi*sqrt(m/(F(1,1)/U(1,1));

% starts analysis
IM(1)=dPGA/g;
sf(1)=(IM(1)*g)/PGA;
ACCdata.sf=sf(1);
ACCdata.tFinal=tFinal;

SDOFDynamicAnalysis(OSpath,ACCdata);
maxDt=load('dsdof.out');
tA=maxDt(end,1);
maxD(1)=max(abs(maxDt(:,2)));
%jure z
DefAll{1}(:,:)=maxDt;
ForceAll{1}(:,:)=load('Fsdof.out');


if round(tA*10)/10>=round(tFinal*10)/10
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
   disp('non-collapse case...continue with analysis');
   %%%%%%%%%%%%%%%%%%%%%%%%
   % finds collpase
   %%%%%%%%%%%%%%%%%%%%%%%%
   disp('-----Starting IDA analysis in order to find capacity point-----')
   collapse(1)=0; i=1;
   while collapse(i)==0
       i=i+1;
       IM(i)=2*IM(i-1);
       sf(i)=(IM(i)*g)/PGA;
       
        % perform analysis
        ACCdata.sf=sf(i);
        SDOFDynamicAnalysis(OSpath,ACCdata);
        maxDt=load('dsdof.out');
        tA(i)=maxDt(end,1);
        maxD(i)=max(abs(maxDt(:,2)));
        %jure z
        DefAll{i}(:,:)=maxDt;
        ForceAll{i}(:,:)=load('Fsdof.out');       
        
        

        if tA(i)>=tFinal & maxD(i)>=Cdrift | tA(i)<tFinal & maxD(i)>=Cdrift
            collapse(i)=1;
        elseif round(tA(i)*10)/10<round(tFinal*10)/10 & maxD(i)<Cdrift
            collapse(i)=2; %non-convergency
        else
            collapse(i)=0;
        end
        disp(['dynamic analysis No.' num2str(i) ', IM=' num2str(IM(i)) ', maxD=' num2str(maxD(i)) ', collapse=' num2str(collapse(i))])
        [tA(i) maxD(i) IM(i) collapse(i)];
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% with bisection finds the capacity point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tolerance=100;
if collapse(i)==1
    disp('------Calculating capacity point 1------')
    tolerance=IM(i)-IM(i-1);
    IMcmin=IM(i);   
    IMnmax=IM(i-1); 
    while tolerance>dPGAc
        i=i+1;
        IM(i)=IMcmin-0.5*(IMcmin-IMnmax);
        sf(i)=(IM(i)*g)/PGA;
        % perform analysis
        ACCdata.sf=sf(i);
        SDOFDynamicAnalysis(OSpath,ACCdata);
        maxDt=load('dsdof.out');
        tA(i)=maxDt(end,1);
        maxD(i)=max(abs(maxDt(:,2)));
        %jure z
        DefAll{i}(:,:)=maxDt;
        ForceAll{i}(:,:)=load('Fsdof.out');
        
      
       % results

       tolerance=(IMcmin-IMnmax)/2;
       disp(['------[IMcmin, IMnmax]=[' num2str(IMcmin) ', ' num2str(IMnmax) '], toleranca=[' num2str(tolerance) '], maxDt=[' num2str(maxD(i)) ']']) 
       
       
       if tA(i)>=tFinal & maxD(i)>=Cdrift | tA(i)<tFinal & maxD(i)>=Cdrift
           IMcmin=IM(i);
           collapse(i)=1; %porušitev
       elseif round(tA(i)*10)/10<round(tFinal*10)/10 & maxD(i)<Cdrift 
           IMcmin=IM(i);
           collapse(i)=2; %nekonvergenca  - v raèunu kapacitete pomeni kot porušitev
       else
           IMnmax=IM(i);
           collapse(i)=0;
       end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine additional points max(dPGA)=dPGAc or max of additional points ==maxPoints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
indexn=find(collapse==0); % index of non-collpase cases
indexc=find(collapse==1); 
SDOFIDAresult.PGAc(1,:)=[max(IM(indexn)) min(IM(indexc))];
            
IMida=IM(indexn);
IMida;
collapse;

% if tolerance<=dPGAc
disp('------Determine additional points------')
%     dIM=diff(IMida);
    %Y=10;
while max(size(IMida))<=maxPoints  %&  Y>=dPGAg
    dIM=diff(IMida);
    % dokler št. toèk ni enako maxPoint ali dolžina intervala <dPGAg % glej spodaj, prekine raèun, èe
    % max(dIM)<=dPGA/g  Ta kriterij ni sprogramiran!!!!
    [Y,I]=max(dIM);                   
    i=i+1;
    IM(i)=IMida(I)+Y/2; %dodatni IM
    sf(i)=(IM(i)*g)/PGA;
    % perform analysis
    ACCdata.sf=sf(i);
    SDOFDynamicAnalysis(OSpath,ACCdata);
    maxDt=load('dsdof.out');
    tA(i)=maxDt(end,1);
    maxD(i)=max(abs(maxDt(:,2)));
    %jure
    DefAll{i}(:,:)=maxDt;
    ForceAll{i}(:,:)=load('Fsdof.out');

%     %marko
%                 SDOFIDAresult.Dsdof{i}=maxDt;
%                 SDOFIDAresult.Fsdof{i}=load('fsdof.out');

    if tA(i)>=tFinal & maxD(i)>=Cdrift | tA(i)<tFinal & maxD(i)>=Cdrift
        collapse(i)=1; %porušitev
%                     IMida(end+1)=IM(i);
%                     IMida=sort(IMida);
%                     dIM=diff(IMida);
    elseif round(tA(i)*10)/10<round(tFinal*10)/10 & maxD(i)<Cdrift 
        collapse(i)=2; %nekonvergenca 
%                     IMida(end+1)=IM(i);
%                     IMida=sort(IMida);
%                     dIM=diff(IMida);
    else
        collapse(i)=0;
        % doda nov IM k IMida
        IMida(end+1)=IM(i);
        IMida=sort(IMida);
        dIM=diff(IMida);
    end
        disp(['dynamic analysis No.' num2str(i) ', IM=' num2str(IM(i)) ',maxD=' num2str(maxD(i)) ', collapse=' num2str(collapse(i))])
        %[tA(i) maxD(i) IM(i) collapse(i)];

        ind0=i;
        if collapse(ind0)==1 || collapse(ind0)==2
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % with bisection finds the capacity point
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            tolerance=100;
            disp('------Calculating lower capacity point------')
            ind11=find(collapse==0 & IM<=IM(i));
            ind1=find(IM(ind11)==max(IM(ind11)));
            tolerance=IM(i)-IM(ind1);
            IMcmin=IM(i);   
            IMnmax=IM(ind1); 
            while tolerance>dPGAc
                i=i+1;
                IM(i)=IMcmin-0.5*(IMcmin-IMnmax);
                sf(i)=(IM(i)*g)/PGA;
                % perform analysis
                ACCdata.sf=sf(i);
                SDOFDynamicAnalysis(OSpath,ACCdata);
                maxDt=load('dsdof.out');
                tA(i)=maxDt(end,1);
                maxD(i)=max(abs(maxDt(:,2)));
                %jure z
                DefAll{i}(:,:)=maxDt;
                ForceAll{i}(:,:)=load('Fsdof.out');

               % results

               tolerance=(IMcmin-IMnmax)/2;
               disp(['------[IMcmin, IMnmax]=[' num2str(IMcmin) ', ' num2str(IMnmax) '], toleranca=[' num2str(tolerance) '], maxDt=[' num2str(maxD(i)) ']']) 


               if tA(i)>=tFinal & maxD(i)>=Cdrift | tA(i)<tFinal & maxD(i)>=Cdrift
                   IMcmin=IM(i);
                   collapse(i)=1; %porušitev
               elseif round(tA(i)*10)/10<round(tFinal*10)/10 & maxD(i)<Cdrift 
                   IMcmin=IM(i);
                   collapse(i)=2; %nekonvergenca  - v raèunu kapacitete pomeni kot porušitev
               else 
                   IMnmax=IM(i);
                   collapse(i)=0;
               end
            end
            indexn=find(collapse==0 & IM<=IM(i)); % index of non-collpase cases
            indexc=find(IM>max(IM(indexn)));
            SDOFIDAresult.PGAc(1,:)=[max(IM(indexn)) min(IM(indexc))];
            IMida=sort(IM(indexn));
            collapse;
            disp('------Determine additional points------')
            ind0=3;
        end
end
% end

% results of analysis
% non-collapse and collapse case
SDOFIDAresult.maxDall=maxD;      
SDOFIDAresult.IMall=IM;      
SDOFIDAresult.collapse=collapse;
SDOFIDAresult.tA=tA;
SDOFIDAresult.sf=sf;
% JURE poiscem minimalno kapaciteto (min IM(kjer je collapse=1 ali collapse=2)
indexporusni=find(collapse~=0);
% JURE V kolikor je veè kapacitet najdem izmed vseh indeksov s porušitvami tistega, kjer je kapaciteta IM najmanjša
minkapaciteta=min(SDOFIDAresult.IMall(indexporusni));
indexn=find(collapse==0 & IM<=minkapaciteta);

[Y sortindexn]=sort(SDOFIDAresult.IMall(indexn)); % sorting index of non-collpase cases

%Jure Z
[Ccc,unisortindexn,ccc]=unique(SDOFIDAresult.IMall(indexn(sortindexn)));


SDOFIDAresult.tA=SDOFIDAresult.tA(indexn(sortindexn(unisortindexn)));
SDOFIDAresult.IM=SDOFIDAresult.IMall(indexn(sortindexn(unisortindexn)));
SDOFIDAresult.sf=SDOFIDAresult.sf(indexn(sortindexn(unisortindexn)));
SDOFIDAresult.maxD=SDOFIDAresult.maxDall(indexn(sortindexn(unisortindexn)));
 %for i=1:size(unisortindexn,2)
 %SDOFIDAresult.Def{i}=SDOFIDAresult.DefAll{indexn(sortindexn(unisortindexn(i)))};
 %SDOFIDAresult.Force{i}=SDOFIDAresult.ForceAll{indexn(sortindexn(unisortindexn(i)))};
 %end





    