%
% @file fit2d.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
file='fit2d.h5';
%
%     Get data from data sets
%
r=hdf5read(file,'/xpt');
t=hdf5read(file,'/ypt');

fcalc=hdf5read(file,'/fcalc');
fexact=hdf5read(file,'/fexact');
errs=hdf5read(file,'/errs');
%
%     Attributes
%
NR=hdf5read(file,'/NX'); NTH=hdf5read(file,'/NY');
NIDBAS1=hdf5read(file,'/NIDBAS1');
NIDBAS2=hdf5read(file,'/NIDBAS2');
MBESS=hdf5read(file,'/MBESS');
LABEL=sprintf('nr = %d, ntheta = %d, nidbas = (%d,%d), mbess = %d', NR, NTH, ...
              NIDBAS1, NIDBAS2, MBESS);

attr=hdf5read(file,'/xpt/title'); x_ann=attr.Data;
attr=hdf5read(file,'/ypt/title'); y_ann=attr.Data;
attr=hdf5read(file,'/fcalc/title'); fcalc_ann=attr.Data;
attr=hdf5read(file,'/fexact/title'); fexact_ann=attr.Data;
attr=hdf5read(file,'/errs/title');errs_ann=attr.Data;

[R,T]=meshgrid(r,t);
x = R.*cos(T); y= R.*sin(T);

figure
subplot(221)
pcolor(double(x),double(y),double(fcalc));
shading interp
xlabel('X'); ylabel('Y')
title(LABEL)
colorbar

subplot(222)
pcolor(double(x),double(y),double(fexact))
shading interp
axis image
xlabel('X'); ylabel('Y')
title('X-Y plane')
colorbar

subplot(223)
surfc(double(x),double(y),double(errs))
xlabel('X'); ylabel('Y');
title(errs_ann)

%% Plot error at theta ~ pi/4
k = max(find(t<pi/4));
err=fcalc(k,:)-fexact(k,:);
figure
plot(r, err)