function Rossby_wave_ray_tracing_ala_Hoskins_Karoly_1981
%% Rossby_wave_ray_tracing_ala_Hoskins_Karoly_1981.m
%% Nathan Arnold, with some errors introduced by Eli
%% For EPS231, spring 2011, for HW-04
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    clear all;

    global Omega a w phijet ujet
    
    Omega=2*pi/86400;   % earth rotation frequency
    a=6.38e6;           % earth radius, meters
    tmax = 86400*10;    % integration length
    dt = 200;           % timestep (s)
    nt = tmax / dt;     % number timesteps
    
    n = 5;              % integer zonal wavenumber (i.e. 1,2,3...)
    
    w=Omega/30;         % constant angular velocity flow
    phijet = 45;        % jet latitude    
    ujet = 0;           % midlatitude jet velocity

    lon_init = 0;       % initial longitude
    lat_init = 5;       % initial latitude
    
    x = nan(nt,1);
    y = nan(nt,1);
    l_hist = nan(nt,1);
    dldt_hist = nan(nt,1);
    x(1) = lon_init*a*pi/180;
    y(1) = phi_to_y(lat_init);
    

    k=n/a;
    l_hist(1) = sqrt(beta_m(y(1))/ubar_m(y(1))-k^2);
    for i=1:nt
                
        Ks = sqrt(beta_m(y(i)) / ubar_m(y(i)));
        
        dldt = k*dbeta_dy(y(i))/Ks^2+(k*w*sech(y(i)/a)*tanh(y(i)/a));
        dldt_hist(i) = dldt;
        
        l_hist(i+1) = l_hist(i) + dldt*dt;
        l_int = (l_hist(i+1)+l_hist(i))/2;    % use "midpoint" l to calculate cg
        
        cg_x = ubar_m(y(i))+beta_m(y(i))*(k^2-l_int^2)/(k^2+l_int^2)^2
        cg_y = 2*beta_m(y(i))*k*l_int/(k^2+l_int^2)^2;

        x(i+1) = x(i) + cg_x*dt;
        y(i+1) = y(i) + cg_y*dt;

    end


    figure(1);
    clf
    %% plot ray in spherical polar projection:
    %subplot(2,2,1)
    axes('position',[.1  .5  .4  .4]);
    hold off;
    h=polar([0,180],[0,80]);
    set(h,'Visible','Off')
    hold on;
    polar(x/a,y_to_phi(y)-90)
    h1=title(sprintf('Ray Path, k=%g/R_e',a*k));
    set([gca,h1],'FontSize',12)
        
    subplot(2,2,2)
    [ax,hl1,hl2]=plotyy([1:dt:tmax]/86400,l_hist(1:nt) ...
                        ,[1:dt:tmax]/86400,dldt_hist(1:nt));
    h1=ylabel(ax(1),'l');
    h2=ylabel(ax(2),'dl/dt');
    set(hl1,'LineWidth',2)
    set(hl2,'LineWidth',2)
    h3=xlabel(ax(1),'Time (days)');
    xlabel(ax(2),' ');
    set(ax(2),'xtick',[])
    set([ax(1),ax(2),h1,h2,h3],'FontSize',12)
    grid on;
    
    
    subplot(2,2,4)
    K_s_plot=sqrt( beta_m(y(1:nt))./ ubar_m(y(1:nt)) );
    [ax,hl1,hl2]=plotyy([1:dt:tmax]/86400,y_to_phi(y(1:nt)) ...
                        ,[1:dt:tmax]/86400,a*K_s_plot );
    h1=ylabel(ax(1),'Ray Latitude');
    h2=ylabel(ax(2),'R_eK_s');
    set(hl1,'LineWidth',2)
    set(hl2,'LineWidth',2)
    h3=xlabel(ax(1),'Time (days)');
    xlabel(ax(2),' ');
    set(ax(2),'xtick',[])
    set([ax(1),ax(2),h1,h2,h3],'FontSize',12)
    grid on;        
    
    %% plot the mean flow:
    subplot(2,2,3)
    yplot=0:1.5*pi*a/2;
    phi=y_to_phi(yplot);
    ubar_plot=ubar_m(yplot).*cos(phi*pi/180);
    hl1=plot(phi,ubar_plot);
    set(hl1,'LineWidth',2)
    h1=ylabel('Ubar');
    h2=xlabel('latitude');
    h3=title(sprintf('$$\\overline{u}(\\phi)$$, with $$u_{jet}=%g$$',ujet),'interpreter','latex');
    set(h3,'units','normalized'); 
    P=get(h3,'Position');
    set(h3,'Position',[P(1) P(2)-0.02 P(3)]);
    set([gca,h1,h2,h3],'FontSize',12)
    axis tight;
    grid on;        
    
    %% save as pdf:
    set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [7 5]);
    set(gcf, 'PaperPosition', [0 0 7 5]); % [left, bottom, width, height];
    saveas(gcf,sprintf('Rossby-rays-a-la-Hoskins-Karoly-1981.pdf'));


function out = phi_to_y(phi)
    global a
    out = a*log((1+sin(phi*pi/180))./cos(phi*pi/180));


function out = y_to_phi(y)
    global a
    out = asin(tanh(y/a))*180/pi;


function out = ubar_m(y)
    global a w phijet ujet
    phi = y_to_phi(y);
    out = ...
        + a*w ...  % constant angular velocity
        + ujet*exp(-(phi-phijet).^2/(15)^2)./cos(phi*pi/180); % jet


function out = beta_m(y)
    global Omega a w ujet phijet
    
    yjet = phi_to_y(phijet);
    dely = 800e3;
    
    angmomflw = 2*(Omega+w)*sech(y/a).^2 / a;

    f = 1/a^2 + 2/dely^2 - 4*(y-yjet).^2/dely^4 - 4*(y-yjet).*tanh(y/a)/(dely^2*a);

    jet = ujet * f .* cosh(y/a) .* exp(-(y-yjet).^2/dely^2);

    out = angmomflw + jet;


function out = dbeta_dy(y)
    global Omega a w ujet phijet
    
    yjet = phi_to_y(phijet);
    dely = 800e3;
    
    angmomflw = -4*(Omega + w) * tanh(y/a).*sech(y/a).^2 / a^2;
    
    f = 1/a^2 + 2/dely^2 - 4*(y-yjet).^2/dely^4 - 4*(y-yjet).*tanh(y/a)/(dely^2*a);
    
    dfdy = -8*(y-yjet)/dely^4 - 4*tanh(y/a)/(dely^2*a) - 4*(y-yjet).*sech(y/a).^2/(a^2*dely^2);
    
    jet = ujet * (dfdy - f.*tanh(y/a)/a - 2*(y-yjet).*f/dely^2) .* cosh(y/a) .* exp(-(y-yjet).^2/dely^2);
    
    out = angmomflw+jet;
    

