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

Changeset 334

Show
Ignore:
Timestamp:
06/14/10 12:19:31
Author:
gdusek
Message:

Total update to the DPWP code. Significant changes made to direspec, specmultiplot, radialtouvw and IMLM code. 6/14/10

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved
  • DPWP/trunk/DPWP/ADCP_splitter/adcpwaves.py

    r315 r334  
    44 
    55NOTE:  Transducer height above bottom is hardwired into this code at .4m. 
    6 Need to update code in the future to be a user selected option""" 
     6Need to update code in the future to be a user selected option.  It is also 
     7hardwired into the Matlab code radialtouvw.m and radial.m""" 
    78 
    89 
  • DPWP/trunk/DPWP/adcp_matlab/radialtouvw.m

    r314 r334  
    1818%orbit=load('orbital data'),sysinfo=load('sysinfo data') and range 
    1919 
    20 %what is the transducer face height off the bottom? 
     20%what is the transducer face height off the bottom 
     21%NOTE: also hardwired into radial.m 
    2122adcpheight=0.4; 
    2223 
     
    251252ID.fs=2; 
    252253%depth 
    253 ID.depth=(0.4+mean(press)); %note changed for work with the FRF ADCP********************** 
     254ID.depth=avgdepth;  
    254255%the spectral matrix structure 
    255256SM.freqs=[0.01:0.01:0.4]; 
    256 SM.dirs=[0:2:360]; 
     257SM.dirs=[0:2:360]; %note that SM.dirs is hardwired into specmultiplot (line 123) 
    257258SM.xaxisdir= 90; 
    258259%the estimation parameter 
    259 EP.method= 'EMEP'; 
    260 EP.iter=100; 
     260EP.method= 'IMLM'; 
     261EP.iter=50; 
     262EP.nfft=256; 
    261263 
    262264 
    263265if data_type == 1 
    264     %For uvw and pressure for bins 2,3,4 
     266    %For uvw and pressure for the highest three bins (3,4,5) 
     267     
    265268    % the datatypes 
    266269    ID.datatypes={'pres' 'velx' 'vely' 'velz' 'velx' 'vely' 'velz' 'velx' 'vely' 'velz'}; 
     
    286289    % the datatypes 
    287290    ID.datatypes={'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial' 'radial'}; 
    288     % the layout using bins 2,3,4 
    289     %ID.layout = [xyzpositions(1,1,2) xyzpositions(1,2,2) xyzpositions(1,3,2) xyzpositions(1,4,2) xyzpositions(1,1,3) xyzpositions(1,2,3) xyzpositions(1,3,3) xyzpositions(1,4,3) xyzpositions(1,1,4) xyzpositions(1,2,4) xyzpositions(1,3,4) xyzpositions(1,4,4); 
    290     %xyzpositions(2,1,2) xyzpositions(2,2,2) xyzpositions(2,3,2) xyzpositions(2,4,2) xyzpositions(2,1,3) xyzpositions(2,2,3) xyzpositions(2,3,3) xyzpositions(2,4,3) xyzpositions(2,1,4) xyzpositions(2,2,4) xyzpositions(2,3,4) xyzpositions(2,4,4);  
    291     %xyzpositions(3,1,2) xyzpositions(3,2,2) xyzpositions(3,3,2) xyzpositions(3,4,2) xyzpositions(3,1,3) xyzpositions(3,2,3) xyzpositions(3,3,3) xyzpositions(3,4,3) xyzpositions(3,1,4) xyzpositions(3,2,4) xyzpositions(3,3,4) xyzpositions(3,4,4)]; 
     291 
     292    % the layout using highest 3 bins depending on the number of bins 
     293    % recorded 
     294     
     295    %highest bin 
     296    bin3=5; 
     297    %second highest 
     298    bin2=4; 
     299    %third highest 
     300    bin1=3; 
     301     
     302    ID.layout = [xyzpositions(1,1,bin3) xyzpositions(1,2,bin3) xyzpositions(1,3,bin3) xyzpositions(1,4,bin3) xyzpositions(1,1,bin2) xyzpositions(1,2,bin2) xyzpositions(1,3,bin2) xyzpositions(1,4,bin2) xyzpositions(1,1,bin1) xyzpositions(1,2,bin1) xyzpositions(1,3,bin1) xyzpositions(1,4,bin1); 
     303    xyzpositions(2,1,bin3) xyzpositions(2,2,bin3) xyzpositions(2,3,bin3) xyzpositions(2,4,bin3) xyzpositions(2,1,bin2) xyzpositions(2,2,bin2) xyzpositions(2,3,bin2) xyzpositions(2,4,bin2) xyzpositions(2,1,bin1) xyzpositions(2,2,bin1) xyzpositions(2,3,bin1) xyzpositions(2,4,bin1);  
     304    xyzpositions(3,1,bin3) xyzpositions(3,2,bin3) xyzpositions(3,3,bin3) xyzpositions(3,4,bin3) xyzpositions(3,1,bin2) xyzpositions(3,2,bin2) xyzpositions(3,3,bin2) xyzpositions(3,4,bin2) xyzpositions(3,1,bin1) xyzpositions(3,2,bin1) xyzpositions(3,3,bin1) xyzpositions(3,4,bin1)]; 
    292305    % the data 
    293     %ID.data = horzcat(orbitnew(:,5),orbitnew(:,6),orbitnew(:,7),orbitnew(:,8),orbitnew(:,9),orbitnew(:,10),orbitnew(:,11),orbitnew(:,12),orbitnew(:,13),orbitnew(:,14),orbitnew(:,15),orbitnew(:,16)); 
    294  
    295     % the layout using bins 3,4,5 
    296     ID.layout = [xyzpositions(1,1,5) xyzpositions(1,2,5) xyzpositions(1,3,5) xyzpositions(1,4,5) xyzpositions(1,1,4) xyzpositions(1,2,4) xyzpositions(1,3,4) xyzpositions(1,4,4) xyzpositions(1,1,3) xyzpositions(1,2,3) xyzpositions(1,3,3) xyzpositions(1,4,3); 
    297     xyzpositions(2,1,5) xyzpositions(2,2,5) xyzpositions(2,3,5) xyzpositions(2,4,5) xyzpositions(2,1,4) xyzpositions(2,2,4) xyzpositions(2,3,4) xyzpositions(2,4,4) xyzpositions(2,1,3) xyzpositions(2,2,3) xyzpositions(2,3,3) xyzpositions(2,4,3);  
    298     xyzpositions(3,1,5) xyzpositions(3,2,5) xyzpositions(3,3,5) xyzpositions(3,4,5) xyzpositions(3,1,4) xyzpositions(3,2,4) xyzpositions(3,3,4) xyzpositions(3,4,4) xyzpositions(3,1,3) xyzpositions(3,2,3) xyzpositions(3,3,3) xyzpositions(3,4,3)]; 
    299     % the data 
    300     ID.data = horzcat(orbitnew(:,17),orbitnew(:,18),orbitnew(:,19),orbitnew(:,20),orbitnew(:,13),orbitnew(:,14),orbitnew(:,15),orbitnew(:,16),orbitnew(:,9),orbitnew(:,10),orbitnew(:,11),orbitnew(:,12)); 
     306    %set up values based on what bins are available 
     307    orb=orbOut; 
     308    ID.data = horzcat(orbitnew(:,orb-3),orbitnew(:,orb-2),orbitnew(:,orb-1),orbitnew(:,orb),orbitnew(:,orb-7),orbitnew(:,orb-6),orbitnew(:,orb-5),orbitnew(:,orb-4),orbitnew(:,orb-11),orbitnew(:,orb-10),orbitnew(:,orb-9),orbitnew(:,orb-8)); 
    301309 
    302310     
  • DPWP/trunk/DPWP/adcp_matlab/specmultiplot.m

    r314 r334  
    1 function [Winfo,spec2d,fd,thetad,td]=specmultiplot(dateinfo1,dateinfo2,timestep);  
     1function [Winfo,spec2d,fd,thetad,td]=specmultiplot(dateinfo1,dateinfo2,timestep) 
    22%where dateinfo1 and 2 are the strings of the first and last date extension 
    33%of the pressure,range,etc files and timestep is the time in hours between 
     
    88%peak period, direction of peak period and dominant dir.  It also generates 
    99%polar plots and 2d freq and directional plots of a variety of different 
    10 %estimation methods, resolutions and the wavesmon output.  
     10%estimation methods.  
     11 
     12%%% IMPT NOTE!!! The ADCP height above bottom is manually set in the 
     13%%% radialtouvw.m and the radial.m code.  In this case the default is set 
     14%%% at .4m  This will be changed in a future release!! 
    1115 
    1216 
     
    3539 
    3640%set up data structure for wave info 
    37 Winfo.setup={'EMEP radial'}; 
     41Winfo.setup={'IMLM radial'}; 
    3842Winfo.hsig=[]; 
    3943Winfo.hconf=[]; 
     
    6468     
    6569    % run diwasp to generate this spectrum with EMEP 
    66     [radialE, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); 
     70    [radialE,~]=dirspec(ID,SM,EP,{'MESSAGE',0,'PLOTTYPE',0}); 
    6771 
    6872    %Use waveplot to plot the polar, freq and dir plots of the spectra 
    6973 
    70     [E_radial_WI,fig1,fig2]=waveplot(radialE,sysinfo);  
     74    [E_radial_WI]=waveplot(radialE);  
    7175     
    72     saveas(fig1,['polar1_' dateinfo '.fig']); 
    73     saveas(fig2,['dirfreq1_' dateinfo '.fig']); 
    74      
    75     movefile(['polar1_' dateinfo '.fig'],'tempdir'); 
    76     movefile(['dirfreq1_' dateinfo '.fig'],'tempdir'); 
    77     close all; 
    78      
    79      
     76    %rename the files 
     77    movefile('tempdir/fig1.fig',['tempdir/polar1_' dateinfo '.fig']); 
     78    movefile('tempdir/fig2.fig',['tempdir/dirfreq1_' dateinfo '.fig']); 
     79   
     80      
    8081    %Fill up the Winfo data structure 
    8182 
     
    9394    Winfo.Ddir=[horzcat(Winfo.Ddir,Ddir)]; 
    9495    Winfo.Spectrum.EMEPradial=[cat(3,Winfo.Spectrum.EMEPradial,radialE.S)]; 
     96     
     97    %save a continuing backup of Winfo in case of an error 
     98    save('backup_Winfo','Winfo'); 
     99    movefile('backup_Winfo.mat','tempdir'); 
    95100end 
    96101 
     
    159164%************************************************************************** 
    160165 
     166%first delete the backup_Winfo file 
     167%delete('tempdir/backup_Winfo.mat'); 
     168 
    161169%now create the time series plots of the wave info 
    162170[fig_sigh,fig_pp,fig_dp,fig_dd]=Winfo_plot(Winfo); 
     
    167175saveas(fig_dd,['domdir_' dateinfo '.fig']); 
    168176 
     177close all; 
     178 
    169179movefile(['sigh_' dateinfo '.fig'],'tempdir'); 
    170180movefile(['peakp_' dateinfo '.fig'],'tempdir'); 
     
    174184movefile('tempdir',strcat('procdata_',dateinfo)); 
    175185 
    176 close all; 
    177186 
    178187 
  • DPWP/trunk/DPWP/adcp_matlab/waveplot.m

    r231 r334  
    1 function [E_radial_WI,fig1,fig2]=waveplot(radialE,sysinfo)  
     1function [E_radial_WI]=waveplot(radialE)  
    22 
    33%make sure to load the EMEP, IMLM and wavesmon samples 
     
    1717scrsz = get(0,'ScreenSize'); 
    1818fig1=figure('Position',[scrsz]); 
    19 subplot(1,1,1); 
    2019subplotspec(radialE,4); 
    2120title('radial velocity data'); 
     21 
     22saveas(fig1,'tempdir/fig1.fig'); 
     23close all force; 
    2224 
    2325 
     
    7173end 
    7274 
    73 legend([h1],'EMEP radial','location','best'); 
     75legend([h1],'IMLM radial','location','best'); 
    7476title('directional wave spectrum integrated over frequency'); 
    7577xlabel('axis angle (degrees true)'); 
     
    8284h2=plot(freqs,EMEPradialfreqUP,'b--'); 
    8385h2a=plot(freqs,EMEPradialfreqLOW,'b--'); 
    84 legend([h1,h2],'EMEP radial','95% confidence limits','location','best'); 
     86legend([h1,h2],'IMLM radial','95% confidence limits','location','best'); 
    8587title('directional wave spectrum integrated over direction'); 
    8688xlabel('frequency in Hz'); 
    8789ylabel('m^2 / hz'); 
     90 
     91saveas(fig2,'tempdir/fig2.fig'); 
     92close all force; 
    8893 
    8994% ______Calculate and display the wave parameters SigH, Tp, Dp,DTp_______ 
     
    110115 
    111116%Display on the screen the SigH,Tp,Dp,DTp 
    112 disp(['EMEP radial']); 
     117disp(['IMLM radial']); 
    113118disp(['SigH (meters): ' num2str(EMEP_radial_Hsig)]); 
    114119disp(['SigH 95% confidence limits: ' num2str(EMEP_radial_HsConf)]); 
  • DPWP/trunk/DPWP/diwasp_1_1GD/dirspec.m

    r314 r334  
    159159%have Ss = the spectral average if the radial velocities are used without a 
    160160%pressure or vertical range datatype first 
     161lags=50; 
    161162compare=strcmp('radial',ID.datatypes(1)); 
    162163if compare == 1 
     
    165166    %the total number of sensors or szd 
    166167    dataUsed=szd; 
    167 end 
    168  
     168     
     169    %also need to now calculate the effective number of windows to 
     170    %accurately calculate the degrees of freedom below. 
     171    Tau=ones(dataUsed); %Tau Matrix 
     172    Nstar=Tau;  %Nstar Matrix 
     173    for i=1:dataUsed 
     174        for k=1:dataUsed 
     175            covi=xcov(data(:,i),lags,'coeff'); 
     176            covk=xcov(data(:,k),lags,'coeff'); 
     177            covik=xcov(data(:,i),data(:,k),lags,'coeff'); 
     178            covki=xcov(data(:,k),data(:,i),lags,'coeff'); 
     179            covicovk=covi.*covk; 
     180            covikcovki=covik.*covki; 
     181            %sign1=sign(covicovk); 
     182            %sign2=sign(covikcovki); 
     183            %part1=(abs(covicovk)).^.5; 
     184            %part2=(abs(covikcovki)).^.5; 
     185            Tau(i,k)=sum(covicovk+covikcovki); 
     186            Nstar=dataLength./Tau; 
     187        end 
     188    end 
     189    %calculate the average value of Nstar for the off-diagnal values 
     190    totalN=sum(sum(Nstar)); 
     191    for i=1:dataUsed 
     192        out(i)=Nstar(i,i); 
     193    end 
     194    autoOut=sum(out); %don't count the auto values 
     195    totalN=totalN-autoOut; 
     196    avgN=totalN/((dataUsed^2)-dataUsed);%find the average N value for all co values 
     197    Nstar=avgN;     
     198else 
     199    %if there is only ONE time series used for the frequency spec 
     200    covi=xcov(data(:,1),lags,'coeff'); 
     201    Tau=sum(covi.*covi); 
     202    Nstar=dataLength./Tau; 
     203end 
     204 
     205%now find the effective number windows for each time series given Nstar 
     206numWindows=Nstar/nfft; 
    169207 
    170208 
     
    172210%See Emery and Thomson, 2004, added on 9/17/08 
    173211windowConstant=2.5164; %for Hamming Window (Matlab default for cpsd) 
    174 degF=dataUsed*2*numWindows*windowConstant; 
     212degF=2*numWindows*dataUsed*windowConstant; 
    175213SM.degF=degF; 
    176214 
     
    180218    % call appropriate estimation function 
    181219disp(['directional spectra using' blanks(1) EP.method ' method']);disp(' ') 
    182 Specout=feval(EP.method,xps(:,:,1:ffs),trm(:,1:ffs,:),kx(:,:,1:ffs,:),Ss(:,1:ffs),pidirs,EP.iter,displ); 
     220 
     221[Specout]=feval(EP.method,xps(:,:,1:ffs),trm(:,1:ffs,:),kx(:,:,1:ffs,:),Ss(:,1:ffs),pidirs,EP.iter,displ); 
    183222 
    184223Specout(find(isnan(Specout)))=0.0; 
     
    196235 
    197236%check Hsig of mapped spectrum and check sufficiently close to original 
    198 Hs2=Hsig(S,SM.freqs(2)-SM.freqs(1),SM.dirs(2)-SM.dirs(1)); 
    199 if (Hs2-Hs)/Hs >0.01 
    200     warning('User defined grid may be too coarse; try increasing resolution of ''SM.freqs'' or ''SM.dirs'''); 
    201 end 
     237%*****As far as I can tell this warning is useless**** Edited out on 3/2/10 
     238 
     239%Hs2=Hsig(S,SM.freqs(2)-SM.freqs(1),SM.dirs(2)-SM.dirs(1)); 
     240%if (Hs2-Hs)/Hs >0.01 
     241%    warning('User defined grid may be too coarse; try increasing resolution of ''SM.freqs'' or ''SM.dirs'''); 
     242%end 
    202243 
    203244%smooth spectrum 
  • DPWP/trunk/DPWP/diwasp_1_1GD/private/IMLM.m

    r168 r334  
    11function [S]=IMLM(xps,trm,kx,Ss,pidirs,miter,displ) 
    22 
    3 gamma=0.1; 
     3%this code was edited significantly on 5/2010 --- see comments for details 
     4 
     5gamma=.1; 
    46beta=1.0; 
    57alpha=0.1; 
     
    3234   end 
    3335    
    34     
    35    invcps=inv(xps(:,:,ff)); 
     36   I=eye(m)*sqrt(eps); 
     37   invcps=pinv(xps(:,:,ff)+I); %changed on 5/5/10 to pinv and +I 
    3638   Sftmp=zeros(ddirs,1); 
    3739   for m=1:szd 
     
    4143                end 
    4244   end 
    43    Eo=(1./Sftmp(:)); 
    44    kappa=1./(ddir*sum(Eo)); 
    45    Eo=kappa*Eo; 
     45   Eo=(1./real(Sftmp(:))); %added 5/5/10 to take just the real part 
     46   [Eo]=normalize(Eo,ddir); %added 5/11/10 to remove negative values before normalizing 
    4647   E=Eo; 
    4748   T=Eo; 
     
    5455                end 
    5556      end 
    56       invcps=inv(ixps); 
     57      invcps=pinv(ixps+I); % changed on 5/5/10 to pinv and +I  
    5758      Sftmp=zeros(ddirs,1); 
    5859      for m=1:szd 
     
    6364        end 
    6465        Told=T; 
    65       T=(1./Sftmp(:)); 
    66       kappa=1./(ddir*sum(T)); 
    67       T=kappa*T; 
     66      T=(1./real(Sftmp(:))); %added on 5/5/10 to take just real part 
     67      [T]=normalize(T,ddir); %added 5/11/10 to remove negative values before normalizing  
    6868       
    6969      %lambda=ones(size(T))-(T./Eo) 
    7070      %ei=gamma*lambda.*E; 
    7171      ei=gamma*((Eo-T)+alpha*(T-Told)); 
     72      %lambda=Eo-T;      %edited on 5/5/10 to calculate ei slightly differently  
     73      %ei=(abs(lambda).^(beta+1))./(gamma.*lambda); 
    7274      E=E+ei; 
    73       kappa=1./(ddir*sum(E)); 
    74         E=kappa*E; 
     75      [E]=normalize(E,ddir);%added 5/11/10 to remove negative values before normalizing 
    7576 
    7677             
  • DPWP/trunk/DPWP/diwasp_1_1GD/private/radial.m

    r231 r334  
    1 function [trm]=radial(ffreqs,dirs,wns,z,depth,xpos,ypos); 
     1function [trm]=radial(ffreqs,dirs,wns,z,depth,xpos,ypos) 
    22 
    33%This transfer function takes the along beam radial velocities from an ADCP 
     
    2121 
    2222% compute the axis angle 
    23 if xpos > 0 & ypos > 0 % 1st quadrant 
     23if xpos > 0 && ypos > 0 % 1st quadrant 
    2424    B= acos(xpos/dist); 
    25 elseif xpos < 0 & ypos < 0 % 3 quadrant 
     25elseif xpos < 0 && ypos < 0 % 3 quadrant 
    2626    B= pi + acos(abs(xpos/dist)); 
    27 elseif xpos < 0 & ypos > 0 % 2nd quadrant  
     27elseif xpos < 0 && ypos > 0 % 2nd quadrant  
    2828    B= pi/2 + asin(abs(xpos/dist)); 
    29 elseif xpos > 0 & ypos < 0 % 4th quadrant 
     29elseif xpos > 0 && ypos < 0 % 4th quadrant 
    3030    B = 3*pi/2 + asin(abs(xpos/dist)); 
    3131end 
     
    6666 
    6767%now put them together to create the transfer matrix 
    68 trm=(FLKZ*(sin(a)*cos(dirs-B))-(i*FLKZi)*A1); 
     68trm=(FLKZ*(sin(a)*cos(dirs-B))-(1i*FLKZi)*A1); 
    6969 
    7070 
  • DPWP/trunk/DPWP/diwasp_1_1GD/private/velx.m

    r168 r334  
    11function trm=velx(ffreqs,dirs,wns,z,depth,xpos,ypos) 
    22 
    3 Kz=cosh(z*wns)./sinh(depth*wns); 
    4 %include a maximum cuttoff for the velocity response function 
    5 Kz(find(Kz<0.1))=0.1; 
    6 Kz(find(isnan(Kz)))=1; 
    7 trm=(ffreqs.*Kz)*cos(dirs); 
     3 
     4%Edits were made on 4/20/10 to fix a bug in the cutoff limits in the KZ 
     5%parts of the transfer function.  
     6 
     7%set the cutoff value (for high frequencies) 
     8cutoff=0.1; 
     9 
     10 
     11KZ=cosh(z*wns)./sinh(depth*wns); 
     12FLKZ=ffreqs.*KZ; 
     13 
     14%find the cutoff and replace the values below it 
     15cutindex=find(FLKZ < cutoff); 
     16FLKZ(cutindex)=cutoff; 
     17 
     18trm=(FLKZ)*cos(dirs); 
  • DPWP/trunk/DPWP/diwasp_1_1GD/private/vely.m

    r168 r334  
    11function trm=vely(ffreqs,dirs,wns,z,depth,xpos,ypos) 
    22 
    3 Kz=cosh(z*wns)./sinh(depth*wns); 
    4 %include a maximum cuttoff for the velocity response function 
    5 Kz(find(Kz<0.1))=0.1; 
    6 Kz(find(isnan(Kz)))=1; 
    7 trm=(ffreqs.*Kz)*sin(dirs); 
     3%Edits were made on 4/20/10 to fix a bug in the cutoff limits in the KZ 
     4%parts of the transfer function.  
     5 
     6%set the cutoff value (for high frequencies) 
     7cutoff=0.1; 
     8 
     9 
     10KZ=cosh(z*wns)./sinh(depth*wns); 
     11FLKZ=ffreqs.*KZ; 
     12 
     13%find the cutoff and replace the values below it 
     14cutindex=find(FLKZ < cutoff); 
     15FLKZ(cutindex)=cutoff; 
     16 
     17 
     18trm=(FLKZ)*sin(dirs); 
  • DPWP/trunk/DPWP/diwasp_1_1GD/private/velz.m

    r168 r334  
    11function trm=velz(ffreqs,dirs,wns,z,depth,xpos,ypos) 
    22 
    3 Kz=sinh(z*wns)./sinh(depth*wns); 
    4 %include a maximum cuttoff for the velocity response function 
    5 Kz(find(Kz<0.1))=0.1; 
    6 Kz(find(isnan(Kz)))=1; 
    7 trm=-i*(ffreqs.*Kz)*ones(size(dirs)); 
     3%Edits were made on 4/20/10 to fix a bug in the cutoff limits in the KZ 
     4%parts of the transfer function.  
     5 
     6%set the cutoff value (for high frequencies) 
     7cutoff=0.1; 
     8 
     9 
     10KZ=sinh(z*wns)./sinh(depth*wns); 
     11FLKZ=ffreqs.*KZ; 
     12 
     13%find the cutoff and replace the values below it 
     14cutindex=find(FLKZ < cutoff); 
     15FLKZ(cutindex)=cutoff; 
     16 
     17 
     18trm=-1i*(FLKZ)*ones(size(dirs)); 
  • DPWP/trunk/DPWP/nortek/nortek_radialtouvw.m

    r314 r334  
    171171%the spectral matrix structure 
    172172SM.freqs=[0.01:0.01:0.4]; 
    173 SM.dirs=[-180:2:180]; 
     173SM.dirs=[0:2:360]; 
    174174SM.xaxisdir= 90; 
    175175%the estimation parameter 
    176176EP.method= 'IMLM'; 
    177 EP.iter=3; 
     177EP.iter=50; 
     178EP.nfft=256; 
    178179 
    179180if data_type == 1 
  • DPWP/trunk/DPWP/nortek/nortek_specmultiplot.m

    r314 r334  
    3838 
    3939%set up data structure for wave info 
    40 Winfo.setup={'EMEP'}'; 
     40Winfo.setup={'IMLM'}'; 
    4141Winfo.hsig=[]; 
    4242Winfo.hconf=[]; 
     
    4444Winfo.dirP=[]; 
    4545Winfo.Ddir=[]; 
    46 Winfo.Spectrum.IMLM=[]; 
     46Winfo.Spectrum=[]; 
    4747 
    4848%load the wave data from the nortek quick wave software 
     
    7070    %For IMLM 
    7171    EP.method='IMLM'; 
    72     EP.iter=3
    73     [I_F01_D2, EPout]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); 
     72    EP.iter=50
     73    [I_F01_D2, ~]=dirspec(ID,SM,EP,{'PLOTTYPE',0,'MESSAGE',0}); 
    7474 
    7575 
     
    7777 
    7878    %[E_WI,I_WI,nortek_WI,fig1,fig2]=nortek_waveplot(E_F01_D2,I_F01_D2,sysinfo,nortekspec); 
    79     [I_WI,fig1,fig2]=nortek_waveplot1(I_F01_D2);  
     79    [I_WI,fig1,fig2]=nortek_waveplot(I_F01_D2);  
    8080     
    8181    saveas(fig1,['polar1_' dateinfo '.fig']); 
     
    100100    Winfo.dirP=[horzcat(Winfo.dirP,dirP)]; 
    101101    Winfo.Ddir=[horzcat(Winfo.Ddir,Ddir)]; 
    102     Winfo.Spectrum.IMLM=[cat(3,Winfo.Spectrum.IMLM,I_F01_D2.S)];     
     102    Winfo.Spectrum=[cat(3,Winfo.Spectrum,I_F01_D2.S)];     
    103103end 
    104104 
  • DPWP/trunk/DPWP/nortek/nortek_waveplot.m

    r314 r334  
    1414freqs=[0.01:0.01:.4]; 
    1515dirres=2; 
    16 dirs=[-180:2:180]; 
     16dirs=[0:2:360]; 
    1717 
    1818% plot the spectrum generated through DIWASP