1 |
function [SM,ID]= makespec(freqlph,theta,spread,weights,Ho,ID,ndat,noise) |
---|
2 |
|
---|
3 |
%DIWASP V1.1 function |
---|
4 |
%makespec: generates a directional wave spectra and pseudo data from the spectra |
---|
5 |
% |
---|
6 |
%[SM,IDout]=makespec(freqlph,theta,spread,weights,Ho,ID,ndat,noise) |
---|
7 |
% |
---|
8 |
%Outputs: |
---|
9 |
%SM A spectral matrix of the generated spectrum |
---|
10 |
%IDout Returns the input ID with data in field ID.data filled |
---|
11 |
% |
---|
12 |
%Inputs: |
---|
13 |
%freqlph 3 component vector [l p h] containing: |
---|
14 |
% the lowest frequency(l),peak, frequency(p) and highest frequency(h) |
---|
15 |
%theta vector with the mean directions of a sea state component |
---|
16 |
%spread vector with the spreading parameters of a sea state component |
---|
17 |
%weights vector with relative weights of sea state components |
---|
18 |
%Ho RMS wave height for generated spectrum |
---|
19 |
%ID Instrument data structure; field ID.data is ignored |
---|
20 |
%ndat length of simulated data |
---|
21 |
%noise level of simulated noise: |
---|
22 |
% Gaussian white noise added with variance of [noise*var(eta)] |
---|
23 |
% |
---|
24 |
% |
---|
25 |
%theta, spread and weights must all be the same length |
---|
26 |
%typical spreading parameter values are 25-100 |
---|
27 |
% |
---|
28 |
%The generated spectrum is based on an TMA spectrum (Bouws et al. 1985 JGR 90 C1,975-986) |
---|
29 |
%with directional spreading calculated with a cosine power function (Mitsuyasu et al.1975 J.Phys Oceanogr.5,750-760) |
---|
30 |
% |
---|
31 |
%"help data_structures" for information on the DIWASP data structures |
---|
32 |
|
---|
33 |
%Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth |
---|
34 |
|
---|
35 |
ID=check_data(ID,1);if isempty(ID) return;end; |
---|
36 |
|
---|
37 |
%setup default values |
---|
38 |
dt=0.5;g=9.806; |
---|
39 |
siga=0.07;sigb=0.09;gamma=2;alpha=0.014; |
---|
40 |
nf=50;nd=60; |
---|
41 |
|
---|
42 |
ncom=size(theta,2); |
---|
43 |
ns=size(ID.layout,2); |
---|
44 |
|
---|
45 |
df=(freqlph(3)-freqlph(1))/nf; |
---|
46 |
ddir=2*pi/nd; |
---|
47 |
|
---|
48 |
ffreqs=[freqlph(1):df:freqlph(3)]'; |
---|
49 |
nf=size(ffreqs,1); |
---|
50 |
fpeak=freqlph(2)*ones(size(ffreqs)); |
---|
51 |
|
---|
52 |
flh=(ffreqs<=fpeak); |
---|
53 |
sigma=siga*flh-sigb*(flh-ones(size(flh))); |
---|
54 |
|
---|
55 |
omgh=2*pi*ffreqs*(ID.depth/g)^0.5; |
---|
56 |
phi=zeros(size(omgh)); |
---|
57 |
phi=phi+0.5*(omgh<=1).*omgh.^2; |
---|
58 |
phi=phi+(omgh>=2); |
---|
59 |
phi=phi+(omgh>1 & omgh <2).*(ones(size(omgh))-0.5*(2*ones(size(omgh))-omgh).^2); |
---|
60 |
|
---|
61 |
Ek=(alpha*g*g*(2*pi)^(-4.0))*phi.*(ffreqs.^(-5)); |
---|
62 |
PhiPM=exp(-(5.0/4.0)*(ffreqs./fpeak).^(-4.0)); |
---|
63 |
PhiJ=exp(log(gamma)*exp(-((ffreqs-fpeak).^2)./((2*sigma.^2).*fpeak.^2))); |
---|
64 |
ETMA=Ek.*PhiPM.*PhiJ; |
---|
65 |
|
---|
66 |
dirs=[-pi:ddir:pi-ddir]; |
---|
67 |
omg=2*pi*ffreqs; |
---|
68 |
|
---|
69 |
for i=1:ncom |
---|
70 |
GS(i,:)=(cos(0.5*(dirs-theta(i)*(pi/180)*ones(size(dirs))))).^(2*spread(i)); |
---|
71 |
sumGS(i)=sum(GS(i,:)); |
---|
72 |
end |
---|
73 |
|
---|
74 |
sumweights=sum(weights); |
---|
75 |
weights=(1/sumweights)*weights; |
---|
76 |
coeff=(1./sumGS).*weights; |
---|
77 |
|
---|
78 |
Gg=zeros(size(dirs)); |
---|
79 |
for i=1:ncom |
---|
80 |
Gg=Gg+coeff(i)*GS(i,:); |
---|
81 |
end |
---|
82 |
|
---|
83 |
spec=(ETMA*Gg); |
---|
84 |
fac=Ho/(sqrt(8*sum(sum(spec))*df*ddir*(180/pi))); |
---|
85 |
spec=fac*fac*spec; |
---|
86 |
|
---|
87 |
SM.freqs=ffreqs;SM.dirs=(180/pi)*dirs;SM.S=spec;SM.xaxisdir=90; |
---|
88 |
|
---|
89 |
disp('plotting spectrum...press any key to make data'); |
---|
90 |
plotspec(SM,1); |
---|
91 |
drawnow; |
---|
92 |
pause |
---|
93 |
|
---|
94 |
disp('writing spectrum matrix to file'); |
---|
95 |
writespec(SM,'specmat.spec'); |
---|
96 |
|
---|
97 |
wns=wavenumber(omg,ID.depth*ones(size(omg))); |
---|
98 |
eamp=sqrt(2*df*ddir*(180/pi)*spec); |
---|
99 |
|
---|
100 |
data=makewavedata(eamp,omg,wns,dirs,ID.layout,ID.datatypes,ID.depth,ID.fs,ndat); |
---|
101 |
|
---|
102 |
for i=1:ns |
---|
103 |
gsnoise=gsamp(0,noise*var(data(:,1)),ndat); |
---|
104 |
data(:,i)=data(:,i)+gsnoise; |
---|
105 |
end |
---|
106 |
|
---|
107 |
ID.data=data; |
---|
108 |
|
---|
109 |
surfout = questdlg('Do you want to see a simulated sea surface?','DIWASP','Yes','No','No'); |
---|
110 |
|
---|
111 |
if strcmp(surfout,'Yes') |
---|
112 |
xx=[1:2:50]; |
---|
113 |
yy=[1:2:100]; |
---|
114 |
surface=makerandomsea(eamp,wns,dirs,xx,yy); |
---|
115 |
|
---|
116 |
[py,px]=meshgrid(yy,xx); |
---|
117 |
|
---|
118 |
surf(px,py,surface); |
---|
119 |
axis equal; |
---|
120 |
|
---|
121 |
end |
---|
122 |
|
---|
123 |
|
---|
124 |
|
---|
125 |
|
---|
126 |
|
---|
127 |
|
---|