% TOPO_EDNO
% simulacija analitickog rjesenja za TOPOGRAFSKE valove,
% generirane u bazenu s eksponencijalnim porastom dubine
% koordinatni sustav; x - okomito na obalu, prema otvorenom moru
%                               y - paralelno s obalom
% primjer (parametri) za juzni Jadran
%
% Mira Pasaric (mpasaric@irb.hr)
% 01.12.2004.
% =======================================
% POZIV POTPROGRAMA:
% sub_nultocka; fja
%
% ------------------------------- fizicki parametri bazena
L_shelf = 70;          % sirina shelfa [km]
Y_shelf = 100;        % duljina shelfa [km], odn. (long. valna duljina) koja se promatra
d00 = 84;               % pocetna dubina dna [m]
b0 = 18.5;              % nagib dna [m]
A = 0.25;                % koeficijent trenja [m^3/s]
% ============== mogucnost interaktivnog unosa parametara 
% L_shelf = input(' Sirina shelfa [km], (npr 70) __?  ');                                  %    sirina shelfa
% Y_shelf = input(' Duljina obale [km] (npr 100) __?  ');                               %    duljina obale koja se crta
% d00 = input(' Dubina dna uz obalu [m] (npr 84) __?  ');                            %    dubina dna uz obalu
% b0 = input(' Koeficijent nagiba dna [m/km] (npr 18.5) __?  ');                   %    koeficijent (eksponencijalnog) nagiba dna
% A = input(' Koeficijent trenja [m^3/s] (npr 0.25) __?  '); 

x = [0:1:L_shelf]';
d0 = d00*1e-3;
b = b0*1e-3;  bb = repmat(b,length(x),1);
y = [0:1:Y_shelf]'; 
gam = 2*pi/Y_shelf; gamm = repmat(gam,length(y),1);                                 % longitudinalni valni broj
f = 1e-4;                                                                                                         % Coriolosov parametar na 45 deg

% % ============== crtanje dubine
dub = [0; repmat(d0,length(x),1).*exp(bb.*x)*1e3]; dub = [dub; dub(end)];
plot([0; x; x(end)*1.2],dub,'LineWidth',4); set(gca,'YDir','reverse');
set(gca,'XLim',[0 x(end)*1.2])
xlabel('Udaljenost od obale (km)'); ylabel('Dubina dna (m)');
title ('Profil dna','FontWeight','bold');
% print('-djpeg90','topo_e_dno')

% ============== odabir moda koji promatramo
k_mod = input(' Broj moda (npr. 1, 2, 3, ...) __?  ');
k = k_mod;

% ==============  odredjivanje transverzalnog valnog broja; aproksimativni i/ili egzaktno odredjivanje

aprox = 0;              %  APROKSIMATIVNO rjesenje           tg(ml) = -m/[abs(gam)+b],     za k-ti mod
if(aprox);        
    ml_res = 10000000;          % rezolucija pri racunanju m
    ml_lim = [(2*k-1)*pi/2 (2*k+1)*pi/2]; dml = diff(ml_lim)/ml_res;
    ml = [ml_lim(1)+dml:dml:ml_lim(2)-dml]; m = ml/L_shelf;
    tgml = tan(ml);
    lin = -m/(gam+b);
    [nula inula] = min(abs(tgml-lin));
    m = m(inula);
end                                 % if(aprox);

[mm] =  sub_nultocka(L_shelf, gam, b);                  %   EGZAKTNO rjesenje za niz modova
om = -2*b*gam./(mm.^2+b^2+gam^2)*f;                   % kruzna frekvencija
per = (2*pi./om)*1/(3600*24); 
plot([1:1:length(mm)], per,'*-','LineWidth',1);
set(gca,'XLim',[0 length(mm)+1]); grid on; 
set(gca,'XTick',[0:1:length(mm)+1]);
XLabel('Mod'); YLabel('Period (dani)');
title('Disperzivna relacija','FontWeight','bold');
print('-djpeg90','topo_e_disprel')

if(aprox ~= 0), 
    fprintf('Razlika izmedju aproksimativnog i egzaktnog rjesenja = %e', abs(mm(k)-m));
end

m = mm(k); m0 = mm(1); clear mm                                              % transverzalni valni broj

%kratki = input('\n\n Zelis li rjesenje uz aproksimaciju kratkih valova (DA = 1) _ ? \n')
kratki = 1;
if(kratki == 1), m = k*pi/L_shelf; m0 = pi/L_shelf; end                        % transverzalni valni broj u aproksimaciji kratkih valova
om = -2*b*gam/(m^2+b^2+gam^2)*f;                   % kruzna frekvencija
om0 = -2*b*gam/(m0^2+b^2+gam^2)*f;            % kruzna frekvencija osnovnog moda
v = om/gam;                                                             % fazna brzina
mm = repmat(m,length(x),1);

N_per = 1;                                                             % Broj perioda simulacije     

ndt = 12;                                                              % vremenska rezolucija (dio perioda; veci broj znaci vecu vrem. rezoluciju)
per = 2*pi/abs(om)                                              %  period
per0 = 2*pi/abs(om0)                                          % period osnovnog moda

% ================ odredjivanje strujne funkcije
sl = 0;
%for it = [0:N_per*ndt];                                                % vremenska petlja
for t = [0:per0/ndt:N_per*per];                                                % vremenska petlja
    sl = sl+1;
    omt = t*om;        
	Fi1 = exp(bb.*x);                                                   % Amplituda; transverzalna ovisnost
	Fi2 = sin(mm.*x);
	Fi = Fi1.*Fi2;
    Psi = cos(gamm.*y-omt);                                        % Longitudinalna i vremenska ovisnost
	for ny = 1:length(y); 
        for nx = 1:length(x);
             Strf(ny,nx) = Fi(nx).*Psi(ny);                         % strujna funkcija za transport
        end                   %  for nx = 1:length(x);
	end                       % 	for ny = 1:length(y); 
	dC = (max(Strf(:))-min(Strf(:)))/10;
    
    if(any(isfinite(dC))),
        figure;
    	contourf(x,y,Strf,[min(Strf(:)):dC:max(Strf(:))]);
        set(gca,'DataAspectRatio',[L_shelf, Y_shelf, 1]);
        rr = rats(t/per); rr1 = deblank(rr(end:-1:1)); rr = deblank(rr1(end:-1:1));
        ttl = [num2str(k) ' mod; t = ' rr ' T'];
        title(ttl)
        cb = colorbar('vert');axes(cb); ylbl = ylabel('Transport [m^3/s]'); set(ylbl,'Rotation',270,'VerticalAlignment','bottom')
        film = 1;      % spremanje *.jpg za animaciju
         if(film);
                 ime_slike = ['Mod_' num2str(k) '_k' sprintf('%2.2i',sl)];
                 print ('-djpeg90 ', ime_slike);    
             end
     else     
         if(it == 0),        fprintf('\n Ne generiraju se topografski valovi \n'), end
    end                     %    if(any(isfinite(dC))),   

end                         %   for it = [0:N_per*ndt]; 