% TLP HydroDYNAMIC response Analysis
% Meg Brogan and katie Wasserman,  Revision 1 : 4/25/03

% Constants 
g=32.17;                                    % gravity  ft/s^2
ro=1.99;                                    % density  slugs/ft^3
E_steel=4.177e9;                            % Youngs MOdulus of Steel lb/ft^2 = 20e10 Pa
nu=35;                                      % Specific Volume ft^3/ton

% Basic Vessel Parameters
c_m=1;
T=85;                                       % Draft of Vessel ft
B_hull=245;                                 % breadth of hull ft
l=2900;
a=1;
d_cais=66.5;                                % Diameter of Caisson ft
l_cais=166;                                 % Full length of Caisson ft
d_tend=2.667;                               % diameter of tendons,  32in
t_tend=.10412;                              % wall thickness of tendons, 1.25in
l_tend=2900;                                % length of tendons ft    
n=12;                                       % number of tendons    
w_pont=35.5;                                % width of ponttons ft
h_pont=23;                                  % hight of pontoons ft
b_pont=B_hull-d_cais;                       % breadth of pontoon ft
b=(B_hull/2)-(d_cais/2);

% Frequency Range
w=(0.001:.001:2);
k=w.*w./g;

% Weights
W_top=20502;                                % total topside weight tons (revised #)
W_hull=12054;                               % hull weight tons (revised #)
W_tend=7500;                                % total weight of all 12 tendons tons
Wtot=W_top+W_hull+0.5*W_tend;               % Total weight of the vessel
Wtot_marg=1.15*Wtot                         % Total weight of the vessel w/ a 15% margin for error
m= Wtot*2200/g;                             % per MIKE T only half the tendon weight is used in these calcs; the 2200 is to get lb from tons

% Displacement, Weight and Tension Determinations

% Caissons

d_sub=T;                                    % submerged depth of caissons ft
A_cais=pi*((d_cais^2)/4);                        % waterplane area per caisson ft^2
V_cais=A_cais*d_sub;                        % submerged volume per caisson ft^3
Vt_cais=4*V_cais ;                          % total submerged volume of caissons ft^3

disp_cais=Vt_cais/nu                       % Total displacement of caisons 


% Pontoons

V_pont=b_pont*w_pont*h_pont;
Vt_pont=4*V_pont;

disp_pont=Vt_pont/nu                       % Total displacement of pontoons

% Tendons

A_tend=pi*(d_tend^2-(d_tend-2*t_tend)^2)/4; %cross sectional area of the tendons
V_tend=l_tend*A_tend;
Vt_tend=12*V_tend;

disp_tend=Vt_tend/nu                       % total displacement of tendons

% Area of the Waterplane

Awp=4*A_cais;

% TOTAL HULL DISPLACEMENT, BOUYANCY

disp_total=disp_cais+disp_pont+disp_tend    % tons

B_tot=disp_total-Wtot_marg                  % Total Bouyant force upwards disp minus weight with fudge factor  

T_tend=B_tot/n                              % Tensions in each tendon is B over # of tendons
T_tend_lb=T_tend*2200                       % Tendon Tension in pounds

%Mass determinations
m= (W_top+W_hull+0.5*W_tend)*2200/g;        % per MIKE T only half the tendon weight is used in these calcs; the 2200 is to get lb from tons

ma_c=0.5*(4/3)*pi*(d_cais/2)^3*ro*4;
ma_p=4*ro*V_pont;
ma_t=ma_c+ma_p;                               % Total added mass of the pontoon and caisson


% FORCES ON THE PLATFORM!

w_nn=sqrt((T_tend*l)/((m+ma_t)*l^2));

Fwx_o=-4*(1+c_m)*(pi/4)*ro*d_cais^2*a.*w.^2.*((1-exp(-k.*T_tend))./k).*cos(k*b); %non time dependant piece

%EOM - (m+ma)l^2*theta_2dot+T*l*theta=Fwx_o*sin(wt)*l


theta_o=(Fwx_o*l)./((-(m+ma_t)*l^2*w.^2)+T_tend*l);

H_w=theta_o./a;  
H_w;
plot(w, abs(H_w));
grid on;
ylim([0 0.003]);
xlabel('w');
ylabel('H_w');


% Transfer function of vessel, theta out from w in

% Sea State Description

zeta=5;                                     % zeta is the significant wave height- H1/3, in m
w_n=.65;                                    % Rad/sec
i=w./w_n;
%x = (w.^5); 
x1 = exp(-1.25 .* (w_n./w).^4); 

S_w= (1.25/4) * zeta^2 * (w_n./w).^4 .* (1./w) .* x1;
figure;
subplot(2,2,1);
plot(w, S_w);
grid on;
ylabel('S(w)');
xlabel('w');
subplot(2,2,2);
plot(i,S_w);
grid on;
ylabel('S(w)');
xlabel('w/w_n');

% Bretschneider Spectrum

% Integration of Vessel Transfer function with input S

S_theta=S_w.*(abs(H_w).^2);

subplot(2,2,3);
plot(w, S_theta);
grid on;
xlabel('w');
ylabel('S(theta)');
subplot(2,2,4);
plot(i,S_theta);
grid on;
xlabel('w/w_n');
ylabel('S(theta)');

rmssqr=trapz(w, S_theta);
Sw_int=trapz(w, S_w)
rms=sqrt(rmssqr);

theta_sig=4*rms
p=theta_sig*2900