clear; clc; close all;

%%
a = 1;
D = 2;
Q = [0,1/4/a];

%%
for ii = 1:3
    switch ii
        case 1
            Q(2) = 1/(6*a);
        case 2
            Q(2) = 1/(4*a);
        case 3
            Q(2) = 1/(2*a);
    end
    if false
        f = figure('Position',[50,50,400,900]); hold on;
        xlim(17.5*D*[-1,1]); ylim(80*D*[0,1]);
    else
        f = figure('Position',[50,50,400,500]); hold on;
        xlim(1.5*D*[-1,1]); ylim(5.5*D*[0,1]);
    end
    axis equal;
    set(gca,'Position',[0,0,1,1]);
    grid on;

    %% REFLEKCOTR
    x = linspace(-10,10,1000);
    
    clear PAR;
    PAR(:,1) = x;
    PAR(:,2) = x.^2*a;
    PAR(abs(x)>D,:)=[];

    plot(PAR(:,1),PAR(:,2),'k-','LIneWidth',3);

    %%

    for alpha = -89:2:89
        t = Q(2);
        m = tand(alpha);

        % strahl y = m*x+t
        % Schnittpunkt mit Parabel y = ax²
        % x = (m +- sqrt(m²+4*a*t))/2*a
        x1 = (m + sqrt(m^2+4*a*t))/2/a;
        y1 = a*x1^2;
        if(abs(x1)<D)
            plot([0,x1],[t,y1],'r-');
            plot([0,-x1],[t,y1],'r-');

            % Steigung am Ort x
            m1 = 2*a*x1;
                % plot(x,m1*(x-x1)+y1,'b--')
                % plot(x,-1/m1*(x-x1)+y1,'b--')
                % axis equal; xlim([-2,2]); ylim([0,4])
            % Reflektierter Strahl aus Vektorgleichung mit Normalvektor n und
            % Richtungsvektor k: k' = k - 2*dot(n,k)*n
            k = [1,m];
            n = [-1,1/m1];
            n = n/norm(n);
            k_ = k - 2*dot(n,k)*n;
                % figure(); hold on; grid on; axis equal;
                % plot([0,k(1)],[0,k(2)]);
                % plot([k(1),k(1)+n(1)],[k(2),k(2)+n(2)]);
                % plot([k(1),k(1)+k_(1)],[k(2),k(2)+k_(2)]);
            plot(x1+1000*[0,k_(1)],y1+1000*[0,k_(2)],'r-');
            plot(-x1-1000*[0,k_(1)],y1+1000*[0,k_(2)],'r-');
        else
            plot(Q(1)+1000*[0,1],Q(2)+1000*[0,m],'r-');
            plot(Q(1)-1000*[0,1],Q(2)+1000*[0,m],'r-');
        end
    end
    plot([0,0],[0,1000],'r-');
    plot(Q(1),Q(2),'kx');
    
    print(f,sprintf('PAR_%d_b.png',ii),'-dpng');
    end