1 |
function [E_UVW_F_WI,E_UVW_D_WI,I_UVW_F_WI,I_UVW_D_WI,fig3,fig4]=waveplot2(EMEP,IMLM,freqE,freqI,dirE,dirI,sysinfo) |
---|
2 |
|
---|
3 |
%make sure to load the EMEP, IMLM and wavesmon samples |
---|
4 |
%EMEP=load('emep'), IMLM=load('imlm'), spec=load('wavesmon') |
---|
5 |
%rangeE=load('EMEP_range'), rangeI=load('IMLM_range') |
---|
6 |
|
---|
7 |
%what is the magnetic variation to nearest degree |
---|
8 |
magvar=10; |
---|
9 |
|
---|
10 |
%what is the start angle |
---|
11 |
heading=sysinfo(18,:); |
---|
12 |
heading=heading/100; |
---|
13 |
sangle=heading+magvar; |
---|
14 |
|
---|
15 |
%what is the frequency and dir resolution for those generated in DIWASP |
---|
16 |
freqres=0.01; |
---|
17 |
freqresH=0.005; |
---|
18 |
freqs=[0.01:0.01:.4]; |
---|
19 |
freqsH=[0.005:0.005:0.4]; |
---|
20 |
dirres=2; |
---|
21 |
dirresH=1; |
---|
22 |
dirs=[-180:2:180]; |
---|
23 |
dirsH=[-180:1:180]; |
---|
24 |
|
---|
25 |
%set up the directions |
---|
26 |
adj_angle=90-sangle+360+180; |
---|
27 |
|
---|
28 |
% plot the spectrum generated through DIWASP |
---|
29 |
scrsz = get(0,'ScreenSize'); |
---|
30 |
fig3=figure('Position',[scrsz]); |
---|
31 |
subplot(2,3,1); |
---|
32 |
subplotspec(EMEP,4); |
---|
33 |
title('EMEP UVW'); |
---|
34 |
|
---|
35 |
subplot(2,3,2); |
---|
36 |
subplotspec(freqE,4); |
---|
37 |
title('EMEP UVW, freq = 0.005'); |
---|
38 |
|
---|
39 |
subplot(2,3,3); |
---|
40 |
subplotspec(dirE,4); |
---|
41 |
title('EMEP UVW, dir res = 1'); |
---|
42 |
|
---|
43 |
subplot(2,3,4); |
---|
44 |
subplotspec(IMLM,4); |
---|
45 |
title('IMLM UVW'); |
---|
46 |
|
---|
47 |
subplot(2,3,5); |
---|
48 |
subplotspec(freqI,4); |
---|
49 |
title('IMLM UVW, freq = 0.005'); |
---|
50 |
|
---|
51 |
subplot(2,3,6); |
---|
52 |
subplotspec(dirI,4); |
---|
53 |
title('IMLM UVW, dir res = 1'); |
---|
54 |
|
---|
55 |
%calculate just the dir energy spectrum on a single graph |
---|
56 |
|
---|
57 |
EMEPdir=sum(EMEP.S)*freqres; |
---|
58 |
IMLMdir=sum(real(IMLM.S))*freqres; |
---|
59 |
freqEdir=sum(freqE.S)*freqresH; |
---|
60 |
dirEdir=sum(dirE.S)*freqres; |
---|
61 |
freqIdir=sum(real(freqI.S))*freqresH; |
---|
62 |
dirIdir=sum(real(dirI.S))*freqres; |
---|
63 |
|
---|
64 |
%calculate just the frequency energy spectrum |
---|
65 |
EMEPfreq=sum(EMEP.S')*dirres; |
---|
66 |
IMLMfreq=sum(real(IMLM.S)')*dirres; |
---|
67 |
freqEfreq=sum(freqE.S')*dirres; |
---|
68 |
dirEfreq=sum(dirE.S')*dirresH; |
---|
69 |
freqIfreq=sum(real(freqI.S)')*dirres; |
---|
70 |
dirIfreq=sum(real(dirI.S)')*dirresH; |
---|
71 |
|
---|
72 |
%plot the directional energy spectrum |
---|
73 |
fig4=figure('Position',[scrsz]); |
---|
74 |
subplot(1,2,1); |
---|
75 |
plot(dirs,EMEPdir,'b'); |
---|
76 |
hold on |
---|
77 |
plot(dirs,freqEdir,'r'); |
---|
78 |
plot(dirsH,dirEdir,'g'); |
---|
79 |
axis(axis); |
---|
80 |
plot(dirs,IMLMdir,'c'); |
---|
81 |
plot(dirs,freqIdir,'m'); |
---|
82 |
plot(dirsH,dirIdir,'k'); |
---|
83 |
|
---|
84 |
|
---|
85 |
legend('EMEP uvw','EMEP uvw freq=0.005','EMEP uvw dir=1','IMLM uvw','IMLM uvw freq=0.005','IMLM uvw dir=1','location','best'); |
---|
86 |
title('directional wave spectrum integrated over frequency'); |
---|
87 |
xlabel('axis angle (degrees true)'); |
---|
88 |
ylabel('m^2 / hz'); |
---|
89 |
|
---|
90 |
%plot the frequency energy spectrum |
---|
91 |
subplot(1,2,2); |
---|
92 |
plot(freqs,EMEPfreq,'b'); |
---|
93 |
hold on |
---|
94 |
plot(freqsH,freqEfreq,'r'); |
---|
95 |
plot(freqs,dirEfreq,'g'); |
---|
96 |
plot(freqs,IMLMfreq,'c'); |
---|
97 |
plot(freqsH,freqIfreq,'m'); |
---|
98 |
plot(freqs,dirIfreq,'k'); |
---|
99 |
legend('EMEP uvw','EMEP uvw freq=0.005','EMEP uvw dir=1','IMLM uvw','IMLM uvw freq=0.005','IMLM uvw dir=1','location','best'); |
---|
100 |
title('directional wave spectrum integrated over direction'); |
---|
101 |
xlabel('frequency in Hz'); |
---|
102 |
ylabel('m^2 / deg'); |
---|
103 |
|
---|
104 |
% ______Calculate and display the wave parameters SigH, Tp, Dp,DTp_______ |
---|
105 |
|
---|
106 |
%For EMEP uvw |
---|
107 |
|
---|
108 |
%calculate the 0,1,2 moments |
---|
109 |
m0=sum(EMEPfreq*freqres); |
---|
110 |
m1=sum(freqs.*EMEPfreq*freqres); |
---|
111 |
m2=sum((freqs.^2).*EMEPfreq*freqres); |
---|
112 |
% Calculate the Sig wave height |
---|
113 |
Hsig=4*sqrt(m0); |
---|
114 |
|
---|
115 |
% Calculate the peak period Tp |
---|
116 |
[P,I]=max(EMEPfreq); |
---|
117 |
Tp=1/(freqs(I)); |
---|
118 |
|
---|
119 |
%Calculate the Direction of the peak period DTp |
---|
120 |
[P,I]=max(real(EMEP.S(I,:))); |
---|
121 |
DTp=dirs(I); |
---|
122 |
|
---|
123 |
%Calculate the Dominant Direction Dp |
---|
124 |
[P,I]=max(EMEPdir); |
---|
125 |
Dp=dirs(I); |
---|
126 |
|
---|
127 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
128 |
disp(['EMEP uvw']); |
---|
129 |
disp(['SigH (meters): ' num2str(Hsig)]); |
---|
130 |
disp(['peak period (seconds): ' num2str(Tp)]); |
---|
131 |
disp(['Dir of peak period: ' num2str(compangle(DTp, EMEP.xaxisdir))]); |
---|
132 |
disp(['Dominant Direction: ' num2str(compangle(Dp, EMEP.xaxisdir))]); |
---|
133 |
disp([' ']); |
---|
134 |
|
---|
135 |
% For IMLM uvw |
---|
136 |
|
---|
137 |
%calculate the 0,1,2 moments |
---|
138 |
m0=sum(IMLMfreq*freqres); |
---|
139 |
m1=sum(freqs.*IMLMfreq*freqres); |
---|
140 |
m2=sum((freqs.^2).*IMLMfreq*freqres); |
---|
141 |
% Calculate the Sig wave height |
---|
142 |
Hsig=4*sqrt(m0); |
---|
143 |
|
---|
144 |
% Calculate the peak period Tp |
---|
145 |
[P,I]=max(IMLMfreq); |
---|
146 |
Tp=1/(freqs(I)); |
---|
147 |
|
---|
148 |
%Calculate the Direction of the peak period DTp |
---|
149 |
[P,I]=max(real(IMLM.S(I,:))); |
---|
150 |
DTp=dirs(I); |
---|
151 |
|
---|
152 |
%Calculate the Dominant Direction Dp |
---|
153 |
[P,I]=max(IMLMdir); |
---|
154 |
Dp=dirs(I); |
---|
155 |
|
---|
156 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
157 |
disp(['IMLM uvw']); |
---|
158 |
disp(['SigH (meters): ' num2str(Hsig)]); |
---|
159 |
disp(['peak period (seconds): ' num2str(Tp)]); |
---|
160 |
disp(['Dir of peak period: ' num2str(compangle(DTp, IMLM.xaxisdir))]); |
---|
161 |
disp(['Dominant Direction: ' num2str(compangle(Dp, IMLM.xaxisdir))]); |
---|
162 |
disp([' ']); |
---|
163 |
|
---|
164 |
% for EMEP freq=0.005 |
---|
165 |
|
---|
166 |
%calculate the 0,1,2 moments |
---|
167 |
m0=sum(freqEfreq*freqresH); |
---|
168 |
m1=sum(freqsH.*freqEfreq*freqresH); |
---|
169 |
m2=sum((freqsH.^2).*freqEfreq*freqresH); |
---|
170 |
% Calculate the Sig wave height |
---|
171 |
E_UVW_F_Hsig=4*sqrt(m0); |
---|
172 |
|
---|
173 |
% Calculate the peak period Tp |
---|
174 |
[P,I]=max(freqEfreq); |
---|
175 |
E_UVW_F_Tp=1/(freqsH(I)); |
---|
176 |
|
---|
177 |
%Calculate the Direction of the peak period DTp |
---|
178 |
[P,I]=max(real(freqE.S(I,:))); |
---|
179 |
E_UVW_F_DTp=dirs(I); |
---|
180 |
|
---|
181 |
%Calculate the Dominant Direction Dp |
---|
182 |
[P,I]=max(freqEdir); |
---|
183 |
E_UVW_F_Dp=dirs(I); |
---|
184 |
|
---|
185 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
186 |
disp(['EMEP UVW freq=0.005']); |
---|
187 |
disp(['SigH (meters): ' num2str(E_UVW_F_Hsig)]); |
---|
188 |
disp(['peak period (seconds): ' num2str(E_UVW_F_Tp)]); |
---|
189 |
disp(['Dir of peak period: ' num2str(compangle(E_UVW_F_DTp, freqE.xaxisdir))]); |
---|
190 |
disp(['Dominant Direction: ' num2str(compangle(E_UVW_F_Dp, freqE.xaxisdir))]); |
---|
191 |
disp([' ']); |
---|
192 |
|
---|
193 |
E_UVW_F_WI.hsig=E_UVW_F_Hsig; |
---|
194 |
E_UVW_F_WI.tp=E_UVW_F_Tp; |
---|
195 |
E_UVW_F_WI.dtp=compangle(E_UVW_F_DTp, freqE.xaxisdir); |
---|
196 |
E_UVW_F_WI.dp=compangle(E_UVW_F_Dp, freqE.xaxisdir); |
---|
197 |
|
---|
198 |
%for EMEP uvw dir=1 |
---|
199 |
|
---|
200 |
%calculate the 0,1,2 moments |
---|
201 |
m0=sum(dirEfreq*freqres); |
---|
202 |
m1=sum(freqs.*dirEfreq*freqres); |
---|
203 |
m2=sum((freqs.^2).*dirEfreq*freqres); |
---|
204 |
% Calculate the Sig wave height |
---|
205 |
E_UVW_D_Hsig=4*sqrt(m0); |
---|
206 |
|
---|
207 |
% Calculate the peak period Tp |
---|
208 |
[P,I]=max(dirEfreq); |
---|
209 |
E_UVW_D_Tp=1/(freqs(I)); |
---|
210 |
|
---|
211 |
%Calculate the Direction of the peak period DTp |
---|
212 |
[P,I]=max(real(dirE.S(I,:))); |
---|
213 |
E_UVW_D_DTp=dirsH(I); |
---|
214 |
|
---|
215 |
%Calculate the Dominant Direction Dp |
---|
216 |
[P,I]=max(dirEdir); |
---|
217 |
E_UVW_D_Dp=dirsH(I); |
---|
218 |
|
---|
219 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
220 |
disp(['EMEP UVW dir=1']); |
---|
221 |
disp(['SigH (meters): ' num2str(E_UVW_D_Hsig)]); |
---|
222 |
disp(['peak period (seconds): ' num2str(E_UVW_D_Tp)]); |
---|
223 |
disp(['Dir of peak period: ' num2str(compangle(E_UVW_D_DTp, dirE.xaxisdir))]); |
---|
224 |
disp(['Dominant Direction: ' num2str(compangle(E_UVW_D_Dp, dirE.xaxisdir))]); |
---|
225 |
disp([' ']); |
---|
226 |
|
---|
227 |
E_UVW_D_WI.hsig=E_UVW_D_Hsig; |
---|
228 |
E_UVW_D_WI.tp=E_UVW_D_Tp; |
---|
229 |
E_UVW_D_WI.dtp=compangle(E_UVW_D_DTp, dirE.xaxisdir); |
---|
230 |
E_UVW_D_WI.dp=compangle(E_UVW_D_Dp, dirE.xaxisdir); |
---|
231 |
|
---|
232 |
% for IMLM freq=0.005 |
---|
233 |
|
---|
234 |
%calculate the 0,1,2 moments |
---|
235 |
m0=sum(freqIfreq*freqresH); |
---|
236 |
m1=sum(freqsH.*freqIfreq*freqresH); |
---|
237 |
m2=sum((freqsH.^2).*freqIfreq*freqresH); |
---|
238 |
% Calculate the Sig wave height |
---|
239 |
I_UVW_F_Hsig=4*sqrt(m0); |
---|
240 |
|
---|
241 |
% Calculate the peak period Tp |
---|
242 |
[P,I]=max(freqIfreq); |
---|
243 |
I_UVW_F_Tp=1/(freqsH(I)); |
---|
244 |
|
---|
245 |
%Calculate the Direction of the peak period DTp |
---|
246 |
[P,I]=max(real(freqI.S(I,:))); |
---|
247 |
I_UVW_F_DTp=dirs(I); |
---|
248 |
|
---|
249 |
%Calculate the Dominant Direction Dp |
---|
250 |
[P,I]=max(freqIdir); |
---|
251 |
I_UVW_F_Dp=dirs(I); |
---|
252 |
|
---|
253 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
254 |
disp(['IMLM UVW freq=0.005']); |
---|
255 |
disp(['SigH (meters): ' num2str(I_UVW_F_Hsig)]); |
---|
256 |
disp(['peak period (seconds): ' num2str(I_UVW_F_Tp)]); |
---|
257 |
disp(['Dir of peak period: ' num2str(compangle(I_UVW_F_DTp, freqI.xaxisdir))]); |
---|
258 |
disp(['Dominant Direction: ' num2str(compangle(I_UVW_F_Dp, freqI.xaxisdir))]); |
---|
259 |
disp([' ']); |
---|
260 |
|
---|
261 |
I_UVW_F_WI.hsig=I_UVW_F_Hsig; |
---|
262 |
I_UVW_F_WI.tp=I_UVW_F_Tp; |
---|
263 |
I_UVW_F_WI.dtp=compangle(I_UVW_F_DTp, freqI.xaxisdir); |
---|
264 |
I_UVW_F_WI.dp=compangle(I_UVW_F_Dp, freqI.xaxisdir); |
---|
265 |
|
---|
266 |
%for IMLM uvw dir=1 |
---|
267 |
|
---|
268 |
%calculate the 0,1,2 moments |
---|
269 |
m0=sum(dirIfreq*freqres); |
---|
270 |
m1=sum(freqs.*dirIfreq*freqres); |
---|
271 |
m2=sum((freqs.^2).*dirIfreq*freqres); |
---|
272 |
% Calculate the Sig wave height |
---|
273 |
I_UVW_D_Hsig=4*sqrt(m0); |
---|
274 |
|
---|
275 |
% Calculate the peak period Tp |
---|
276 |
[P,I]=max(dirIfreq); |
---|
277 |
I_UVW_D_Tp=1/(freqs(I)); |
---|
278 |
|
---|
279 |
%Calculate the Direction of the peak period DTp |
---|
280 |
[P,I]=max(real(dirI.S(I,:))); |
---|
281 |
I_UVW_D_DTp=dirsH(I); |
---|
282 |
|
---|
283 |
%Calculate the Dominant Direction Dp |
---|
284 |
[P,I]=max(dirIdir); |
---|
285 |
I_UVW_D_Dp=dirsH(I); |
---|
286 |
|
---|
287 |
%Display on the screen the SigH,Tp,Dp,DTp |
---|
288 |
disp(['IMLM UVW dir=1']); |
---|
289 |
disp(['SigH (meters): ' num2str(I_UVW_D_Hsig)]); |
---|
290 |
disp(['peak period (seconds): ' num2str(I_UVW_D_Tp)]); |
---|
291 |
disp(['Dir of peak period: ' num2str(compangle(I_UVW_D_DTp, dirI.xaxisdir))]); |
---|
292 |
disp(['Dominant Direction: ' num2str(compangle(I_UVW_D_Dp, dirI.xaxisdir))]); |
---|
293 |
disp([' ']); |
---|
294 |
|
---|
295 |
I_UVW_D_WI.hsig=I_UVW_D_Hsig; |
---|
296 |
I_UVW_D_WI.tp=I_UVW_D_Tp; |
---|
297 |
I_UVW_D_WI.dtp=compangle(I_UVW_D_DTp, freqI.xaxisdir); |
---|
298 |
I_UVW_D_WI.dp=compangle(I_UVW_D_Dp, freqI.xaxisdir); |
---|
299 |
|
---|
300 |
%function to change from axis angles to compass bearings |
---|
301 |
|
---|
302 |
function angle=compangle(angle,xaxisdir) |
---|
303 |
angle=xaxisdir*ones(size(angle))-angle; |
---|
304 |
angle=angle+360*(angle<0); |
---|
305 |
angle=angle-360*(angle>360); |
---|
306 |
|
---|
307 |
|
---|
308 |
|
---|
309 |
|
---|
310 |
|
---|