% TOPO_PDNO
% simulacija analitickog rjesenja za TOPOGRAFSKE valove,
% generirane u bazenu oblika rotacijskog paraboloida;
% koordinatni sustav; r - radijalno
%                               th - kutno
% primjer (parametri) za juzni Jadran
%
% crtanje topografije bazena
% crtanje disperzivne relacije
% crtanje u vremenu strujne funkcije transporta, za odredjeni mod (broj cvornih dijametara te broj cvornih kruznica)
%
% Mira Pasaric (mpasaric@irb.hr)
% 01.12.2004.
% =======================================
% POZIV funkcije: hypergeom (Symbolic Math Toolbox)
%
% ------------------------------- fizicki parametri bazena
            R = 157;                   % polumjer bazena (km)
			f = 1e-4;                   % Coriolisov parametar
			
			dR = R/50;                 % prostrorni korak, radijalni
			dTh = 2*pi/50;           % prostrorni korak, kutni
			
			r = [0:dR:R]';              % prostorna mreza
			th = [0:dTh:2*pi];
			
			[Th,rr] = meshgrid(th,r);
			rrR = rr/R;
			[X,Y] = pol2cart(Th,rr); 
            
			% --------------------------- crtanje topografije bazena
crta_topo = 1;            
if(crta_topo)       
            D0 = 100;                                                        % dubina dna (m) u sredistu bazena   
            D = D0*(1-rrR.^2);
            
            vwD = [-30 15];                                              % rotacija 3d slike; azimut i elevacija;  
            surf(X,Y,-D,'FaceColor','interp', 'EdgeColor','none'); daspect([1.5 1.5 1]); axis tight             % 3-d slika dna
            xlb = xlabel('km'); set(xlb,'VerticalAlignment','top','HorizontalAlignment','center');
            ylb = ylabel('km'); set(ylb,'VerticalAlignment','bottom','HorizontalAlignment','center');
            zlabel('Dubina (m)'); 
            title('Bazen kruznog oblika s parabolickim dnom','FontWeight','bold');
            view(vwD);
end            
			% --------------------------- crtanje disperzivne relacije

crta_disp = 1;            
if(crta_disp);            
			kmax = 5;                           % maksimalni broj cvornih dijametara
			jmax = 5;                           % maksimalni broj cvornih kruznica
			k = [1:kmax]';
			j = [0:jmax]';
			
			[J K] = meshgrid(j, k);
            om = f.*K./(2.*(J+2).*(K+J+1)-K);            % kruzna frekvencija    om = om(k,j)
			per = 2*pi./om; per = per/(3600*24);        % period(dani)

            figure; contourf(K, J, per);
            set(gca,'XTick',k,'YTick',j); grid on; set(gca,'layer','top');
            XLabel('Broj cvornih dijametara'); YLabel('Broj cvornih kruznica')
            title('Topografski valovi (parabolicno dno); Disperzivna relacija','FontSize',12);
            cb= colorbar('horiz'); axes(cb); xlabel('Period (dani)');
            clear k j K J om
end			

% ============== odabir moda
 k = input('Modovi: broj cvornih dijametara (n > 1, npr. 1) __ ? \n')
 j = input('Modovi: broj cvornih kruznica (n >= 0, npr. 1) __ ? \n')
n = 2*j+k+6;

alfa = (n+k-4)/2;
beta = (2+k-n)/2;
delta = k+1;
fprintf(' Pricekajte, racuna se! \n')
Z = hypergeom([alfa,beta],delta,rrR.^2);                        % suma Gaussovog hipergeometrijskog reda
Fi = double(Z.*(rrR.^k));                                                 % radijalna ovisnost 

om = f*k/(2*(j+2)*(k+j+1)-k);                                  % kruzna frekvencija
per = 2*pi/om;                                                           % period
fprintf(' Period, T = %4.1f  (dani) \n\n', per/(3600*24));
fprintf(' Brzina sirenja vala duz oboda bazena, v = %4.1f  (m/s) \n\n', om*R*1e3);

rper = 8;                                                                      % vremenski korak; dio perioda         
sl = 0;
for t = [0:per/rper:per];                                               % vremenska petlja
            sl = sl+1;
            omt = om*t;
            Lmb = cos(k*Th-omt);                                    % kutna i vremenska ovisnost
            Psi = Fi.*Lmb;                                                  % strujna funkcija za transport
            
            figure
            contourf(X,Y,Psi,10); hold on;                          % crtanje izolinija strujne funkcije
            set(gca,'XTickLabel',[],'yTickLabel',[],'DataAspectRatio',[1 1 1]); grid on;
            xlm = get(gca,'XLim'); ylm = get(gca,'YLim');

            tt = rats(t/per); tt1 = deblank(tt(end:-1:1)); tt = deblank(tt1(end:-1:1));
            ttl = ['t = ' tt ' T'];
            title(ttl,'FontWeight','bold');            
            txt = ['Broj cvornih dijametara = ' num2str(k) ', broj cvornih kruznica = ' num2str(j)];
            text(xlm(2), ylm(1), txt,'VerticalAlignment','top','HorizontalAlignment','right');

            cb = colorbar('vert');axes(cb); ylbl = ylabel('Transport [m^3/s]'); set(ylbl,'Rotation',270,'VerticalAlignment','bottom')

        film = 0;                                                               % spremanje *.jpg za animaciju
         if(film);
                 ime_slike = ['Topo_p_' num2str(k) '_' num2str(j) '_' sprintf('%2.2i',sl)];
                 print ('-djpeg90 ', ime_slike);    
         end
             
end                         % for t = [0:per/rper:per];
