% Script to illustrate the thickness estimate method by
%
% Farinotti, D., Huss, M., Bauder, A., Funk, M., & Truffer, M. (2009).
% A method to estimate the ice volume and ice-thickness distribution
% of alpine glaciers. Journal of Glaciology, 55(191), 422-430.
%
% This code calculates the along-flow bed profile of Rhonegletscher
% given 2007 data of glacier outline, surface and flowline coordinates.
% For Rhonegletscher, absence of debris cover is assumed.
%
% Wherever possible without too much extra code, variables are initiated
% with zeros to speed calculation of loops.
clear
close all

% Define file directory ***(adjust accordingly)***
% and data file names.
dpath='./';
fname_outl='glacier_outline_rhone_2007.txt';
fname_flowl='flowline_rhone_2007.txt';
fname_surfdem='surface_dem_rhone_2007.asc';

% Read in data.
[outl header_outl]=read_asc_rhone([dpath,fname_outl],0);
[flowl header_flowl]=read_asc_rhone([dpath,fname_flowl],1);
[surfdem header_surfdem]=read_asc_rhone([dpath,fname_surfdem],6);

% Pick glacier outlines of "connected" catchment.
outl=outl(1:954,:);

% Plot glacier outlines and flowline.
figure(1); hold on; box on
plot(outl(:,1),outl(:,2),'k')
plot(flowl(:,1),flowl(:,2),'r.')
axis equal
xlabel('Easting (m in CH1903/LV03)')
ylabel('Norhthing (m in CH1903/LV03)')

% Calculate curvlinear coordinate, which starts
% at glacier head and grows towards terminus.
flowl_coord=zeros(size(flowl,1),1);
flowl_coord(1)=0;
for i=2:size(flowl,1)
   flowl_coord(i,1)=flowl_coord(i-1,1)+sqrt( (flowl(i,1)-flowl(i-1,1))^2 + (flowl(i,2)-flowl(i-1,2))^2);
end

% Plot glacier bed and surface along
% flow coordinate. Glacier bed is 
% known from radar and considered as
% "ground truth" for our purposes.
figure(2); hold on; box on
plot(flowl_coord,flowl(:,3),'k--')
plot(flowl_coord,flowl(:,4),'k')
legend({'Glacier surface';'Glacier bed (from radar)'})
xlabel('Along-flow coordinate (m)')
ylabel('Elevation above sea level (m)')

% Calculate easting (x) and northing (y)
% coordinate grid and create meshgrid.
x=671237.50:25:671237.50+(167-1)*25;
y=157387.50:25:157387.50+(405-1)*25;
[X,Y] = meshgrid(x,y);
Z=flipud(surfdem);

% Downsample (e.g. by a factor of 10) to
% reduce computational time.
X=X(1:10:end,1:10:end);
Y=Y(1:10:end,1:10:end);
Z=Z(1:10:end,1:10:end);

% Plot glacier in 3D. Add flow line and
% outline (above surface for illustration).
figure(3); hold on; box on
s=surf(X,Y,Z);
plot3(outl(:,1),outl(:,2),4000*ones(size(outl(:,2))),'k*')
plot3(flowl(:,1),flowl(:,2),flowl(:,3),'k.')
xlabel('Easting (m in CH1903/LV03)')
ylabel('Northing (m in CH1903/LV03)')
zlabel('Elevation above sea level (m)')
s.EdgeColor='none';
axis equal

% Reshape mesh grid to vectors.
X_vec=reshape(X,[numel(X) 1]);
Y_vec=reshape(Y,[numel(Y) 1]);
Z_vec=reshape(Z,[numel(Z) 1]);

% Find grid points, which locate within
% glacier boundaries.
in = inpolygon(X_vec,Y_vec,outl(:,1),outl(:,2));
X_vec_in=X_vec(in);
Y_vec_in=Y_vec(in);
Z_vec_in=Z_vec(in);

% Plot grid points within glacier boundary with
% altitude color coded. This avoids creation
% of a new mesh grid.
figure(4); hold on; box on
scatter(X_vec_in,Y_vec_in,10,Z_vec_in,'filled')
axis equal
xlabel('Easting (m in CH1903/LV03)')
ylabel('Northing (m in CH1903/LV03)')
h_cb=colorbar;
set(get(h_cb,'Label'),'String','Elevation above sea level (m)')

% We need to divide the surface into individual 
% elements, for which we can later calculate the
% specific mass balance. Since the surface is regularly
% gridded, this is straightforward. However, in order to
% better fit the glacier surface within the specified outlines,
% we divide the surface into triangles.
TRI = delaunay(X_vec_in,Y_vec_in); % Returns indeces for each triangle.
A=zeros(size(TRI,1),1); %
maxdim=zeros(size(TRI,1),1);
for i=1:length(TRI)
    % Triangle area
    A(i,1) = polyarea(X_vec_in(TRI(i,:)),Y_vec_in(TRI(i,:)));
    % Largest side of triangle
    maxdim(i,1)=max(max(abs([diff(X_vec_in(TRI(i,[1 2 3 1]))) diff(Y_vec_in(TRI(i,[1 2 3 1])))])));
end

% Only keep small triangles.
TRI_small=TRI(A<50000 & maxdim<600,:);
A_small=A(A<50000 & maxdim<600);
maxdim_small=maxdim(A<50000 & maxdim<600);

% Plot Delauney mesh and 
% flow line coordinates.
figure(5); hold on; box on
for i=1:length(TRI_small)
    % Each of the following lines plots one side of the Delauney triangle.
    plot([X_vec_in(TRI_small(i,1)) X_vec_in(TRI_small(i,2))],[Y_vec_in(TRI_small(i,1)) Y_vec_in(TRI_small(i,2))],'k')
    plot([X_vec_in(TRI_small(i,2)) X_vec_in(TRI_small(i,3))],[Y_vec_in(TRI_small(i,2)) Y_vec_in(TRI_small(i,3))],'k')
    plot([X_vec_in(TRI_small(i,3)) X_vec_in(TRI_small(i,1))],[Y_vec_in(TRI_small(i,3)) Y_vec_in(TRI_small(i,1))],'k')
    % Add flow line coordinates on top.
    plot(flowl(:,1),flowl(:,2),'r.')
end
axis equal
xlabel('Easting (m in CH1903/LV03)')
ylabel('Northing (m in CH1903/LV03)')

% For each Delauney triangle calculate a
% center coordinate and elevation. This 
% results in the variables labeled *_f, 
% which constitute the final set of 
% variables to work with.
X_f=zeros(1,size(TRI_small,1));
Y_f=zeros(1,size(TRI_small,1));
Z_f=zeros(1,size(TRI_small,1));
for i=1:length(TRI_small)
    X_f(i)=mean(X_vec_in(TRI_small(i,:)));
    Y_f(i)=mean(Y_vec_in(TRI_small(i,:)));
    Z_f(i)=mean(Z_vec_in(TRI_small(i,:)));
end
A_f=A_small;
    
% *** We need to determine the ELA of the apparent
% mass balance. Z_ELA=2900 is a good choice. Write
% a short piece of code to support this value. 
% Hint: you may use a graphical approach and assume
% the values below for the apparent mass balance 
% gradients gacc and gabl.***
gacc=1.0*0.5*10^(-2);
gabl=1.0*0.9*10^(-2);
Z_ELA=2900;

figure(6); hold on; box on 
title('Placeholder for your figure.')


% Calculate volume flux. This requires integrating
% the apparent mass balance over the contributing
% glacier surface upstream of the coordinate. Identifying
% the contributing surface is difficult to automate. Here we
% simply assume that the contributing surface is at higher
% elevations than the flowline coordinate. *** Can you think
% of a different (better) way? To see how reasonable our
% assumption is, mark the contributing surface for 
% various flow line coordinates (change check_flow). ***
check_flow=18;
figure(7); hold on; box on
scatter(X_f,Y_f,5,Z_f,'filled')
plot(flowl(:,1),flowl(:,2),'r.')
for i=1:length(flowl)
    text(flowl(i,1)+20,flowl(i,2),num2str(i))
end
is_above=Y_f>flowl(check_flow,2);

h_s1=scatter(X_f(is_above),Y_f(is_above),30,'ko');
h_s2=scatter(flowl(check_flow,1),flowl(check_flow,2),30,'go','filled');
xlabel('Easting (m in CH1903/LV03)')
ylabel('Northing (m in CH1903/LV03)')
h_cb=colorbar;
set(get(h_cb,'Label'),'String','Elevation above sea level (m)')
axis equal
h_l=legend([h_s1,h_s2],{'Picked flow line';'Contributing glacier surface'});
set(h_l,'Location','NorthWest')

% Calculate slope angle for each flowline grid point using the
% grid points behind and in front of it.
alpha_deg=zeros([size(flowl,1),1]);
for i=2:length(alpha_deg)-1
    alpha_deg(i,1)=abs(atand( (flowl(i+1,3)-flowl(i-1,3)) / sqrt( (flowl(i+1,1)-flowl(i-1,1))^2 + (flowl(i+1,2)-flowl(i-1,2))^2) ));
end

% Plot surface elevation and slope.
figure(8); hold on; box on
yyaxis left
plot(flowl_coord,flowl(:,3),'k')
set(gca,'YColor','k')
ylabel('Surface elevation (m)')
yyaxis right
plot(flowl_coord,alpha_deg,'r')
xlabel('Along-flow coordinate (m)')
ylabel('Surface slope ({\circ})')
set(gca,'YColor','r')

% Calculate thickness. Note that glacier width w
% was read off the map manually. *** Feel free
% to improve. *** We go only to the 34th flow
% coordinate beyond which there is not glacier.
% *** Change the value of C to investigate the
% sensitivity of the thickness calculation.
w=[0 1057 1533 1587 1788 1980 1980 1903 2036 2205 2732 2874 2507 2421 2027 1328 1447 1119 886 516 654 712 902 511 787 660 660 494 312 332 272 236 130 164 106 0 0 0 0 0 0 0];
A=2.4*10^(-15);
n=3;
C=0.25;
%C=0.38;
%C=1.5*C;
rho=900;
g=9.81;
% Loop over all flow coordinates.
B=zeros(length(flowl),1);
h=zeros(length(flowl),1);
qbar=zeros(length(flowl),1);
for i=1:34%length(flowl)

    % Find surface above flow coordinate elevation.
    is_above=Z_f>flowl(i,3);
    
    % Pick elevation and areas of all surface elements
    % above flow coordinate elevation.
    Z_ab=Z_f(is_above);
    A_ab=A_f(is_above);
    % Integrate mass balance of all surface elements
    % above flow coordinate elevation.
    b=zeros(length(Z_ab),1);
    for j=1:length(Z_ab)
        if Z_ab(j)<Z_ELA
            b(j,1)=(Z_ab(j)-Z_ELA) * gabl * A_ab(j);
        elseif Z_ab(j)>Z_ELA
            b(j,1)=(Z_ab(j)-Z_ELA) * gacc * A_ab(j);
        end
    end
    B(i)=sum(b);
    % Ice volume flux
    qbar(i)=B(i)/w(i);
    % Equation (7) in Farinotti et al. (2009)
    h(i)=( qbar(i)/(2*A) * (n+2)/(C*rho*g*abs(sind(alpha_deg(i))))^n )^(1/(n+2));
    
end

% Plot glacier surface with bed from
% radar and from calculation.
figure(9); hold on; box on
plot(flowl_coord,flowl(:,3),'k--')
plot(flowl_coord,flowl(:,4),'k')
plot(flowl_coord(1:34),flowl(1:34,3)-h(1:34),'co-')
legend({'Glacier surface';'Glacier bed (radar)';'Glacier bed (calculated)'})
xlabel('Along-flow coordinate (m)')
ylabel('Elevation (m)')









