1 |
function EP=testspec(ID,theta,spread,weights,EP); |
---|
2 |
|
---|
3 |
%DIWASP V1.1 function |
---|
4 |
%testspec: Testing routine for directional wave spectrum estimation methods |
---|
5 |
% |
---|
6 |
%[EPout] = testspec(ID,theta,spread,weights,EP) |
---|
7 |
% |
---|
8 |
%Outputs: |
---|
9 |
%EPout The estimation parameters structure used in the test. |
---|
10 |
% |
---|
11 |
%Inputs: |
---|
12 |
%ID An instrument data structure containing the measured data. The ID.data field is ignored. |
---|
13 |
%theta vector with the mean directions of a sea state component |
---|
14 |
%spread vector with the spreading parameters of a sea state component |
---|
15 |
%weights vector with relative weights of sea state components |
---|
16 |
%EP The estimation parameters structure with the values under test used. Default settings are used where not specified. |
---|
17 |
% |
---|
18 |
%All inputs are required |
---|
19 |
% |
---|
20 |
%Testspec details: |
---|
21 |
% |
---|
22 |
%The fields ID.layout and ID.datatypes and ID.depth are used to specify the arrangement of the imaginary sensors. |
---|
23 |
% |
---|
24 |
%The function outputs a plot of the specified spreading function (solid line) and the estimated spreading shape (dotted line). |
---|
25 |
% |
---|
26 |
%The calculation is carried out for a frequency of 0.2 Hz. |
---|
27 |
%The inputs theta, spread and weights determine the shape of the directional spreading function. |
---|
28 |
%Each of these inputs is a vector of length n where n is the number of sea state components. |
---|
29 |
%Each sea state component has a mean direction and a spreading parameter. |
---|
30 |
%The directional spreading is calculated with a cosine power function (Mitsuyasu et al.1975 J.Phys Oceanogr.5,750-760) |
---|
31 |
% |
---|
32 |
%"help data_structures" for information on the DIWASP data structures |
---|
33 |
% |
---|
34 |
%All of the implemented calculation algorithms are as described by: |
---|
35 |
%Hashimoto,N. 1997 "Analysis of the directional wave spectrum from field data" |
---|
36 |
%In: Advances in Coastal Engineering Vol.3. Ed:Liu,P.L-F. Pub:World Scientific,Singapore |
---|
37 |
% |
---|
38 |
%Copyright (C) 2002 Coastal Oceanography Group, CWR, UWA, Perth |
---|
39 |
|
---|
40 |
S=1;F=0.2;ddir=2*pi/EP.dres;ncom=size(theta,2); |
---|
41 |
|
---|
42 |
ID=check_data(ID,1);if isempty(ID) return;end; |
---|
43 |
EP=check_data(EP,3);if isempty(EP) return;end; |
---|
44 |
|
---|
45 |
szd=size(ID.layout,2); |
---|
46 |
szdt=size(ID.datatypes,2); |
---|
47 |
|
---|
48 |
dirs=[-180:1:179]; |
---|
49 |
|
---|
50 |
disp(' ');disp('calculating.....');disp(' '); |
---|
51 |
|
---|
52 |
disp('wavenumbers') |
---|
53 |
wns=wavenumber(2*pi*F,ID.depth); |
---|
54 |
|
---|
55 |
dres=EP.dres; |
---|
56 |
pidirs=[-pi:(2*pi/dres):pi-(2*pi/dres)]; |
---|
57 |
|
---|
58 |
%generate spreading function |
---|
59 |
for ni=1:ncom |
---|
60 |
GS(ni,:)=(cos(0.5*(pidirs-theta(ni)*(pi/180)*ones(size(pidirs))))).^(2*spread(ni)); |
---|
61 |
sumGS(ni)=sum(GS(ni,:)); |
---|
62 |
end |
---|
63 |
|
---|
64 |
sumweights=sum(weights); |
---|
65 |
weights=(1/sumweights)*weights; |
---|
66 |
coeff=(1./(sumGS*ddir)).*weights; |
---|
67 |
|
---|
68 |
Gg=zeros(size(pidirs)); |
---|
69 |
for ni=1:ncom |
---|
70 |
Gg=Gg+coeff(ni)*GS(ni,:); |
---|
71 |
end |
---|
72 |
|
---|
73 |
disp('transfer parameters'); |
---|
74 |
for m=1:szd |
---|
75 |
trm(m,:,:)=feval(ID.datatypes{m},2*pi*F,pidirs,wns,ID.layout(3,m),ID.depth); |
---|
76 |
for n=1:szd |
---|
77 |
kx(m,n,:,:)=wns*((ID.layout(1,n)-ID.layout(1,m))*cos(pidirs)+(ID.layout(2,n)-ID.layout(2,m))*sin(pidirs)); |
---|
78 |
end |
---|
79 |
end |
---|
80 |
|
---|
81 |
disp('cross power spectra');disp(' '); |
---|
82 |
for m=1:szd |
---|
83 |
Ss(m,1)=S; |
---|
84 |
for n=1:szd |
---|
85 |
expx(1:dres)=exp(-i*kx(m,n,1,1:dres)); |
---|
86 |
Hh(1:dres)=trm(m,1,1:dres); |
---|
87 |
Hhs(1:dres)=conj(trm(n,1,1:dres)); |
---|
88 |
Htemp=(Hh.*Hhs.*expx); |
---|
89 |
expG=Htemp.*Gg; |
---|
90 |
xps(m,n,1)=sum(expG)*ddir; |
---|
91 |
end |
---|
92 |
end |
---|
93 |
|
---|
94 |
|
---|
95 |
% call appropriate estimation function |
---|
96 |
disp(['directional spectra using' blanks(1) EP.method ' method']); |
---|
97 |
disp(' ') |
---|
98 |
Specout=feval(EP.method,xps,trm,kx,Ss,pidirs,200,2); |
---|
99 |
|
---|
100 |
% map spectrum onto user defined direction vector |
---|
101 |
Dd=pidirs*(180/pi); |
---|
102 |
|
---|
103 |
disp('finished...plotting spectrum'); |
---|
104 |
fig=figure; |
---|
105 |
plot(Dd,real(Specout),'b.',Dd,Gg,'k'); |
---|
106 |
|
---|
107 |
nse=sum((real(Specout)-Gg).^2)./sum(Gg.^2); |
---|
108 |
|
---|
109 |
T=['Directional spreading estimated using ' blanks(1) EP.method ' method' blanks(2) 'NSE:' num2str(nse,3)]; |
---|
110 |
title(T); |
---|
111 |
ylabel(''); |
---|
112 |
xlabel('direction [degrees]'); |
---|
113 |
zlabel('m^2'); |
---|
114 |
|
---|
115 |
|
---|