| harmonicseries(AMP,PER,PHA,LEN,DIR,T,X,Y)
|
function Z = harmonicseries(AMP,PER,PHA,LEN,DIR,T,X,Y)
%HARMONICSERIES Sum of harmonic space-time series.
%
% Syntax:
% Z = harmonicseries(AMP,PER,PHA,LEN,DIR,T,X,Y);
%
% Input:
% AMP - Vector of length NH defining each harmonic amplitudes.
% Default: ones.
% PER - Vector of length NH defining each harmonic periods.
% Default: Infs (zero temporal frequency)
% PHA - Vector of length NH defining each harmonic initial phases.
% Default: zeros.
% LEN - Vector of length NH defining each harmonic wavelengths
% Default: Infs (zero spatial frequency)
% DIR - Vector of length NH defining each harmonic direction.
% Default: Infs (zero spatial frequency)
% T - Vector of length NT defining the time space.
% Default: zero.
% X - Matrix of size NYxNX defining the X cartesian coordinates.
% Default: zero.
% Y - Matrix of size NYxNX defining the Y cartesian coordinates.
% Default: zero.
%
% Output:
% Z - 3-dimensional array as the sum of the harmonics, time in the
% rows and space (Y,X) en the columns and third dimension
% respectively. That is: size(Z) = NT x NY x NX.
%
% Desciption:
% This program performs the sum of spatial-time harmonic series defined
% as:
% z(t,y,x) = sum [ a*cos(k*x + l*y - w*t + theta) ]
% where:
% a = AMP - amplitude
% k = K*cos(rdir) - x-component of wavenumber vector
% l = K*sin(rdir) - y-component of wavenumber vector
% K = sqrt(k^2+l^2) - magnitude of wavenumber vector
% = (2*pi)/LEN
% rdir = DIR*(pi/180) - direction of wavenumber vector in
% radians
% w = (2*pi)/PER - Temporal angular frequency
% theta = PHA*(pi/180) - Initial phase in radians.
%
% Notes:
% * Any input can be an empty argument. In this case, defaults are
% used.
% * If the user wants to change the output dimensions, for example,
% space in the rows and columns, he can use the PERMUTE function:
% Z = permute(Z,[2 3 1]); % Z(T,Y,X) => Z(Y,X,T).
% Or the SQUEEZE funtion when T is empty for example:
% Z = squeeze(Z); % Z([],Y,X) => Z(Y,X).
% See the examples.
% * The directions should be in degrees and starting from the x-axis.
% Positive rotation is anti-clockwise.
% * The resultant Z will have the same units as the amplitudes.
% * About the initial phase, for a temporal harmonic series (X,Y
% empty), the initial phase is related with the lapse time
% neccesary to reach the harmonic crest:
% A*cos(w*(t-t0))==A when t = t0 := THETA/w.
% Note that others authors used the complement of this initial
% phase, that is, the spent time from the last maxima:
% A*cos(wt-PHA) (this) vs A*cos(wt+PHA2) (others)
% PHA+PHA2 = 2*pi!
% * NH is the number of harmonics, NT the temporal series length, and
% NYxNX the number of spatial points. Time and space do not need to
% be homogeneous.
%
% Example:
% % Parameters:
% tini = 0; tfin = 100; dt = 1;
% xini = -100; xfin = 100; dx = 4;
% yini = -50; yfin = 50; dy = 5;
% amp = [ 10 11 9];
% per = [200 70 30];
% pha = [ 0 90 180];
% len = [ 10 20 50];
% dir = [ 0 45 150];
% % Grid:
% t = tini:dt:tfin;
% x = xini:dx:xfin;
% y = yini:dy:yfin;
% [x,y] = meshgrid(x,fliplr(y)); % y increase upwards
% % Temporal series:
% zt = harmonicseries(amp,per,pha,[],[],t);
% subplot(221), plot(t,zt)
% title('Temporal series')
% % Espatial 2D series:
% zxy = harmonicseries(amp,[],pha,len,dir,[],x,y);
% zxy = squeeze(zxy);
% subplot(222), surf(x,y,zxy), shading interp, colorbar, view(2)
% title('Espatial series')
% % Waves 1 and 2D (movie):
% ztx = harmonicseries(amp,per,pha,len,[],t,x(1,:));
% ztx = squeeze(ztx);
% ztxy = harmonicseries(amp,per,pha,len,dir,t,x,y);
% zxyt = permute(ztxy,[2 3 1]);
% for k = 1:length(t)/2
% subplot(223), plot(x(1,:),ztx(k,:)), axis tight,
% set(gca,'YLim',sum(amp)*[-1 1])
% xlabel x, title(['1D wave, t: ' num2str(t(k)) ' s'])
% subplot(224), surf(x,y,zxyt(:,:,k)), shading interp, axis tight
% view(2), colorbar, set(gca,'CLim',sum(amp)*[-1 1])
% xlabel x, ylabel y, title(['2D Wave, t: ' num2str(t(k)) ' s'])
% drawnow
% end
% Copyright 2008 Carlos Adrian Vargas Aguilera
% $Revision: 2.0 $ $Date: 2008/06/21 11:30:00 $
% Written by
% M.S. Carlos Adrian Vargas Aguilera
% Physical Oceanography PhD candidate
% CICESE
% Mexico, 2004-2008
% nubeobscura@hotmail.com
%
% Download from:
% http://www.mathworks.com/matlabcentral/fileexchange/loadAuthor.do?objec
% tType=author&objectId=1093874
% Check inputs:
if nargin<6 || nargin>8
error('Harmonicseries:IncorrectNumberOfInputs',...
'At least 6 input arguments are needed (although empties) and less than 8' )
end
if isempty(AMP), AMP = 1; end
if isempty(PER), PER = Inf; end
if isempty(PHA), PHA = 0; end
if isempty(LEN), LEN = Inf; end
if isempty(DIR), DIR = 0; end
if isempty(T), T = 0; end
if nargin<7 || isempty(X), X = 0; end
if nargin<8 || isempty(Y), Y = 0; end
% Check harmonic sizes:
NHs = [length(AMP) length(PER) length(PHA) length(LEN) length(DIR)];
NH = max(NHs);
ind = (NHs<NH);
if any(ind & (NHs~=1))
error('Hamonicseries:IncorrectInputSize',...
'Harmonic parameters should have the same length or a single value for each one.')
end
if ind(1), AMP = repmat(AMP(1),1,NH); end
if ind(2), PER = repmat(PER(1),1,NH); end
if ind(3), PHA = repmat(PHA(1),1,NH); end
if ind(4), LEN = repmat(LEN(1),1,NH); end
if ind(5), DIR = repmat(DIR(1),1,NH); end
% Check spatial sizes:
[NYx,NXx] = size(X);
[NYy,NXy] = size(Y);
NXYx = NYx*NXx;
NXYy = NYy*NXy;
if NXYx~=NXYy
if NXYx==1
X = repmat(X,NYy,NXy);
elseif NXYy==1
Y = repmat(Y,NYx,NXx);
else
error('Harmonicseries:IncorrectInputSize',...
'X and Y must be of the same size, or one of them a single value.')
end
end
[NY,NX] = size(X);
% Time size:
NT = length(T);
% Convertions:
rdir = DIR*(pi/180);
theta = PHA*(pi/180);
w = (2*pi)./PER;
K = (2*pi)./LEN;
k = K.*cos(rdir);
l = K.*sin(rdir);
% Pre-allocation:
Z = zeros(NT,NY,NX);
% 3D Harmonic space-time series:
for ny = 1:NY
for nx = 1:NX
for nt = 1:NT
Z(nt,ny,nx) = sum( AMP .* cos( k*X(ny,nx) + l*Y(ny,nx) - w*T(nt) + theta ) );
end
end
end |
|
