function Transform2hcwd(s,TRmat)
% function Transform2hcwd(s,TRmat)
%
% This function nearshore wave time-series on the basis of offshore time-series and a 
% transformation matrix that was made with Delft3D-waevs / SWAN for the Holland coast and
% Wadden coast (only North sea coast!). The offshore time-series are derived from DONAR 
% (i.e. from openDAP server at Deltares) and supplemented with MATROOS time-series for the
% recent data. The SWAN model (HCWD) uses 391 wave conditions (10° wave direction bins and 
% 0.5m wave height bins). Four wave stations are used, which are Europlatform, IJmuiden 
% munitie platform, Eierlandse gat and Schiermonnikoog Noord. Space varying wind fields 
% are applied for each condition (i.e the wind conditions at the five stations are used).
% 
% INPUT:
%         s         Structure with information on the nearshore location
%                    .StartDate     Datenum of the start date
%                    .EndDate       Datenum of the end date
%                    .Loc.X         X-location of nearshore point [m RD]
%                    .Loc.Y         Y-location of nearshore point [m RD]
%                    .OutputFile    String with output filename of .ASC file
%         TRmat     String with filename of transformation matrix OR a structure with
%                    the fields 'Grid','Input','Trans','WCmat','classes' and 'Depth'
%
% OUTPUT:
%         ASC file with nearshore wave time-series
%         Plot of the nearshore wave data

    fprintf('Transform wave conditions for X=%1.0fm and Y=%1.0fm\n',s.Loc.X,s.Loc.Y);
    jjREF=2;
    Locations    = {'EURPFM.nc','IJMDMNTSPS.nc','EIELSGT.nc','SCHIERMNOND.nc','SCHIERMNOND.nc'};  
    %Locations   = {'SCHOUWBACSH.nc','EURPFM.nc','IJMDMNTSPS.nc','EIELSGT.nc','SCHIERMNOND.nc','SCHIERMNOND.nc'};
    Parameters   = {'sea_surface_wave_significant_height','sea_surface_wave_from_direction','sea_surface_wind_wave_tm02'};
    fieldnms     = {'EUR','YM6','ELD','SON','SON'}; 
    %fieldnms    = {'SCH','EUR','YM6','ELD','SON','SON'};
    debugmode    = 0;


    [pthnm,fn,en] = fileparts(which('Transform2hcwd.m'));
    if ~strcmpi(pthnm(end),filesep)
        pthnm=[pthnm,filesep];
    end;
    addpath([pthnm,'functions\']);
    if nargin<2;
        TRmat = 'TransMatrix_InputParameters_hcwd.mat';
    end

    %-----------------------------------------------------------------------------------------
    %% LOADING TRANSFORMATION MATRIX
    %-----------------------------------------------------------------------------------------
    fprintf('  - Loading transformation matrix\n');
    if isstr(TRmat)
        load([pthnm,TRmat]);                                              % load input matrix for conversion
    else
        Grid        =  TRmat.Grid;
        Input       =  TRmat.Input;
        Trans       =  TRmat.Trans;
        WCmat       =  TRmat.WCmat;
        classes     =  TRmat.classes;
        Depth       =  TRmat.Depth;
    end
    % Grid, Input, Trans, WCmat, alpha, classes;
    NRstations = length(fieldnms);
    % waveconditions from SWAN simulations
    for jj=1:NRstations
        WC.(fieldnms{jj}) = Input.data((2+jj):(NRstations+4):length(Input.data),:);
    end

    %-----------------------------------------------------------------------------------------
    %% LOADING DONAR AND MATROOS DATA
    %-----------------------------------------------------------------------------------------
    fprintf('  - Loading offshore data from DONAR and MATROOS\n');
    D=getDONARdata_hcwd(Locations,Parameters);

    % ------MATROOS--------------
    IDfs=findstr(pthnm,'\');
    urlChar = [pthnm(1:IDfs(end-1)),'WaveTable\'];
    if ~exist(urlChar,'dir')
        urlChar=['p:',filesep,'webtools',filesep,'tools',filesep,'WaveTable',filesep];        %urlChar = 'd:\repositories\wavetransformationtool\WaveTable\'; %'p:/webtools/tools/WaveTable/';
    end
    for jj=1:length(fieldnms)
        ASCmatroosfiles{jj} = [urlChar fieldnms{jj},'_matroos.asc'];
        ASCdatafiles{jj} = [urlChar fieldnms{jj},'_data.asc'];
    end
    if exist([urlChar,'ASCdata.mat'],'file');
        load([urlChar,'ASCdata.mat']);
    else
        [M]=readASCwavefiles(ASCmatroosfiles);
        [B]=readASCwavefiles(ASCdatafiles);
    end
    
    %-----------------------------------------------------------------------------------------
    %% Fill small gaps in data with interpolation
    %-----------------------------------------------------------------------------------------
    fields1 = {'hm0','th0','tm02','dspr','u10','udir','wl','wlastro','wlsurge'};
    Ddt = (D(1).hm0(1).time(2)-D(1).hm0(1).time(1))*24; % time step in hours
    Mdt = (M(1).time(2)-M(1).time(1))*24; % time step in hours
    Bdt = (B(1).time(2)-B(1).time(1))*24; % time step in hours
    [D]=fillGAPSindata(D,fields1,round(6/Ddt));
    [M]=fillGAPSindata(M,fields1,round(6/Mdt));
    [B]=fillGAPSindata(B,fields1,round(6/Bdt));
    
    %% Combine DONAR (D) and MATROOS (M)
    timestep = 3; % hours in-between timesteps
    [W]=addDATA(D,M,B,fields1,timestep);

    clear D M B id
    %% PLOT DATA
    if debugmode
        colours={'b-','r-','g-','m-','c-','y-','b--','r--','g--','m--','c--'};
        for kk=1:length(fields1)
            figure(kk);set(gcf,'Position',[140,70,1100,900])
            for jj=1:NRstations
                subplot(NRstations,1,jj);
                plot(W(jj).time(kk).datenum,W(jj).(fields1{kk}).data,colours{kk});datetick
                title([fields1{kk},' at station ',num2str(jj)]);
            end
        end
    end


    %-----------------------------------------------------------------------------------------
    %% DETERMINE M AND N INDEX OF SPECIFIED XY COORDINATES
    %-----------------------------------------------------------------------------------------
    fprintf('  - Determine M and N index of specified XY-coordinates\n');
    aa                  = sqrt(((Grid.X - s.Loc.X).^2 + (Grid.Y - s.Loc.Y).^2));
    [s.Loc.M,s.Loc.N]   = find(aa==min(min(aa)));

    %-----------------------------------------------------------------------------------------
    %% Determine depth at grid cell
    %-----------------------------------------------------------------------------------------
    DEPloc = Depth(s.Loc.M,s.Loc.N);

    %-----------------------------------------------------------------------------------------
    %% convert timeseries
    %-----------------------------------------------------------------------------------------
    % make new timeline and fill the data to this timeline
    %   3-hourly date up to 01-Jan-1987
    %   hourly data up to 01-Jul-2007
    %   10 minutes interval up to now
    fprintf('  - Setting up offshore time series\n');
    tstart = [];
    for jj=1:NRstations
        tmp=[];
        for kk=1:length(W(jj).time)
            tmp =[tmp,W(jj).time(kk).datenum(1)];
        end
        tstart =  [tstart,tmp];
    end

    dt       = 3;  % hour between output timesteps
    Tnew     = unique([min(tstart):dt/24:W(jj).time(1).datenum(end)]');
    %fields1 = {'hm0','th0','tm02'};
    for jj=1:NRstations
        for i=1:length(W(jj).time)
            % IJMUIDEN
            % remove doubles
            [dum,idt]                = unique(W(jj).time(i).datenum);
            W(jj).time(i).datenum    = W(jj).time(i).datenum(idt);
            W(jj).(fields1{i}).data  = W(jj).(fields1{i}).data(idt);

            % line up timeseries to Tnew
            dumT                = nan(length(Tnew),1);
            dumV                = nan(length(Tnew),1);
            a                   = ismember(round(24*6*Tnew),round(24*6*W(jj).time(i).datenum));
            b                   = ismember(round(24*6*W(jj).time(i).datenum),round(24*6*Tnew));
            dumT(a)             = W(jj).time(i).datenum(b); 
            W(jj).time(i).datenum = Tnew;
            dumV(a)= W(jj).(fields1{i}).data(b);  W(jj).(fields1{i}).data = dumV;
            clear dumT dumV a b dum idt;
        end
    end

    % %% if offshore direction 30 < dir < 200, then no transformation
    % % because these data are not classified in the matrix
    % for jj=1:NRstations
    % nid                 = find(W(jj).th0.data >= 30 & W(jj).th0.data <= 200);
    % W(jj).hm0.data(nid)   = NaN;      W(jj).th0.data(nid)   = NaN;      W(jj).tm02.data(nid)  = NaN;      % W(jj).surge.data(nid) = NaN;
    % end

    % transformation of tm02 to tp (because no realible results for Tm, only for tp)
    for jj=1:NRstations
        W(jj).tp              = W(jj).tm02;
        W(jj).tp.data         = W(jj).tp.data * 1.28; 
    end

    %-----------------------------------------------------------------------------------------
    %% Determine to which class each timepoint belongs
    %-----------------------------------------------------------------------------------------
    fprintf('  - Determine the class-bin for all timepoints\n');
    IDii=zeros(size(W(jj).hm0.data));
    IDjj=zeros(size(W(jj).hm0.data));
    for i=1:length(classes{1})
        ID1 = (W(jjREF).hm0.data >  classes{1}(i,1));
        ID2 = (W(jjREF).hm0.data <=  classes{1}(i,2));
        IDii = IDii+ID1.*ID2.*i;
    end
    for j=1:length(classes{2})
        ID3 = (W(jjREF).th0.data >  classes{2}(j,1));
        ID4 = (W(jjREF).th0.data <=  classes{2}(j,2));
        IDjj = IDjj+ID3.*ID4.*j;
    end
    ID = find(W(jj).hm0.data==0);
    IDii(ID)=1;
    ID = find(W(jj).th0.data==0);
    IDjj(ID)=1;

    %% All station data
    fields2 = {'hm0','th0','tm02','tp','dspr','u10','udir','wl','wlastro','wlsurge'};
    for kk=1:length(fields2)
        stations.(fields2{kk}) = [];
        for jj=1:NRstations        
            stations.(fields2{kk}) = [stations.(fields2{kk}),W(jj).(fields2{kk}).data];
        end
    end

    %-----------------------------------------------------------------------------------------
    %% PLOT correlation between stations
    %-----------------------------------------------------------------------------------------
    if debugmode
    for kk=1:length(fields2)
    figure(11+kk);clf;set(gcf,'Position',[140,70,1100,900])
    for jj=1:NRstations;
        for ii=1:NRstations;
            subplot(NRstations,NRstations,(jj-1)*NRstations+ii);
            plot(stations.(fields2{kk})(:,jj),stations.(fields2{kk})(:,ii),'.');
            if ii==1;ylabel([fields2{kk},' station ',num2str(jj)]);end
            if jj==NRstations;xlabel([fields2{kk},' station ',num2str(ii)]);end
        end;
    end
    end
    end

    %-----------------------------------------------------------------------------------------
    %% Compute nearshore wave conditions from transformation matrix and offshore wave conditions
    %-----------------------------------------------------------------------------------------
    fprintf('  - Transforming offshore data to nearshore [');
    % find right transformation factor from the matrix
    % step 1: classify data
    % step 2: classify SWAN runs
    % step 3: find right factor
    vid                 = find(Tnew>=s.StartDate & Tnew<s.EndDate);
    output.hm0.data     = repmat(NaN,size(vid)); 
    output.hdir.data    = repmat(NaN,size(vid)); 
    output.tp.data      = repmat(NaN,size(vid));
    output.dspr.data    = repmat(NaN,size(vid)); 
    output.u10.data     = repmat(NaN,size(vid)); 
    output.udir.data    = repmat(NaN,size(vid));
    output.wl.data      = repmat(NaN,size(vid)); 
    output.wlastro.data = repmat(NaN,size(vid)); 
    output.wlsurge.data = repmat(NaN,size(vid));
    %hh = waitbar(0,'Transforming offshore data to nearshore, Please wait .......');
    for n=1:length(vid)
        %waitbar(n/length(vid),hh)
        % classify values, so that the transformation matrix can be used
        ii=IDii(vid(n));
        jj=IDjj(vid(n));

        if (ii~=0 && jj~=0)
            nc           = WCmat(ii,jj);

            output.hm0.data(n)      = nan;
            output.hdir.data(n)     = nan;
            output.tp.data(n)       = nan;
            output.dspr.data(n)     = nan;
            output.u10.data(n)      = nan;
            output.udir.data(n)     = nan;
            output.wl.data(n)       = nan;
            output.wlastro.data(n)  = nan;
            output.wlsurge.data(n)  = nan;
            output.depthavg.data(n) = mean(DEPloc(:));
            if ~isnan(nc)
                IDstat(nc)   = Trans(nc).IDstation(s.Loc.N);
                MAT_hs       = Trans(nc).Hs(s.Loc.M,s.Loc.N);
                MAT_hdir     = Trans(nc).Hdir(s.Loc.M,s.Loc.N);
                MAT_tp       = Trans(nc).Tp(s.Loc.M,s.Loc.N);
                output.hm0.data(n)      = MAT_hs     .* stations.hm0(vid(n),IDstat(nc));
                output.hdir.data(n)     = MAT_hdir   +  stations.th0(vid(n),IDstat(nc));  % directions, so summation
                output.tp.data(n)       = MAT_tp     .* stations.tp(vid(n),IDstat(nc));

                % for other parameters DONAR and MATROOS data are needed.
                MAT_dspr     = Trans(nc).Dspr(s.Loc.M,s.Loc.N);
                MAT_u10      = Trans(nc).U10(s.Loc.M,s.Loc.N);
                MAT_udir     = Trans(nc).Udir(s.Loc.M,s.Loc.N);
                output.dspr.data(n)    = MAT_dspr   .* stations.dspr(vid(n),IDstat(nc));  % directions, so summation
                output.u10.data(n)     = MAT_u10    .* stations.u10(vid(n),IDstat(nc));  % directions, so summation
                output.udir.data(n)    = MAT_udir   +  stations.udir(vid(n),IDstat(nc));  % directions, so summation

                % extra information (from nearest station)
                MAT_surge    = Trans(nc).Surge(s.Loc.M,s.Loc.N);
                output.wl.data(n)      = stations.wl(vid(n),IDstat(nc));
                output.wlastro.data(n) = stations.wlastro(vid(n),IDstat(nc));
                output.wlsurge.data(n) = MAT_surge  +  stations.wlsurge(vid(n),IDstat(nc));
            end
        end
        if mod(n,floor(length(vid)/25))==0;fprintf('#');end
    end
    %close(hh);
    fprintf(']\n');

    % nan the periods without waves
    idhm0_0                        = find(output.hm0.data==0);
    output.hm0.data(idhm0_0)       = NaN; 
    output.hdir.data(idhm0_0)      = NaN; 
    output.tp.data(idhm0_0)        = NaN; 
    output.dspr.data(idhm0_0)      = NaN; 
    output.u10.data(idhm0_0)       = NaN;
    output.udir.data(idhm0_0)      = NaN;
    output.wl.data(idhm0_0)        = NaN;
    output.wlastro.data(idhm0_0)   = NaN;
    output.wlsurge.data(idhm0_0)   = NaN;

    % make directions from 0-360 deg.
    output.hdir.data            = mod(output.hdir.data,360);
    output.udir.data            = mod(output.udir.data,360);
    if ~isfield(output,'depthavg');output.depthavg.data=nan;end
    
    %-----------------------------------------------------------------------------------------
    %% OUTPUT DATA STRUCTURE
    %-----------------------------------------------------------------------------------------   
    Output.file   = s.OutputFile;
    Output.X      = s.Loc.X;
    Output.Y      = s.Loc.Y;
    Output.Data   = [output.hm0.data(:) output.hdir.data(:) output.tp.data(:) output.dspr.data(:) output.wlsurge.data(:) output.wlastro.data(:) output.wl.data(:) output.u10.data(:) output.udir.data(:)];
    Output.Time   = {datestr(Tnew(vid),'yyyymmdd') datestr(Tnew(vid),'HHMMSS')};
    Output.Depavg = output.depthavg.data(1);
    
    %-----------------------------------------------------------------------------------------
    %% EXPORT DATA TO ASC FILE
    %-----------------------------------------------------------------------------------------
    fprintf('  - Export data to ''%s''\n',s.OutputFile);
    exportASCfile(Output);

    %-----------------------------------------------------------------------------------------
    %% PLOTTING
    %-----------------------------------------------------------------------------------------
    fprintf('  - Plotting data\n\n');
    plotwavedata(Tnew,vid,Output,s);

end



%-----------------------------------------------------------------------------------------
%% READ ASCwavefiles
%-----------------------------------------------------------------------------------------
function [M]=readASCwavefiles(ASCwavefiles)
    M=struct;
    NRstations = length(ASCwavefiles);
    for jj=1:NRstations
        if sum(strcmpi(ASCwavefiles(1:jj),ASCwavefiles{jj}))>1
            % same station is used for two locations
            ID = find(strcmpi(ASCwavefiles(1:jj),ASCwavefiles{jj}));
            M(jj) = M(ID(1));
        else
            % data  = fileread([urlChar fieldnms{jj},'_matroos.asc']);
            % id1 = regexp(data,'\n');
            % data  = str2num(reshape(data,id1(1),(length(data)/id1(1)))');
            data  = load(ASCwavefiles{jj});
            %     time             hs           xdir           tm02         ';
            M(jj).file = ASCwavefiles{jj};
            M(jj).time = data(:,1);
            M(jj).hm0 = data(:,2);
            M(jj).th0 = data(:,3);
            M(jj).tm02 = data(:,4);
            if size(data,2)>4
                %     time             hs           xdir           tm02           tp         dirspr       windspeed      winddir          wl         wlastro        wlsurge';
                M(jj).dspr = data(:,6);
                M(jj).u10 = data(:,7);
                M(jj).udir = data(:,8);
                M(jj).wl = data(:,9);
                M(jj).wlastro = data(:,10);
                M(jj).wlsurge = data(:,11);
            else
                M(jj).dspr = nan(size(data(:,1)));
                M(jj).u10 = nan(size(data(:,1)));
                M(jj).udir = nan(size(data(:,1)));
                M(jj).wl = nan(size(data(:,1)));
                M(jj).wlastro = nan(size(data(:,1)));
                M(jj).wlsurge = nan(size(data(:,1)));
            end
        end
    end
end
    
%-----------------------------------------------------------------------------------------
%% COMBINE DATA
%-----------------------------------------------------------------------------------------
function [W]=addDATA(D,M,B,fields1,timestep)

    if nargin<5
        timestep=6; % hours
    end
    W = struct;
    NRstations = length(D);
    fprintf('  - adding data of different sources [');
    
    Tmin=datenum('2000-01-01','yyyy-mm-dd');
    Tmax=datenum('2000-01-01','yyyy-mm-dd');
    for jj=1:NRstations
        Dtime = D(jj).time(1).datenum(:);
        Mtime = M(jj).time(1).datenum(:);
        Btime = B(jj).time(1).datenum(:);
        Tmin = min(Tmin,round(24*timestep*min([Dtime;Mtime;Btime]))/24/timestep);
        Tmax = max(Tmax,round(24*timestep*max([Dtime;Mtime;Btime]))/24/timestep);
    end
    for jj=1:NRstations
        fprintf('#')
        for kk=1:length(fields1)
            
            % timelimits
            Dtime = D(jj).time(kk).datenum(:);
            Mtime = M(jj).time(kk).datenum(:);
            Btime = B(jj).time(kk).datenum(:);
            % Tmin = round(24*timestep*min([Dtime;Mtime;Btime]))/24/timestep;
            % Tmax = round(24*timestep*max([Dtime;Mtime;Btime]))/24/timestep;
            timenum  = [Tmin:timestep/24:Tmax]';
            
            Ddata = D(jj).(fields1{kk}).data(:);
            Mdata = M(jj).(fields1{kk}).data(:);
            Bdata = B(jj).(fields1{kk}).data(:);

            % combine DONAR and MATROOS
            [DMdata]=combineDATA(Dtime,Ddata,Mtime,Mdata,timenum);
            [DMBdata]=combineDATA(timenum,DMdata,Btime,Bdata,timenum);
            W(jj).time(kk).datenum     = timenum;
            W(jj).(fields1{kk}).data   = DMBdata;
        end
    end
    fprintf('#]\n')
end

function [data12]=combineDATA(time1,data1,time2,data2,time12)
% combine data with different time signal

    % find intersection between time12 and time1
    %  -> apply data2 to corresponding time12 instances
    if nargin<5
        time12                 = time1;
        data12                 = data1;
    else
        [a,ID1,ID2] = intersect(floor(time12*10^3)/10^3,floor(time1*10^3)/10^3);
        data12                = nan(size(time12));
        data12(ID1)           = data1(ID2);
        %datatmp(datatmp==-999) = nan;
    end
    % find intersection between time12 and time2
    %  -> apply data2 to corresponding time12 instances
    %  -> remove possible missingvalues (-999 to nan)
    [a,ID1,ID2] = intersect(floor(time12*10^3)/10^3,floor(time2*10^3)/10^3);
    datatmp                = nan(size(time12));
    datatmp(ID1)           = data2(ID2);
    %datatmp(datatmp==-999) = nan;
    % fill new data structure
    IDnan                  = find(isnan(data12));
    data12(IDnan)          = datatmp(IDnan);
end


%-----------------------------------------------------------------------------------------
%% EXPORT ASC FILE
%-----------------------------------------------------------------------------------------
function exportASCfile(Output)
    fid = fopen(Output.file,'w');
    fprintf(fid,'%s\n','Transformed wave data from the offshore platforms ELD, YM6 and EUR to nearshore');
    fprintf(fid,'%s\n','NaN values mean that the wave direction is between 30-200 deg, so no transformation is done');
    fprintf(fid,'%s\n',['* Output Location X: ' num2str(Output.X) ' [m; RD]']);
    fprintf(fid,'%s\n',['* Output Location Y: ' num2str(Output.Y) ' [m; RD]']);
    fprintf(fid,'%s\n',['* Depth without surge or tide: ' num2str(Output.Depavg,'%3.1f') ' [m MSL]']);
    fprintf(fid,'%s\n','* column 1 = Date');
    fprintf(fid,'%s\n','* column 2 = Time');
    fprintf(fid,'%s\n','* column 3 = hm0 (m)');
    fprintf(fid,'%s\n','* column 4 = hdir (deg N)');
    fprintf(fid,'%s\n','* column 5 = tp (s)');
    fprintf(fid,'%s\n','* column 6 = directional spreading (deg)');
    fprintf(fid,'%s\n','* column 7 = Surge (m)');
    fprintf(fid,'%s\n','* column 8 = Tide (astronomical) (m MSL)');
    fprintf(fid,'%s\n','* column 9 = Total water level (m MSL)');
    fprintf(fid,'%s\n','* column 10 = Wind velocity (m/s)');
    fprintf(fid,'%s\n','* column 11 = Wind direction (deg N)');
    for i=1:size(Output.Data,1)
        fprintf(fid,'%s',[Output.Time{1}(i,:) '  ' Output.Time{2}(i,:)]);
        fprintf(fid,[repmat('%8.3f   ',1,size(Output.Data,2)) '\n'],Output.Data(i,:)');
    end
    fclose(fid);
end

%-----------------------------------------------------------------------------------------
%% PLOT wave data
%-----------------------------------------------------------------------------------------
function plotwavedata(Tnew,vid,Output,s)
    time = Tnew(vid);
    xx   = (time(end)-time(1))/5;
    if xx > 0.8 
        xx   = round((time(end)-time(1))/5);
        xxx  = time(1):xx:time(end);
        xxxx = datestr(xxx,1);
        xl   = 'date [dd-mmm-yyyy]';
    elseif xx <= 0.8
        xx   = round(24*((time(end)-time(1))/5))/24;
        xxx  = time(1):xx:time(end);
        xxxx = datestr(xxx,16);  
        xl   = 'date [HH:MM]';
    end

    % -------------------------------------------------------------------------
    % get large nan sections (larger than 5)
    id_nan = [find(~isnan(Output.Data(:,1))) ; size(Output.Data,1)]; 

    k=0;
    if length(id_nan)>1;
        for i=1:length(id_nan)
            if i==1
                t_diff(i,:) = [1 id_nan(i)-1];
            else
                t_diff(i,:) = [id_nan(i-1)+1 id_nan(i)-1];
            end
        end
        for i=1:length(t_diff)
            if t_diff(i,2)-t_diff(i,1) > 5;
                k=k+1;
                tt_diff(k,:) = t_diff(i,:);
            end
        end
    end
    % -------------------------------------------------------------------------

    %% determine IDnan's
    IDnan=[];
    if exist('tt_diff')
        for i=1:size(tt_diff,1)
            IDnan = [IDnan;[tt_diff(i,1):tt_diff(i,2)]'];
        end
    end
    Output.Data(:,1) = interpNANs(Output.Data(:,1));
    Output.Data(:,2) = interpNANs(Output.Data(:,2));
    Output.Data(:,3) = interpNANs(Output.Data(:,3));
    Output.Data(IDnan,:)=nan;

    % -------------------------------------------------------------------------

    ColorNaN = [1 1 0.88];

    figure('PaperType','A4','PaperOrientation','landscape','PaperPositionMode','manual',...
        'PaperUnits','centimeters','Paperposition',[0.5 0.5 28.7 20],'Position',[100 100 1124 768]);
    as_1(1) = axes('Position',[.005 .005 .99 .99],'Xtick',[],'Ytick',[],'box','on','handlevisibility','off');
    as_1(2) = axes('Position',[0.25 0.68 0.72 0.24]);

    if ~isempty(find(~isnan(Output.Data)))
    if exist('tt_diff')
        for i=1:size(tt_diff,1)
            t = [[time(tt_diff(i,1));time(tt_diff(i,2))]]; t=[t;flipud(t)];
            v = [0; 0; max(Output.Data(:,1)); max(Output.Data(:,1))];
            fill(t,v,ColorNaN,'EdgeColor',ColorNaN); hold on;
        end
    end
    ll = legend('Wave Direction between 30-200°');
    set(ll,'Position',[0.83 0.93 0.1 0.05]);
    plot(time,Output.Data(:,1),'b');
    set(gca,'XTick',xxx,'XTickLabel',[])
    xlim([time(1) time(end)])
    yy = [0:0.5:max(Output.Data(:,1))];
    set(gca,'YTick',yy,'YTickLabel',yy,'FontSize',8)
    ylim([yy(1) yy(end)])
    ylabel('hm0 [m]')
    title('hm0 [m]','FontWeight','bold');
    grid on; set(gca,'Layer','top');
    box on;

    as_1(3) = axes('Position',[0.25 0.38 0.72 0.24]);
    if exist('tt_diff')
        for i=1:size(tt_diff,1)
            t = [[time(tt_diff(i,1));time(tt_diff(i,2))]]; t=[t;flipud(t)];
            v = [0; 0; 360; 360];
            fill(t,v,ColorNaN,'EdgeColor',ColorNaN); hold on;
        end
    end
    plot(time,Output.Data(:,2),'b')
    set(gca,'XTick',xxx,'XTickLabel',[])
    xlim([time(1) time(end)])
    yy = [0:30:360];
    set(gca,'YTick',yy,'YTickLabel',yy,'FontSize',8)
    ylabel('direction [°N]')
    ylim([yy(1) yy(end)])
    title('Wave direction [°N]','FontWeight','bold');
    grid on; set(gca,'Layer','top');
    box on;

    as_1(4) = axes('Position',[0.25 0.08 0.72 0.24]);
    if exist('tt_diff')
        for i=1:size(tt_diff,1)
            t = [[time(tt_diff(i,1));time(tt_diff(i,2))]]; t=[t;flipud(t)];
            v = [0; 0; 10; 10];
            fill(t,v,ColorNaN,'EdgeColor',ColorNaN); hold on;
        end
    end
    plot(time,Output.Data(:,3),'b')
    set(gca,'XTick',xxx,'XTickLabel',xxxx,'FontSize',8)
    xlim([time(1) time(end)])
    yy = [0:1:10];
    set(gca,'YTick',yy,'YTickLabel',yy,'FontSize',8)
    ylim([yy(1) yy(end)])
    xlabel(xl);
    ylabel('tp [s]')
    title('tp [s]','FontWeight','bold');
    grid on; set(gca,'Layer','top');
    box on;

    as_1(5) = axes('Position',[0.02 0.02 0.15 0.95]);
    textline{1} = ['Computed wave conditions using the Wave Transformation Table'];
    textline{2} = ['at location X:' num2str([s.Loc.X]) ' Y:' num2str([s.Loc.Y])]; 
    textline{3} = ['from ' datestr(time(1)) ' to ' datestr(time(end))];
    textline{4} = ['developed for Building with Nature HK 3.2'];
    for j=1:length(textline)
        k=length(textline)+1-j;
        text(j/5,0.99,textline{k},'rotation',270,'FontSize',10)
    end
    set(gca,'XTick',[],'YTick',[])
    text(0.09,0.01,'Deltares','rotation',0,'FontSize',6)
    box on;

    PapSize = get(gcf,'PaperSize');             % dimension figure * dimension A4 to get the [width/length] correct
    as_1(6) = axes('Position',[0.03 0.04 (0.25*(746/1033)*(PapSize(2)/PapSize(1))) 0.25]);
    set(gca,'XTick',[],'YTick',[]); 
    box on; set(gca,'Layer','top');
    end

    %%
    %print(gcf,'-dpng','-painters',[s.OutputFile(1:end-4)]);
    print(gcf,'-dpdf','-painters',[s.OutputFile(1:end-4)]);
    close all; 
    fclose all;
end


