% Laboratory in Oceanography - Data and Methods
% SMAST, UMass Dartmouth
%
% For understanding rotary spectra
%
% Written by Miles A. Sundermeyer, 5/5/09

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% start by creating a u,v time series with a single frequency, but varying amplitudes
T = 10*86400;				% (s) length of time array in days
N = 1048;				% number of points in time array (power of 2)
dt = T/N;				% (s)

t = (1:N)*dt;				% (s)
omega = 1*1/86400;			% (cps)

if(0)
  % u
  A = 1*[1 1e-15 1 1];					% (m/s)
  B = 1*[1 1 1e-15 1];					% (m/s)

  % v
  C = [1 1 1e-15 -1];					% (m/s)
  D = [1 1e-15 1 -1];					% (m/s)
else
  % u
  A = [1 0 1];					% (m/s)
  B = [1 1 0];					% (m/s)

  % v
  C = [1 1 0];					% (m/s)
  D = [1 0 1];					% (m/s)
end

for m=1:length(A)
  u = A(m)*cos(2*pi*omega*t) + B(m)*sin(2*pi*omega*t);		% (m/s)
  v = C(m)*cos(2*pi*omega*t) + D(m)*sin(2*pi*omega*t);		% (m/s)

  U = u + i*v;				% (m/s) complex velocity

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(m)
  clf

  subplot(3,1,1)
  plot(t/86400,u,'b-',t/86400,v,'r--','linewidth',2)
  title({	['u=',num2str(A(m)),'cos(\omegat)+',num2str(B(m)),'sin(\omegat)']; ...
	['v=',num2str(C(m)),'cos(\omegat)+',num2str(D(m)),'sin(\omegat)']});
  xlabel('t (days)')
  ylabel('u,v (m/s)')
  grid on

  subplot(3,2,3)
  plot(U,'linewidth',2);
  xlabel('u (m/s)')
  ylabel('v (m/s)')
  title('U = u + iv')
  axis equal
  axis square
  axis(2*[-1 1 -1 1])
  grid on

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % now compute rotary components, cw & ccw
  Rp = 0.5*[A(m)+D(m) + i*(C(m)-B(m))];
  Rm = 0.5*[A(m)-D(m) + i*(C(m)+B(m))];

  % magnitudes
  absRp = abs(Rp);
  absRm = abs(Rm);

  % major, minor axes of ellipse
  major = sqrt(abs(absRp + absRm));
  minor = sqrt(abs(absRp - absRm));
  %major = sqrt(absRp.^2 + absRm.^2);
  %minor = sqrt(absRp.^2 - absRm.^2);

  % orientation
  epsilonp = atan((C(m)-B(m))/(A(m)+D(m)));
  epsilonm = atan((C(m)+B(m))/(A(m)-D(m)));
  theta = 0.5*(epsilonp + epsilonm);
  phase = 0.5*(epsilonp - epsilonm);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  subplot(3,2,4)
  axis([0 1 0 1])
  axis off
  text(0.5,0.5, ...
	{['|R^+|=',num2str(absRp),', |R^-|=',num2str(absRm)]; ...
%	['\epsilon^+=',num2str(epsilonp),', \epsilon^-=',num2str(epsilonm)]; ...
	['major axis=',num2str(major),', minor axis=',num2str(minor)]; ...
	[' \theta=',num2str(180/pi*theta),', phase@t=0=',num2str(180/pi*phase)]}, ...
	'horizontalAlignment','center','verticalalignment','middle', ...
	'horizontalAlignment','center','verticalalignment','middle' ...
      ); 

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % create and plot the rotary pieces
  Rpt = Rp*(cos(2*pi*omega*t) + i*sin(2*pi*omega*t));
  Rmt = Rm*(cos(2*pi*omega*t) - i*sin(2*pi*omega*t));

  subplot(3,2,5)
  plot(Rpt,'linewidth',2);
  xlabel('u^+ (m/s) (CCW)')
  ylabel('v^+ (m/s) (CCW)')
  axis equal
  axis square
  axis(2*[-1 1 -1 1])
  grid on

  subplot(3,2,6)
  plot(Rmt,'linewidth',2);
  xlabel('u^- (m/s) (CW)')
  ylabel('v^- (m/s) (CW)')
  axis equal
  axis square
  axis(2*[-1 1 -1 1])
  grid on

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Now compute the rotary auto spectrum
  df = 1/T;
  f = df*(-N/2:(N/2-1))*86400;

  fRpt = fft(Rpt)/N;
  fRmt = fft(Rmt)/N;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure(100+m)
  clf
  subplot(2,2,1)
  plot(f, (fftshift(fRpt.*conj(fRpt))), 'b');
  hold on
  plot(f, (fftshift(fRmt.*conj(fRmt))), 'r');
  xlabel('frequency (cpd)');
  ylabel('Pxx (m^2/s^2)');
  title('Rotary Spectra (b=''+''(CCW), r=''-''(CW)');

  set(gca,'xlim',10*[-1 1]);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if(1)
    figure(m)
    eval(['print -djpeg rotary_spectra_',num2str(m)]);
    figure(100+m)
    eval(['print -djpeg rotary_spectra_10',num2str(m)]);
  end

end