NCCOOS Trac Projects: Top | Web | Platforms | Processing | Viz | Sprints | Sandbox | (Wind)

Changeset 225

Show
Ignore:
Timestamp:
02/12/09 11:31:55
Author:
gdusek
Message:

Re-organized directory structure and cleaned up some code.

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved
  • DPWavesProc/trunk/DPWavesProc/adcp_matlab/Winfo_plot.m

    r168 r225  
    1 function [fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot(Winfo); 
     1function [fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot1(Winfo); 
    22 
    33%This function takes the wave info (sigH, peak period, direction of peak and dominant direction 
     
    99fig_sigh=figure('Position',[scrsz]); 
    1010h1=plot(Winfo.time,Winfo.hsig(1,:),'b.-'); 
    11 hold on 
    12 h2=plot(Winfo.time,Winfo.hsig(2,:),'r.-'); 
    13 h3=plot(Winfo.time,Winfo.hsig(3,:),'g.-'); 
    14 h4=plot(Winfo.time,Winfo.hsig(4,:),'c.-'); 
    15 h5=plot(Winfo.time,Winfo.hsig(5,:),'m.-'); 
    16 h6=plot(Winfo.time,Winfo.hsig(7,:),'k.-'); 
    1711 
    1812datetick('x',15,'keepticks'); 
    19 legend([h1, h2, h3, h4, h5, h6],'EMEP uvw','IMLM uvw','EMEP range','IMLM range','EMEP radial','wavesmon','location','best'); 
     13legend([h1],'EMEP radial','location','best'); 
    2014title(['Significant Wave Height starting on ',datestr(Winfo.time(:,1),0)]); 
    2115xlabel('time'); 
     
    2519fig_pp=figure('Position',[scrsz]); 
    2620h1=plot(Winfo.time,Winfo.peakP(1,:),'b.-'); 
    27 hold on 
    28 h2=plot(Winfo.time,Winfo.peakP(2,:),'r.-'); 
    29 h3=plot(Winfo.time,Winfo.peakP(3,:),'g.-'); 
    30 h4=plot(Winfo.time,Winfo.peakP(4,:),'c.-'); 
    31 h5=plot(Winfo.time,Winfo.peakP(5,:),'m.-'); 
    32 h6=plot(Winfo.time,Winfo.peakP(7,:),'k.-'); 
    3321 
    3422datetick('x',15,'keepticks'); 
    35 legend([h1, h2, h3, h4, h5, h6],'EMEP uvw','IMLM uvw','EMEP range','IMLM range','EMEP radial','wavesmon','location','best'); 
     23legend([h1],'EMEP radial','location','best'); 
    3624title(['Peak Period starting on ',datestr(Winfo.time(:,1),0)]); 
    3725xlabel('time'); 
     
    4129fig_dp=figure('Position',[scrsz]); 
    4230h1=plot(Winfo.time,Winfo.dirP(1,:),'b.-'); 
    43 hold on 
    44 h2=plot(Winfo.time,Winfo.dirP(2,:),'r.-'); 
    45 h3=plot(Winfo.time,Winfo.dirP(3,:),'g.-'); 
    46 h4=plot(Winfo.time,Winfo.dirP(4,:),'c.-'); 
    47 h5=plot(Winfo.time,Winfo.dirP(5,:),'m.-'); 
    48 h6=plot(Winfo.time,Winfo.dirP(7,:),'k.-'); 
    4931 
    5032datetick('x',15,'keepticks'); 
    51 legend([h1, h2, h3, h4, h5, h6],'EMEP uvw','IMLM uvw','EMEP range','IMLM range','EMEP radial','wavesmon','location','best'); 
     33legend([h1],'EMEP radial','location','best'); 
    5234title(['Direction of Peak Period starting on ',datestr(Winfo.time(:,1),0)]); 
    5335xlabel('time'); 
     
    5840fig_dd=figure('Position',[scrsz]); 
    5941h1=plot(Winfo.time,Winfo.Ddir(1,:),'b.-'); 
    60 hold on 
    61 h2=plot(Winfo.time,Winfo.Ddir(2,:),'r.-'); 
    62 h3=plot(Winfo.time,Winfo.Ddir(3,:),'g.-'); 
    63 h4=plot(Winfo.time,Winfo.Ddir(4,:),'c.-'); 
    64 h5=plot(Winfo.time,Winfo.Ddir(5,:),'m.-'); 
    65 h6=plot(Winfo.time,Winfo.Ddir(7,:),'k.-'); 
    6642 
    6743datetick('x',15,'keepticks'); 
    68 legend([h1, h2, h3, h4, h5, h6],'EMEP uvw','IMLM uvw','EMEP range','IMLM range','EMEP radial','wavesmon','location','best'); 
     44legend([h1],'EMEP radial','location','best'); 
    6945title(['Dominant Direction starting on ',datestr(Winfo.time(:,1),0)] ); 
    7046xlabel('time'); 
  • DPWavesProc/trunk/DPWavesProc/adcp_matlab/specmultiplot.m

    r222 r225  
    6767    %Use waveplot to plot the polar, freq and dir plots of the spectra 
    6868 
    69     [E_radial_WI,fig1,fig2]=waveplottest1(radialE,sysinfo);  
     69    [E_radial_WI,fig1,fig2]=waveplot(radialE,sysinfo);  
    7070     
    7171    saveas(fig1,['polar1_' dateinfo '.fig']); 
     
    151151 
    152152%now create the time series plots of the wave info 
    153 [fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot1(Winfo); 
     153[fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot(Winfo); 
    154154 
    155155saveas(fig_sigh,['sigh_' dateinfo '.fig']); 
  • DPWavesProc/trunk/DPWavesProc/adcp_matlab/waveplot.m

    r224 r225  
    1 function [E_UVW_WI,I_UVW_WI,E_range_WI,I_range_WI,wmon_WI,fig1,fig2]=waveplot(EMEP,IMLM,rangeE,rangeI,spec,sysinfo)  
     1function [E_radial_WI,fig1,fig2]=waveplottest1(radialE,sysinfo)  
    22 
    33%make sure to load the EMEP, IMLM and wavesmon samples 
     
    88 
    99%what is the magnetic variation to nearest degree 
    10 magvar=-10; 
    11  
    12 %first change to m^2/Hz/deg 
    13 wmon.S=spec/(360*1000*1000); 
     10magvar=10; 
    1411 
    1512%what is the start angle 
     
    2623%set up the directions 
    2724adj_angle=90-sangle+360+180; 
    28 wmon.dirs=[adj_angle:-4:-(356-adj_angle)];  
    29 wmondirres=4; 
    30 wmon.xaxisdir=90; 
    31 wmon.freqs=[0.00781250:0.00781250:1]; 
    32 wmonfreqres=0.00781250; 
    3325 
    3426% plot the spectrum generated through DIWASP  
     
    3628scrsz = get(0,'ScreenSize'); 
    3729fig1=figure('Position',[scrsz]); 
    38 subplot(2,3,1); 
    39 subplotspec(EMEP,4); 
    40 title('EMEP uvw'); 
     30subplot(1,1,1); 
     31subplotspec(radialE,4); 
     32title('radial velocity data'); 
    4133 
    42 subplot(2,3,2); 
    43 subplotspec(IMLM,4); 
    44 title('IMLM uvw'); 
    45  
    46 subplot(2,3,3); 
    47 subplotspec(wmon,4); 
    48 title('Wavesmon output'); 
    49  
    50 subplot(2,3,4); 
    51 subplotspec(rangeE,4); 
    52 title('EMEP Range'); 
    53  
    54 subplot(2,3,5); 
    55 subplotspec(rangeI,4); 
    56 title('IMLM Range'); 
    5734 
    5835%calculate just the dir energy spectrum on a single graph 
    5936 
    60 EMEPdir=sum(EMEP.S)*freqres; 
    61 IMLMdir=sum(real(IMLM.S))*freqres; 
    62 EMEPrangedir=sum(rangeE.S)*freqres; 
    63 IMLMrangedir=sum(real(rangeI.S))*freqres; 
    64 wmondir=sum(wmon.S)*wmonfreqres; 
     37 
     38EMEPradialdir=sum(radialE.S)*freqres; 
    6539 
    6640%calculate just the frequency energy spectrum 
    67 EMEPfreq=sum(EMEP.S')*dirres; 
    68 IMLMfreq=sum(real(IMLM.S)')*dirres; 
    69 EMEPrangefreq=sum(rangeE.S')*dirres; 
    70 IMLMrangefreq=sum(real(rangeI.S)')*dirres; 
    71 wmonfreq=sum(wmon.S')*wmondirres; 
    72  
    73 %Compute the coefficient for the upper and lower error bounds for the power 
    74 %spectrum assuming 95% confidence.  Added on 9/17/08 
    75 degF=EMEP.degF; 
    76 chiUp=chi2inv(.975,degF); 
    77 chiLow=chi2inv(.025,degF); 
    78 coeffUp=degF/chiLow; 
    79 coeffLow=degF/chiUp; 
     41EMEPradialfreq=sum(radialE.S')*dirres; 
    8042 
    8143%Find the maximum for the directional spectrum so we can set up the proper 
    8244%x-axis 
    83 [maxvalue,maxindex] = max(EMEPdir); 
     45[maxvalue,maxindex] = max(EMEPradialdir); 
    8446maxdir=dirs(maxindex); 
    8547% set up the x-axis for all of the spectra depending on the max 
     
    8951    index2=find(dirs > -1); 
    9052    dirs(index1)=dirs(index1) +360; 
    91     %for wavesmon 
    92     Aindex=find(wmon.dirs < 0);                
    93     Bindex=find((-1 < wmon.dirs) & (wmon.dirs < 361)); 
    94     Bindex2=find(wmon.dirs > 360); 
    95     wmon.dirs(Bindex2)=wmon.dirs(Bindex2)-360; 
    96     wmon.dirs(Aindex)=wmon.dirs(Aindex)+360; 
    9753    %plot the directional energy spectrum 
    9854    fig2=figure('Position',[scrsz]); 
    9955    subplot(1,2,1); 
    100     h1 = plot(dirs(index2),EMEPdir(index2),'b'); 
     56    h1 = plot(dirs(index2),EMEPradialdir(index2),'b'); 
    10157    hold on 
    102     h1a= plot(dirs(index1),EMEPdir(index1),'b'); 
    103     h2 = plot(dirs(index2),EMEPrangedir(index2),'r'); 
    104     h2a= plot(dirs(index1),EMEPrangedir(index1),'r'); 
    105     h3 = plot(dirs(index2),IMLMrangedir(index2),'g'); 
    106     h3a= plot(dirs(index1),IMLMrangedir(index1),'g'); 
    107     %plot the wavesmon data 
    108     h4 = plot(wmon.dirs(Bindex2),wmondir(Bindex2),'k'); 
    109     h4a = plot(wmon.dirs(Bindex),wmondir(Bindex),'k'); 
    110     h4b= plot(wmon.dirs(Aindex),wmondir(Aindex),'k'); 
    111     axis(axis); 
    112     h5 = plot(dirs(index2),IMLMdir(index2),'c'); 
    113     h5a= plot(dirs(index1),IMLMdir(index1),'c'); 
    114  
     58    h1a= plot(dirs(index1),EMEPradialdir(index1),'b'); 
    11559else 
    11660    %for diwasp spectra do nothing 
    11761     
    118     %for wavesmon 
    119     Aindex=find(wmon.dirs > 180); 
    120     Bindex=find(wmon.dirs < 181); 
    121     wmon.dirs(Aindex)=wmon.dirs(Aindex)-360; 
    12262    %plot the directional energy spectrum 
    12363    fig2=figure('Position',[scrsz]); 
    12464    subplot(1,2,1); 
    125     h1 = plot(dirs,EMEPdir,'b'); 
    126     hold on 
    127     h2 = plot(dirs,EMEPrangedir,'r'); 
    128     h3 = plot(dirs,IMLMrangedir,'g'); 
    129     %plot the wavesmon data 
    130     h4 = plot(wmon.dirs(Aindex),wmondir(Aindex),'k'); 
    131     h4a = plot(wmon.dirs(Bindex),wmondir(Bindex),'k'); 
    132     axis(axis); 
    133     h5 = plot(dirs,IMLMdir,'c'); 
     65    h1 = plot(dirs,EMEPradialdir,'b'); 
    13466end 
    13567 
    136 legend([h1, h2, h3, h4, h5],'EMEP uvw','EMEP range','IMLM range','wavesmon','IMLM uvw','location','best'); 
     68legend([h1],'EMEP radial','location','best'); 
    13769title('directional wave spectrum integrated over frequency'); 
    13870xlabel('axis angle (degrees true)'); 
     
    14173%plot the frequency energy spectrum 
    14274subplot(1,2,2); 
    143 plot(freqs,EMEPfreq,'b'); 
    144 hold on 
    145 plot(freqs,EMEPrangefreq,'r'); 
    146 plot(freqs,IMLMrangefreq,'g'); 
    147 plot(wmon.freqs,wmonfreq,'k'); 
    148 axis(axis); 
    149 plot(freqs,IMLMfreq,'c'); 
    150 legend('EMEP uvw','EMEP range','IMLM range','wavesmon','IMLM uvw','location','best'); 
     75plot(freqs,EMEPradialfreq,'b'); 
     76legend('EMEP radial','location','best'); 
    15177title('directional wave spectrum integrated over direction'); 
    15278xlabel('frequency in Hz'); 
     
    15581% ______Calculate and display the wave parameters SigH, Tp, Dp,DTp_______ 
    15682 
    157 %For EMEP uvw 
     83 
     84% for Radial of EMEP 
     85 
    15886% Calculate the Sig wave height 
    159 EMEP_UVW_Hsig=Hsig(EMEP); 
     87EMEP_radial_Hsig=Hsig(radialE); 
    16088%Use the function HsigConf.m to calculate the sigH confidence limits 
    161 EMEP_HsConf=HsigConf(EMEP); 
     89EMEP_radial_HsConf=HsigConf(radialE); 
    16290 
    16391% Calculate the peak period Tp 
    164 [P,I]=max(EMEPfreq); 
    165 EMEP_UVW_Tp=1/(freqs(I)); 
     92[P,I]=max(EMEPradialfreq); 
     93EMEP_radial_Tp=1/(freqs(I)); 
    16694 
    16795%Calculate the Direction of the peak period DTp 
    168 [P,I]=max(real(EMEP.S(I,:))); 
    169 EMEP_UVW_DTp=dirs(I); 
     96[P,I]=max(real(radialE.S(I,:))); 
     97EMEP_radial_DTp=dirs(I); 
    17098 
    17199%Calculate the Dominant Direction Dp 
    172 [P,I]=max(EMEPdir); 
    173 EMEP_UVW_Dp=dirs(I); 
     100[P,I]=max(EMEPradialdir); 
     101EMEP_radial_Dp=dirs(I); 
    174102 
    175103%Display on the screen the SigH,Tp,Dp,DTp 
    176 disp(['EMEP uvw']); 
    177 disp(['SigH (meters): ' num2str(EMEP_UVW_Hsig)]); 
    178 disp(['SigH 95% confidence limits: ' num2str(EMEP_HsConf)]); 
    179 disp(['peak period (seconds): ' num2str(EMEP_UVW_Tp)]); 
    180 disp(['Dir of peak period: ' num2str(compangle(EMEP_UVW_DTp, EMEP.xaxisdir))]); 
    181 disp(['Dominant Direction: ' num2str(compangle(EMEP_UVW_Dp, EMEP.xaxisdir))]); 
     104disp(['EMEP radial']); 
     105disp(['SigH (meters): ' num2str(EMEP_radial_Hsig)]); 
     106disp(['SigH 95% confidence limits: ' num2str(EMEP_radial_HsConf)]); 
     107disp(['peak period (seconds): ' num2str(EMEP_radial_Tp)]); 
     108disp(['Dir of peak period: ' num2str(compangle(EMEP_radial_DTp, radialE.xaxisdir))]); 
     109disp(['Dominant Direction: ' num2str(compangle(EMEP_radial_Dp, radialE.xaxisdir))]); 
    182110disp(['  ']); 
    183111 
    184 E_UVW_WI.hsig=EMEP_UVW_Hsig; 
    185 E_UVW_WI.hconf=EMEP_HsConf; 
    186 E_UVW_WI.tp=EMEP_UVW_Tp; 
    187 E_UVW_WI.dtp=compangle(EMEP_UVW_DTp, EMEP.xaxisdir); 
    188 E_UVW_WI.dp=compangle(EMEP_UVW_Dp, EMEP.xaxisdir); 
     112E_radial_WI.hsig=EMEP_radial_Hsig; 
     113E_radial_WI.hconf=EMEP_radial_HsConf; 
     114E_radial_WI.tp=EMEP_radial_Tp; 
     115E_radial_WI.dtp=compangle(EMEP_radial_DTp, radialE.xaxisdir); 
     116E_radial_WI.dp=compangle(EMEP_radial_Dp, radialE.xaxisdir); 
    189117 
    190  
    191 % For IMLM uvw 
    192  
    193 % Calculate the Sig wave height 
    194 IMLM_UVW_Hsig=Hsig(IMLM); 
    195 %Use the function HsigConf.m to calculate the sigH confidence limits 
    196 IMLM_HsConf=HsigConf(IMLM); 
    197  
    198 % Calculate the peak period Tp 
    199 [P,I]=max(IMLMfreq); 
    200 IMLM_UVW_Tp=1/(freqs(I)); 
    201  
    202 %Calculate the Direction of the peak period DTp 
    203 [P,I]=max(real(IMLM.S(I,:))); 
    204 IMLM_UVW_DTp=dirs(I); 
    205  
    206 %Calculate the Dominant Direction Dp 
    207 [P,I]=max(IMLMdir); 
    208 IMLM_UVW_Dp=dirs(I); 
    209  
    210 %Display on the screen the SigH,Tp,Dp,DTp 
    211 disp(['IMLM uvw']); 
    212 disp(['SigH (meters): ' num2str(IMLM_UVW_Hsig)]); 
    213 disp(['SigH 95% confidence limits: ' num2str(IMLM_HsConf)]); 
    214 disp(['peak period (seconds): ' num2str(IMLM_UVW_Tp)]); 
    215 disp(['Dir of peak period: ' num2str(compangle(IMLM_UVW_DTp, IMLM.xaxisdir))]); 
    216 disp(['Dominant Direction: ' num2str(compangle(IMLM_UVW_Dp, IMLM.xaxisdir))]); 
    217 disp(['  ']); 
    218  
    219 I_UVW_WI.hsig=IMLM_UVW_Hsig; 
    220 I_UVW_WI.hconf=IMLM_HsConf; 
    221 I_UVW_WI.tp=IMLM_UVW_Tp; 
    222 I_UVW_WI.dtp=compangle(IMLM_UVW_DTp, IMLM.xaxisdir); 
    223 I_UVW_WI.dp=compangle(IMLM_UVW_Dp, IMLM.xaxisdir); 
    224  
    225 % for Range of EMEP 
    226  
    227 % Calculate the Sig wave height 
    228 EMEP_range_Hsig=Hsig(rangeE); 
    229 %Use the function HsigConf.m to calculate the sigH confidence limits 
    230 EMEP_range_HsConf=HsigConf(rangeE); 
    231  
    232 % Calculate the peak period Tp 
    233 [P,I]=max(EMEPrangefreq); 
    234 EMEP_range_Tp=1/(freqs(I)); 
    235  
    236 %Calculate the Direction of the peak period DTp 
    237 [P,I]=max(real(rangeE.S(I,:))); 
    238 EMEP_range_DTp=dirs(I); 
    239  
    240 %Calculate the Dominant Direction Dp 
    241 [P,I]=max(EMEPrangedir); 
    242 EMEP_range_Dp=dirs(I); 
    243  
    244 %Display on the screen the SigH,Tp,Dp,DTp 
    245 disp(['EMEP range']); 
    246 disp(['SigH (meters): ' num2str(EMEP_range_Hsig)]); 
    247 disp(['SigH 95% confidence limits: ' num2str(EMEP_range_HsConf)]); 
    248 disp(['peak period (seconds): ' num2str(EMEP_range_Tp)]); 
    249 disp(['Dir of peak period: ' num2str(compangle(EMEP_range_DTp, rangeE.xaxisdir))]); 
    250 disp(['Dominant Direction: ' num2str(compangle(EMEP_range_Dp, rangeE.xaxisdir))]); 
    251 disp(['  ']); 
    252  
    253 E_range_WI.hsig=EMEP_range_Hsig; 
    254 E_range_WI.hconf=EMEP_range_HsConf; 
    255 E_range_WI.tp=EMEP_range_Tp; 
    256 E_range_WI.dtp=compangle(EMEP_range_DTp, rangeE.xaxisdir); 
    257 E_range_WI.dp=compangle(EMEP_range_Dp, rangeE.xaxisdir); 
    258  
    259  
    260 %for Range of IMLM 
    261  
    262 % Calculate the Sig wave height 
    263 IMLM_range_Hsig=Hsig(rangeI); 
    264 %Use the function HsigConf.m to calculate the sigH confidence limits 
    265 IMLM_range_HsConf=HsigConf(rangeI); 
    266  
    267 % Calculate the peak period Tp 
    268 [P,I]=max(IMLMrangefreq); 
    269 IMLM_range_Tp=1/(freqs(I)); 
    270  
    271 %Calculate the Direction of the peak period DTp 
    272 [P,I]=max(real(rangeI.S(I,:))); 
    273 IMLM_range_DTp=dirs(I); 
    274  
    275 %Calculate the Dominant Direction Dp 
    276 [P,I]=max(IMLMrangedir); 
    277 IMLM_range_Dp=dirs(I); 
    278  
    279 %Display on the screen the SigH,Tp,Dp,DTp 
    280 disp(['IMLM range']); 
    281 disp(['SigH (meters): ' num2str(IMLM_range_Hsig)]); 
    282 disp(['SigH 95% confidence limits: ' num2str(IMLM_range_HsConf)]); 
    283 disp(['peak period (seconds): ' num2str(IMLM_range_Tp)]); 
    284 disp(['Dir of peak period: ' num2str(compangle(IMLM_range_DTp, rangeI.xaxisdir))]); 
    285 disp(['Dominant Direction: ' num2str(compangle(IMLM_range_Dp, rangeI.xaxisdir))]); 
    286 disp(['  ']); 
    287  
    288 I_range_WI.hsig=IMLM_range_Hsig; 
    289 I_range_WI.hconf=IMLM_range_HsConf; 
    290 I_range_WI.tp=IMLM_range_Tp; 
    291 I_range_WI.dtp=compangle(IMLM_range_DTp, rangeI.xaxisdir); 
    292 I_range_WI.dp=compangle(IMLM_range_Dp, rangeI.xaxisdir); 
    293  
    294 % for wavesmon 
    295  
    296 %calculate the 0,1,2 moments 
    297 m0=sum(wmonfreq*wmonfreqres);                               
    298 m1=sum(wmon.freqs.*wmonfreq*wmonfreqres); 
    299 m2=sum((wmon.freqs.^2).*wmonfreq*wmonfreqres); 
    300 % Calculate the Sig wave height 
    301 wmon_Hsig=4*sqrt(m0); 
    302  
    303 % Calculate the peak period Tp 
    304 [P,I]=max(wmonfreq); 
    305 wmon_Tp=1/(wmon.freqs(I)); 
    306  
    307 %Calculate the Direction of the peak period DTp 
    308 [P,I]=max(real(wmon.S(I,:))); 
    309 wmon_DTp=wmon.dirs(I); 
    310  
    311 %Calculate the Dominant Direction Dp 
    312 [P,I]=max(wmondir); 
    313 wmon_Dp=wmon.dirs(I); 
    314  
    315 %Display on the screen the SigH,Tp,Dp,DTp 
    316 disp(['Wavesmon output']); 
    317 disp(['SigH (meters): ' num2str(wmon_Hsig)]); 
    318 disp(['peak period (seconds): ' num2str(wmon_Tp)]); 
    319 disp(['Dir of peak period: ' num2str(compangle(wmon_DTp, wmon.xaxisdir))]); 
    320 disp(['Dominant Direction: ' num2str(compangle(wmon_Dp, wmon.xaxisdir))]); 
    321 disp(['  ']); 
    322  
    323 wmon_WI.hsig=wmon_Hsig; 
    324 wmon_WI.tp=wmon_Tp; 
    325 wmon_WI.dtp=compangle(wmon_DTp, wmon.xaxisdir); 
    326 wmon_WI.dp=compangle(wmon_Dp, wmon.xaxisdir); 
    327118 
    328119