001 function Snew=dir2enc(S,d,v)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
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
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;
095 w2 = [-flipud(w(2:end));w];
096
097
098
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
107
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)));
112
113
114
115 in=zeros(length(th),length(thw));
116 for j=1:length(w)
117
118
119
120
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);
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
176
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
192 Snew.S(im) = 0.;
193 L2 = dw*sum(Snew.S);
194
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