Home > wafo > spec > private > dir2k1d.m

dir2k1d

PURPOSE ^

Transforms directional spectrum to wavenumber spectrum.

SYNOPSIS ^

Snew=dir2k1d(S);

DESCRIPTION ^

  DIR2K1D Transforms directional spectrum to wavenumber spectrum.

CROSS-REFERENCE INFORMATION ^

This function calls: 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