1 |
function [Winfo]=specmultiplot(dateinfo1,dateinfo2,timestep); |
---|
2 |
%where dateinfo1 and 2 are the strings of the first and last date extension |
---|
3 |
%of the pressure,range,etc files and timestep is the time in hours between |
---|
4 |
%each burst |
---|
5 |
|
---|
6 |
%This function takes the pressure, range, orbit and sysinfo files created |
---|
7 |
%in python and outputs a Winfo data structure with the sig wave height, |
---|
8 |
%peak period, direction of peak period and dominant dir. It also generates |
---|
9 |
%polar plots and 2d freq and directional plots of a variety of different |
---|
10 |
%estimation methods, resolutions and the wavesmon output. |
---|
11 |
|
---|
12 |
|
---|
13 |
%first set up the time vector going from dateinfo1 to dateinfo2 |
---|
14 |
Winfo.time=[]; |
---|
15 |
|
---|
16 |
date1=datevec(dateinfo1, 'yymmddHHMM'); |
---|
17 |
date2=datevec(dateinfo2,'yymmddHHMM'); |
---|
18 |
|
---|
19 |
if date1 == date2 |
---|
20 |
Winfo.time=[datenum(date1)]; |
---|
21 |
else |
---|
22 |
Winfo.time=[datenum(date1)]; |
---|
23 |
newdate=date1; |
---|
24 |
while datenum(newdate) ~= datenum(date2) |
---|
25 |
newdate=newdate+[0 0 0 timestep 0 0]; |
---|
26 |
Winfo.time=[horzcat(Winfo.time,datenum(newdate))]; |
---|
27 |
end |
---|
28 |
end |
---|
29 |
|
---|
30 |
|
---|
31 |
%set up time, which is Winfo.time as date str in the yymmddHHMM format |
---|
32 |
time=datestr(Winfo.time, 'yymmddHHMM'); |
---|
33 |
|
---|
34 |
%set up data structure for wave info |
---|
35 |
Winfo.setup={'EMEP UVW','EMEP UVW F005','EMEP UVW D01','IMLM UVW','IMLM UVW F005','IMLM UVW D01','EMEP range','IMLM range','EMEP radial','IMLM radial','wavesmon'}'; |
---|
36 |
Winfo.hsig=[]; |
---|
37 |
Winfo.peakP=[]; |
---|
38 |
Winfo.dirP=[]; |
---|
39 |
Winfo.Ddir=[]; |
---|
40 |
|
---|
41 |
%Load the data and run the script |
---|
42 |
for i=1:length(time(:,1)) |
---|
43 |
dateinfo=time(i,:); |
---|
44 |
|
---|
45 |
pressure=strcat('pressure_',dateinfo,'.txt'); |
---|
46 |
range=strcat('range_',dateinfo,'.txt'); |
---|
47 |
orbit=strcat('orbit_',dateinfo,'.txt'); |
---|
48 |
sysinfo=strcat('sysinfo_',dateinfo,'.txt'); |
---|
49 |
spec=strcat('DSpec',dateinfo,'.txt'); |
---|
50 |
|
---|
51 |
pressure=load(pressure); |
---|
52 |
range=load(range); |
---|
53 |
orbit=load(orbit); |
---|
54 |
sysinfo=load(sysinfo); |
---|
55 |
spec=load(spec); |
---|
56 |
|
---|
57 |
%set up data with uvw and pressure, freq at 0.01 and dir at 2 |
---|
58 |
[ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,1); |
---|
59 |
|
---|
60 |
%run diwasp to generate this spectrum |
---|
61 |
[E_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
62 |
|
---|
63 |
%now for IMLM |
---|
64 |
EP.method='IMLM'; |
---|
65 |
EP.iter=3; |
---|
66 |
[I_uvw_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
67 |
|
---|
68 |
%now for IMLM with increased Freq resolution |
---|
69 |
SM.freqs=[0.005:0.005:0.4]; |
---|
70 |
[I_uvw_F005_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
71 |
|
---|
72 |
%now for IMLM with increased Dir resolution |
---|
73 |
SM.freqs=[0.01:0.01:0.4]; |
---|
74 |
SM.dirs=[-180:1:180]; |
---|
75 |
[I_uvw_F01_D1, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
76 |
|
---|
77 |
%now for EMEP with increased Freq res and then Dir res |
---|
78 |
EP.method='EMEP'; |
---|
79 |
EP.iter=100; |
---|
80 |
SM.freqs=[0.005:0.005:0.4]; |
---|
81 |
SM.dirs=[-180:2:180]; |
---|
82 |
[E_uvw_F005_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
83 |
|
---|
84 |
SM.freqs=[0.01:0.01:0.4]; |
---|
85 |
SM.dirs=[-180:1:180]; |
---|
86 |
[E_uvw_F01_D1, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
87 |
|
---|
88 |
|
---|
89 |
%set up data with ranges, freq and dir at default |
---|
90 |
[ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,2); |
---|
91 |
|
---|
92 |
% run diwasp to generate this spectrum with EMEP |
---|
93 |
[E_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
94 |
|
---|
95 |
%now with IMLM |
---|
96 |
EP.method='IMLM'; |
---|
97 |
EP.iter=3; |
---|
98 |
[I_range_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
99 |
|
---|
100 |
%set up data with radial velocities, freq and dir at default |
---|
101 |
[ID,SM,EP]=radialtouvw(pressure,range,orbit,sysinfo,3); |
---|
102 |
|
---|
103 |
% run diwasp to generate this spectrum with EMEP |
---|
104 |
[E_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
105 |
|
---|
106 |
%now with IMLM |
---|
107 |
EP.method='IMLM'; |
---|
108 |
EP.iter=3; |
---|
109 |
[I_radial_F01_D2, EPout]=dirspectest(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); |
---|
110 |
|
---|
111 |
%Use waveplot to plot the polar, freq and dir plots of the spectra |
---|
112 |
|
---|
113 |
[E_UVW_WI,I_UVW_WI,E_range_WI,I_range_WI,wmon_WI,fig1,fig2]=waveplot(E_uvw_F01_D2,I_uvw_F01_D2,E_range_F01_D2,I_range_F01_D2,spec,sysinfo); |
---|
114 |
|
---|
115 |
saveas(fig1,['polar1_' dateinfo '.fig']); |
---|
116 |
saveas(fig2,['dirfreq1_' dateinfo '.fig']); |
---|
117 |
|
---|
118 |
[E_UVW_F_WI,E_UVW_D_WI,I_UVW_F_WI,I_UVW_D_WI,fig3,fig4]=waveplot2(E_uvw_F01_D2,I_uvw_F01_D2,E_uvw_F005_D2,I_uvw_F005_D2,E_uvw_F01_D1,I_uvw_F01_D1,sysinfo); |
---|
119 |
|
---|
120 |
saveas(fig3,['polar2_' dateinfo '.fig']); |
---|
121 |
saveas(fig4,['dirfreq2_' dateinfo '.fig']); |
---|
122 |
|
---|
123 |
[E_radial_WI,I_radial_WI,fig5,fig6]=radialwaveplot(E_radial_F01_D2,I_radial_F01_D2,spec,sysinfo); |
---|
124 |
|
---|
125 |
saveas(fig5,['polar3_' dateinfo '.fig']); |
---|
126 |
saveas(fig6,['dirfreq3_' dateinfo '.fig']); |
---|
127 |
|
---|
128 |
close all; |
---|
129 |
|
---|
130 |
%Fill up the Winfo data structure |
---|
131 |
|
---|
132 |
hsig=[vertcat(E_UVW_WI.hsig,E_UVW_F_WI.hsig,E_UVW_D_WI.hsig,I_UVW_WI.hsig,I_UVW_F_WI.hsig,I_UVW_D_WI.hsig,E_range_WI.hsig,I_range_WI.hsig,E_radial_WI.hsig,I_radial_WI.hsig,wmon_WI.hsig)]; |
---|
133 |
peakP=[vertcat(E_UVW_WI.tp,E_UVW_F_WI.tp,E_UVW_D_WI.tp,I_UVW_WI.tp,I_UVW_F_WI.tp,I_UVW_D_WI.tp,E_range_WI.tp,I_range_WI.tp,E_radial_WI.tp,I_radial_WI.tp,wmon_WI.tp)]; |
---|
134 |
dirP=[vertcat(E_UVW_WI.dtp,E_UVW_F_WI.dtp,E_UVW_D_WI.dtp,I_UVW_WI.dtp,I_UVW_F_WI.dtp,I_UVW_D_WI.dtp,E_range_WI.dtp,I_range_WI.dtp,E_radial_WI.dtp,I_radial_WI.dtp,wmon_WI.dtp)]; |
---|
135 |
Ddir=[vertcat(E_UVW_WI.dp,E_UVW_F_WI.dp,E_UVW_D_WI.dp,I_UVW_WI.dp,I_UVW_F_WI.dp,I_UVW_D_WI.dp,E_range_WI.dp,I_range_WI.dp,E_radial_WI.dp,I_radial_WI.dp,wmon_WI.dp)]; |
---|
136 |
|
---|
137 |
Winfo.hsig=[horzcat(Winfo.hsig,hsig)]; |
---|
138 |
Winfo.peakP=[horzcat(Winfo.peakP,peakP)]; |
---|
139 |
Winfo.dirP=[horzcat(Winfo.dirP,dirP)]; |
---|
140 |
Winfo.Ddir=[horzcat(Winfo.Ddir,Ddir)]; |
---|
141 |
|
---|
142 |
end |
---|
143 |
|
---|
144 |
% change Winfo.time back to a matlab datenum |
---|
145 |
|
---|
146 |
%Winfo.time=datenum(Winfo.time); |
---|
147 |
|
---|
148 |
|
---|