DIR2K1D Transforms directional spectrum to wavenumber spectrum.
returns the constant acceleration of gravity | |
Numerical integration with the Simpson method | |
Calculates spectral moments from spectrum | |
Plot a spectral density | |
Difference and approximate derivative. | |
Hold current graph. | |
2-D interpolation (table lookup). | |
Linearly spaced vector. | |
Average or mean value. | |
Linear plot. | |
Remove fields from a structure array. |
001 function Snew=dir2k1d(S); 002 % DIR2K1D Transforms directional spectrum to wavenumber spectrum. 003 004 Snew=S; 005 g=gravity; 006 % Commence with adding S(.,theta)+S(.,theta+pi) for -pi/2<theta<pi/2 007 % since instantaneously there is no difference between waves coming 008 % from direction -pi/4 and 3*pi/4 for example. 009 dth=diff(S.theta(1:2)); 010 th0=min(S.theta(S.theta>=0)); 011 th=[fliplr(th0-dth:-dth:-pi/2) th0:dth:pi/2]'; 012 [r2,c]=size(S.S(S.theta>pi/2+1e-13,:)); 013 [r3,c]=size(S.S(S.theta<-pi/2-1e-13,:)); 014 015 m1=length(th); 016 if r2==0 % no thetas over pi/2, which is the largest theta? 017 m1=find(abs(S.theta(end)-th)<=1e-8); 018 end 019 m2=1; 020 if r3==0 % no thetas under -pi/2, which is the smallest theta? 021 m2=find(abs(S.theta(1)-th)<=1e-8); 022 end 023 S0=zeros(length(th),length(S.w)); 024 S0(m2:m1,:)=S.S(r3+1:end-r2,:); 025 S0(1:r2,:)=S0(1:r2,:)+S.S(end-r2+1:end,:); 026 S0(end-r3+1:end,:)=S0(end-r3+1:end,:)+S.S(1:r3,:); 027 028 % Divide the integral in two parts, abs(th)<=a0 and abs(th)>a0 029 k=linspace(0,S.w(end)^2/g,length(S.w))'; 030 a0=pi/4; 031 % the case abs(th)<=a0 032 th1=th(th>=-a0 &th<=a0); 033 S1=S0(th>=-a0 &th<=a0,:); 034 k1=k(2:end,ones(1,length(th1))); 035 th1=th1(:,ones(1,length(S.w)-1))'; 036 wmax=max(max(sqrt(g*k1./cos(th1)))); 037 S.w=[S.w; S.w(end)+diff(S.w(end-1:end));wmax*1.1]; 038 S1=[S1 zeros(size(S1,1),2)]; % add zeros to prevent NaN in interpolation 039 Sip=interp2(S.w,th1(1,:),S1,sqrt(g*k1./cos(th1)),th1); 040 Snew.S=[0; sqrt(g./k(2:end)).*simpson(th1(1,:),Sip'./sqrt(cos(th1))',1).'/2]; 041 042 % the case abs(th)>a0 043 th1=th; 044 S1=S0; 045 wmax=sqrt(g*k(end)/cos(a0)); 046 n=2; 047 if wmax*1.1>S.w(end) 048 S.w=[S.w; wmax*1.1]; 049 n=3; 050 end 051 S1=[S1 zeros(size(S1,1),n)]; % add zeros to prevent NaN in interpolation 052 S1=[S1(end,:); S1; S1(1,:)]; % extend angles to prevent NaN in interpolation 053 th1=[th1(1)-dth;th1;th1(end)+dth]; 054 Sk=zeros(size(k)); 055 for j=2:length(k) 056 w=S.w(S.w>=sqrt(g*k(j)/cos(a0))); 057 thw=acos(g*k(j)./w.^2); 058 Sip=(interp2(S.w,th1,S1,w,thw)+(interp2(S.w,th1,S1,w,-thw)))./sqrt(w.^4-(g*k(j))^2); 059 if length(w(Sip>0))==2 060 Sk(j)=mean(Sip(Sip>0))*diff(w(Sip>0)); 061 elseif length(w(Sip>0))>2 062 Sk(j)=g*simpson(w(Sip>0),Sip(Sip>0),1); 063 % not perfect here, somtimes change of integration limit in simpson 064 end 065 end 066 hold off 067 plot(k,Snew.S,k,Sk) 068 069 Snew.S=Snew.S+Sk; 070 Snew=rmfield(Snew,'w'); 071 Snew=rmfield(Snew,'theta'); 072 Snew.k=k; 073 Snew.type='rotk1d'; 074 spec2mom(Snew) 075 hold on 076 wspecplot(Snew,1,'r')
Comments or corrections to the WAFO group