%% omni_read_2months.m
% OMNI2 data reading: January–February, solar wind parameters + Kp only

clear; clc; close all;

%% === 1) File ===
fname = "datanl.txt";

%% === 2) Read table ===
opts = detectImportOptions(fname, 'FileType','text', 'CommentStyle','#');
T = readtable(fname, opts);

% Time (OMNI2 format)
timeStr = string(T{:,1}) + " " + string(T{:,2});
t = datetime(timeStr, 'InputFormat','dd-MM-yyyy HH:mm:ss', ...
             'Format','dd-MM-yyyy HH:mm');

%% === 3) Time interval: 2 months starting from January ===
t_start = datetime(2025,1,1);
t_end   = datetime(2025,2,28,23,59,59);

mask_time = (t >= t_start) & (t <= t_end);
t = t(mask_time);

%% === 4) Select ONLY required parameters ===
Bmag = T{mask_time, 3};      % |B| (nT)
Bz   = T{mask_time, 6};      % Bz (nT)
Nsw  = T{mask_time, 8};      % plasma density (cm^-3)
Vsw  = T{mask_time, 9};      % solar wind speed (km/s)
Pdyn = T{mask_time,10};      % dynamic pressure (nPa)
Kp   = T{mask_time,11} / 10; % Kp index (scaled)

% figure;
% scatter(Bmag, Bz, 8, 'filled');
% grid on;
% set(gcf,'Color','w','MenuBar','none','ToolBar','none');
% xlabel('|B| (nT)');
% ylabel('Bz (nT)');
% title('IMF: |B| vs Bz');

%% === DATA CLEANING ===

bad = ...
    Bmag <= 0   | Bmag > 100   | ...   % |B| (nT)
    abs(Bz) > 50               | ...   % Bz (nT)
    Vsw  < 200  | Vsw  > 2000  | ...   % solar wind speed (km/s)
    Pdyn <= 0   | Pdyn > 50    | ...   % dynamic pressure (nPa)
    Nsw  <= 0   | Nsw  > 200   | ...   % plasma density (cm^-3)
    isnan(Bmag) | isnan(Bz) | isnan(Vsw) | isnan(Pdyn) | isnan(Nsw);

% Apply mask
t    = t(~bad);
Bmag = Bmag(~bad);
Bz   = Bz(~bad);
Vsw  = Vsw(~bad);
Nsw  = Nsw(~bad);
Pdyn = Pdyn(~bad);
Kp   = Kp(~bad);

%% === Data for clustering (Bz, V) ===
X  = [Bz, Vsw];   % 2D feature space
Xz = zscore(X);

%% === 6) PLOTS ===
figure(1);
scatter(Bz, Vsw, 8, 'filled');
grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('Solar wind states distribution (Bz vs V)');

figure(2);
scatter(Pdyn, Vsw, 8, 'filled');
grid on;
set(gcf,'Color','w');
xlabel('P_{dyn} (nPa)');
ylabel('V (km/s)');
title('Solar wind states distribution');

figure(3);
scatter(Bmag, Bz, 8, 'filled');
grid on;
set(gcf,'Color','w');
xlabel('|B| (nT)');
ylabel('Bz (nT)');
title('IMF: |B| vs Bz');

figure(4);
scatter(Bz, Vsw, 10, Kp, 'filled');
grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('Bz vs V with Kp color coding');
cb = colorbar;
cb.Label.String = 'Kp';

%% === 7) Save prepared data ===
save("omni_JanFeb_raw.mat", ...
     "t","Bmag","Bz","Vsw","Nsw","Pdyn","Kp");

disp("Saved: omni_JanFeb_raw.mat");

%% === K-means clustering (K = 3) ===
K = 3;
rng(1);   % reproducibility

[idx, C] = kmeans(Xz, K, ...
    'Replicates', 30, ...
    'MaxIter', 1000);

figure(5);
gscatter(Bz, Vsw, idx, 'bry', '.', 10);
grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('K-means clustering of solar wind states (Bz vs V, K = 3)');
legend('Cluster 1','Cluster 2','Cluster 3','Location','best');

fprintf('Median Kp by clusters:\n');
for k = 1:K
    fprintf('Cluster %d: median Kp = %.2f, N = %d\n', ...
        k, median(Kp(idx==k)), sum(idx==k));
end
%% === Cluster centers in physical units (K = 3) ===
mu = mean(X, 1);
sigma = std(X, 0, 1);

C_phys = C .* sigma + mu;   % centers in (Bz, V)
meanKp3 = arrayfun(@(k) mean(Kp(idx==k), 'omitnan'), 1:K);
figure;
hold on;

colors = lines(K);
for k = 1:K
    scatter(Bz(idx==k), Vsw(idx==k), 10, colors(k,:), 'filled');
end

% cluster centers as filled circles
scatter(C_phys(:,1), C_phys(:,2), 120, colors, ...
        'filled', 'MarkerEdgeColor','k');

grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('K-means clustering of solar wind states (Bz vs V, K = 3)');

legend(arrayfun(@(k) ...
    sprintf('Cluster %d (mean Kp = %.2f)', k, meanKp3(k)), ...
    1:K, 'UniformOutput', false), ...
    'Location','best');



%% === K-means clustering (K = 5) ===
% Assumes you already have:
%   Bz, Vsw  (cleaned)
%   Xz = zscore([Bz, Vsw])

K = 5;
rng(1);  % reproducibility


[idx5, C5] = kmeans(Xz, K, ...
    'Replicates', 50, ...
    'MaxIter', 2000);

figure(7);
gscatter(Bz, Vsw, idx5, lines(K), '.', 10);
grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('K-means clustering of solar wind states (Bz vs V, K = 5)');
legend(arrayfun(@(k) sprintf('Cluster %d', k), 1:K, 'UniformOutput', false), ...
       'Location','best');
C5_phys = C5 .* sigma + mu;
meanKp5 = arrayfun(@(k) mean(Kp(idx5==k), 'omitnan'), 1:5);
figure(8);
hold on;

colors = lines(5);
for k = 1:5
    scatter(Bz(idx5==k), Vsw(idx5==k), 10, colors(k,:), 'filled');
end

scatter(C5_phys(:,1), C5_phys(:,2), 120, colors, ...
        'filled', 'MarkerEdgeColor','k');

grid on;
set(gcf,'Color','w');
xlabel('Bz (nT)');
ylabel('V (km/s)');
title('K-means clustering of solar wind states (Bz vs V, K = 5)');

legend(arrayfun(@(k) ...
    sprintf('Cluster %d (mean Kp = %.2f)', k, meanKp5(k)), ...
    1:5, 'UniformOutput', false), ...
    'Location','best');

hold off;

%% === Cluster sizes for K = 5 (counts + percentages) ===
K5 = 5;

counts5 = accumarray(idx5, 1, [K5 1]);
perc5 = 100 * counts5 / sum(counts5);

figure(9);
bar(1:K5, counts5);
grid on;
set(gcf,'Color','w');
xlabel('Cluster number');
ylabel('Number of points');
title('Number of data points in each cluster (K = 5)');

% add percentage labels above bars
ymax = max(counts5);
for k = 1:K5
    text(k, counts5(k) + 0.02*ymax, sprintf('%.1f%%', perc5(k)), ...
        'HorizontalAlignment','center', 'FontWeight','bold');
end

%% === Elbow method (within-cluster sum of squares vs K) ===
Kmax = 10;
wcss = zeros(Kmax,1);

rng(1); % reproducibility
for K = 1:Kmax
    % Sum of distances within clusters (kmeans returns it as 'sumd')
    [~, ~, sumd] = kmeans(Xz, K, 'Replicates', 30, 'MaxIter', 2000);
    wcss(K) = sum(sumd);
end

figure(10);
plot(1:Kmax, wcss, '-o', 'LineWidth',1.5, 'MarkerSize',6);
grid on;
set(gcf,'Color','w');
xlabel('Number of clusters (K)');
ylabel('Total within-cluster sum of squares');
title('Elbow method for selecting K');



%% === 11) Export cleaned data to NetCDF (time series) ===

out_nc = "omni_JanFeb_cleaned.nc";

% If file exists, delete to avoid conflicts
if isfile(out_nc)
    delete(out_nc);
end

% --- Time encoding (CF-style): hours since Unix epoch ---
t0 = datetime(1970,1,1,0,0,0,'TimeZone','UTC');
t_utc = t; 
t_utc.TimeZone = 'UTC';
time_hours = hours(t_utc - t0);   % numeric vector

nT = numel(time_hours);

%% Create dimension "time"
nccreate(out_nc, 'time', ...
    'Dimensions', {'time', nT}, ...
    'Datatype', 'double');

ncwrite(out_nc, 'time', time_hours);

ncwriteatt(out_nc, 'time', 'standard_name', 'time');
ncwriteatt(out_nc, 'time', 'long_name', 'Time');
ncwriteatt(out_nc, 'time', 'units', 'hours since 1970-01-01 00:00:00 UTC');
ncwriteatt(out_nc, 'time', 'calendar', 'gregorian');

%% Helper to create 1D time variables
create_1d = @(vname, dtype) nccreate(out_nc, vname, ...
    'Dimensions', {'time', nT}, 'Datatype', dtype);

%% Main solar wind / IMF variables
create_1d('Bmag', 'single');  ncwrite(out_nc,'Bmag', single(Bmag));
ncwriteatt(out_nc,'Bmag','long_name','IMF magnitude');
ncwriteatt(out_nc,'Bmag','units','nT');

create_1d('Bz', 'single');    ncwrite(out_nc,'Bz', single(Bz));
ncwriteatt(out_nc,'Bz','long_name','IMF south-north component (GSM/GSE as provided by OMNI)');
ncwriteatt(out_nc,'Bz','units','nT');

create_1d('Vsw', 'single');   ncwrite(out_nc,'Vsw', single(Vsw));
ncwriteatt(out_nc,'Vsw','long_name','Solar wind speed');
ncwriteatt(out_nc,'Vsw','units','km s-1');

create_1d('Nsw', 'single');   ncwrite(out_nc,'Nsw', single(Nsw));
ncwriteatt(out_nc,'Nsw','long_name','Solar wind proton number density');
ncwriteatt(out_nc,'Nsw','units','cm-3');

create_1d('Pdyn', 'single');  ncwrite(out_nc,'Pdyn', single(Pdyn));
ncwriteatt(out_nc,'Pdyn','long_name','Solar wind dynamic pressure');
ncwriteatt(out_nc,'Pdyn','units','nPa');

create_1d('Kp', 'single');    ncwrite(out_nc,'Kp', single(Kp));
ncwriteatt(out_nc,'Kp','long_name','Planetary geomagnetic index Kp');
ncwriteatt(out_nc,'Kp','units','1');
ncwriteatt(out_nc,'Kp','comment','Kp provided as standard 0..9 scale');

%% Cluster labels (optional but useful)
% K=3 labels
if exist('idx','var')
    create_1d('cluster_K3', 'int16');
    ncwrite(out_nc,'cluster_K3', int16(idx));
    ncwriteatt(out_nc,'cluster_K3','long_name','K-means cluster label in (Bz,V) space, K=3');
    ncwriteatt(out_nc,'cluster_K3','units','1');
end

% K=5 labels
if exist('idx5','var')
    create_1d('cluster_K5', 'int16');
    ncwrite(out_nc,'cluster_K5', int16(idx5));
    ncwriteatt(out_nc,'cluster_K5','long_name','K-means cluster label in (Bz,V) space, K=5');
    ncwriteatt(out_nc,'cluster_K5','units','1');
end

%% === Global attributes (publication-style metadata) ===
% --- REQUIRED/COMMON CF-style metadata ---
ncwriteatt(out_nc, '/', 'Conventions', 'CF-1.10');
ncwriteatt(out_nc, '/', 'title', 'Cleaned OMNI hourly solar wind parameters and Kp (Jan–Feb 2025)');
ncwriteatt(out_nc, '/', 'summary', ['Hourly OMNI-derived solar wind/IMF parameters with Kp index for Jan–Feb 2025. ' ...
    'Data were cleaned by removing missing and non-physical values and prepared for exploratory clustering analysis in the Bz–V space.']);
ncwriteatt(out_nc, '/', 'keywords', 'solar wind, IMF, OMNI, geomagnetic activity, Kp, clustering, space weather');
ncwriteatt(out_nc, '/', 'keywords_vocabulary', 'NASA/GCMD Science Keywords (suggested)');

% --- Provenance / source ---
ncwriteatt(out_nc, '/', 'source', 'NASA OMNI Hourly Data, OMNIWeb (GSFC)');
ncwriteatt(out_nc, '/', 'references', 'https://omniweb.gsfc.nasa.gov/');
ncwriteatt(out_nc, '/', 'processing_level', 'Derived subset + cleaning + optional k-means labels');
ncwriteatt(out_nc, '/', 'history', sprintf('Created on %s (local time) using MATLAB.', datestr(datetime('now'))));

% --- Creator / contact  ---
ncwriteatt(out_nc, '/', 'creator_name', 'Daniil Dementjev');
ncwriteatt(out_nc, '/', 'creator_institution', 'TalTech');
ncwriteatt(out_nc, '/', 'creator_email', 'dademe@taltech.ee');
ncwriteatt(out_nc, '/', 'project', 'Solar wind regime clustering and relation to Kp');
ncwriteatt(out_nc, '/', 'acknowledgment', 'NASA OMNIWeb and contributing space physics data providers.');

% --- Repository / identifier (EDIT THIS) ---
% Put your university repository landing page URL or DOI if assigned later
ncwriteatt(out_nc, '/', 'identifier', 'To be added...');

% --- License / usage ---
ncwriteatt(out_nc, '/', 'license', 'CC BY 4.0 (recommended) or follow university repository requirements');
ncwriteatt(out_nc, '/', 'comment', ['Time coordinate is given as hours since 1970-01-01 00:00:00 UTC. ' ...
    'Cluster labels (cluster_K3, cluster_K5) are included if computed and represent k-means results in normalized (Bz,V) space.']);

% --- Time coverage (CF/ACDD-style) ---
t_utc = t; t_utc.TimeZone = 'UTC';
ncwriteatt(out_nc, '/', 'time_coverage_start', datestr(min(t_utc), 'yyyy-mm-ddTHH:MM:SSZ'));
ncwriteatt(out_nc, '/', 'time_coverage_end',   datestr(max(t_utc), 'yyyy-mm-ddTHH:MM:SSZ'));

% --- Geospatial coverage (optional; OMNI is near-Earth / L1 context) ---
ncwriteatt(out_nc, '/', 'geospatial_bounds', 'Near-Earth solar wind conditions (OMNI); location varies by upstream monitors mapped to Earth');
ncwriteatt(out_nc, '/', 'geospatial_bounds_crs', 'EPSG:4326');

disp("Saved NetCDF: " + out_nc);

