% KELVIN
% simulacija Kelvinovih valova;
% x-os se proteze duz obale,  y-os okomito na obalu
%
% crtanje disperzivne relacije
% crtanje  u vremenu denivelacije i vektora struja
%
% Mira Pasaric (mpasaric@irb.hr)
% 01.12.2004.
% =======================================

XX = 1000*1e3;     % max. duljina obale koju promatramo
YY = 500*1e3;       % max. udaljenost od obale
f = 1e-4;                % Coriolosov parametar na 45 deg
g = 9.81;               % ubrzanje sile teze
H = input('Dubina dna [m]; (npr 50) __ ? \n')       % dubina dna

 x = [0:10e3:XX]';          % koordinate paralelno s obalom
 y = [0:10e3:YY]';          %       "          okomito na obalu

R = sqrt(g*H)/f;                                       % Rossby radius deformacije 
fprintf('\n Rossby-ev radius deformacije, R = %i  km\n\n', round(R*1e-3));
c = sqrt(g*H);                                          % fazna brzina
fprintf(' Fazna brzina, c = %3.1f  m/s \n\n', c);

 % ======================= DISPERZIVNA RELACIJA
if(0)           % crtanje disp. relacije
	Hmin = 10;
	RR = sqrt(9.81.*Hmin)/f;           % odgovarajuci Rossby barotropni radijusi deformacije
	kmin = 0; kmax = 1/RR; dk = (kmax-kmin)/10;        % raspon valnih brojeva
	k = [0:dk:kmax]';
    om = k*c;                                                   % kruzna frekvencija
			plot(k*1e3,om*3600,'-','LineWidth',1); hold on; 
			plot(0,f*3600,'xr','MarkerSize',12); plot(1/R*1e3,0,'rx','MarkerSize',12); grid on;
			xlabel('Valni broj [km^-^1]');
			ylabel('Kruzna frekvencija [h^-^1]');
            text(1/R*1e3,0,'1/R','FontWeight','bold','VerticalALignment','top','FontSize',12);
            text(1/R*1e3, 0, ['(H=' num2str(H) ' m)'],'VerticalALignment','bottom');
            text(0,f*3600,'f_4_5  ','FontWeight','bold','HorizontalALignment','right','FontSize',12);
            title('Kelvinovi valovi; disperzivna relacija','FontSize',12);

	clear RR k* om      
end

alamb = input(' Valna duljina [m] = a*Rossby r.d.;     a = (npr. 2)   __ ? \n');        % valna duljina

k = 2*pi/(alamb*R);                                     % valni broj    
om = k*c;                                                    % kruzna frekvencija
per = 2*pi/om;                                             % period
fprintf(' Period = %4.1f  h \n', per/3600);

[X Y] = meshgrid(x, y);                                 % generiranje 2d mreze                                     

vw = [-30 30];                                              % rotacija 3d slike; azimut i elevacija;  
sl = 0;
rper = 4;                                                       % vremenski korak (dio perioda)
for t = [0:per/rper:per];                                % vremenska petlja
     sl = sl+1;
     omt = om*t;
     zeta = exp(-Y./R).*cos(k.*X - omt);            % denivelacija (m)
     u = sqrt(g/H) * zeta;                                   % brzina u x smjeru
     figure
     %  ---------------------------------     crtanje denivelacije
     mesh(X*1e-3,Y*1e-3,zeta); view(vw); hold on;           % crtanje 3-d slike
     set(gca,'XLim',[0 XX]*1e-3);            set(gca,'YLim',[0 YY]*1e-3);
     xlb = xLabel('Udaljenost duz obale [km]'); set(xlb,'VerticalAlignment','top','HorizontalAlignment','center')
     ylb = yLabel('Udaljenost od obale [km]'); set(ylb,'VerticalAlignment','bottom','HorizontalAlignment','center')
     %  ---------------------------------    crtanje vektora struja
    v = u; v(:,:) = 0;
    Xq = X(1:10:end,1:10:end)*1e-3;    Yq = Y(1:10:end,1:10:end)*1e-3;
    uq = u(1:10:end,1:10:end)*1e-3;    vq = v(1:10:end,1:10:end)*1e-3;
    q1 = quiver(Xq,Yq, uq, vq) ;   hidden,'off';
    set(q1,'LineWidth',2,'Color','black');
    view(vw)
    
    tt = rats(t/per); tt1 = deblank(tt(end:-1:1)); tt = deblank(tt1(end:-1:1));
    ttl = ['t = ' tt '*T; H = ' num2str(H) ' m'];
    title(ttl)
    %       generiranje jpg slika za animaciju
    film = 0;       
    if(film);
          ime_slike = ['Kelvin_h' sprintf('%i',H) '_' sprintf('%2.2i',sl)];
          print ('-djpeg90 ', ime_slike);    
    end

end                  % for t = [0:per/16:per];            % vremenska petlja
