% created by Nils Hase 2016 (nilshase@math.uni-bremen.de)

% Plot results and show information specially adapted for the Barnett
% testcase

function show_barnett(methodname,f_opt,alpha_opt,b_index)

load_sol = ['f',strrep(lower(methodname(3:end)),' ','_')];

% Some tricky translation from methodname to savename
if strcmp(load_sol,'frnett') == 1;
    load_sol = 'f_EDGAR';
elseif strcmp(load_sol,'fandard_inversion')
    load_sol = 'f_st_inv';
elseif strcmp(load_sol,'fansform_inversion')
    load_sol = 'f_tr_inv';
elseif strcmp(load_sol,'fgrange_multiplier')
    load_sol = 'f_lm';
elseif strcmp(load_sol,'fbbs_sampler')  
    load_sol = 'f_gibbs';
end

load('synthetic_methane_data.mat','area')
load('inverse_methods_solutions.mat',load_sol)
eval(['f_old = ',load_sol,';'])

barnett_total_est = sum((f_opt(b_index) - f_old(b_index)).*area(b_index));
total_dif = sum((f_opt - f_old).*area);

fprintf('Warning: assumed loss is 1.5 percent!\n')
fprintf('   %s\n',methodname)
fprintf('   %2.2f %% of total barnett emissions\n',barnett_total_est*100/(1.5*438.66));
fprintf('   barnett total: %3.4f\n',barnett_total_est)
fprintf('   total dif: %3.4f\n',total_dif)
fprintf('   chosen regularization parameter: %3.4f\n',alpha_opt)

% Show the change in surface flux induced by additional sources in the
% Barnett (see below)
map_result_barnett(f_opt - f_old,methodname);

end




% function to create a map of surface fluxes specially for the Barnett
% testcase

function h = map_result_barnett(f,title_string)

if nargin == 2;
    root_plot = 1;
elseif nargin == 1;
    root_plot = 1;
    title_string = '';
end

% Dimensions for US methane setup
dim_x = 94; dim_y = 31;

lon = (-145:-52) + 0.5;
lat = (25:55) + 0.5;

% Specify a colormap and linecolor
[my_colormap,linecolor] = get_colors();
c_min = -0.011/4; c_max = 0.08/4;
my_colormap = my_colormap([2,3:2:end],:);

% Load shapes for borders and land grids
load('shapes.mat')

% Make a scaled image plot
my_pos = [500, 300, 600/2, 330/2];

h = figure('position',my_pos);
hold on
f_temp = zeros(dim_y,dim_x);
f_temp(land_grids) = f;
imagesc(lon,lat,nthroot(f_temp,root_plot));
set(gca,'YDir','normal');
for k = 1:length(US_STATE.ncst)
    plot(US_STATE.ncst{k}(:,1),US_STATE.ncst{k}(:,2),'color',linecolor,'linewidth',1);
end
for k = 1:length(COASTLINE.ncst)
    lon_points = COASTLINE.ncst{k}(:,1);
    lat_points = COASTLINE.ncst{k}(:,2);
    if sum(lon_points > lon(1)) > 0 && ...
       sum(lon_points < lon(end)) > 0 && ...
       sum(lat_points > lat(1)) > 0 && ...
       sum(lat_points < lat(end))
   
        plot(lon_points,lat_points,'color',linecolor,'linewidth',0.5);
    end
end

plot([-99,-99,-96,-96,-99],[32,34,34,32,32],'r','linewidth',2)
hold off
xlim([-106,-89])
ylim([28,36])
set(gca,'xtick',[-105:5:-89])
set(gca,'ytick',[29:2:35])
colormap(my_colormap);
caxis([c_min,c_max])
title(title_string)
box on
drawnow
fprintf('\n');

end
