 function [out] = interpolate_surface(year,month,type)
%%% Script needs:
% Niet meer gebruiken, plot interpolation is beter !
% oetsettings
% Vlugtenburg scripts folder added to path!


%Read Lines from vluchtenburg.nc
[survey] = extract_surveys(year,month,type);


if 1 %% rotate yes(1) or no(0)
    x0=69141.4936;      %% coordinates northern Jarkus pole (exracted from the 5 Aug 2009 survey)
    y0=447448.7293;
    theta0 =-131.65534968173412+90;  %% angle ref line (fence) in add_surveygrid!

    Xrect = -survey.X_rd.*sind(theta0) + survey.Y_rd.*cosd(theta0);
    Yrect =  -survey.X_rd*cosd(theta0) - survey.Y_rd*sind(theta0);
end


if 1
    figure

    scatter(Xrect(:),Yrect(:),3,survey.altitude(:))
    hold on
    %add_surveygrid
    axis image
    % xlim([6.69e4 6.946e4])
    % ylim([4.459e5 4.484e5])
    xlabel('RD coordinates [10^3 m]');
    ylabel('RD coordinates [10^3 m]');
    colorbar
    title('Survey Grid');
end
%interpolate
if 1
    [out.xi,out.yi] = meshgrid(3.797e5:(3.803e5-3.797e5)/100:3.803e5,...
        2.455e5:(2.460e5-2.455e5)/100:2.460e5);
    out.zi = griddata(Xrect,Yrect,survey.altitude,out.xi,out.yi);

    figure;pcolor(out.xi,out.yi,out.zi)
    shading interp
end
