function katastrophbsp

% Simulation und Visualisierung der Zeeman Katastrophenmaschine
% Mit rechtem Mausklick in Figure 1 Fadenkreuz aktivieren und
% Kontrollparameter (a,b) waehlen

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% 24.10.2014
% Version 2

% Kontrollparameter a,b --> Zahl der Extrema des Potentials
a = -2:0.05:0;
b = sqrt(-4/27*a.^3);

% Abbildungsumgebung oeffnen und platzieren
fh = figure;

% Abbildung Kontrollparameter
ah(1) = subplot(2,1,1);
set(ah(1),'Units','normalized')
plot(b,a),hold on, plot(-b,a)
hpl = plot(0,0,'ro','markersize',8);
ylim([min(a),-min(a)/4])
grid on
xlabel('b'),ylabel('a')

% Abbildung des korrespondierenden Potentials
ah(2) = subplot(2,1,2);

x=-1.5:0.05:1.5;
V = potential(x,0,0);
% figure, plot(x,V)

ax=kuspe;
% Callback Funktion zur graphischen Auswahl der Kontrollparameter
set(fh,'windowbuttondownfcn',{@neuwert,ah(2),hpl,ax});

potplot(x,V,ah(2))



    function v = potential(x,b,a)

    % Berechnung des Potentials
        v = 1/4*x.^4 + 1/2*a.*x.^2 + b.*x;
    end

    function potplot(x,v,hax)

        %Plotten des Potentials
    plot(hax,x,v)
    xlabel(hax,'x'),ylabel(hax,'Vpot')
    end

    function neuwert(hobject,eventdata,hax,hpl,ax)
    % Callbackfunktion
    
    % Auswahl der Kontrollparameter
    [b,a] = ginput(1);

    set(hpl,'Xdata',b,'Ydata',a)

    % Polynom Ableitung des Potentials
    px = [1 0 a b];
    
    % Berechnung und Auswahl der rellen Nullstellen
    nust = roots(px);
    index = find([isreal(nust(1)),isreal(nust(2)),isreal(nust(3))]);
    nust = nust(index);
    % Festlegung der Plotgrenzen
    xmin=min(nust)-0.75;
    xmax=max(nust)+0.75;
    
    % Plotten des Potentials
    x = linspace(xmin,xmax);
    V = potential(x,b,a);
    potplot(x,V,hax)

        %Plotten eines Schnitts (festes u) durch den Kontrollbereich 
    x2=-1.5:0.1:1.5;
    v2=-a*x2 - x2.^3;
    plot(ax,v2,x2)
    xlabel(ax,'b'),ylabel(ax,'x extremal')
    end

    function ax=kuspe
        % Visualisierung der Extrema des Potentials
        ak = linspace(-1.5,0.,25);
        xk = linspace(-1.75,1.75,25);
        [Ak,Xk] = meshgrid(ak,xk);
        Vk = - (Xk.^3 + Ak.*Xk);
        figure, subplot(2,1,1)
        surfc(Vk,Ak,Xk), shading interp, axis tight
        xlabel('b'),ylabel('a'),zlabel('x extremal') 
        ax = subplot(2,1,2);
    end
        

end


