001 function [m,mtext] = spec2mom(S,nr,vari,even)
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
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054 if nargin<2|isempty(nr)
055 nr=2;
056 end
057 if nargin<4|isempty(even)
058 even=1;
059 end
060 even=logical(even);
061
062 if strcmpi(S.type,'freq')| strcmpi(S.type,'enc')|...
063 strcmpi(S.type,'k1d')
064 if isfield(S,'w')
065 f=S.w(:);
066 vari='t';
067 elseif isfield(S,'k')
068 f=S.k(:);
069 vari='x';
070 else
071 f=2*pi*S.f(:);
072 S.S=S.S/2/pi;
073 vari='t';
074 end
075 S1=abs(S.S(:));
076 m=simpson(f,S1);
077 mtext={'m0'};
078 if nr>0
079 if ~even
080 m=[m simpson(f, f.*S1)];
081 mtext(end+1)={'mi'};
082 end
083 if nr>1
084 m=[m simpson(f,f.^2.*S1)];
085 mtext(end+1)={'mii'};
086 if nr>2
087 if ~even
088 m=[m simpson(f,f.^3.*S1)];
089 mtext(end+1)={'miii'};
090 end
091 if nr>3
092 m=[m simpson(f,f.^4.*S1)];
093 mtext(end+1)={'miiii'};
094 end
095 end
096 end
097 mtext(end+1)={strcat(' where i=',vari)};
098 end
099
100
101
102
103
104
105 else
106 if (nargin<3|isempty(vari))&nr<=1;
107 vari='x';
108 Nv=1;
109 elseif nargin<3|isempty(vari);
110 vari='xt';
111 Nv=2;
112 else% secure the mutual order ('xyt')
113 vari=sort(lower(vari));
114 Nv=length(vari);
115
116 if ( (vari(1)=='t') & (Nv>1)),
117 vari=[vari(2:Nv), 't'];
118 end
119 end
120 S1=S;
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137 if ~strcmpi(S.type(end-2:end),'dir')
138 S1=spec2spec(S,strcat(S.type(1:end-3),'dir'));
139 end
140 if ~isfield(S,'phi') | isempty(S.phi), S1.phi=0;
141
142 end
143 th=S1.theta(:)-S1.phi;
144 Sw=simpson(th,S1.S).';
145 m=simpson(S1.w,Sw);
146 mtext={'m0'};
147
148 if nr>0
149 w=S1.w(:);
150 nw=length(w);
151 if strcmpi(vari(1),'x')
152 Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
153
154 end
155 if strcmpi(vari(1),'y')
156 Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).';
157
158 if strcmpi(vari(1),'x')
159 Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
160 end
161 end
162 if ~isfield(S1,'g')
163 S1.g=gravity;
164 end
165 kx=w.^2/S1.g(1);
166 ky=w.^2/S1.g(end);
167
168 if Nv>=1
169 switch vari
170 case 'x'
171 vec=[kx.*Sc];
172 mtext(end+1)={'mx'};
173 case 'y'
174 vec=[ky.*Ss];
175 mtext(end+1)={'my'};
176 case 't'
177 vec=[w.*Sw];
178 mtext(end+1)={'mt'};
179 end
180 else
181 vec=[kx.*Sc ky.*Ss w*Sw];
182 mtext(end+(1:3))={'mx', 'my', 'mt'};
183 end
184 if nr>1
185 if strcmpi(vari(1),'x')
186 Sc=simpson(th,S1.S.*(cos(th)*ones(1,nw))).';
187
188 Sc2=simpson(th,S1.S.*(cos(th).^2*ones(1,nw))).';
189
190 end
191 if strcmpi(vari(1),'y')|strcmpi(vari(2),'y')
192 Ss=simpson(th,S1.S.*(sin(th)*ones(1,nw))).';
193
194 Ss2=simpson(th,S1.S.*(sin(th).^2*ones(1,nw))).';
195
196 if strcmpi(vari(1),'x')
197 Scs=simpson(th,S1.S.*((cos(th).*sin(th))*ones(1,nw))).';
198
199 end
200 end
201 if ~isfield(S1,'g')
202 S1.g=gravity;
203 end
204
205 if Nv==2
206 switch vari
207 case 'xy'
208 vec=[kx.*Sc ky.*Ss kx.^2.*Sc2 ky.^2.*Ss2 kx.*ky.*Scs];
209 mtext(end+(1:5))={'mx','my','mxx', 'myy', 'mxy'};
210 case 'xt'
211 vec=[kx.*Sc w.*Sw kx.^2.*Sc2 w.^2.*Sw kx.*w.*Sc];
212 mtext(end+(1:5))={'mx','mt','mxx', 'mtt', 'mxt'};
213 case 'yt'
214 vec=[ky.*Ss w.*Sw ky.^2.*Ss2 w.^2.*Sw ky.*w.*Ss];
215 mtext(end+(1:5))={'my','mt','myy', 'mtt', 'myt'};
216 end
217 else
218 vec=[kx.*Sc ky.*Ss w.*Sw kx.^2.*Sc2 ky.^2.*Ss2 w.^2.*Sw kx.*ky.*Scs kx.*w.*Sc ky.*w.*Ss];
219 mtext(end+(1:9))={'mx','my','mt','mxx', 'myy', 'mtt', 'mxy', 'mxt', 'myt'};
220 end
221 if nr>3
222 if strcmpi(vari(1),'x')
223 Sc3=simpson(th,S1.S.*(cos(th).^3*ones(1,nw))).';
224
225 Sc4=simpson(th,S1.S.*(cos(th).^4*ones(1,nw))).';
226
227 end
228 if strcmpi(vari(1),'y')|strcmpi(vari(2),'y')
229 Ss3=simpson(th,S1.S.*(sin(th).^3*ones(1,nw))).';
230
231 Ss4=simpson(th,S1.S.*(sin(th).^4*ones(1,nw))).';
232
233 if strcmpi(vari(1),'x')
234 Sc2s=simpson(th,S1.S.*((cos(th).^2.*sin(th))*ones(1,nw))).';
235
236 Sc3s=simpson(th,S1.S.*((cos(th).^3.*sin(th))*ones(1,nw))).';
237
238 Scs2=simpson(th,S1.S.*((cos(th).*sin(th).^2)*ones(1,nw))).';
239
240 Scs3=simpson(th,S1.S.*((cos(th).*sin(th).^3)*ones(1,nw))).';
241
242 Sc2s2=simpson(th,S1.S.*((cos(th).^2.*sin(th).^2)*ones(1,nw))).';
243
244 end
245 end
246 if Nv==2
247 switch vari
248 case 'xy'
249 vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 kx.^3.*ky.*Sc3s ...
250 kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3];
251 mtext(end+(1:5))={'mxxxx','myyyy','mxxxy','mxxyy','mxyyy'};
252 case 'xt'
253 vec=[vec kx.^4.*Sc4 w.^4.*Sw kx.^3.*w.*Sc3 ...
254 kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc];
255 mtext(end+(1:5))={'mxxxx','mtttt','mxxxt','mxxtt','mxttt'};
256 case 'yt'
257 vec=[vec ky.^4.*Ss4 w.^4.*Sw ky.^3.*w.*Ss3 ...
258 ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss];
259 mtext(end+(1:5))={'myyyy','mtttt','myyyt','myytt','myttt'};
260 end
261 else
262 vec=[vec kx.^4.*Sc4 ky.^4.*Ss4 w.^4.*Sw kx.^3.*ky.*Sc3s ...
263 kx.^2.*ky.^2.*Sc2s2 kx.*ky.^3.*Scs3 kx.^3.*w.*Sc3 ...
264 kx.^2.*w.^2.*Sc2 kx.*w.^3.*Sc ky.^3.*w.*Ss3 ...
265 ky.^2.*w.^2.*Ss2 ky.*w.^3.*Ss kx.^2.*ky.*w.*Sc2s ...
266 kx.*ky.^2.*w.*Scs2 kx.*ky.*w.^2.*Scs];
267 mtext(end+(1:15))={'mxxxx','myyyy','mtttt','mxxxy','mxxyy',...
268 'mxyyy','mxxxt','mxxtt','mxttt','myyyt','myytt','myttt','mxxyt','mxyyt','mxytt'};
269
270 end
271 end
272 end
273 m=[m simpson(w,vec)];
274 end
275
276 end
277