Home > wafo > spec > private > dir2k1d.m

# dir2k1d

## PURPOSE

Transforms directional spectrum to wavenumber spectrum.

Snew=dir2k1d(S);

## DESCRIPTION

`  DIR2K1D Transforms directional spectrum to wavenumber spectrum.`

## CROSS-REFERENCE INFORMATION

This function calls:
 gravity returns the constant acceleration of gravity simpson Numerical integration with the Simpson method spec2mom Calculates spectral moments from spectrum wspecplot Plot a spectral density diff Difference and approximate derivative. hold Hold current graph. interp2 2-D interpolation (table lookup). linspace Linearly spaced vector. mean Average or mean value. plot Linear plot. rmfield Remove fields from a structure array.
This function is called by:

## SOURCE CODE

```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')```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003