clear;clc;
filename = 'gaussian_M140_1000Hz_50_90_jmax10000_imax200_G100.h5';
%r0 = 9.1003;
r0 = 0.091;
rmax = 25;
zmax = 800;
M = 140;
%R = 1;
R = 0.1;
jmax = 10000;
imax = 200;
pm = 6.4;
fid = H5F.open(filename);
j = 1;
idx = 1;
eps = 0.1256;
% sigma = 0.02;
sigma = 0.08;

num_of_j = jmax/10;
Lps = zeros(num_of_j,imax);
Lpsz = zeros(num_of_j,imax);
Ws = zeros(num_of_j,imax);
Wsz = zeros(num_of_j,imax);
tau = linspace(-2*pi,2*pi,1000);
p0 = 2e-005;
while j <= jmax-9
    idx
    g_name = strcat('/G/j',num2str(j));
    h_name = strcat('/H/j',num2str(j));
    gid = H5D.open(fid,g_name);
    hid = H5D.open(fid,h_name);
    G = H5D.read(gid);
    H = H5D.read(hid);
    for i=1:imax
        f = 0;
        fz = 0;
        for n=1:M
            f = f + G(n,i)*sin(n*tau)+H(n,i)*cos(n*tau);
            if n<21
                fz = fz + G(n,i)*sin(n*tau)+H(n,i)*cos(n*tau);
            end
        end
        Ws(idx,i) = max(f)*pm;
        Wsz(idx,i) = max(fz)*pm;
        Lps(idx,i) = 20*log10(max((Ws(idx,i)/p0),1));
        Lpsz(idx,i) = 20*log10(max((Wsz(idx,i)/p0),1));
    end
    idx = idx+1;
    j = j+10;
end
Lp = [flip(Lps,2), Lps];
Lpz = [flip(Lpsz,2), Lpsz];
x = linspace(-rmax*R,rmax*R,2*imax);
y = linspace(0,zmax*r0,num_of_j);
%[X,Y] = meshgrid(x,y);
figure(1)

pcolor(x,y,Lp);
ylabel('z [m]')
xlabel('r [m]')
shading interp
colorbar
colormap(hot)
figure(2)
surf(x,y,Lp,'FaceAlpha',0.7)
ylabel('z [m]')
xlabel('r [m]')
zlabel('L_p [dB]')
shading interp
colorbar
figure(3)
pcolor(x,y,Lpz);
ylabel('z [m]')
xlabel('r [m]')
shading interp
colorbar
colormap(hot)
figure(4)
surf(x,y,Lpz,'FaceAlpha',0.7)
ylabel('z [m]')
xlabel('r [m]')
zlabel('L_p [dB]')
shading interp
colorbar