clc;
clear;
close all;

data = readtable('eopc01.1846-now.csv');

mjd = data.MJD;
xp = data.x_pole;
yp = data.y_pole;

t0 = mjd(1);
t(:) = mjd(:) - t0;

% 2D plot
figure();
plot(xp,yp,'k-');

xlabel('x-Achse in ["]');
ylabel('y-Achse in ["]');
title('Polkoordinaten');

grid on;

% 3D-Plot
figure();
plot3(xp,yp,t,'k-');

xlabel('x-Achse in ["]');
ylabel('y-Achse in ["]');
zlabel('Time in [days] nach 1846')
title('Polkoordinaten');

grid on;

% Fouriertransformation
dt = mean(diff(t));
fs = 1/dt;                              % [1/Tag]

% Transformation der Polarkoordinaten
N = length(xp);

% FFT
XP = fft(xp - mean(xp));
YP = fft(yp - mean(yp));

% Frequenzachse:
f = (0:N-1)*(fs/N);

% Amplitudenspektrum (nur positive Frequenzen)
X_amp = abs(XP(1:floor(N/2)));
Y_amp = abs(YP(1:floor(N/2)));
f_pos = f(1:floor(N/2));
peri = 1./f_pos;


% Fouriertransformation mit Frequenz
figure;
plot(f_pos, X_amp, 'r', 'DisplayName', 'x\_pole');
hold on;
plot(f_pos, Y_amp, 'b', 'DisplayName', 'y\_pole');
xlabel('Periode [Tag]');
ylabel('Amplitude');
legend;
title('Fouriertransformierte der Polbewegung');
grid on;
xlim([0, 450]);  % z. B. nur niedrige Frequenzen zeigen

% Fouriertransformation plotten mit Periode
peri = 1 ./ f_pos;
figure;
plot(peri, X_amp, 'r', peri, Y_amp, 'b');
xlabel('Periode [Tage]');
ylabel('Amplitude');
set(gca, 'XDir','reverse');  % damit große Perioden links stehen
grid on;
xlim([100, 1000]);  % z. B. nur Chandler und Jahreswellen
title('Spektrum der Polbewegung nach FFT');
legend('x\_pole','y\_pole');


%% Polpewegungen analysieren mittels der Lomb-Scargle Periodogramm 
% Matlab benötigt für PWelch eine gleichmäßige Abtastrate. Wenn diese nicht vorliegt werden die Ergebnisse physikalisch verzerrt.
% Das liegt hier vor. PWelch benutzt FFT, die die Zeit in gleichmäßige
% Intervalle einteilt.

% Die benötigten Variablen:
% xp    -   x_pol
% yp    -   y_pol
% t     -   Zeit ab 1845

XP = xp - mean(xp);
YP = yp - mean(yp);

% Lomb-Scargle-Spektrum
[amp_x, f] = plomb(xp, t);
amp_y = plomb(yp, t);

% Umrechnung in Perioden
peri = 1 ./ f;

% Plot
figure;
plot(peri, amp_x, 'r', 'DisplayName', 'xp');
hold on;
plot(peri, amp_y, 'b', 'DisplayName', 'yp');
set(gca, 'XDir', 'reverse');
xlabel('Periode [Tage]');
ylabel('Leistungsdichte');
legend;
title('Spektralanalyse der Polbewegung mit Lomb-Scargle');
xlim([100, 1000]);  % Chandler und Jahresgang sichtbar
grid on;


% Berechnung der Hauptperiode
% Vorsortierung der Daten, im Bereich zwischen 100 und 1000
filt = (peri >= 100) & (peri <= 1000);
amp_x = amp_x(filt);
amp_y = amp_y(filt);
peri_lim = peri(filt);

[peaks_x,loc_x] = findpeaks(amp_x, 'SortStr', 'descend');
[peaks_y,loc_y] = findpeaks(amp_y, 'SortStr', 'descend');

% Auswahl der höchsten Peaks
top_peri_x = peri_lim(loc_x(1:3));      % Umrechnen in Perioden
top_amp_x = peaks_x(1:3);             % Amplituden der Peaks
top_peri_y = peri_lim(loc_y(1:3));      % Umrechnen in Perioden
top_amp_y = peaks_y(1:3);             % Amplituden der Peaks

fprintf("Top 3 Perioden im Bereich 100–1000 Tage:\n");
fprintf("Für den x-Pol:\n")
for i = 1:min(3, length(peaks))
    fprintf("  %d: Periode = %.2f Tage, Amplitude = %.4f\n", ...
        i, top_peri_x(i), peaks_x(i));
end
fprintf("Für den y-Pol:\n")
for i = 1:min(3, length(peaks))
    fprintf("  %d: Periode = %.2f Tage, Amplitude = %.4f\n", ...
        i, top_peri_y(i), peaks_y(i));
end



%% Sekundäre Polbewegung
% Zeit in Jahre umwandeln:
t_years = t/365.25;

% Linearität
p_x = polyfit(t_years, xp, 1);
p_y = polyfit(t_years, yp, 1);
px_trend = polyval(p_x, t_years);
py_trend = polyval(p_y, t_years);

% Geschwindigkeit und Richtung
vx = p_x(1);
vy = p_y(1);

vel = hypot(vx, vy);
dir = atan2(vy, vx) * 180/pi;

fprintf("Driftgeschwindigkeit: %.2f mas/Tag\n", vel);
fprintf("Richtung der Bewegung: %.1f° (Nord = 0°)\n", dir);

% Visualisieren
figure;
quiver(0, 0, vx, vy, 0, 'k', 'LineWidth', 2, 'MaxHeadSize', 2);
hold on;
plot(0, 0, 'ko', 'MarkerFaceColor', 'k');  % Ursprung
plot(xp, yp, '.', 'DisplayName', 'Polbewegung');
axis equal;
grid on;
xlabel('Trend in x_p [mas/Jahr]');
ylabel('Trend in y_p [mas/Jahr]');
title('Säkulare Polbewegung als Trendvektor');
xlim([-1 1]*abs(vx)*2);  % automatische Skalierung
ylim([-1 1]*abs(vy)*2);

% Linearität zeigen
figure;
subplot(2,1,1);
plot(t, xp, '.', t, px_trend, '-', 'LineWidth', 2);
ylabel('x_p [mas]'); grid on;
legend('Messung','Trend');
title('Lineare Drift in x_p');

subplot(2,1,2);
plot(t, yp, '.', t, py_trend, '-', 'LineWidth', 2);
ylabel('y_p [mas]'); xlabel('Jahr'); grid on;
legend('Messung','Trend');
title('Lineare Drift in y_p');

% 3D-Plot
figure;
hold on;
grid on;
box on;

% Originaldaten (Messpunkte)
plot3(xp, yp, t_years, '.', 'DisplayName', 'Messung');

% Trendlinie (als Gerade in 3D)
plot3(px_trend, py_trend, t_years, 'r-', 'LineWidth', 2, 'DisplayName', 'Trend');

zlabel('Zeit [Jahre]');
xlabel('x_p [mas]');
ylabel('y_p [mas]');
title('3D-Visualisierung der Polbewegung mit Trend');
legend;
view(45, 20);  % Ansicht perspektivisch (frei anpassbar)




