% MATLAB file natural
% Parcel Trajectory for static geopotential field.
% Diagnoses force and acceleration terms in natural coordinates.
% But solves momentum equation in x-y coordinates following parcel.
% If the initial condition is not a steady state solution,
% the parcel tends to wobble as it circles the pressure center
%

%%%%%%% things to tinker with %%%%%%%%%%%%%

% the default system is a closed low (swap=0)
% change swap to 1 if you want a high instead
% but if you want to attempt to find gradient wind
% balance around a high, you must also scale phi so the 
% gradient wind is not an imaginary number (Eq 3.15)
% otherwise there is no steady state (but the model 
% will still run)
swap=0; % 0= low, 1= high, or for sinusoidal 1 swaps H & L
scale=10;  % scale of phi in m/s 10 is good for a low
%scale=3;  % scale of phi in m/s 3 is good for a high

% You might want to 
% change the initial parcel velocity to try to 
% find the gradient wind solution (in other words, 
% to make partial phi / partial s as small as possible)
% This will reduce the wobbles
% You might also want to try to make an anomalous high or low
% solution
% Or possibly you want to try something far out 
v=zeros(2,1);
 v(2)=-6.75; % about right for regular low
% v(2)=32;    % should make an anomalous low
% v(2)=2.95; %a regular high when you have swap=1;
% v(2)=22.1; %an anomalous high when you have swap=1;
% v(1)=-5; v(2)=5; % what do you think this will do?
 
% you can change abs(R) here, where x=-abs(R)+500
% initial parcel position  
X=zeros(2,1);
X(1)=250*1e3;  

%runtime = input('integration time in days  ')
runtime = 2;
%runtime = 10;  % makes nice spirograph looking plot

%%%%% next are things to leave alone for the most part %%%%%%
%%%%% unless you are a compulsive tinkerer, like me    %%%%%%

sinusoidal=0; % 0 makes phi a gaussian, 1 makes phi sinusoidal

% numerical solution
% Best choice is probably the 3rd order Adam's Bashforth. 
% Alternatives are trapezoidal & forwards differencing 
% for acceleration and
% Coriolis terms, but not PGF. 
% The wobble is damped with a forward solution 
% which is timestep dependent (more damping for longer step).
%
trap=1; % use trapezoidal differencing, otherwise use forward
ab3=1; % use Adam's Bashforth 3rd order differencing, trumps trap

% probably should not increase above 360, unless you are interested
% in exploring errors in the numerics
% When the timestep is too big, the
% solution strays away from forming concentric circles even 
% when started with the correct gradient wind speed
%dt = input('timestep in sec  ')
%dt = 1800;  
dt = 360;

cor=1.e-4;
k=2*pi/2.e6;  m=pi/1.e6; % wavenumber if sinusoidal phi
L=.35e6;  % Width of gaussian 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% below here there is nothing fun to change  %%%%%%%%%%%

% First define the grid points on which fields are computed:
% (Distances are in units of 1000 km).
xx=linspace(-1000,1000,30);         % 30 gridpoints in x   
%xx=linspace(-100,1000,30);         % 30 gridpoints in x   
yy=linspace( -500,500,20);          % 20 gridpoints in y
[x,y]=meshgrid(xx*1000,yy*1000);    % plotting matrix for grid system
                                    % in x and y in meters
if swap, scale=-scale; end
if sinusoidal,
  subplot(2,1,1);
  phi = -cor*scale/k*sin(k*x).*cos(m*y);
else, % gaussian shaped high/low
  hi=0; % turn on/off the high for gaussian
  subplot(2,2,1);
  xo=.5e6;
  phi = cor*scale*L*(-exp(-((x-xo).^2.+y.^2)/L.^2)+...
     +hi*exp(-((x+xo).^2.+y.^2)/L.^2));
end
top = .9*max(max(abs(phi)))/9.8;
v1= [0:top/4:top];
contour(x/1000,y/1000,phi/9.8,v1,'k-');
set(gca,'fontsize',18)
hold on
contour(x/1000,y/1000,phi/9.8, [0 0],'r-');
v2=[-top:top/4:0];
contour(x/1000,y/1000,phi/9.8,v2,'k--');
hold off
if swap
 text(-500,0,'L');text(500,0,'H')
else
 text(500,0,'L');text(-500,0,'H')
end
xlabel('Distance along Latitude Circles')
ylabel('Distance along Meridians')
title('Contours of Geopotential')

time = runtime*24.*3600.;
N=time/dt;
r=zeros(2,1); dphi_dx=zeros(2,1);
vprimn= zeros(2,1); vprim1=vprimn; vprim2=vprimn; dt12 = dt/12;
xprimn= zeros(2,1); xprim1=xprimn; xprim2=xprimn; 
XX=zeros(N,1); YY=zeros(N,1);
dphi_dnat=zeros(N,2); 
hold; plot(X(1)/1000,X(2)/1000,'b.','markersize',30);hold

Mat=ones(2,2); Mat(2,1)=dt*cor; 
if trap,  Mat(2,1)=dt*cor/2; end
Mat(1,2)=-Mat(2,1); invM=inv(Mat);

for n= 1:N
  if sinusoidal,
    dphi_dx(1) = -cor*scale*cos(k*X(1)).*cos(m*X(2));
    dphi_dx(2) = cor*scale*m/k*sin(k*X(1)).*sin(m*X(2));
  else
    dphi_dx(1) = -cor*scale*2/L*(-(X(1)-xo)*exp(-((X(1)-xo).^2.+X(2).^2)/L.^2)+...
     +hi*(X(1)+xo)*exp(-((X(1)+xo).^2.+X(2).^2)/L.^2));
    dphi_dx(2) = -cor*scale*2/L*(-(X(2))*exp(-((X(1)-xo).^2.+X(2).^2)/L.^2)+...
     +(X(2))*hi*exp(-((X(1)+xo).^2.+X(2).^2)/L.^2));
  end

  if ab3, % third order AB differencing
    vprim2=vprim1;  vprim1=vprimn;
    vprimn(1)=-dphi_dx(1)+cor*v(2);  vprimn(2)=-dphi_dx(2)-cor*v(1); 
    xprim2=xprim1;  xprim1=xprimn;  xprimn=v;
    v=v+dt12*(23*vprimn-16*vprim1+5*vprim2);
    X=X+dt12*(23*xprimn-16*xprim1+5*xprim2);
  else, % either trapezoidal or forward
    r=v-dt*dphi_dx;
    if trap
      r(1)=r(1)+dt*cor/2*v(2);
      r(2)=r(2)-dt*cor/2*v(1);
    end
    vold=v;
    v=invM*r;

    if trap
      X=X+(v+vold)/2*dt;
    else
      X=X+v*dt;
    end
  end

% apply the periodic BC
  if X(1) > 3.e6
     X(1)= -6.e6+X(1);
  end
  if X(1) < -3.e6
    X(1) = 6.e6+X(1);
  end

  XX(n)=X(1);    YY(n)=X(2);

  % diagnose the motion in natural coordinates
  alpha(n) = atan(v(2)/v(1));
  if (v(1)<0), alpha(n)=alpha(n)+pi; end
  rotate = cos(alpha(n))*ones(2,2);
  rotate(1,2) = sin(alpha(n)); rotate(2,1) = -rotate(1,2);
  
  dphi_dnat(n,:) = rotate*dphi_dx;
  V(n)=sqrt(sum(v.^2));
  Co(n)=cor*V(n);
  Ce(n)=-dphi_dnat(n,2)-Co(n);
  usav(n)=v(1);  vsav(n)=v(2);
end

skip=round(3600/dt);  
n=1:skip:N;
hold; plot(XX(n)/1000,YY(n)/1000,'r-');hold
hold; plot(XX(n)/1000,YY(n)/1000,'r.');hold
if sinusoidal==0, axis([0 1000 -500 500]); end

%fid=fopen('natural.txt','w');
%fprintf(fid,'%8.2f  %8.2f %8.2f %8.2f %8.2f \n',[XX/1000; YY/1000; V; usav; vsav]);


t=1:dt:time; t=t/3600;  n=1:1:N;
subplot(2,2,3);
plot(t(n),-dphi_dnat(n,1),   'linewidth',2);
hn=legend('-\partial \phi  / \partial s');set(hn,'fontsize',12)
set(gca,'fontsize',18)
set(hn,'fontsize',12)
xlabel('hour')
xlim([0 inf])
subplot(2,2,4);
plot(t(n),-dphi_dnat(n,2),t(n),Co(n),'--',t(n),Ce(n),'-.',...
   'linewidth',2);
hn=legend('-\partial \phi  / \partial n','Co','Ce');
set(gca,'fontsize',18)
set(hn,'fontsize',12)
xlabel('hour')
xlim([0 inf])
mean(-dphi_dnat)