function writeMDA(mda_filename_tmp, baseline_ldb, varargin)
%writeMDA: Writes a unibest mda-file based on a baseline (ldb). Can also
%          include a shoreline (ldb), it will then automatically determine
%          the cross-shore (perpendicular) distance between the baseline
%          and shoreline. You can also specify a minimum grid resolution
%          or cut-up your baseline with a certain resolution. Finally, the
%          output can automatically be visualized using a keyword:
%
%  Syntax:
%     writeMDA(mda_filename, baseline_ldb, <keyword,value>);
% 
%  Input:
%    mda_filename       Required: String with output filename of mda-file
%    baseline_ldb       Required: String with a filename of polygon/ldb for 
%                                 the baseline (reference line) OR a [Nx2]
%                                 matrix containing the X,Y coordinates of
%                                 the baseline (can also be [2xN] matrix).
%    <coastline_ldb>    Optional: String with a filename of polygon/ldb for 
%                                 the coastline OR a [Nx2] matrix 
%                                 containing the X,Y coordinates of the 
%                                 coastline (can also be [2xN] matrix). You
%                                 are advised to not let N be too large
%                                 (N > 10000) as the function polyintersect
%                                 requires some time to be computed. Use
%                                 ldbTool/RGFGRID/QUICKIN to change the ldb.
%                                 Default: coastline_ldb = baseline_ldb.
%    <min_resolution>   Optional: Single numerical value [1x1] specifying
%                                 the minimal resolution of the final
%                                 baseline grid as used in Unibest.
%                                 Default: min_resolution = [] (baseline
%                                 resolution is used).
%    <baseline_dx>      Optional: Single numerical value [1x1] specifying
%                                 the minimal distance of baseline
%                                 sections. Cuts-up the baseline in finer
%                                 equidistant sections, but only if 
%                                 required (does not coarsen it).
%                                 Default: baseline_dx = [] (original
%                                 baseline resolution is used).
%    <plot_figure>      Optional: Boolean (true/false, 1/0) to switch on
%                                 function output plotting to a figure. Can
%                                 be used to check the function and easily
%                                 detect unintended behavior. Old figures
%                                 generated by this functions are closed.
%                                 Default: plot_figure = 0 (no figure)
%    <save_to_mda>      Optional: Boolean (true/false, 1/0) to switch off
%                                 writing the output to the *.mda file.
%                                 This can be used to first inspect the
%                                 behaviour of the function. When setting
%                                 save_to_mda to false (0), plot_figure is
%                                 automatically set to true (1).
%                                 Default: save_to_mda = 1 (saving to file)
%    <check_other_side> Optional: Boolean (true/false, 1/0) to switch off
%                                 the search for a coastline on the
%                                 negative (y < 0) side of the baseline as
%                                 well. When switched on, if no coastline
%                                 is found for y >= 0, it will also search
%                                 for y < 0. This is particulary useful for
%                                 a case where the coastline is located
%                                 behind the baseline (see example 2).
%                                 Default: check_other_side = 1 (negative
%                                 y values are allowed)
%
%  Call writeMDA() to get an overview of the default <keyword,value> pairs 
% 
%  Output:
%     *.mda file as specified by mda_filename (or figure if save_to_mda = 0)
%
%  Example 1:
%     baseline  = [0:100:1000; 0:100:1000]';
%     coastline = [-500:1000; 500+100*sin(0.01*[-500:1000])+[-500:1000]]';
%     writeMDA('test.MDA',baseline,...
%              'coastline_ldb',coastline,...
%              'baseline_dx',50,...
%              'min_resolution',25,...
%              'plot_figure',true,...
%              'save_to_mda',false);
%
%   See also readMDA landboundary add_equidist_points polyintersect

%% Copyright notice
%   --------------------------------------------------------------------
%
%   Copyright (C) 2014 Deltares
%       Freek Scheel
%       freek.scheel@deltares.nl	
%
%       Deltares
%       Rotterdamseweg 185
%       PO Box Postbus 177
%       2600MH Delft
%       The Netherlands
%
%   This library 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 2.1 of the License, or (at your option) any later version.
%
%   This library 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
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and 
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute 
% your own tools.

    %% Check and handle the function input:

    keywords.coastline_ldb    = [];
    keywords.min_resolution   = [];
    keywords.baseline_dx      = [];
    keywords.plot_figure      = 0;
    keywords.save_to_mda      = 1;
    keywords.check_other_side = 1;
    keywords.interp_missing_ldb = 1;

    if nargin == 0
        disp('Default keywords:');
        disp(' ');
        disp(keywords)
        return
    end

    if odd(length(varargin)) % Check for input consistency
        error('Odd number of input parameters for <keyword,value> pairs supplied, please check this');
    end

    % Set user defined changes to default keywords:

    keywords = setproperty(keywords,varargin);

    % baseline_ldb_tmp checks:
    if ~isempty(baseline_ldb)
        if isstruct(baseline_ldb)
            baseline_ldb = cell2mat(baseline_ldb);
        end
        if isstr(baseline_ldb)
            baseline_ldb = landboundary('read',baseline_ldb);
        end
        if isnumeric(baseline_ldb)
            if length(size(baseline_ldb)) == 2
                ldb_size = size(baseline_ldb);
                if ~isempty(find(ldb_size==2))
                    if length(find(ldb_size==2))==1
                        if ldb_size(1) == 2
                            baseline_ldb = baseline_ldb';
                        else
                            % Just a normal ldb is found, good
                        end
                    else
                        % [2x2] is found, bit small, but alright, we'll use it..
                    end
                else
                    error(['Unexpected format of baseline ldb, should be [Nx2] or [2xN], it is [' num2str(size(baseline_ldb)) ']']);
                end
            else
                error(['Unexpected dimension of baseline ldb, should be [Nx2] or [2xN], it is [' num2str(size(baseline_ldb)) ']']);
            end
            % Check for NaN's in the baseline:
            if sum(sum(isnan(baseline_ldb)))>0
                % Remove all preceding and trailing NaN's if they exist
                while sum(isnan(baseline_ldb(end,:)))>0
                    baseline_ldb(end,:) = [];
                end
                while sum(isnan(baseline_ldb(1,:)))>0
                    baseline_ldb(1,:) = [];
                end
                if sum(sum(isnan(baseline_ldb)))>0
                    plot(baseline_ldb(:,1),baseline_ldb(:,2),'k'); hold on; axis equal; grid on; title('Your baseline ldb including NaN connected points (red circles)');
                    plot(baseline_ldb([find(sum(isnan(baseline_ldb),2)./sum(isnan(baseline_ldb),2)==1)+1; find(sum(isnan(baseline_ldb),2)./sum(isnan(baseline_ldb),2)==1)-1],1),baseline_ldb([find(sum(isnan(baseline_ldb),2)./sum(isnan(baseline_ldb),2)==1)+1; find(sum(isnan(baseline_ldb),2)./sum(isnan(baseline_ldb),2)==1)-1],2),'ro');
                    error('NaN''s are found cutting-up your baseline (see figure), this is not allowed, please use the ldbTool or RGFGRID/QUICKIN to change this');
                end
            end
        else
            error(['Unexpected baseline ldb input, should be a [Nx2] or [2xN] vector or ldb filename']);
        end
    else
        error('empty baseline_ldb specified, please check this')
    end
    baseline     = baseline_ldb;
    baseline_ori = [];

    % mda_filename_tmp checks:
    if ~isempty(mda_filename_tmp)
        if isstruct(mda_filename_tmp)
            if length(mda_filename_tmp)==1
                mda_filename_tmp = cell2mat(mda_filename_tmp);
            else
                error('too many mda_filename names specified, please check this')
            end
        end
        if isstr(mda_filename_tmp)
            if length(mda_filename_tmp)<5
                mda_filename_tmp = [mda_filename_tmp '.MDA'];
            else
                if ~(strcmp(mda_filename_tmp(1,end-3:end),'.mda') | strcmp(mda_filename_tmp(1,end-3:end),'.MDA'))
                    mda_filename_tmp = [mda_filename_tmp '.MDA'];
                end
            end
        else
            error('mda_filename should be specified as a string, please check this')
        end
    else
        error('empty mda_filename specified, please check this')
    end
    mda_filename = mda_filename_tmp;

    % baseline_dx checks:
    if ~isempty(keywords.baseline_dx)
        if isstruct(keywords.baseline_dx)
            keywords.baseline_dx = cell2mat(keywords.baseline_dx);
        end
        if ~isnumeric(keywords.baseline_dx)
            error('Please specify a number for baseline_dx');
        else
            if length(keywords.baseline_dx)>1
                error('Please specify one number for baseline_dx');
            end
        end
    end

    % min_resolution checks:
    if ~isempty(keywords.min_resolution)
        if isstruct(keywords.min_resolution)
            keywords.min_resolution = cell2mat(keywords.min_resolution);
        end
        if ~isnumeric(keywords.min_resolution)
            error('Please specify a number for min_resolution');
        else
            if length(keywords.min_resolution)>1
                error('Please specify one number for min_resolution');
            end
        end
    end

    % plot_figure checks:
    if ~isempty(keywords.plot_figure)
        if isstruct(keywords.plot_figure)
            keywords.plot_figure = cell2mat(keywords.plot_figure);
        end
        if isnumeric(keywords.plot_figure) | islogical(keywords.plot_figure)
            if length(keywords.plot_figure)>1
                error('Please specify true/false (1/0) for plot_figure');
            else
                if keywords.plot_figure ~= 0
                    keywords.plot_figure = 1;
                else
                    keywords.plot_figure = 0;
                end
            end
        else
            error('Please specify true/false (1/0) for plot_figure');
        end
    else
        error('Please specify true/false (1/0) for plot_figure');
    end
   
    % save_to_mda checks:
    if ~isempty(keywords.save_to_mda)
        if isstruct(keywords.save_to_mda)
            keywords.save_to_mda = cell2mat(keywords.save_to_mda);
        end
        if isnumeric(keywords.save_to_mda) | islogical(keywords.save_to_mda)
            if length(keywords.save_to_mda)>1
                error('Please specify true/false (1/0) for save_to_mda');
            else
                if keywords.save_to_mda ~= 0
                    keywords.save_to_mda = 1;
                else
                    keywords.save_to_mda = 0;
                    if keywords.plot_figure == 0
                        disp('plot_figure is set to true, as MDA saving is set to false');
                    end
                    % keywords.save_to_mda = 0 so set figure plotting on:
                    keywords.plot_figure = 1;
                end
            end
        else
            error('Please specify true/false (1/0) for save_to_mda');
        end
    else
        error('Please specify true/false (1/0) for save_to_mda');
    end

    % coastline_ldb checks:
    if ~isempty(keywords.coastline_ldb)
        coast_is_baseline = 0;
        if isstruct(keywords.coastline_ldb)
            keywords.plot_figure = cell2mat(keywords.coastline_ldb);
        end
        if isstr(keywords.coastline_ldb)
            keywords.coastline_ldb = landboundary('read',keywords.coastline_ldb);
        end
        if isnumeric(keywords.coastline_ldb)
            if length(size(keywords.coastline_ldb)) == 2
                ldb_size = size(keywords.coastline_ldb);
                if ~isempty(find(ldb_size==2))
                    if length(find(ldb_size==2))==1
                        if ldb_size(1) == 2
                            keywords.coastline_ldb = keywords.coastline_ldb'
                        else
                            % Just a normal ldb is found, good
                        end
                    else
                        % [2x2] is found, bit small, but alright, we'll use it..
                    end
                else
                    error(['Unexpected format of coastline ldb, should be [Nx2] or [2xN], it is [' num2str(size(keywords.coastline_ldb)) ']']);
                end
            else
                error(['Unexpected dimension of coastline ldb, should be [Nx2] or [2xN], it is [' num2str(size(keywords.coastline_ldb)) ']']);
            end
        end
        coastline = keywords.coastline_ldb;
    else
        coast_is_baseline = 1;
        coastline = baseline;
    end

    %% Do some data manipulation

    if ~isempty(keywords.baseline_dx)
        baseline_ori = baseline;
        baseline     = add_equidist_points(keywords.baseline_dx,baseline,'equi');
        baseline     = baseline(2:end-1,:);
    end

    % loop through points of the baseline and find for each point the
    % perpendicular distance to the coastline if coast ~= baseline
    leg_tel = 0;
    if coast_is_baseline
        y = zeros(size(baseline,1),1);
        if keywords.plot_figure
            close(findobj('type','figure','tag','writeMDA')); % closes old figures made by this function, handy when looping
            fig = figure; set(fig,'tag','writeMDA'); hold on; grid on; box on; axis equal;
            set(fig,'outerPosition',get(0,'ScreenSize').*[0 0 1 1],'color','w');
            xlabel('X-coordinate'); ylabel('Y-coordinate');
            title('This baseline equals the coastline position');
            if ~isempty(baseline_ori)
                leg_tel = leg_tel + 1;
                leg_text{leg_tel} = ['Original baseline points'];
                l(leg_tel) = plot(baseline_ori(:,1),baseline_ori(:,2),'x','linewidth',2,'markersize',14,'color','k'); % [139 69 19]/255
                leg_tel  = leg_tel + 1;
                leg_text{leg_tel} = ['Added baseline points (''baseline_dx'')'];
                l(leg_tel) = plot(baseline(:,1),baseline(:,2),'.','linewidth',2,'markersize',16,'color','b');
                plot(baseline(:,1),baseline(:,2),'-','linewidth',2,'markersize',16,'color','k');
            else
                leg_tel = leg_tel + 1;
                leg_text{leg_tel} = ['Baseline points (equals coastline position)'];
                l(leg_tel) = plot(baseline(:,1),baseline(:,2),'x-','linewidth',2,'markersize',14,'color','k');
            end
            try
                legend(l,leg_text,'interpreter','none','location','SouthOutside','AutoUpdate','off');
            catch
                legend(l,leg_text,'interpreter','none','location','SouthOutside');
            end
        end
    else
        max_dist_to_coast = ceil(sqrt(sum((diff([min([min(coastline(:,1)) min(baseline(:,1))]) min([min(coastline(:,2)) min(baseline(:,2))]); max([max(coastline(:,1)) max(baseline(:,1))]) max([max(coastline(:,2)) max(baseline(:,2))])])).^2)));
        if size(coastline,1)>10000 & size(baseline,1)>5
            disp(' ');
            disp(['You''ve supplied quite a large coastline landboundary ([Nx2] with N=' num2str(size(coastline,1)) ')']);
            disp('Determining the cross-shore distance from the baseline to the coastline may take some time');
            disp('You are advised to cut up the coastline to parts really needed for this analysis (if possible, see the help)');
        end
        disp(' ');
        fprintf(1,'\nComputing cross-sectional distances from baseline to coastline:      ')
        try
            for ii=1:size(baseline,1)
                fprintf(1,[repmat('\b',1,(length(num2str(100*((ii-1)/size(baseline,1)),'%9.0f'))+2)) num2str(100*(ii/size(baseline,1)),'%9.0f') ' %%']);
                if ii==1
                    baseline_angle(ii,1) = mod((xy2degN(baseline(ii,1),baseline(ii,2),baseline(ii+1,1),baseline(ii+1,2)) - 90),360);
                elseif ii==size(baseline,1)
                    baseline_angle(ii,1) = mod((xy2degN(baseline(ii-1,1),baseline(ii-1,2),baseline(ii,1),baseline(ii,2)) - 90),360);
                else
                    baseline_angle(ii,1) = mod((xy2degN(baseline(ii-1,1),baseline(ii-1,2),baseline(ii+1,1),baseline(ii+1,2)) - 90),360);
                end
                [X_crs Y_crs] = polyintersect([baseline(ii,1) baseline(ii,1)+(sind(baseline_angle(ii,1))*max_dist_to_coast)],[baseline(ii,2) baseline(ii,2)+(cosd(baseline_angle(ii,1))*max_dist_to_coast)],coastline(:,1),coastline(:,2));
                flip_factor   = 1;
                if keywords.plot_figure
                    if ii == 1
                        close(findobj('type','figure','tag','writeMDA')); % closes old figures made by this function, handy when looping
                        fig = figure; set(fig,'tag','writeMDA'); hold on; grid on; box on; axis equal;
                        set(fig,'outerPosition',get(0,'ScreenSize').*[0 0 1 1],'color','w');
                        xlabel('X-coordinate'); ylabel('Y-coordinate');
                        ht=title(['Overview of actions undertaken by writeMDA for ' mda_filename]);set(ht,'interpreter','none');
                        if ~isempty(baseline_ori)
                            leg_tel = leg_tel + 1;
                            leg_text{leg_tel} = ['Original baseline points'];
                            l(leg_tel) = plot(baseline_ori(:,1),baseline_ori(:,2),'x','linewidth',2,'markersize',14,'color','k'); % [139 69 19]/255
                            leg_tel  = leg_tel + 1;
                            leg_text{leg_tel} = ['Final baseline points (''baseline_dx'')'];
                            l(leg_tel) = plot(baseline(:,1),baseline(:,2),'.','linewidth',2,'markersize',16,'color','b');
                            plot(baseline(:,1),baseline(:,2),'-','linewidth',2,'markersize',16,'color','k');
                        else
                            baseline_ori = baseline;
                            leg_tel = leg_tel + 1;
                            leg_text{leg_tel} = ['Baseline points'];
                            l(leg_tel) = plot(baseline(:,1),baseline(:,2),'x-','linewidth',2,'markersize',14,'color','k');
                        end
                        leg_tel = leg_tel + 1;
                        leg_text{leg_tel} = ['Coastline'];
                        l(leg_tel) = plot(coastline(:,1),coastline(:,2),'-','linewidth',3,'color',[238 201 0]/255);
                        try
                            legend(l,leg_text,'interpreter','none','location','SouthOutside','AutoUpdate','off'); drawnow;
                        catch
                            legend(l,leg_text,'interpreter','none','location','SouthOutside'); drawnow;
                        end
                    end
                end
                if size(X_crs,1)==0 && keywords.check_other_side
                    [X_crs Y_crs] = polyintersect([baseline(ii,1) baseline(ii,1)-(sind(baseline_angle(ii,1))*max_dist_to_coast)],[baseline(ii,2) baseline(ii,2)-(cosd(baseline_angle(ii,1))*max_dist_to_coast)],coastline(:,1),coastline(:,2));
                    flip_factor  = -1;
                end
                if size(X_crs,1)==0
                    y(ii,1) = NaN;
                else
                    if size(X_crs,1)>1
                        X_crs = X_crs(find(sqrt(((X_crs - baseline(ii,1)).^2) + ((Y_crs - baseline(ii,2)).^2)) == min(sqrt(((X_crs - baseline(ii,1)).^2) + ((Y_crs - baseline(ii,2)).^2))),1,'first'),1);
                        Y_crs = Y_crs(find(sqrt(((X_crs - baseline(ii,1)).^2) + ((Y_crs - baseline(ii,2)).^2)) == min(sqrt(((X_crs - baseline(ii,1)).^2) + ((Y_crs - baseline(ii,2)).^2))),1,'first'),1);
                    end
                    y(ii,1) = flip_factor * sqrt(((X_crs - baseline(ii,1)).^2) + ((Y_crs - baseline(ii,2)).^2));
                    if keywords.plot_figure
                        if flip_factor == 1
                            plot([baseline(ii,1) X_crs],[baseline(ii,2) Y_crs],'k--','linewidth',2);
                        else
                            plot([baseline(ii,1) X_crs],[baseline(ii,2) Y_crs],'k:','linewidth',2);
                        end
                        text(mean([baseline(ii,1) X_crs]),mean([baseline(ii,2) Y_crs]),['y_0 = ' num2str(y(ii,1),'%9.1f') ' mtr.'],'horizontalalignment','center','clipping','on'); drawnow;
                    end
                end
                if ii==size(baseline,1)
                    fprintf(1,[repmat('\b',1,(length(num2str(100*(ii/size(baseline,1)),'%9.0f'))+2)) num2str(100*(ii/size(baseline,1)),'%9.0f') ' %%']);
                    fprintf(1,'\n');
                    disp(' ');
                end
            end
        catch err
            fprintf(1,'\n');
            disp(' ');
            error(err.message);
        end

        if sum(isnan(y)) ~= size(baseline,1)
            % Remove preceeding and trailing NaN's:
            num_ini_nans_removed = 0;
            num_end_nans_removed = 0;
            while isnan(y(end,1))
                if keywords.plot_figure
                    rem_handle = plot(baseline(end,1),baseline(end,2),'rx','markersize',16,'linewidth',3);
                end
                y(end,:)              = [];
                baseline(end,:)       = [];
                baseline_angle(end,:) = [];
                num_end_nans_removed = num_end_nans_removed + 1;
            end
            while isnan(y(1,1))
                if keywords.plot_figure
                    rem_handle = plot(baseline(1,1),baseline(1,2),'rx','markersize',16,'linewidth',3);
                end
                y(1,:)              = [];
                baseline(1,:)       = [];
                baseline_angle(1,:) = [];
                num_ini_nans_removed = num_ini_nans_removed + 1;
            end
        else
            num_ini_nans_removed = 0;
            num_end_nans_removed = 0;
        end
        
        % remove NaNs
        if keywords.interp_missing_ldb
            y=interpNANs(y);
        end
        
        if ((num_ini_nans_removed + num_end_nans_removed)>0) & keywords.plot_figure
            leg_tel = leg_tel + 1;
            leg_text{leg_tel} = ['Removed (outer) baseline points'];
            l(leg_tel) = rem_handle;
            try
                legend(l,leg_text,'interpreter','none','location','SouthOutside','AutoUpdate','off');
            catch
                legend(l,leg_text,'interpreter','none','location','SouthOutside');
            end
        end

        if sum(isnan(y))>0
            if keywords.plot_figure
                for ii = find(isnan(y)==1)'
                    CS_handle = plot([baseline(ii,1)-(sind(baseline_angle(ii,1))*max_dist_to_coast) baseline(ii,1)+(sind(baseline_angle(ii,1))*max_dist_to_coast)],[baseline(ii,2)-(cosd(baseline_angle(ii,1))*max_dist_to_coast) baseline(ii,2)+(cosd(baseline_angle(ii,1))*max_dist_to_coast)],'r--','linewidth',3); drawnow;
                end
                leg_tel = leg_tel + 1;
                leg_text{leg_tel} = ['Cross-sections without perpendicular coastline data'];
                l(leg_tel) = CS_handle;
                try
                    legend(l,leg_text,'interpreter','none','location','SouthOutside','AutoUpdate','off');
                catch
                    legend(l,leg_text,'interpreter','none','location','SouthOutside');
                end
                xlim([min([baseline_ori(:,1); coastline(:,1)]) max([baseline_ori(:,1); coastline(:,1)])]);
                ylim([min([baseline_ori(:,2); coastline(:,2)]) max([baseline_ori(:,2); coastline(:,2)])]);
                title('The red line(s) indicate the cross-section(s) without perpendicular coastline data');
                drawnow;
                if sum(isnan(y)) == size(baseline,1)
                    title({'The red line(s) indicate the cross-section(s) without perpendicular coastline data';'None of them feature coastline information...'});
                    error('Not a single line crosses the coastline (see figure)');
                else
                    fprintf(['Warning : Not all internal baseline points have perpendicular coastline data (you can check the red lines in the figure for help)\n']);
                    error(['Some internal baseline points have no perpendicular coastline data (you can check the red lines in the figure for help)']);
                end
            end
            if sum(isnan(y)) == size(baseline,1)
                error('Not a single line crosses the coastline (you can set plot_figure = true (1) to help you show the problem in a figure)');
            else
                fprintf(['Warning : Not all internal baseline points have perpendicular coastline data (you can set plot_figure = true (1) to help you show the problem in a figure)\n']);
                error(['Some internal baseline points have no perpendicular coastline data (you can set plot_figure = true (1) to help you show the problem in a figure)']);
            end
        end

        if sum(y<0) == size(baseline,1)
            disp(['ATTENTION: all coastline information resembles negative y-values, you are advised to flip the orientation of the baseline']);
            disp(['You can do this by calling baseline = flipud(baseline);']);
            disp(' ');
        elseif sum(y<0) > 0
            disp(['ATTENTION: ' num2str(sum(y<0)) ' baseline points feature negative y-values']);
        end

        if num_ini_nans_removed>0
            disp(['ATTENTION: ' num2str(num_ini_nans_removed) ' initial baseline points were removed as no perpendicular coastline data was found']);
        end
        if num_end_nans_removed>0
            disp(['ATTENTION: ' num2str(num_end_nans_removed) ' trailing baseline points were removed as no perpendicular coastline data was found']);
        end
        if ((sum(y<0) > 0) | ((num_ini_nans_removed+num_end_nans_removed)>0)) & (keywords.plot_figure == 0)
            disp('If you want, you can inspect this behaviour by setting plot_figure to true (1)');
        elseif ((sum(y<0) > 0) | ((num_ini_nans_removed+num_end_nans_removed)>0)) & (keywords.plot_figure == 1)
            disp('Please check the figure to see if you agree with this behaviour');
        end

    end

    % Add additional gridpoints in between the baselinepoints if keywords.min_resolution requires this:   
    if ~isempty(keywords.min_resolution)
        N = [0; ceil(diff(pathdistance(baseline(:,1),baseline(:,2)))/keywords.min_resolution)];
        N = min(99,N); % limiter
        N(find(N(2:end)<1)+1) = 1; % Make sure each value (except the first) is above 0 
    else
        N = [0; ones(size(baseline,1)-1,1)];
    end
    Ray = [1:size(baseline,1)]';

    %% Write everything to the MDA file:
    if keywords.save_to_mda
        fid=fopen(mda_filename,'wt');
        fprintf(fid,'%s\n',' BASISPOINTS');
        fprintf(fid,'%4.0f\n',length(N));
        fprintf(fid,'%s\n','     Xw             Yw             Y              N              Ray');
        fprintf(fid,'%18.4f   %18.4f %12.4f %11.0f %11.0f\n',[baseline y(:) N Ray]');
        fclose(fid);
        if keywords.plot_figure
            disp([mda_filename ' was saved succesfully']);
        end
    else
        disp(' ');
        disp('Script has completed succesfully, see the figure for information');
    end
end


function [data2]=interpNANs(data,useextrapolation)
    % function [data2]=interpNANs(data)
    %
    % INPUT:
    %     data                Array [Nx1] with data
    %     useextrapolation    switch to use extrapolation (use 'nearest' or 'linear' for extrapolation)
    %                         this does not affect the interpolation which will be 'linear' 
    %
    % data = [10, 9, nan, 8, nan, 7.5, nan, nan]

    if nargin<2
        useextrapolation='nearest';  %<- use extrapolation if it is set to 1
    end

    rotate=0;
    if max(size(data,1),2)<size(data,2)
        rotate=1;
        data=data';
    end

    if size(data,2)==1
        x = [1:length(data)]';
        y = data;
    else
        x = data(:,1);
        y = data(:,2);
    end

    %% remove nans
    id = find(~isnan(y));
    x2 = x(id);
    y2 = y(id);

    if length(id)==length(x)
         ynew = y2;
    else
        if length(x2)>1
            if strcmpi(useextrapolation,'linear')
                if ~isempty(find(isnan(x)))
                    xv  = [1:length(x)];
                    idn = find(~isnan(x));
                    x   = interp1(xv(idn),x(idn),xv,'linear','extrap');
                    x   = x(:);
                end
                idn = find(~isnan(x2));
                ynew = interp1(x2(idn),y2(idn),x,'linear','extrap');
            elseif strcmpi(useextrapolation,'nointerp') || ~isempty(findstr(useextrapolation,'no'))
                y3 = interp1(x2,y2,x,'linear');
                id = find(~isnan(y3));
                x2 = x(id);
                y2 = y3(id);
                ynew = interp1(x2,y2,x,'linear');
            else
                y3 = interp1(x2,y2,x,'linear');
                id = find(~isnan(y3));
                x2 = x(id);
                y2 = y3(id);
                ynew = interp1(x2,y2,x,'nearest','extrap');
            end
        elseif length(x2)==1
            ynew=repmat(y2,size(x));
        else
            %fprintf('Need more data points than 1\n')
            ynew=nan(size(x));
        end
    end

    if size(data,2)==2
        data2 = [x,ynew];
    else
        data2 = [ynew];
    end

    if rotate==1
        data2 = data2';
    end
end
