% Monte Carlo Latin Hypercube fatige study program
% Must use SimulationCreator.m to create Latin Hypercube samples before using
% Program created for use on Condor High Throughput Network

% Sara Roberts
% Graduate Student
% Clemson University
% sarar@g.clemson.edu
% Date Created:  June 2010
% Enviro = Corrosion environment choice (1, 2 or 3), refer to SimulationCreator.m
% LPchoice = Shear key Degradation choice (1, 2, 3, 4 or 5) refer to SimulationCreator.m
% simStart = starting simulation number
% stride = number of simulations done in one process
% process = process number, used in Condor program to split simulations and send to 
% multiple computers. If not using this through Condor, make process 1


function[] = HCfatigueLHCRt(Enviro, LPchoice, simStart, stride, process)

sims = 10000;

Dps = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 0 (sims-1) 0]);
b = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 1 (sims-1) 1]);
%d = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 2 (sims-1) 2]);
fc = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 3 (sims-1) 3]);
Ag = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 4 (sims-1) 4]);
S = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 5 (sims-1) 5]);
Md = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 6 (sims-1) 6]);
n = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 7 (sims-1) 7]);
e = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 8 (sims-1) 8]);
Ss = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 9 (sims-1) 9]);
Les = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 10 (sims-1) 10]);
Lc = dlmread(['Sims', num2str(Enviro), num2str(LPchoice), '.txt'], '\t', [0 11 (sims-1) 11]);
Ls = dlmread(['Sims', num2str(Enviro), num2str(LPchoice), '.txt'], '\t', [0 12 (sims-1) 12]);
T1 = dlmread(['Sims', num2str(Enviro), num2str(LPchoice), '.txt'], '\t', [0 13 (sims-1) 13]);
rcorr = dlmread(['Sims', num2str(Enviro), num2str(LPchoice),'.txt'], '\t', [0 14 (sims-1) 14]);

Year = zeros([sims 1]); DIc = zeros([sims 1]); DIs = zeros([sims 1]);
fcf = zeros([sims 1]);

%Simulation Range
start = simStart + (process*stride);
finish = simStart + (process + 1)*stride - 1;

% Bridge Span
L = 50;  %ft

% Ultimate Prestress 
fpu = 270;  %ksi

% Tractor Load: lognormal distribution
MUpt1 = 5; %kips
Vpt1 = 0.12;
%SIGpt1 = MUpt1*Vpt1;
SIGLNpt1 = sqrt(log(Vpt1^2 + 1));
MULNpt1 = log(MUpt1) - 0.5*SIGLNpt1^2;

% Trailer Load: Wiebull Extreme Smallest Distribution
MUpt2 = 28;  %kips
Vpt2 = 0.20;
SIGpt2 = MUpt2*Vpt2;
wpt2 = 12;  %kips, lower bound
ka = 0.1; kb = 30; err = 1;
while err > 0.01
    ua = (MUpt2 - wpt2 + wpt2*gamma(1 + 1/ka))/gamma(1 + 1/ka);
    ub = (MUpt2 - wpt2 + wpt2*gamma(1 + 1/kb))/gamma(1 + 1/kb);
    SIGa = sqrt((ua - wpt2)^2*(gamma(1 + 2/ka) - (gamma(1 + 1/ka))^2));
    SIGb = sqrt((ub - wpt2)^2*(gamma(1 + 2/kb) - (gamma(1 + 1/kb))^2));
    kx = (ka + kb)/2;
    ux = (MUpt2 - wpt2 + wpt2*gamma(1 + 1/kx))/gamma(1 + 1/kx);
    SIGx = sqrt((ux - wpt2)^2*(gamma(1 + 2/kx) - (gamma(1 + 1/kx))^2));
    if (SIGa - SIGpt2)*(SIGx - SIGpt2) < 0
        kb = kx;
    elseif (SIGb - SIGpt2)*(SIGx - SIGpt2) < 0
        ka = kx;
    end
    err = abs(SIGx - SIGpt2);
end
kpt2 = kx; upt2 = ux;

% Average trucks per hour
AADT = 3000;  %vehicles/day
Ptrucks = 10;  %Percent trucks
VPH = AADT/24;  %vehicles/hr
TPH = ceil(VPH*Ptrucks/100);   %trucks/hr

% Average number of trucks per year
Its = TPH*8760;

for k = start:finish
    fprintf('%g\n',k)    
    tc=0;
    counter = 0;
    while tc < (438000 - 168) && DIc(k) < 1 && DIs(k) < 1 && fcf(k) ~= 1
        IT = exprnd(1/TPH, [its 1]);  % Interarrival time generation
        tv = cumsum(IT);  %cumulative time
        fprintf('tc = %g\n', tc/8760)
        % Cycle Vector
        Cyc = (1:its)' + counter*its;
        counter = counter + 1;
        
        % Diameter of wires (in)
        if (tc/8760 + .05) <= T1(k)
            Dpw = Dps(k);
        elseif (tc/8760 + .05) > T1(k) && (tc/8760 + .05) < (T1(k) + Dps(k)/(rcorr(k)/1000))
            Dpw = Dps(k) - (rcorr(k)/1000)*(tc/8760 - T1(k));
        else
            Dpw = 0;
        end
        
        % Area of prestress steel (18 – 7 wire strands)
        Aps = 18*7*pi/4*Dpw^2; %in^2
        
        % Relaxation Loss
        fpi = 600*fpu; %psi
        if tc < 1
            Lr = 0;
        else
            Lr = fpi*log10(tc)/10*(fpi/(fpu*900) - 0.55); %psi
        end
        
        % Prestress Loss
        pl = Les(k) + Lc(k) + Ls(k) + Lr; %psi
        
        % Prestress Force
        P = Aps*(.75*fpu - pl/1000); %kips
        
        % Dead stress in top of section
        ftD = P*(1/Ag(k) - e(k)/S(k)) + Md(k)/S(k);  %ksi

        %Truck Load Simulation
        Pt1 = lognrnd(MULNpt1, SIGLNpt1, its, 1);
        Fx = rand(its, 1);
        Pt2 = wpt2 + (upt2 - wpt2).*(-log(1 - Fx)).^(1/kpt2);
        impact = unifrnd(0.10, 0.15, its, 1);
        Pt1 = Pt1.*(1 + impact);
        Pt2 = Pt2.*(1 + impact);
        
        % Truck Load Moment
        Ml = ((27.3.*(L.*(Pt1 + 0.7436.*Pt2) - 41.3.*(Pt1 + 0.4915.*Pt2)))./L).*12; %kip-in
        
        % Cracking Moment
        Mcr = S(k)*(P/Ag(k) + P*e(k)/S(k) + 7.5*sqrt(fc(k))/1000); %kip-in
        
        % Live stress in top of section (max stress)
        ftL = zeros([its 1]);
        dfbL = zeros([its 1]);
        
        if LPchoice == 1
            Mle = 1/3.*Ml;
        elseif LPchoice == 2
            Mle = Ml;
        elseif LPchoice == 3
            if Cyc(its) < 20000000 
                LP = (1/30000000).*Cyc + 1/3;
            elseif Cyc(1) < 20000000 
                ind = 20000000 - Cyc(1);
                LP(1:ind - 1) = (1/30000000).*Cyc(1:ind - 1) + 1/3;
                LP(ind:its) = 1;
            else
                LP = 1;
            end
            Mle = LP.*Ml;
        elseif LPchoice == 4
            if Cyc(its) < 35000000 
                LP = (1/52500000).*Cyc + 1/3;
            elseif Cyc(1) < 35000000 
                ind = 35000000 - Cyc(1);
                LP(1:ind - 1) = (1/52500000).*Cyc(1:ind - 1)+1/3;
                LP(ind:its) = 1;
            else
                LP = 1;
            end
            Mle = LP.*Ml;
        elseif LPchoice == 5
            if Cyc(its) < 50000000 
                LP = (1/75000000).*Cyc + 1/3;
            elseif Cyc(1) < 50000000 
                ind = 50000000 - Cyc(1);
                LP(1:ind - 1) = (1/75000000).*Cyc(1:ind - 1) + 1/3;
                LP(ind:its) = 1;
            else
                LP = 1;
            end
            Mle = LP.*Ml;
        end
        Mle = sort(Mle);
        MT = Mle + Md(k);
        if MT(its) < Mcr 
            ftL = ftD + Mle./S(k); %ksi
            dfbL = Mle./Ss(k); %ksi
        else
            crack = find(MT>Mcr,1,'first'); 
            ftL(1:(crack-1)) = ftD+Mle(1:(crack-1))./S(k); %ksi
            dfbL(1:(crack - 1)) = Mle(1:(crack - 1))./Ss(k); %ksi
            
            Atps = n(k)*Aps;
            MTi = (floor(MT(crack)):ceil(MT(its)))';
            sMTi = length(MTi);
            ftLi = zeros(sMTi);
            dfbLi = zeros(sMTi);
            for l = 1:sMTi
                Ca = 0.1; cb = 40; err = 1;
                while err > 0.5
                    cx = (ca + cb)/2;
                    ta = ca - 4.5;
                    tb = cb - 4.5;
                    tx = cx - 4.5;
                    if ta <= 0
                        pcirca = 0;
                        ycirca = 0;
                        Ica = 0;
                    elseif ta > 0 && ta < 12
                        pcirca = -0.000638*ta^3 + 0.011489*ta^2 + 0.038572*ta - 0.007160;
                        ycirca = -0.011*ta^2 + 0.658*ta + 4.491;
                        Ica = -0.115*ta^4 + 3.204*ta^3 - 17.28*ta^2 + 34.15*ta - 15.9;
                    else
                        pcirca = 1;
                        ycirca = 10.5;
                        Ica = pi*6^4/4;  %(6 refers to the radius of the void)
                    end
                    if tb <= 0
                        pcircb = 0;
                        ycircb = 0;
                        Icb = 0;
                    elseif tb > 0 && tb < 12
                        pcircb = -0.000638*tb^3 + 0.011489*tb^2 + 0.038572*tb - 0.007160;
                        ycircb = -0.011*tb^2 + 0.658*tb + 4.491;
                        Icb = -0.115*tb^4 + 3.204*tb^3 - 17.28*tb^2 + 34.15*tb - 15.9;
                    else
                        pcircb = 1;
                        ycircb = 10.5;
                        Icb = pi*6^4/4;  %(6 refers to the radius of the void)
                    end
                    if tx <= 0
                        pcircx = 0;
                        ycircx = 0;
                        Icx = 0;
                    elseif tx > 0 && tx < 12
                        pcircx = -0.000638*tx^3 + 0.011489*tx^2 + 0.038572*tx - 0.007160;
                        ycircx = -0.011*tx^2 + 0.658*tx + 4.491;
                        Icx = -0.115*tx^4 + 3.204*tx^3 - 17.28*tx^2 + 34.15*tx - 15.9;
                    else
                        pcircx = 1;
                        ycircx = 10.5;
                        Icx = pi*6^4/4;  %(6 refers to the radius of the void)
                    end
                    Ata = Atps + b(k)*ca - 2*pcirca*pi*6^2;  %(6 refers to the radius of the void, next 9 lines)
                    Atb = Atps + b(k)*cb - 2*pcircb*pi*6^2;
                    Atx = Atps + b(k)*cx - 2*pcircx*pi*6^2;
                    yta = (b(k)*ca*ca/2 - 2*pcirca*pi*6^2*ycirca + Atps*(21/2 + e(k)))/Ata;
                    ytb = (b(k)*cb*cb/2 - 2*pcircb*pi*6^2*ycircb + Atps*(21/2 + e(k)))/Atb;
                    ytx = (b(k)*cx*cx/2 - 2*pcircx*pi*6^2*ycircx + Atps*(21/2 + e(k)))/Atx;
                    Ita = b(k)*ca^3/12 + b(k)*ca*(yta - ca/2)^2 + Atps*((21/2 + e(k)) - yta)^2 - 2*Ica - 2*pcirca*pi*6^2*(yta - ycirca);
                    Itb = b(k)*cb^3/12 + b(k)*cb*(ytb - cb/2)^2 + Atps*((21/2 + e(k)) - ytb)^2 -  2*Icb - 2*pcircb*pi*6^2*(ytb - ycircb);
                    Itx = b(k)*cx^3/12 + b(k)*cx*(ytx - cx/2)^2 + Atps*((21/2 + e(k)) - ytx)^2 - 2*Icx - 2*pcircx*pi*6^2*(ytx - ycircx);
                    ynaa = ca - yta;
                    ynab = cb - ytb;
                    ynax = cx - ytx;
                    Minta = MTi(l) - P*((21/2 + e(k)) - yta);
                    Mintb = MTi(l) - P*((21/2 + e(k)) - ytb);
                    Mintx = MTi(l) - P*((21/2 + e(k)) - ytx);
                    PAta = P/Ata;
                    PAtb = P/Atb;
                    PAtx = P/Atx;
                    MIya = Minta/(Ita/ynaa);
                    MIyb = Mintb/(Itb/ynab);
                    MIyx = Mintx/(Itx/ynax);
                    diffa = PAta - MIya;
                    diffb = PAtb - MIyb;
                    diffx = PAtx - MIyx;
                    err = abs(ca - cb);
                    if diffa == 0
                        Atx = Ata;
                        Mintx = Minta;
                        Itx = Ita;
                        ynax = ynaa;
                        ytx = yta;
                        cb = ca;
                    elseif diffb == 0
                        Atx = Atb;
                        Mintx = Mintb;
                        Itx = Itb;
                        ynax = ynab;
                        ytx = ytb;
                        ca = cb;
                    elseif diffa*diffx < 0
                        cb = cx;
                    elseif diffb*diffx < 0
                        ca = cx;
                    end
                end
                ftLi(l) = P/Atx + Mintx/(Itx/ynax);
                dfbLi(l) = (Mintx/(Itx/((21/2 + e(k)) - ytx)))*n(k);
            end
            for j = crack:its
                index = find(MTi > MT(j), 1, 'first');
                ftL(j) = ftLi(index);
                dfbL(j) = dfbLi(index);
            end
        end
        
        % Compression limit check
        if max(ftL) > fc(k) 
            i = find(ftL > fc(k), 1); 
            tf = tc + sum(tv(1:i)); 
            fcf(k) = 1;
        end
        
        % # of cycles to failure for concrete in compression
        mp = 1./(0.069.*(1-ftD./ftL)); 
        Dc = 1./(10.^(mp).*(10.^((ftL/fc(k)).*-mp)));

        % # of cycles to failure for prestressing steel
        Ds = 1./(10^11.*(dfbL).^(-3.5));

        % Damage Index
        if fcf(k) == 1
            tc = tf;
            Dcc = sum(Dc(1:i)); Dsc = sum(Ds(1:i)); 
        elseif DIc(k) + sum(Dc)>1 
            CUMc = cumsum(Dc) + DIc(k);
            l = find(CUMc>1,1);
            tc = tc + sum(IT(1:l));
            Dcc = sum(Dc(1:l)); Dsc = sum(Ds(1:l));
        elseif DIs(k) + sum(Ds) > 1 
            CUMs = cumsum(Ds) + DIs(k);
            l = find(CUMs>1,1);
            tc = tc + sum(IT(1:l));
            Dcc = sum(Dc(1:l)); Dsc = sum(Ds(1:l));
        else
            tc = sum(IT) + tc;
            Dcc = sum(Dc); Dsc = sum(Ds);
        end
        DIc(k) = DIc(k) + Dcc; DIs(k) = DIs(k) + Dsc;
    end
    Year(k) = tc/8760; %years
    fid = fopen(['Sims', num2str(Enviro), num2str(LPchoice), '-', num2str(start,'%05d'), '-', num2str(finish,'%05d'), '.txt'], 'a');
    fprintf(fid,'%g\t%g\t%g\t%g\n', Year(k), DIc(k), DIs(k), fcf(k));
    fclose(fid);
end