function moveSun(y,m,d,ut1_vec,time_loc,Phi,Lambda,loc,season)

% Linenstyle
line_style = '-r';

% Initialisierung von Azimut und Zenitwinkel
A_vec = zeros(1, length(ut1_vec));
z_vec = zeros(1, length(ut1_vec));

for i = 1:length(ut1_vec)
    % Berechnung der Zeit nach dem Julianischen Jahrhundert
    ut1 = ut1_vec(i);
    [t,mjd,jd] = julianjh(y,m,d,ut1);

    % Bestimmung der Sonne:
    [alpha,delta,rsun] = solarpos(t);
    % Nutations- und Präzessionswinkel:
    [eps0,deleps,delpsi] = nutwink(t);
    [za,thetaa,zetaa] = prezwink(t);
    % Bestimmung von GAST:
    [gast,gmst] = jul2gast(ut1,t);
    
    % Test:
    % alpha = -90;
    % delta = 90;

    % Nutation:
    RM1 = R1(-eps0-deleps);
    RM3 = R3(-delpsi);
    RM11 = R1(eps0);

    N = RM1*RM3*RM11;
    % N = [1 0 0; 0 1 0; 0 0 1];

    % Präzession:
    RM3 = R3(-za);
    RM2 = R2(thetaa);
    RM31 = R3(-zetaa);

    P = RM3*RM2*RM31;
    % P = [1 0 0; 0 1 0; 0 0 1];

    % Berechnung von GAST
    GAST = 15 * gast;                   % in [deg]
    RGAST = R3(GAST);
    % RGAST = [1 0 0; 0 1 0; 0 0 1];


    % Vektor im inertialen System:
    ri = [cosd(delta)*cosd(alpha); cosd(delta)*sind(alpha); sind(delta)];

    % S1;
    S1 = [-1 0 0;
        0 1 0;
        0 0 1];

    % Transformation vom globalen ins locale Koordinatensystem
    RM2 = R2(90-Phi);
    RM3 = R3(Lambda);

    % Vollständige Transformation
    rg = S1 * RM2 * RM3 * RGAST * N * P * ri;
    % rg = [5;-5;-1];

    % Rückrechnung der lokalen kartesischen Koodinaten zu Azimut und
    % Zenitwinkel:
    norm = sqrt(rg(1)^2+rg(2)^2+rg(3)^2);
    z = acos(rg(3)/norm);
    z = rad2deg(z);
    A = mod(atan2(rg(2),rg(1)), 2*pi);
    A = rad2deg(A);

    A_vec(i) = A;
    z_vec(i) = z;
end

% Benötigen Elevation:
elev = 90 - z_vec;

% Alles in einer Matrix (Zeit, Elevation, Azimut):
All(1,:) = time_loc(:);
All(2,:) = A_vec(:);
All(3,:) = elev(:);

% Skyplot
figure();
hsky = skyplot(A_vec,elev,line_style);
% title(sprintf('Sonnenverlauf über %s im Jahr %d', loc, y))
title(sprintf('Sonnenverlauf am %s im Jahr %d über %s', season, y, loc));

hsky = hsky;

%% Duration Tag bestimmen
sgn = sign(elev);
start_end = find(diff(sgn)~=0);

dday = time_loc(start_end(2) - start_end(1));

fprintf('Die Tageslänge für %s im Jahr %d beträgt: %.2f Stunden\n', season, y, dday);



end


%%RotationsMatrizen:

% R1 um die x-Achse:
function R = R1(theta)
R = [1 0 0;
    0 cosd(theta) sind(theta);
    0 -sind(theta) cosd(theta)];
end

%R2 um die y-Achse
function R = R2(theta)
R = [cosd(theta) 0 -sind(theta);
    0 1 0;
    sind(theta) 0 cosd(theta)];
end

%R3 um die z-Achse
function R = R3(theta)
R = [cosd(theta) sind(theta) 0;
    -sind(theta) cosd(theta) 0;
    0 0 1];
end