%calculating eddy diffusivity
clear
time_series = h5read('scalers/scalers_s1.h5','/scales/sim_time');
psi1 = h5read('scalers/scalers_s1.h5','/tasks/psi1');
psi2 = h5read('scalers/scalers_s1.h5','/tasks/psi2');
ky_read = h5read('scalers/scalers_s1.h5','/scales/ky');
%ky_total = [ky_read(1:256);2*ky_read(256)-ky_read(255);ky_read(257:end)];
ky_total = [ky_read(1:256);0;ky_read(257:end)];
kx_total = ky_total';

D_time_series = zeros(1000,1);

for j = 1:1000
    barotropic_stream_function = 0.5*(psi1(:,:,j) + psi2(:,:,j));
    barotropic_velocity_hat = 1j*kx_total.*fft2(barotropic_stream_function);
    barotropic_velocity = real(ifft2(barotropic_velocity_hat));
    baroclinic_stream_function = 0.5*(psi1(:,:,j) - psi2(:,:,j));
    D = 2*mean(baroclinic_stream_function.*barotropic_velocity,[1,2]);
    D_time_series(j) = D;
end

    
for i =2:50
    disp(i)
    filename = sprintf("scalers/scalers_s%d.h5",i);
    %{
    eke_up_layer_read = h5read(filename,'/tasks/eke_up_layer');
    eke_up_layer_append = squeeze(eke_up_layer_read);
    eke_bottom_layer_read = h5read(filename,'/tasks/eke_bottom_layer');
    eke_bottom_layer_append = squeeze(eke_bottom_layer_read);
    time_series_append = h5read(filename,'/scales/sim_time');
    eke_up_layer = [eke_up_layer;eke_up_layer_append];
    eke_bottom_layer = [eke_bottom_layer;eke_bottom_layer_append];
    time_series = [time_series;time_series_append];
    %}
    time_series_append = h5read(filename,'/scales/sim_time');
    psi1 = h5read(filename,'/tasks/psi1');
    psi2 = h5read(filename,'/tasks/psi2');
    D_time_series_append = zeros(1000,1);
    for j = 1:1000
        barotropic_stream_function = 0.5*(psi1(:,:,j) + psi2(:,:,j));
        barotropic_velocity_hat = 1j*kx_total.*fft2(barotropic_stream_function);
        barotropic_velocity = real(ifft2(barotropic_velocity_hat));
        baroclinic_stream_function = 0.5*(psi1(:,:,j) - psi2(:,:,j));
        D = 2*mean(baroclinic_stream_function.*barotropic_velocity,[1,2]);
        D_time_series_append(j) = D;
    end
    time_series = [time_series;time_series_append];
    D_time_series = [D_time_series;D_time_series_append];
end
