% Simulation Creator for Monte Carlo Latin Hypercube Fatigue Study
% Makes a set number of simulated bridges for use in HCfatigueLHCRt.m

% Sara Roberts
% Graduate Student
% Clemson University
% sarar@g.clemson.edu
% Date Created:  June 2010


clear
clc
sims = 10000;
I = linspace(1, sims, sims)';
L = 50;  %ft  (must also be changed in fatigue program)
fpu = 270;  %ksi  (must also be changed in fatigue program)

% Initial Diameter of Prestress Wires, Lognormal
MUDps = 0.16682;  %in
VDps = 0.0125;
SIGDps = MUDps*VDps;
SIGLNDps = sqrt(log(VDps^2 + 1));
MULNDps = log(MUDps) - 0.5*SIGLNDps^2;
a = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
Dps1 = logninv(a1, MULNDps, SIGLNDps);
x = randperm(sims);
Dps = zeros([sims 1]);
for j = 1:sims
   Dps(j,1) = Dps1(x(j));
end

% Beam Width, Normal
MUb = 36 + 5/32;  %in  (see Table 4.1 for mean and coefficient of variation equations)
Vb = 0.25*(36 + 5/32)^(-1);
SIGb = MUb*Vb;
a  = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
b1 = norminv(a1, MUb, SIGb);
x = randperm(sims);
b = zeros([sims 1]);
for j = 1:sims
   b(j,1) = b1(x(j));
end

% Beam Depth
Mud = 21;  %in  (see Table 4.1 for mean and coefficient of variation equations)
Vd = 1/(4*21);  
SIGd = MUd*Vd;
a  = unifrnd(0,1,sims,1);
a1 = (I - 1)./sims + a./sims;
d1 = norminv(a1, MUd, SIGd);
x = randperm(sims);
d = zeros([sims 1]);
for j = 1:sims
   d(j,1) = d1(x(j));
end

% Compressive strength of concrete, Lognormal
MUfc = 0.96*(1.205 + 0.108*1)*6;  %ksi (6 refers to design strength of concrete in ksi)
Vfc = 0.15;  %(see Table 4.1 for mean and coefficient of variation equations)
SIGfc = MUfc*Vfc;
SIGLNfc = sqrt(log(Vfc^2 + 1));
MULNfc = log(MUfc) - 0.5*SIGLNfc^2;
a  = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
fc1 = logninv(a1, MULNfc, SIGLNfc);
x = randperm(sims);
fc = zeros([sims 1]);
for j = 1:sims
   fc(j,1) = fc1(x(j));
end

% Gross Area of Section
Ag = b.*d - 2*pi*6^2;  %in^2  (6 refers to the radius of the void)

% Elastic Section Modulus for top & bottom
S = b.*d.^2./6 - 2*pi*12^3/32;  %in^3  (6 refers to the radius of the void)

% Moment of Inertia of Hollow Core
ybH = d./2;  %Neutral Axis, in
IgH = b.*d.^3./12 - 2*pi*6^4/4; %in^4  (6 refers to the radius of the void)

% Weight of Concrete, Normal
MUW = 150;  %lb/ft^3  (see Table 4.1 for mean and coefficient of variation equations)
VW = 0.10;
SIGW = MUW*VW;
a  = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
w1 = norminv(a1, MUW, SIGW);
x = randperm(sims);
w = zeros([sims 1]);
for j = 1:sims
   w(j,1) = w1(x(j));
end
W = w./(1000*12^3);  %kip/in^3
% Weight of section
WS = W.*Ag;  %kip/in

% Dead Load Moment
Md = WS.*(L*12)^2./8;  %kip-in

% Modular Ratio
n = 28500000./(w.^1.5.*33.*sqrt(fc*1000));

% Prestress steel eccentricity, Normal
MUe = 7.5 + 1/8;  %in  (see Table 4.1 for mean and coefficient of variation equations)
Ve = (11/32)*(7.5 + 1/8)^(-1);
SIGe = MUe*Ve;
a  = unifrnd(0,1,sims,1);
a1 = (I - 1)./sims + a./sims;
e1 = norminv(a1, MUe, SIGe);
x = randperm(sims);
e = zeros([sims 1]);
for j = 1:sims
   e(j,1) = e1(x(j));
end

% Elastic Section Modulus with respect to steel centroid
Ss = IgH./e;

% Concrete stress at cg of strands
fcsi = (-0.75.*fpu.*18.*7.*pi.*(Dps./2).^2.*(1./Ag + e./Ss) + Md./Ss)*1000;  %psi

% Elastic Shortening loss
Les = abs(n.*fcsi);  %psi

% Creep Loss
Lc = abs(n.*2.0.*fcsi);  %psi

% Shrinkage Loss
Ls = 8.2.*10.^(-6).*1.0.*28.5.*10.^6.*(1 - 0.06.*Ag./(2.*b + 2.*d))*(100 - 75);

% Corrosion Parameters
% Choose Environment Type
% 1: Low corrosion; almost no deicing salt use, not near salt water
% 2: Medium corrosion; moderate deicing salt use, near salt water without
% direct spray
Enviro = 2;

% Corrosion Initiation Time (years)
if Enviro == 1
    MUT1 = 25;
    VT1 = 0.2;
elseif Enviro == 2
    MUT1 = 10;
    VT1 = 0.6;
end
SIGT1 = MUT1*VT1;
SIGLNT1 = sqrt(log(VT1^2 + 1));
MULNT1 = log(MUT1) - 0.5*SIGLNT1^2;
a  = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
T11 = logninv(a1, MULNT1, SIGLNT1);
X = randperm(sims);
T1 = zeros([sims 1]);
for j = 1:sims
   T1(j,1) = T11(x(j));
end

% Rate of Corrosion (mil/yr)
if Enviro == 1
    MUrcorr = 0.5*.75;
    Vrcorr = 0.1;
elseif Enviro == 2
    MUrcorr = 3.0*.75;
    Vrcorr = 0.3;
end
SIGrcorr = MUrcorr*Vrcorr;
SIGLNrcorr = sqrt(log(Vrcorr^2 + 1));
MULNrcorr = log(MUrcorr) - 0.5*SIGLNrcorr^2;
a  = unifrnd(0, 1, sims, 1);
a1 = (I - 1)./sims + a./sims;
rcorr1 = logninv(a1, MULNrcorr, SIGLNrcorr);
x = randperm(sims);
rcorr = zeros([sims 1]);
for j = 1:sims
   rcorr(j,1) = rcorr1(x(j));
end

% Shear Key Degradation Choice (1, 2, 3, 4 or 5)
% 1: 1/3 load factor for lifetime (upper bound)
% 2: Deteriorate from 1/3 to 1 in 50 million cycles (linear)
% 3: Deteriorate from 1/3 to 1 in 35 million cycles (linear)
% 4: Deteriorate from 1/3 to 1 in 20 million cycles (linear)
% 5: 1 load factor for lifetime (lower bound)

LPchoice = 5;


% Writing Files
xlswrite(['Sims', num2str(Enviro), num2str(LPchoice), '.xls'],...
    [Dps b d fc Ag S Md n e Ss Les Lc Ls T1 rcorr])