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:
 gravity returns the constant acceleration of gravity simpson Numerical integration with the Simpson method spec2spec Transforms between different types of spectra error Display message and abort function. interp1 1-D interpolation (table lookup) interp2 2-D interpolation (table lookup). isfield True if field is in structure array. linspace Linearly spaced vector. meshgrid X and Y arrays for 3-D plots. rmfield Remove fields from a structure array. strcmpi Compare strings ignoring case.
This function is called by:
 spec2spec Transforms between different types of spectra

## 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