Home > wafo > spec > private > dir2enc.m

dir2enc

PURPOSE ^

Transform a dir. spectrum to an encountered. (Used in spec2spec)

SYNOPSIS ^

Snew=dir2enc(S,d,v)

DESCRIPTION ^

 DIR2ENC Transform a dir. spectrum to an encountered. (Used in spec2spec)
 
  CALL:  Snew=dir2enc(S,d,v)
 
       Snew = encountered spectrum, directional or frequency type
       S    = directional spectrum
       d    = dimension of result, 1=frequency type, 2=directional type
                                                  (default 1)
       v    = speed of ship     (default 0)
 
  Evaluates the spectrum of the seaway encountered by a ship travelling 
  along the x-axis (rotated by in a constant angle (S.phi)) 
  with constant speed (v) from the directional spectrum of the
  seaway.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function Snew=dir2enc(S,d,v)
002 %DIR2ENC Transform a dir. spectrum to an encountered. (Used in spec2spec)
003 %
004 % CALL:  Snew=dir2enc(S,d,v)
005 %
006 %      Snew = encountered spectrum, directional or frequency type
007 %      S    = directional spectrum
008 %      d    = dimension of result, 1=frequency type, 2=directional type
009 %                                                 (default 1)
010 %      v    = speed of ship     (default 0)
011 %
012 % Evaluates the spectrum of the seaway encountered by a ship travelling 
013 % along the x-axis (rotated by in a constant angle (S.phi)) 
014 % with constant speed (v) from the directional spectrum of the
015 % seaway.
016 %
017 
018 % References: Podgorski et al. (2000)
019 %             Exact distributions for apparent waves in irregular seas
020 %             Ocean Engng. 27, p. 979-1016.
021 
022 % Tested on MATLAB 5.3  
023 % History: 
024 % revised:
025 % pab 21.09.2004: Added todo statements  
026 %  ir 26.06.04 BUG's removals plus major changes.
027 %  ir 04.04.01 BUG's removals plus major changes.
028 %  es 05.06.00 simpson in stead of trapz, clean up 
029 %  by es 18.08.1999
030 
031 % TODO % The numerical accuracy of the method is pure. 
032 % TODO % Should remove singularity by e.g. splitting the integral.
033 
034 error(nargchk(1,3,nargin))  
035 if nargin<0|isempty(S)
036   error('Needs an input spectrum')
037 end
038 
039 if ~any(strcmpi(S.type,{'dir','freq'}))
040   error('Spectrum must be of type dir or freq ')
041 end
042 
043 if (nargin<2|isempty(d))
044   d=1;
045 end
046 if ( nargin<3|isempty(v))
047   v=0;
048 end
049 
050 Snew=S;
051 if isfield(Snew,'phi')
052   phi=Snew.phi;
053 else
054   phi=0;
055   Snew.phi=phi;
056 end
057 
058 if ~isfield(Snew,'theta')
059   Snew.theta=0;
060 end
061 
062 Snew.v=v;
063 if v<0
064   disp('Message: Negative speed, interpreted as opposite direction')
065   v=abs(v);
066   phi=phi+pi;
067 end
068 
069 
070 if (v==0)
071   if d==1
072     Snew=spec2spec(Snew,'freq');
073     Snew.type='enc';
074   else
075     Snew.type='encdir';
076   end
077   return
078 end
079 
080 
081 th=Snew.theta;
082 if isfield(Snew,'f')
083   w=2*pi*Snew.f(:);
084   Sf=Snew.S/2/pi;
085 else
086   w=Snew.w(:);
087   Sf=Snew.S;
088 end
089 
090 if length(th)>1
091   % IR changes 26 june 2004, formula (6) in Podgorski et al. 
092   [w1,th1] = meshgrid(w,th);
093   Sfneg    = interp2(w,S.theta,Sf,w1,mod(th1+2*pi,2*pi)-pi);
094   Sf       = [fliplr(Sfneg(:,2:end)) Sf]/2; % unitary spectum - symmetric in w
095   w2       = [-flipud(w(2:end));w];
096 
097   % Now we need to make the transformation that 
098   %        Sf1(theta,omega)<= Sf(theta+phi,omega).
099 
100   [w1,th1] = meshgrid(w2,th);
101   
102   Sf1      = interp2(w2,S.theta,Sf,w1,mod(th1+phi+pi,2*pi)-pi);
103 
104   Snew.S = Sf1;
105   g      = gravity;
106   % Since the encountered spectrum has usually higher frequencies we
107   % need to use longer vector of frequencies omega.
108 
109   w        = linspace(0,2*max(w),2*length(w));
110   thw      = pi*ones(size(w));
111   thw(w>0) = acos(max(-1,-g/4/v./w(w>0))); % This computes integration 
112   % limits in  theta
113   %plot(w,thw)
114 
115   in=zeros(length(th),length(thw));
116   for j=1:length(w)                     
117     % The limits depends on omega
118     % and hence for each omega we
119     % check which theta satisfies
120     % the integration limits.
121     in(abs(th)<thw(j),j)=1;
122   end
123 
124   in = logical(in);
125   s1 = zeros(length(th),length(w2));
126   s2 = s1;
127   for j=1:length(th)
128     ix = find(in(j,:));
129     if any(ix)
130       h = g/(2*v*cos(th(j)));
131       a = sqrt(1+2*w(ix)/h);
132       s1(j,ix) = (interp1(w2,Sf1(j,:),h*(-1+a)))./a;
133       s2(j,ix) = (interp1(w2,Sf1(j,:),h*(-1-a)))./a;
134     end
135   end
136   s1(isnan(s1)) = 0;
137   s2(isnan(s2)) = 0;
138   Snew.S        = 2*(s1+s2);
139   Snew.w        = w;
140   Snew.type     = 'encdir';
141   Snew.phi      = 0.;
142 
143   if d==1
144     Snew.S=simpson(Snew.theta,Snew.S,1); %integrate out angle
145   end
146 
147 else 
148   if abs(phi-pi/2)>0.00001
149 
150     g  = gravity;
151     c  = 4*v/g;
152     w  = linspace(0,2*max(S.w),2*length(S.w));
153     Dpl   = zeros(size(w));
154     Dmin  = Dpl;
155     kapa  = phi;      
156     h1    = 2./c./cos(phi);
157     delta = 1+c*w*cos(phi);
158     ind   = find(delta>=0);
159     y1    = sqrt(delta(ind));
160     lambda1 = h1.*(-1+y1);
161     lambda2 = h1.*(-1-y1);
162     D_1   = interp1(S.w,S.S,lambda1,'spline',0)./y1;
163     D_2   = interp1(S.w,S.S,lambda2,'spline',0)./y1; 
164     Dpl(ind) = D_1+D_2; 
165     
166     h1    = 2./c./cos(phi);
167     delta = 1-c*w*cos(phi);
168     ind   = find(delta>=0);
169     y1    = sqrt(delta(ind));
170     lambda1 = h1.*(-1+y1);
171     lambda2 = h1.*(-1-y1);
172     D_1 = interp1(S.w,S.S,lambda1,'spline',0)./y1;
173     D_2 = interp1(S.w,S.S,lambda2,'spline',0)./y1; 
174     %
175     % Note that actually the encountered spectrum is always a directional
176     % spectrum. 
177     %
178     Dmin(ind) = D_1+D_2;
179     Snew.w    = [fliplr(-w) w];
180     Snew.S    = [fliplr(Dmin) Dpl];
181     Snew.type = 'encdir';
182     if d==1
183        Snew.type = 'enc';
184        Snew.S    = (Dmin+Dpl)';
185        Snew.w    = w';
186        if abs(phi)>pi/4
187           [mm,im] = max(Snew.S);
188           dw0     = abs(S.w(end)-S.w(end-1));
189           dw      = abs(w(end)-w(end-1));
190           L0      = dw0*sum(S.S);
191           %L1=dw*sum(Snew.S)-0.5*mm*dw
192           Snew.S(im) = 0.;
193           L2         = dw*sum(Snew.S);
194           %[L0 L1 2*(L0-L2)/dw]
195           Snew.S(im)=max(0,2*(L0-L2)/dw);
196        end
197     end
198   end
199   
200 end
201 
202 if d==1
203   Snew = rmfield(Snew,'theta');
204   Snew.type = 'enc';
205 end
206

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