clear;           % clear all variables from memory
close all;       % close all open windows
clc;

tic

for K =1:32
%         K=11;
    if K == 1
        fc0 = 20;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 2
        fc0 = 20;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 3
        fc0 = 20;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 4
        fc0 = 20;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 5
        fc0 = 20;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 6
        fc0 = 20;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 7
        fc0 = 20;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 8
        fc0 = 20;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 9
        fc0 = 30;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 10
        fc0 = 30;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 11
        fc0 = 30;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 12
        fc0 = 30;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 13
        fc0 = 30;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 14
        fc0 = 30;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 15
        fc0 = 30;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 16
        fc0 = 30;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 17
        fc0 = 50;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 18
        fc0 = 50;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 19
        fc0 = 50;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 20
        fc0 = 50;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 21
        fc0 = 50;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 22
        fc0 = 50;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 23
        fc0 = 50;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 24
        fc0 = 50;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 25
        fc0 = 70;
        D = 150;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 26
        fc0 = 70;
        D = 150;
        fultf = 600;
        Ef = 2.5e4;
        dst = 5;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 27
        fc0 = 70;
        D = 300;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 28
        fc0 = 70;
        D = 300;
        fultf = 600;
        Ef = 2.5e4;
        dst = 8;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 29
        fc0 = 70;
        D = 600;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 30
        fc0 = 70;
        D = 600;
        fultf = 600;
        Ef = 2.5e4;
        dst = 13;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 31
        fc0 = 70;
        D = 1200;
        fultf = 1.2e3;
        Ef = 1e5;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
    
    if K == 32
        fc0 = 70;
        D = 1200;
        fultf = 600;
        Ef = 2.5e4;
        dst = 15;
        run_analises(fc0,D,fultf,Ef,dst,K);
    end
end
toc

function parsave(fname,P)
save(fname,'P')
end

function [] = run_analises(fc0,D,fultf,Ef,dst,K)
cov = 25;
Re = D/2;
Rc = Re - cov;

if D == 150; nb = 4;
elseif D == 300; nb = 8;
elseif D == 600; nb = 16;
elseif D == 1200; nb = 32;end

rho_l = 0.01;
dsl = sqrt(rho_l*D^2/nb);
Esl = 2e5;
fyl = 450;
b = 0.005;

rho_s = [0.0000000000000000000000000000001 0.5 1 2 3 4];
Est = 2e5;
fyt = 450;
eults = 0.1;

if mod(K, 2) == 0
    rho_f = [0 0.50 1.0 2 4 6];
else
    rho_f = [0 0.25 0.5 1 2 3];
end

s = 4*(0.25*pi*dst^2)./(0.01*rho_s*(D-2*cov));
tf = D*rho_f/400;
if abs(fc0) <= 40
    eff = (2*tf*Ef/D)*(-0.00003)+0.6806;
else
    eff = (2*tf*Ef/D)*(-0.00003)+0.6378;
end

% eff = 0.004/(fultf/Ef)*ones(1,6);

ec0 = 0.00076+((0.626*fc0-4.33)*10^-7)^0.5;
Ec = 5000*fc0.^0.5;

opt_steel = 2; %1 for spirals and 2 for ties

ke = (1-0.5*(s-dst*ones(1,6))./(D-2*cov));
ke(ke<0)=0;
ke(ke>0)=1;

ke = ke.*(1-0.5*(s-dst*ones(1,6))./(D-2*cov)).^opt_steel/(1-nb*dsl^2/D^2);
ke(ke>1)=1;

fls = 0.5*ke.*rho_s*450/100;
fls = fls.*[1 1 1 1 1 1]';
flf = 0.5*rho_f*fultf.*eff/100;
flf = flf'.*[1 1 1 1 1 1];
Kf = 0.5*rho_f*Ef;
Kf = Kf'.*[1 1 1 1 1 1];
Ks = 0.5*Est*rho_s.*ke;
Ks = Ks.*[1 1 1 1 1 1]';
flf_fc = flf/fc0;

% Cf = (Ks./Kf')*((D-2*cov)^2/D^2);

cf = 100*(fls./flf)*((D-2*cov)^2/D^2);
Cs = 100*100*(fls./(0.5*rho_f'*Ef/100))*((D-2*cov)^2/D^2);

ecc = [0.05 0.15 0.25]*D;

% load(sprintf('P_%d_complete.mat',K));

for ff = 1:6
    for ss = 1:6
        for ee = 1:3
                                                    ff = 1;
                                                    ss = 4;
                                                    ee = 1;
            tf1 = tf(ff);
            eff1 = eff(ff);
            s1 = s(ss);
            ecc1 = ecc(ee);
            
            if ff <= 3 && ss <= 3 && ee==3
                eccu = 0.004;
            elseif ff <= 3 && ss <= 3
                eccu = 0.011;
            elseif ff > 3 && ss > 3 && ee==3
                eccu = 0.005;
            else
                eccu = 0.022;
            end
            
            if ee == 3
                numSteps = 1000;
            else
                numSteps = 500;
            end

            eccu = 0.01;
            
            filename = 'Run_ResistanceModel.tcl';
            exename = 'OpenSees.exe';
            cmd = sprintf('"%s" < "%s"', exename, filename);
            
            write_file(fc0,ec0,Ec,fyl,Esl,b,dsl,fyt,Est,eults,dst,Ef,fultf,eff1,D,s1,cov,tf1,eccu,nb,ecc1,opt_steel,numSteps);
            system(cmd);
            
            
            fileID = fopen('P_disp.out');
            D2 = textscan(fileID,'%f %f %f %f');
            fclose(fileID);
            
            p1 = D2{1,2};
            [S,I] = max(p1);

            m1 =D2{1,4};
            
            p(ff,ss,ee) = S;
            m(ff,ss,ee) = m1(I);
     
        end
    end
end
P = p;
P(:,:,:,2) = m;
P(:,:,:,3)=m./(p*D);

parsave(sprintf('P_%d_complete.mat',K),P);
end

function [] = write_file(fc0,ec0,Ec,fyl,Esl,b,dsl,fyt,Est,eults,dst,Ef,fultf,eff,D,s,c,tf,eccu,nbars,ecc,opt_steel,numSteps)
fid=fopen('Run_ResistanceModel.tcl','w');
fprintf(fid,'wipe;\n');
fprintf(fid,'set fc0 %0.8e;\n', -1*fc0);
fprintf(fid,'set ec0 %0.8e;\n', -1*ec0);
fprintf(fid,'set Ec %0.8e;\n', Ec);
fprintf(fid,'set fyl %0.8e;\n', fyl);
fprintf(fid,'set Esl %0.8e;\n', Esl);
fprintf(fid,'set b %0.8e;\n', b);
fprintf(fid,'set dsl %0.8e;\n', dsl);
fprintf(fid,'set fyt %0.8e;\n', fyt);
%fprintf(fid,'set fyt %0.8e;\n', 0);
fprintf(fid,'set Est %0.8e;\n', Est);
fprintf(fid,'set eults %0.8e;\n', eults);
fprintf(fid,'set dst %0.8e;\n', dst);
fprintf(fid,'set Ef %0.8e;\n', Ef);
fprintf(fid,'set fultf %0.8e;\n', fultf);
fprintf(fid,'set eff %0.8e;\n', eff);
fprintf(fid,'set D %0.8e;\n', D);
fprintf(fid,'set c %0.8e;\n', c);
fprintf(fid,'set s %0.8e;\n', s);
fprintf(fid,'set tf %0.8e;\n', tf);
fprintf(fid,'set eccu %0.8e;\n', -1*eccu);
fprintf(fid,'set nbars %d;\n', nbars);
fprintf(fid,'set ecc %0.8e;\n', ecc);
fprintf(fid,'set opt_steel %0.8e;\n', opt_steel);
fprintf(fid,'set numSteps %d;\n', numSteps);

fprintf(fid,'source Model_crossectional_compression_controlled.tcl;\n');
% fprintf(fid,'source Model_crossectional_tension_controlled.tcl;\n');
fclose(fid);
end