clc
clear all
close all
myData = importdata('bathymetry.dat');
x = myData(:,1);
z = myData(:,2);
figure(1)
plot(x,z,'b')
grid on
axis tight
xlabel ('Longitude {\circ}')
ylabel ('Water depth (m)')
title('Bathymetry Profile')
load('FDS.mat')
longitude = FDS(:,1);
depth = -FDS(:,2);
figure(2);
plot(longitude,-depth);
grid on
axis tight
xlabel ('Longitude {\circ}')
ylabel ('Water depth (m)')
title('Fully-Developed Sea Profile')
gravity = 9.81;
EarthRadius = 6371000;
Tmin = 3;
Tmax = 30;
nFreq = 30;
fmin = 1/Tmax;
fmax = 1/Tmin;
f = linspace(fmin,fmax,nFreq);
Hs = 5;
Tp = 20;
gam = 3;
S = JONSWAP_spectrum(f,Hs,Tp,gam);
figure(3);
plot(f,S)
xlabel('Frequency (Hz)')
ylabel('Wave spectrum in (m^2/s)')
title('The JONSWAP Spectrum')
nWaves = numel(f);
df = f(2) - f(1);
for i = 1:nWaves
a(i,1) = sqrt(2*S(i)*df);
omega(i,1) = 2*pi*f(i);
k(i,1) = kfromw(omega(i,1), depth(1), gravity);
cg(i,1) = 0.5*omega(i,1)/k(i);
xg(i,1) = 0;
dt(i,1) = 0;
end
ndepth = size(FDS,1);
for it = 2:ndepth
for i = 1:nWaves
a(i,it) = a(i,1);
omega(i,it) = omega(i,1);
k(i,it) = kfromw(omega(i,it), depth(it),gravity);
cg(i,it) = 0.5*omega(i,1)/k(i,it);
dtheta = (longitude(it)-longitude(it-1))/180*pi;
dx = EarthRadius*dtheta;
dt(i,it) = dx/cg(i,it);
end
end
0 Comments
Sign in to comment.