function stress_calc()

%Function to calculate stress in column 11
%CEE 744 final project. Nasser M. Abbasi

%FRAME FORCES
%frame(%s)  station(%f) out(%s) case(%s) step(%s) stepNum(%f) P(%f) V2(%f) V3(%f)
%T(%f) M2(%f)  M3(%f)  S11(%f)  ptS11Max(%s) x2S11Max(%f) x3S11Max(%f)
%S11Min(%f) ptS11Min(%s)  x2S11Min(%f) x3S11Min(%f) frameElem(%s) elemStation(%f)

clear all;
close all;

TIME=6;
P=7;
M2=11;
M3=12;


cd('g:/nabbasi/data/nabbasi_web_Page/my_courses/univ_wisconson_madison/spring_2013/CEE_744_advanced_structural_dynamics/HWs/HW4');
fid = fopen('final_station_zero_forces_no_header.txt');
format_string = '%s%f%s%s%s%f%f%f%f%f%f%f%f%s%f%f%f%s%f%f%s%f';
C=textscan(fid,format_string);
figure();
N=length(C{TIME});

plot(C{TIME}(1:800),C{P}(1:800))
title('P vs. time, column 3, station 0');
xlabel('time (sec)');
ylabel('P (kips)');
grid

figure();
plot(C{TIME}(1:800),C{M2}(1:800))
title('M22 vs. time, column 3, station 0');
xlabel('time (sec)');
ylabel('M22 (kip-ft)');
grid

figure();
plot(C{TIME}(1:800),C{M3}(1:800))
title('M33 vs. time, column 3, station 0');
xlabel('time (sec)');
ylabel('M33 (kip-ft)');
grid

stress_=zeros(N,1);

for i=1:N
    m2      = abs(C{M2}(i));
    m3      = abs(C{M3}(i));
    p       = C{P}(i);
    
    stress_(i) = find_max_stress(p,m2,m3);
end

P_DEAD = -101.64;
M2_DEAD = -19.253;
M3_DEAD = 26.24;
stress_due_dead = find_max_stress(P_DEAD,abs(M2_DEAD),abs(M3_DEAD));
stress_due_dead/144


P_LIVE = -40.210;
M2_LIVE = 2.164;
M3_LIVE = 14.0489;
stress_due_live = find_max_stress(P_LIVE,abs(M2_LIVE),abs(M3_LIVE));
stress_due_live/144

stress_ = stress_ + stress_due_dead + stress_due_live;

figure();
plot(C{TIME},stress_)
title({sprintf('stress vs. time, column 3, station 0.');
    sprintf('Maximum tensile stress %6.3f kips/in^2, maximum compressive stress %6.3f kips/in^2 ',max(stress_)/144,min(stress_)/144)});

xlabel('time (sec)');
ylabel('stress slugs/ft^2');
grid

end
%-----------------------
function stress1 = find_max_stress(p,m2,m3)
A=0.8594;
C2=0.537;
C3=0.586;

stress1 = p/A + m2/C2 + m3/C3;
stress2 = p/A + m2/C2 - m3/C3;

if abs(stress2)>abs(stress1)
    stress1 = stress2;
end;

stress3 = p/A - m2/C2 + m3/C3;
if abs(stress3)>abs(stress1)
    stress1 = stress3;
end;

stress4 = p/A - m2/C2 - m3/C3;
if abs(stress4) > abs(stress1)
    stress1 = stress4;
end;

end