%% Modified from readin_ram_mat.m 
%% Using this code to generate the TL map 

clear
close all

%% Input -Level 1


date_str1 = '11Aug00z';
date_str2 = '11Aug 000Z';

%ssp_dir  = '/projects/plusinp/RAM/Aug07_Z000/ocean/';
%ram_dir  = '/projects/plusinp/RAM/Aug07_Z000/ramout/';
ssp_dir  = '/projects/plusinp/RAM/Aug11_Z000/ocean/';
ram_dir  = '/projects/plusinp/RAM/Aug11_Z000/ram_600hz/ramout/';

%% Input -Level 2
%grid_dir = '/home/jinshan/Matlab/work/plusinp/grid_ver2/';
grid_dir = '/projects/plusinp/RAM/grid/grid_ver5/';


%% Input -Level 3
rcv_max_depth = 500;  % meter

%% Load in the grid: 
load([grid_dir,'acoustic_grid'], 'up_lon','up_lat')

name_i=1;

for ix = 1:size(up_lon,1)
    for iy =1:size(up_lon,2)
        name_idx(name_i)=name_i;
        name_i=name_i+1;
    end
end

%% Load in sound speed sections

%% 
name_i = 78;  % center point
fname = ['PlusINP_',num2str(name_i)];
load([ssp_dir,'sound_',fname,'.mat'],'sections','nsec')

%sections.coordinate =[up_lat(ix,iy) up_lon(ix,iy)]; 
%sections.range = 10.1; % km
%sections.angle = 30;  % degree
%sections.src_lat = sections.coordinate(1);
%sections.src_lon = sections.coordinate(2);
%sections.ang_range = 360;


DEPTH = 5:5:rcv_max_depth;
RANGE = 10:10:round(sections.range)*1000;

z_size = length(DEPTH);
r_size = length(RANGE);

%% TL lon, lat, bearing angle, depth, range

if(nsec == (sections.ang_range/sections.angle))
display('Angle section is right')
else
display('Check Bearing Angle range')
break
end

TL = NaN(size(up_lon,1),size(up_lon,2),nsec,z_size,r_size);

name_i=0;

for ix = 1:size(up_lon,1)
    ix
    for iy =1:size(up_lon,2)
        name_i=name_i+1;
        
        fname = ['PlusINP_',num2str(name_i)];
        
        load([ssp_dir,'sound_',fname,'.mat'],'sections','nsec')
        
        for isec =1:nsec
            
            load([ram_dir,'TL_grid_',num2str(name_i),'ang_',num2str(isec)])
            
            [dep_y,rng_x] = size(tl);
            
            if dep_y > z_size                 
            %temp = NaN(1,1,1,100,700);            
                %TL(ix,iy,isec,1:100,1:rng_x)= tl(1:100,1:rng_x);
                TL(ix,iy,isec,1:z_size,1:r_size)= tl(1:z_size,1:r_size);

            else                
                TL(ix,iy,isec,1:dep_y,1:r_size)= tl(:,1:r_size);
            end
                                    
            
        end
    end
end

%% Plot 1 
Zmin =  0;  % Receiver Depth - minimum 
Zmax = 101; % Receiver Depth - maximum 

TL_threshold_min = 60;  % TL minimum
TL_threshold_max = 70;  % TL maximum 

depth_idx = find(DEPTH>Zmin & DEPTH<Zmax);

TL_no_bearing = squeeze(nanmean(TL,3)); % average along the bearing angles
TL_no_bearing_dep_avg = squeeze(nanmean(TL_no_bearing(:,:,depth_idx,:),3));

TL_dep_avg = squeeze(nanmean(TL(:,:,:,depth_idx,:),4));
TL_dep_avg_no_bearing = squeeze(nanmean(TL_dep_avg,3));

%% Bearing across the shore (40 degree) -> 2, maybe and 3(avg): and (210 degree - 8) 
TL_bearing_0 = squeeze(TL(:,:,2,:,:));

TL_bearing_0_dep_avg = squeeze(nanmean(TL_bearing_0(:,:,depth_idx,:),3));

%% Bearing along the shore (130 degree ( 5 - 120 degree)) and  11 - 300
%% degree
TL_bearing_90 = squeeze(TL(:,:,5,:,:));

TL_bearing_90_dep_avg = squeeze(nanmean(TL_bearing_90(:,:,depth_idx,:),3));

%TL 4
depth_idx_mix = find(DEPTH>15 & DEPTH< 30);
TL_25_depth= squeeze(nanmean(TL(:,:,:,depth_idx_mix,:),4));
TL_no_bearing_dep_avg_mixing_layer = squeeze(nanmean(TL_no_bearing(:,:,depth_idx_mix,:),3));

rng_idx = find(RANGE>5000 & RANGE<7000);

TL_5_intruder = nanmean(TL_no_bearing_dep_avg(:,:,rng_idx),3);
TL_5_intruder_sigma = nanstd(TL_no_bearing_dep_avg(:,:,rng_idx),0,3);

temp = squeeze(nanmean(TL_dep_avg(:,:,:,rng_idx),4));
TL_5_intruder_sigma_bearing = nanstd(temp,0,3);

%temp = squeeze(nanmean(TL_no_bearing(:,:,:,rng_idx),4));
temp = squeeze(nanmean(TL_no_bearing(:,:,depth_idx,rng_idx),4));

TL_5_intruder_sigma_depth = nanstd(temp,0,3);


%%
TL3_range_depth_sigma=NaN(size(up_lon));

for ix = 1:size(up_lon,1)
    % ix
    for iy =1:size(up_lon,2)
%% TL 1        
        idx = find(TL_no_bearing_dep_avg(ix,iy,:)> TL_threshold_min & ...
            TL_no_bearing_dep_avg(ix,iy,:) <TL_threshold_max);
        
        if isempty(idx)
            TL1_range(ix,iy) = NaN;
            TL1_range_sigma(ix,iy)= NaN;
        else
            
        TL1_range(ix,iy) = nanmean(RANGE(idx));
        TL1_range_sigma(ix,iy)= nanstd(RANGE(idx));
        
        end
        
        % Depth sigma
        
        for idep = 1:length(depth_idx)
            idx = find(TL_no_bearing(ix,iy,depth_idx(idep),:)> TL_threshold_min & ...
                TL_no_bearing(ix,iy,depth_idx(idep),:)< TL_threshold_max);
 
            if isempty(idx)
                TL1_range_depth_sigma(ix,iy) = NaN;                
            else
                TL1_range_depth_sigma(ix,iy)= nanstd(RANGE(idx));
            end
        end
        
        % Bearing sigma
        
        for ibearing = 1:nsec
            idx = find(TL_dep_avg(ix,iy,ibearing,:)> TL_threshold_min & ...
                TL_dep_avg(ix,iy,ibearing,:)< TL_threshold_max);
            
            if isempty(idx)
                TL1_range_bearing_sigma(ix,iy) = NaN;                
            else
                TL1_range_bearing_sigma(ix,iy)= nanstd(RANGE(idx));
            end
            
        end
            
%% TL 2
        idx = find(TL_bearing_0_dep_avg(ix,iy,:)> TL_threshold_min & ...
            TL_bearing_0_dep_avg(ix,iy,:) <TL_threshold_max);
        
        if isempty(idx)
            TL2_range(ix,iy) = NaN;
            TL2_range_sigma(ix,iy)= NaN;
        else
            
        TL2_range(ix,iy) = nanmean(RANGE(idx));
        TL2_range_sigma(ix,iy) = nanstd(RANGE(idx));
        end
        
        temp = NaN(length(depth_idx),length(RANGE));

        for idep = 1:length(depth_idx)
            idx = find(TL_bearing_0(ix,iy,depth_idx(idep),:)> TL_threshold_min & ...
                TL_bearing_0(ix,iy,depth_idx(idep),:)< TL_threshold_max);
 
            if isempty(idx)
%                TL2_range_depth_sigma(ix,iy) = NaN;                
            else
                temp(idep,(1:length(idx)))=RANGE(idx);
%                TL2_range_depth_sigma(ix,iy)= nanstd(RANGE(idx));
            end
        end
        
        temp = reshape(temp,1,numel(temp));
        temp(isnan(temp))=[];
        
        TL2_range_depth_sigma(ix,iy)=std(temp);
%        TL2_range(ix,iy) =  mean(temp);        
        
        
%% TL 3
        idx = find(TL_bearing_90_dep_avg(ix,iy,:)> TL_threshold_min & ...
            TL_bearing_90_dep_avg(ix,iy,:) <TL_threshold_max);
        
        if isempty(idx)
            TL3_range(ix,iy) = NaN;
            TL3_range_sigma(ix,iy) = NaN;
        else
            
        TL3_range(ix,iy) = nanmean(RANGE(idx));
        TL3_range_sigma(ix,iy) = nanstd(RANGE(idx));
        
        end
        
        temp = NaN(length(depth_idx),length(RANGE));
        
        for idep = 1:length(depth_idx)
            idx = find(TL_bearing_90(ix,iy,depth_idx(idep),:)> TL_threshold_min & ...
                TL_bearing_90(ix,iy,depth_idx(idep),:)< TL_threshold_max);
 
            if isempty(idx)
%                TL3_range_depth_sigma(ix,iy) = NaN;                
            else
                temp(idep,(1:length(idx)))=RANGE(idx);
%                TL3_range_depth_sigma(ix,iy)= nanstd(RANGE(idx));
            end
        end        
        
         temp = reshape(temp,1,numel(temp));
        temp(isnan(temp))=[];
       TL3_range_depth_sigma(ix,iy)=std(temp);
        
%% TL 4
        idx = find(TL_no_bearing_dep_avg_mixing_layer(ix,iy,:)> TL_threshold_min & ...
            TL_no_bearing_dep_avg_mixing_layer(ix,iy,:) <TL_threshold_max);
        
        if isempty(idx)
            TL4_range(ix,iy) = NaN;
            TL4_range_sigma(ix,iy) = NaN;
        else
            
        TL4_range(ix,iy) = nanmean(RANGE(idx));
        TL4_range_sigma(ix,iy) = nanstd(RANGE(idx));
        
        end
        
        
        temp = NaN(nsec,length(RANGE));
        
        for ibearing = 1:nsec
            idx = find(TL_25_depth(ix,iy,ibearing,:)> TL_threshold_min & ...
                TL_25_depth(ix,iy,ibearing,:)< TL_threshold_max);
            
            if isempty(idx)
%                TL4_range_bearing_sigma(ix,iy) = NaN;                
            else
                
            temp(ibearing,(1:length(idx)))=RANGE(idx);        
                
%                TL4_range_bearing_sigma(ix,iy)= nanstd(RANGE(idx));
            end
            
        end
        
        temp = reshape(temp,1,numel(temp));
        temp(isnan(temp))=[];
        
        TL4_range_bearing_sigma(ix,iy) = std(temp);
        
    end
end


figure 
depth_old  = depth;
load /home/jinshan/Matlab/work/plusinp/kauai_bath.mat
z2D = depth;
depth= depth_old; 
ind=find(z2D>0);
%%
z2D(ind)=NaN;  % needed for later interpolation
V=[0 -5 -10 -50 -500 -1000 -2000 -2500];
[c,h]=contour(lon,lat,z2D,V);
clabel(c,h)

xlim([-159.85 -159.6])
ylim([21.75 22.0])

xbox1=[-159.7561 -159.7967 -159.6583 -159.6044 -159.7561];
ybox1=[  21.9789   21.9383   21.8333   21.8928   21.9789];

xbox2=[-159.7561 -159.855 -159.7192 -159.6044 -159.7561];
ybox2=[  21.9789   21.8725   21.7667   21.8928   21.9789];

%[lat2, lon2]=track(ybox2,xbox2); 
%[lat1, lon1]=track(ybox1,xbox1);
%hold on
%plot(lon2,lat2,'r','linewidth',4)
%plot(lon1,lat1,'k--','linewidth',4)

saveas(gcf,'bath_con','fig')
% shading flat

figure, 
pcolor(up_lon,up_lat,TL1_range/1000), 
clim = caxis;
colorbar
%shading flat
title('Range of Day (km):  TL in range of [60 70] dB, Depth [1 100]m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

open('bath_con.fig')
hold on
pcolor(up_lon,up_lat,TL1_range/1000), 
caxis(clim);
colorbar
%shading flat
% title('Range of Day (km)')
title('Mean Expected Range of Transmission (km) for source at 600 Hz, 40 m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

caxis([1.5 5.2])

saveas(gcf, ['mean_range_trans_plot1_',date_str1],'fig')
saveas(gcf, ['mean_range_trans_plot1_',date_str1],'png')

% % Sigma PLot1
% figure, 
% pcolor(up_lon,up_lat,TL1_range_sigma/1000), 
% clim = caxis;
% close
% 
% open('bath_con.fig')
% hold on
% pcolor(up_lon,up_lat,TL1_range_sigma/1000), 
% caxis(clim);
% colorbar
% title('STD of Expected Range of Transmission (km) for source at 600 Hz, 40 m')
% xlabel(['                     Longitude                  ',date_str2])
% ylabel('Latitude')
% 
% saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma'],'fig')
% saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma'],'png')

% Depth Sigma 
figure, 
pcolor(up_lon,up_lat,TL1_range_depth_sigma/1000), 
clim = caxis;
close

open('bath_con.fig')
hold on
pcolor(up_lon,up_lat,TL1_range_depth_sigma/1000), 
caxis(clim);
colorbar
%shading flat
% title('Range of Day (km)')
title('STD in Depths of Expected Range of Transmission (km) for source at 600 Hz, 40 m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

%caxis([0.75 2.75])
caxis([0.35 3])
saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma_in_depth'],'fig')
saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma_in_depth'],'png')

% Bearing Sigma 
figure, 
pcolor(up_lon,up_lat,TL1_range_bearing_sigma/1000), 
clim = caxis;
close

open('bath_con.fig')
hold on
pcolor(up_lon,up_lat,TL1_range_bearing_sigma/1000), 
caxis(clim);
colorbar
%shading flat
% title('Range of Day (km)')
title('STD in Bearing angles of Expected Range of Transmission (km) for source at 600 Hz, 40 m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

caxis([0 2.95])
saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma_in_bearing'],'fig')
saveas(gcf, ['mean_range_trans_plot1_',date_str1,'_sigma_in_bearing'],'png')

%% PLOT 2
%%

figure, 
pcolor(up_lon,up_lat,TL2_range/1000)
clim = caxis;
close

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL2_range/1000)
caxis(clim);
colorbar
%shading flat
title('Mean Expected Range of Transmission (km), cross-slope, for source at 600 Hz, 40m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')
%caxis([0.5 5.5])
caxis([0.5 6.])
saveas(gcf, ['mean_range_trans_plot2_',date_str1],'fig')
saveas(gcf, ['mean_range_trans_plot2_',date_str1],'png')

% % Sigma of plot2
% figure, 
% pcolor(up_lon,up_lat,TL2_range_sigma/1000)
% clim = caxis;
% close
% 
% open('bath_con.fig')
% hold on
% pcolor(up_lon,up_lat,TL2_range_sigma/1000)
% caxis(clim);
% colorbar
% title('STD of Expected Range of Transmission (km), cross-slope, for source at 600 Hz, 40m')
% xlabel(['                     Longitude                  ',date_str2])
% ylabel('Latitude')
% saveas(gcf, ['mean_range_trans_plot2_',date_str1,'_sigma'],'fig')
% saveas(gcf, ['mean_range_trans_plot2_',date_str1,'_sigma'],'png')

%% Sigma in Depth of plot2
figure, 
pcolor(up_lon,up_lat,TL2_range_depth_sigma/1000)
clim = caxis;
close

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL2_range_depth_sigma/1000)
caxis(clim);
colorbar

title('STD in Depths of Expected Range of Transmission (km), cross-slope, for source at 600 Hz, 40m')
xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

caxis([0.35 3])

saveas(gcf, ['mean_range_trans_plot2_',date_str1,'_sigma_in_depth'],'fig')
saveas(gcf, ['mean_range_trans_plot2_',date_str1,'_sigma_in_depth'],'png')

%% Plot 3

figure, 
pcolor(up_lon,up_lat,TL3_range/1000)
clim = caxis;
close

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL3_range/1000)
caxis(clim);
colorbar
%shading flat
title('Mean Expected Range of Transmission (km), along-slope, for source at 600 Hz,40m ')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')
caxis([0.5 6.])
saveas(gcf, ['mean_range_trans_plot3_',date_str1],'fig')
saveas(gcf, ['mean_range_trans_plot3_',date_str1],'png')

% % Sigma
% figure, 
% pcolor(up_lon,up_lat,TL3_range_sigma/1000)
% clim = caxis;
% close
% 
% open('bath_con.fig')
% hold on
% 
% pcolor(up_lon,up_lat,TL3_range_sigma/1000)
% caxis(clim);
% colorbar
% title('STD of Expected Range of Transmission (km), along-slope, for source at 600 Hz,40m ')
% 
% xlabel(['                     Longitude                  ',date_str2])
% ylabel('Latitude')
% saveas(gcf, ['mean_range_trans_plot3_',date_str1,'_sigma'],'fig')
% saveas(gcf, ['mean_range_trans_plot3_',date_str1,'_sigma'],'png')

% Sigma in Depth of PLOT3
figure, 
pcolor(up_lon,up_lat,TL3_range_depth_sigma/1000)
clim = caxis;
close

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL3_range_depth_sigma/1000)
caxis(clim);
colorbar
%shading flat
title('STD in Depths of Expected Range of Transmission (km), along-slope, for source at 600 Hz,40m ')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')
%caxis([0.35 3])
caxis([1.0 3])
saveas(gcf, ['mean_range_trans_plot3_',date_str1,'_sigma_in_depth'],'fig')
saveas(gcf, ['mean_range_trans_plot3_',date_str1,'_sigma_in_depth'],'png')

%% Plot 4
figure, 
pcolor(up_lon,up_lat,TL4_range/1000)
clim = caxis;
close 

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL4_range/1000)
caxis(clim);
colorbar
%shading flat
title('Mean Expected Range of Transmission (km) for source at 600 Hz,40m,receiver at 20 m ')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

caxis([1.5 5.2])

saveas(gcf, ['mean_range_trans_plot4_',date_str1],'fig')
saveas(gcf, ['mean_range_trans_plot4_',date_str1],'png')

% %Sigma
% figure, 
% pcolor(up_lon,up_lat,TL4_range_sigma/1000)
% clim = caxis;
% close 
% 
% open('bath_con.fig')
% hold on
% 
% pcolor(up_lon,up_lat,TL4_range_sigma/1000)
% caxis(clim);
% colorbar
% title('STD of Expected Range of Transmission (km) for source at 600 Hz,40m,receiver at 20 m ')
% 
% xlabel(['                     Longitude                  ',date_str2])
% ylabel('Latitude')
% saveas(gcf, ['mean_range_trans_plot4_',date_str1,'_sigma'],'fig')
% saveas(gcf, ['mean_range_trans_plot4_',date_str1,'_sigma'],'png')

% Sigma in Bearing 
figure, 
pcolor(up_lon,up_lat,TL4_range_bearing_sigma/1000)
clim = caxis;
close 

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL4_range_bearing_sigma/1000)
caxis(clim);
colorbar
title('STD in Bearing angles of Expected Range of Transmission (km) for source at 600 Hz,40m,receiver at 20 m ')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')
%caxis([0 2.95])
caxis([1.5 2.75])
saveas(gcf, ['mean_range_trans_plot4_',date_str1,'_sigma_in_bearing'],'fig')
saveas(gcf, ['mean_range_trans_plot4_',date_str1,'_sigma_in_bearing'],'png')

%% Plot 5


figure, 
pcolor(up_lon,up_lat,-TL_5_intruder)
clim = caxis;
close 

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,-TL_5_intruder)
caxis(clim);
colorbar
%shading flat
title('Mean TL(dB) for source at 600 Hz, 40m (mean over bearings, 5-7km ranges and 1-100 m depths )')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

caxis([-76 -66])

saveas(gcf, ['TL_5_7km_plot5_',date_str1],'fig')
saveas(gcf, ['TL_5_7km_plot5_',date_str1],'png')

% % Sigma plot 5
% figure, 
% pcolor(up_lon,up_lat,-TL_5_intruder_sigma)
% clim = caxis;
% close 
% 
% open('bath_con.fig')
% hold on
% 
% pcolor(up_lon,up_lat,-TL_5_intruder_sigma)
% caxis(clim);
% colorbar
% title('STD of TL(dB) for source at 600 Hz, 40m (mean over bearings, 5-7km ranges and 1-100 m depths )')
% 
% xlabel(['                     Longitude                  ',date_str2])
% ylabel('Latitude')
% 
% saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma'],'fig')
% saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma'],'png')
% 

% Sigma in depth of plot 5
figure, 
pcolor(up_lon,up_lat,TL_5_intruder_sigma_depth)
clim = caxis;
close 

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL_5_intruder_sigma_depth)
caxis(clim);
colorbar
%shading flat
title('STD in Depths of TL(dB) for source at 600 Hz, 40m (mean over bearings, 5-7km ranges and 1-100 m depths )')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')

%caxis([-35 -0.05])
caxis([0.01 4.5])
saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma_in_depth'],'fig')
saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma_in_depth'],'png')


% Sigma in bearing of plot 5
figure, 
pcolor(up_lon,up_lat,TL_5_intruder_sigma_bearing)
clim = caxis;
close 

open('bath_con.fig')
hold on

pcolor(up_lon,up_lat,TL_5_intruder_sigma_bearing)
caxis(clim);
colorbar
%shading flat
title('STD in Bearing Angles of TL(dB) for source at 600 Hz, 40m (mean over bearings, 5-7km ranges and 1-100 m depths )')

xlabel(['                     Longitude                  ',date_str2])
ylabel('Latitude')
%caxis([-12 -0.5])
caxis([0.5 20])
saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma_in_bearing'],'fig')
saveas(gcf, ['TL_5_7km_plot5_',date_str1,'_sigma_in_bearing'],'png')




