%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                            Michael Pokojovy                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

modus = 2; % 2 entspricht besserer Aufloesung

[pde_fig,ax] = pdeinit;
pdetool('appl_cb', 1);
pdetool('snapon', 'on');
set(ax, 'DataAspectRatio', [1 1 1]);
set(ax, 'PlotBoxAspectRatio', [1.5 1 1]);
set(ax, 'XTickMode', 'auto');
set(ax, 'YTickMode', 'auto');
pdetool('gridon', 'on');

if (modus == 1)
    % Geometrie des Gebietes Omega:
    X = [22, 40, 40, 51, 64, 72, 91, 95, 75, 72, 67, 56, 34, 24, 09, 20, 13, 17, 35, 39, 37];
    Y = [54, 45, 42, 34, 36, 24, 21, 13, 12, 07, 09, 23, 31, 32, 29, 36, 42, 45, 32, 34, 42];

    set(ax, 'XLim', [min(X) max(X)]);
    set(ax, 'YLim', [min(Y) max(Y)]);
    
    pdepoly(X, Y, 'Omega');

    % Mengenformel (set formula):
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'), 'String', 'Omega')

    % Randbedingungen : Dirichletsche Randbedingungen
    pdetool('changemode', 0);
    rb = '0';
    for k = 1:21
         pdesetbd(k, 'dir', 1, '1', rb);     
    end
else
    % Geometrie des Gebietes Omega:
    X = [270, 250, 234, 227, 230, 208, 208, 175, 137, 129, 125, 118, 112, ...
         101,  91,  89,  63,  48,  40,  27,  17,  11,   2,  15,  30,  48, ...
          57,  54,  33,  32,  42,  58,  48,  52,  94, 104, 109, 117, 128, ...
         137, 140, 137, 124, 121, 110, 104,  81,  62,  66,  75, 107, 140, ...
         146, 167, 195, 210, 223, 237, 250, 258, 263, 263, 281, 290, 297, ...
         309, 313, 326, 330, 338, 343, 354, 353, 332, 327, 322, 316, 311, ...
         311, 299, 292, 281, 276, 274];
    Y = [  1,   0,  11,  13,  19,  34,  47,  66,  84,  86,  94,  96,  93, ...
          94,  94,  98,  96,  92,  86,  80,  75,  80,  83,  83,  87,  99, ...
         109, 113, 123, 129, 132, 124, 134, 135, 116, 100,  98, 101,  98, ...
          96,  99, 109, 117, 134, 143, 141, 154, 171, 179, 177, 155, 142, ...
         125, 110,  99, 104, 100,  95,  96,  87,  74,  68,  59,  66,  55, ...
          51,  46,  47,  44,  46,  35,  29,  20,  19,  10,  16,  10,  12, ...
          19,   8,   6,  12,  12,   6];
      
    mX = mean(X);
    mY = mean(Y);
    
    X = X - mX;
    Y = Y - mY;
      
    diam = sqrt((X(73) - X(49))^2 + (Y(73) - Y(79))^2);

    
    X = X/diam*63000;
    Y = Y/diam*63000;
    
    set(ax, 'XLim', [min(X) max(X)]);
    set(ax, 'YLim', [min(Y) max(Y)]);

    U = [ 70,  81,  89,  83,  72];
    V = [117, 114, 106, 104, 108];
    
    U = U - mX;
    V = V - mY;
    
    U = U/diam*63000;
    V = V/diam*63000;
    
    pdepoly(X, Y, 'Omega1');
    pdepoly(U, V, 'Omega2');
    % Mengenformel (set formula):
    set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'), 'String', 'Omega1-Omega2')
    
    % Randbedingungen:
    pdetool('changemode', 0);
    rb = '0';
    for k = 1:89
        %pdesetbd(k, 'dir', '1', '0', rb)
        pdesetbd(k, 'neu', '1', '0', rb)
    end
end

% Erstellung der Triangulierung:
refine_count = 2;
setappdata(pde_fig, 'Hgrad', 1.6);
setappdata(pde_fig, 'refinemethod', 'regular');
pdetool('initmesh')
for k = 1:refine_count
    pdetool('refine')
end

u = get(findobj(pde_fig,'Tag','PDEPlotMenu'),'UserData');
h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu');
hp=findobj(get(h,'Children'),'flat','Tag','PDEInitMesh');
he=findobj(get(h,'Children'),'flat','Tag','PDERefine');
ht=findobj(get(h,'Children'),'flat','Tag','PDEMeshParam');
p=get(hp,'UserData'); 
e=get(he,'UserData');
t=get(ht,'UserData');

% Vorgabe der Parameter
T = 40;

N = 500;
tlist = linspace(0, T, N);

b = @pdebound;
d = '1';
a = '0';
f = '0';

c = '1463^2';

u0 = '0';
u1 = '((x - 2.90e4).^2 + (y+0.90e4).^2 <= 10000)*1E6';
u1 = '((x - 3.90e4).^2 + (y+1.20e4).^2 <= 10000)*1E6';

u = hyperbolic(u0, u1, tlist, b, p, e, t, c, a, f, d);

% Plotten
figure(1);

set(gcf, 'PaperUnits', 'centimeters');
xSize = 24; ySize = 14;
xLeft = (21 - xSize)/2; yTop = (30 - ySize)/2;
set(gcf,'PaperPosition', [xLeft yTop xSize ySize]);
set(gcf,'Position',[0 0 xSize*50 ySize*50]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
hold on;

h = title('Bodensee zentriert um sein Baryzentrum', 'interpreter', 'latex');
set(h, 'FontSize', 16);   

xlabel('x_{1}');
ylabel('x_{2}');
axis([min(X) max(X) min(Y) max(Y) -1E5 1E5]);

pdemesh(p, e, t);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);
hold on;

h = title('L\"osung der Wellengleichung $u_{tt} - c^{2} \triangle u = 0$', 'interpreter', 'latex');
set(h, 'FontSize', 16);

k = 0;
for i = 1:3
for j = 1:3
    k = k + 1;
    subplot(3, 3, k);
    
    tt = min(1 + floor((k - 1)*N/8), N);
    
    hold on;
    
    pdeplot(p, e, t, 'xydata', u(:, tt), 'zdata', u(:, tt), 'colormap', 'jet');
    
    caxis([-10000 10000]);
    axis([min(X) max(X) min(Y) max(Y) -1E5 1E5]);
    
    h = title(['$t = $ ' num2str(tlist(tt)) ' [s]'], 'interpreter', 'latex');
    set(h, 'FontSize', 16);
    view(0, 90);
    
    xlabel('x_{1}');
    ylabel('x_{2}');
    zlabel('u(t, x_{1}, x_{2})');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
hold on;

tt = 320;
h = title(['Akustischer Druck zum Zeitpunkt $t = $' num2str(tlist(tt)) ' [s]'], 'interpreter', 'latex');
set(h, 'FontSize', 16);   

pdeplot(p, e, t, 'xydata', u(:, tt), 'zdata', u(:, tt), 'colormap', 'jet');
    
caxis([-10000 10000]);
axis([min(X) max(X) min(Y) max(Y) -1E5 1E5]);

h = title(['Akustischer Druck zum Zeitpunkt $t = $ ' num2str(tlist(tt)) ' [s]'], 'interpreter', 'latex');
set(h, 'FontSize', 16);   
view(0, 90);

xlabel('x_{1}');
ylabel('x_{2}');
zlabel('u(t, x_{1}, x_{2})');