%In this file, there is 1 line that should be changed at every time
clear
clc

%change the next line
date_='8_24_06';
load ([date_, '_TL_forecast.mat'])
Start=[1 8 14 20];%[1 20]%this is to divide realizations in 3 parts
jj=1;
for bearing=1:8
    for f=[400, 100]
        for division=1:3%1:1
            figure(jj)
            k=1;
            for sz=[5,40,80]
                x=eval(['TL_forecast.B', int2str(bearing), '.f', num2str(f), '.sz', num2str(5), '.TL_realizations']);
                dr=eval(['TL_forecast.B', int2str(bearing), '.f', num2str(f), '.sz', num2str(5), '.dr']);
                ndr=eval(['TL_forecast.B', int2str(bearing), '.f', num2str(f), '.sz', num2str(5), '.ndr']);
                dz=eval(['TL_forecast.B', int2str(bearing), '.f', num2str(f), '.sz', num2str(5), '.dz']);
                ndz=eval(['TL_forecast.B', int2str(bearing), '.f', num2str(f), '.sz', num2str(5), '.ndz']);
                [m p n]=size(x);
                rz_i=round(linspace(2,p,3));
                for i=1:3 %rz
                    subplot(3,3,k)
                    plot([1:n]*dr*ndr,squeeze(x(Start(division):Start(division+1)-1,rz_i(i),:)))
                    view(0,270), grid on
                    
                    title(['D', int2str(division),' B', int2str(bearing), ' f=', num2str(f), 'Hz sz=', num2str(sz), 'm rz=', num2str((rz_i(i)-1)*dz*ndz), 'm']);
                    k=k+1;
                end
            end
            jj=jj+1;
            file_name=[date_, '_D', int2str(division), '_B', int2str(bearing), 'f', num2str(f)];
            print('-dpng', file_name)            
        end
    end
    
    water_grid_r=eval(['SVP_PE.B', int2str(bearing), '.water_grid_r']);
    water_grid_z=eval(['SVP_PE.B', int2str(bearing), '.water_grid_z']);
    water_c_rl_1=eval(['SVP_PE.B', int2str(bearing), '.SVP_PE']);
    rb=eval(['SVP_PE.B', int2str(bearing), '.rb']);
    zb=eval(['SVP_PE.B', int2str(bearing), '.zb']);
    cb=eval(['SVP_PE.B', int2str(bearing), '.cb']);
    seabed_speed_grid_z=eval(['SVP_PE.B', int2str(bearing), '.seabed_speed_grid_z']);
    rhob=eval(['SVP_PE.B', int2str(bearing), '.rhob']);
    seabed_density_grid_z=eval(['SVP_PE.B', int2str(bearing), '.seabed_density_grid_z']);
    attn=eval(['SVP_PE.B', int2str(bearing), '.attn']);
    seabed_attn_grid_z=eval(['SVP_PE.B', int2str(bearing), '.seabed_attn_grid_z']);
    
    figure(jj)
    subplot(2,2,1)
    contourf(water_grid_r /1000, water_grid_z, water_c_rl_1, 70),colorbar('vert');
    hold on;
    shading flat, colormap jet, caxis([1488 1507]), colorbar;
    
    area(rb/1000, zb,1500,'FaceColor',[.5 .5 .5]);
    view(0,270)
    xlabel('(km)'); ylabel('(m)'); title('Principal Estimate of Sound Velocity Profile (m/s)')
    
    subplot(2,2,2)   
    [m n]=size(cb);
    r=[];z=[];x=[];
    for j=1:n
        for i=1:m
            r=[r, water_grid_r(j)];
            z=[z, zb(j)+seabed_speed_grid_z(i)];
            x=[x, cb(i,j)];
        end
    end
    zi=[0, min(zb):1.2*water_grid_z(length(water_grid_z))];
    [RI ZI]=meshgrid(water_grid_r, zi);
    XI=griddata(r, z, x, RI, ZI);
    pcolor(RI/1000, ZI, XI);
    view(0,270), shading flat, colormap jet, caxis([min(min(cb))*0.9 max(max(cb))*1.1]), colorbar, hold on;
    area(rb/1000, zb, 0,'FaceColor',[1 1 1]);
    xlabel('(km)'); ylabel('(m)'); title('Seabed Sound Velocity Profile (m/s)')    
    
    subplot(2,2,3)
    [m n]=size(rhob);
    r=[];z=[];x=[];
    for j=1:n
        for i=1:m
            r=[r, water_grid_r(j)];
            z=[z, zb(j)+seabed_density_grid_z(i)];
            x=[x, rhob(i,j)];
        end
    end
    zi=[0, min(zb):1.2*water_grid_z(length(water_grid_z))];
    [RI ZI]=meshgrid(water_grid_r, zi);
    XI=griddata(r, z, x, RI, ZI);
    pcolor(RI/1000, ZI, XI);
    view(0,270), shading flat, colormap jet, caxis([min(min(rhob))*0.9 max(max(rhob))*1.1]), colorbar, hold on;
    area(rb/1000, zb, 0,'FaceColor',[1 1 1]);
    xlabel('(km)'); ylabel('(m)'); title('Seabed Density Profile (g/cc)') 
    
    subplot(2,2,4)
    [m n]=size(attn);
    r=[];z=[];x=[];
    for j=1:n
        for i=1:m
            r=[r, water_grid_r(j)];
            z=[z, zb(j)+seabed_attn_grid_z(i)];
            x=[x, attn(i,j)];
        end
    end
    zi=[0, min(zb):1.2*water_grid_z(length(water_grid_z))];
    [RI ZI]=meshgrid(water_grid_r, zi);
    XI=griddata(r, z, x, RI, ZI);
    pcolor(RI/1000, ZI, XI);
    view(0,270), shading flat, colormap jet, caxis([min(min(attn))*0.9 max(max(attn(1:m-1,:)))*1.1]), colorbar, hold on;
    area(rb/1000, zb, 0,'FaceColor',[1 1 1]);
    xlabel('(km)'); ylabel('(m)'); title('Seabed Attenuation Coefficient Profile (dB/wavelength)') 
    
    file_name=[date_, '_B', int2str(bearing), '_SVP_PE'];
    print('-dpng', file_name)
    jj=jj+1;
end