% Calculate difference between distance from center of channel and WSE

clear

max_dist = 27000; % km
max_dist_xs = max_dist/250;

%% Load observation xs and rcl ids

load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'obs_rcl_ids.mat']);
load(['/Users/Colin/Documents/MATLAB/DFlow/influences_chap2/'...
    'obs_xs_ids.mat']);

[obs_num,~] = size(obs_rcl_ids);

%% Names of input flows (180 and 300 cms)

floc_root = ['/Users/Colin/Documents/MATLAB/DFlow/abq_models/'...
    'abq2010add2012bath/'];

f180fp060name = strcat(floc_root,['run3049_180cmsbank060/'...
    'abq10a12b180cmsrst02_0000_20160104_140000_his.nc']);
f300fp060name = strcat(floc_root,['run3050_300cmsbank060/'...
    'abq10a12bf300cmsfp060rst02_0000_20160103_200000_his.nc']);
f180fp040name = strcat(floc_root,['run3083_180cmsbank040/'...
    'abq10a12bf180cmsbank040rst02_0000_20160104_160000_his.nc']);
f300fp040name = strcat(floc_root,['run3084_300cmsbank040/'...
    'abq10a12bf300cmsbank040rst02_0000_20160103_210000_his.nc']);
f180fp080name = strcat(floc_root,['run3085_180cmsbank080/'...
    'abq10a12bf180cmsbank080rst02_0000_20160104_140000_his.nc']);
f300fp080name = strcat(floc_root,['run3086_300cmsbank080/'...
    'abq10a12bf300cmsbank080rst02_0000_20160103_210000_his.nc']);
f180fp100name = strcat(floc_root,['run3087_180cmsbank100/'...
    'abq10a12bf180cmsbank100rst02_0000_20160104_130000_his.nc']);
f300fp100name = strcat(floc_root,['run3088_300cmsbank100/'...
    'abq10a12bf300cmsbank100rst02_0000_20160103_200000_his.nc']);

flow_cell = {f180fp060name f300fp060name f180fp040name f300fp040name ...
    f180fp080name f300fp080name f180fp100name f300fp100name};

%% Load channel center x- and y-coordinates

x_c = ncread(f180fp060name,'station_x_coordinate');
y_c = ncread(f180fp060name,'station_y_coordinate');

for i = 1:max_dist_xs
    xs_center_xy(i,1) = mean(x_c(obs_xs_ids==i&obs_rcl_ids==2,1));
    xs_center_xy(i,2) = mean(y_c(obs_xs_ids==i&obs_rcl_ids==2,1));
end

%% Load flow data for 180 and 300 cms and retrieve

[~,c] = size(flow_cell);
for i = 1:c
    fname_i = flow_cell{1,i};
    % Depth and WSE information for the steady flow run
    wse_i = ncread(fname_i,'waterlevel');
    depth_i = ncread(fname_i,'Waterdepth');
    wse_i_steady = wse_i(:,end);
    depth_i_steady = depth_i(:,end);
    % Calculate mean channel WSE
    for j = 1:max_dist_xs
        xs_chan_z(j,1) = mean(...
            wse_i_steady(obs_rcl_ids(:,1)==2 &...
            obs_xs_ids(:,1)==j & depth_i_steady(:,1)>0,1));
    end
    flow_cell{2,i} = xs_chan_z;
    % Coordinate information
    xc_i = ncread(fname_i,'station_x_coordinate');
    yc_i = ncread(fname_i,'station_y_coordinate');
    a = 0;
    for j = 1:obs_num
        if obs_xs_ids(j,1) <= max_dist_xs && depth_i(j,1) > 0 && ...
                obs_rcl_ids(j,1) == 1
            a = a + 1;
            wse_diff_i(a,1) = obs_xs_ids(j,1);
            wse_diff_i(a,2) = obs_rcl_ids(j,1);
            wse_diff_i(a,3) = wse_i_steady(j,1) - ...
                xs_chan_z(obs_xs_ids(j,1),1);
            wse_diff_i(a,4) = sqrt(...
                (xs_center_xy(obs_xs_ids(j,1),1) - xc_i(j,1))^2 + ...
                (xs_center_xy(obs_xs_ids(j,1),2) - yc_i(j,1))^2);
        elseif obs_xs_ids(j,1) <= max_dist_xs && depth_i(j,1) > 0 && ...
                obs_rcl_ids(j,1) == 3
            a = a + 1;
            wse_diff_i(a,1) = obs_xs_ids(j,1);
            wse_diff_i(a,2) = obs_rcl_ids(j,1);
            wse_diff_i(a,3) = wse_i_steady(j,1) - ...
                xs_chan_z(obs_xs_ids(j,1),1);
            wse_diff_i(a,4) = -sqrt(...
                (xs_center_xy(obs_xs_ids(j,1),1) - xc_i(j,1))^2 + ...
                (xs_center_xy(obs_xs_ids(j,1),2) - yc_i(j,1))^2);
        end
    end
    flow_cell{3,i} = wse_diff_i;
    clear wse_diff_i
end

%% Plot differences
figure(1)
set(gcf,'color','white');
for i = [1 3 5 7]
    wse_diff_i = flow_cell{3,i};
    scatter(wse_diff_i(wse_diff_i(:,1)==30,4),...
        wse_diff_i(wse_diff_i(:,1)==30,3));
    hold on
end
xlabel('Lateral difference');
ylabel('FP - CH elevation');
legend('060', '040', '080', '100');
%ylim([-1 1]);