% Generate an EISCAT HDF5-file from mat-files generated in a Guisdap analysis

function [Storepath,EISCAThdf5file] = mat2hdf5_vel(matfile_vel,datapath,addfigs,addnotes) 

global path_GUP hdf5ver %result_path 
% if nargin<4, addnotes = []; else addnotes = 1; end 
% if nargin<3, addfigs = []; else addfigs = 1; end 
% if nargin==1, error('Not enough input parameters, path to matfiles folder and path to datastore folder needed'); end
% if nargin<1
%     matpath = result_path;
%     datapath = result_path;
% end
% if isstring(datapath)
%     datapath = char(datapath);    % need to be char class
% end

software = 'https://git.eiscat.se/cvs/guisdap9';
level2_link = '';

load(matfile_vel)

GuisdapParFile = fullfile(path_GUP,'matfiles','Guisdap_Parameters.xlsx'); % path to the .xlsx file
[~,text] = xlsread(GuisdapParFile);
parameters_list = text(:,1);

if exist('V_area','var'), V_area = col(V_area); end

vdate = [];
vleap = [];
vpos  = [];
vg    = [];
vgv   = [];
varea = [];
lrec  = [];

% Sorting time
tolsec = 1e-2;
tol = tolsec/max(Vdate(1,:)*24*3600);       % tolerence = 0,01 s
t1all=uniquetol(sort(Vdate(1,:)*24*3600),tol);
for t1=t1all
    d=find(Vdate(1,:)*24*3600>(t1-tolsec) & Vdate(1,:)*24*3600<(t1+tolsec));
    tol = tolsec/max(Vdate(2,d)*24*3600);   % tolerence = 0,01 s
    t2all=uniquetol(sort(Vdate(2,d)*24*3600),tol);
    for t2=t2all
        d1=d(find(Vdate(2,d)*24*3600>(t2-tolsec) & Vdate(2,d)*24*3600<(t2+tolsec)));
        lrec  = [lrec; length(d1)];
        vdate = [vdate;Vdate(:,d1(1))'];
        if exist('Vleap','var')
            vleap = [vleap;Vleap(:,d1(1))'];
        end
        vpos  = [vpos;Vpos(d1,:)]; 
        vg    = [vg;Vg(d1,:)];
        vgv   = [vgv;Vgv(d1,:)];
        if exist('V_area','var')
            varea = [varea;V_area(d1)];
        end
    end
end
nrecs = length(lrec);

a = find(strcmp('h_nrec',parameters_list)==1);
[~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
info(1) = {info{1}(3:end)};
info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};        
if length(unique(lrec)) == 1
    matfile.data.par0d = lrec(1);
    matfile.metadata.par0d = info';
else
    matfile.data.par1d = lrec;
    matfile.metadata.par1d = info';
end

if exist('Vdate','var')
    parameters_vdate = {'time1' 'time2'};
    if ~isempty(vleap)
        [matfile.data.utime(:,1),leaps(:,1)] = timeconv([vdate(:,1) vleap(:,1)],'mat2unx');
        [matfile.data.utime(:,2),leaps(:,2)] = timeconv([vdate(:,2) vleap(:,2)],'mat2unx');
    else
        matfile.data.utime(:,1) = timeconv(vdate(:,1),'mat2unx');
        matfile.data.utime(:,2) = timeconv(vdate(:,2),'mat2unx');
    end
    for ii = 1:2
        a = find(strcmp(char(parameters_vdate(ii)),parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(4) = {'Vdate'};
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        matfile.metadata.utime(:,ii) = info';
    end
end

l = 0; m = 0; n = 0; 

if exist('Vpos','var')
    parameters_vpos = {'lat' 'lon' 'h_h'};
    vpos(:,3) = vpos(:,3)*1000;   % km --> m
    
    for ii = 1:3
        
        a = find(strcmp(char(parameters_vpos(ii)),parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        if parameters_vpos{ii} == 'h_h' 
            info{1} = 'h'; end
        info(4) = {'Vpos'};
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
    
        start = 0;
        par_1d = [];
        for kk = 1:nrecs
            par_rec = vpos(start+1:start+lrec(kk),ii);
            if length(unique(par_rec)) == 1
                par_1d = [par_1d; par_rec(1)];
            end
            start = start + lrec(kk);
        end
        
        if length(par_1d) == nrecs
            if length(unique(par_1d)) == 1
                l = l + 1; 
                matfile.data.par0d(:,l) = par_1d(1);
                matfile.metadata.par0d(:,l) = info';
            else
                m = m + 1; 
                matfile.data.par1d(:,m) = par_1d;
                matfile.metadata.par1d(:,m) = info';
            end
        else
            n = n + 1; 
            matfile.data.par2d(:,n) = vpos(:,ii);
            matfile.metadata.par2d(:,n) = info';
        end
    end
end

if exist('Vg','var')
    parameters_vg = {'vi_east' 'vi_north' 'vi_up'};
    
    for ii = 1:3
        a = find(strcmp(char(parameters_vg(ii)),parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        
        start = 0;
        par_1d = [];
        for kk = 1:nrecs
            par_rec = vg(start+1:start+lrec(kk),ii);
            if length(unique(par_rec)) == 1
                par_1d = [par_1d; par_rec(1)];
            end
            start = start + lrec(kk);
        end

        if length(par_1d) == nrecs
            m = m + 1; 
            matfile.data.par1d(:,m) = par_1d;
            matfile.metadata.par1d(:,m) = info';
        else
            n = n + 1; 
            matfile.data.par2d(:,n) = vg(:,ii);
            matfile.metadata.par2d(:,n) = info';
        end
    end
end

if exist('Vgv','var')
    parameters_vgv = {'var_vi_east' 'var_vi_north' 'var_vi_up' 'vi_crossvar_12' 'vi_crossvar_23' 'vi_crossvar_13'};
    for ii = 1:6
        a = find(strcmp(char(parameters_vgv(ii)),parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        
        start = 0;
        par_1d = [];
        for kk = 1:nrecs
            par_rec = vgv(start+1:start+lrec(kk),ii);
            if length(unique(par_rec)) == 1
                par_1d = [par_1d; par_rec(1)];
            end
            start = start + lrec(kk);
        end

        if length(par_1d) == nrecs
            m = m + 1; 
            matfile.data.par1d(:,m) = par_1d;
            matfile.metadata.par1d(:,m) = info';
        else
            n = n + 1;
            matfile.data.par2d(:,n) = vgv(:,ii);
            matfile.metadata.par2d(:,n) = info';
        end
    end
end

if exist('V_area','var')

    parameters_varea = {'vi_solidangle'};
    a = find(strcmp(char(parameters_varea),parameters_list)==1);
    [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
    info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        
    start = 0;
    par_1d = [];
    for kk = 1:nrecs
        par_rec = varea(start+1:start+lrec(kk));
        if length(unique(par_rec)) == 1
            par_1d = [par_1d; par_rec(1)];
        end
        start = start + lrec(kk);
    end
    
    circ_area = pi*(pi/2)^2;        % half sphere, projected to circle
    sph_sr = 2*pi;                  % half sphere, sr
    
    if length(par_1d) == nrecs
        m = m + 1; 
        matfile.data.par1d(:,m) = par_1d/circ_area*sph_sr;
        matfile.metadata.par1d(:,m) = info';
    else
        n = n + 1; 
        matfile.data.par2d(:,n) = varea/circ_area*sph_sr;
        matfile.metadata.par2d(:,n) = info';
    end
end

% store the Vinput content
if exist('Vinputs','var')
    Vinputs_Fields = fieldnames(Vinputs);
    for jj = 1:length(Vinputs)
        for field = Vinputs_Fields.'
            metapar = [];
            if ischar(Vinputs(jj).(char(field)))
                matfile.metadata.Vinputs.(char(field))(jj) = {Vinputs(jj).(char(field))}; 
            else
                if strcmp(field,'UpConstriant')
                    matfile.metadata.Vinputs.UpConstraint(jj) = {num2str(Vinputs(jj).UpConstriant)};
                else
                    matfile.metadata.Vinputs.(char(field))(jj) = {num2str(Vinputs(jj).(char(field)))};
                end
            end
        end      
    end
    for field = Vinputs_Fields.'
        if strcmp(field,'UpConstriant')
            field = {'UpConstraint'};
        end
        matfile.metadata.Vinputs.(char(field)) = matfile.metadata.Vinputs.(char(field))';
    end
end
    
starttime = datestr(t1,'yyyy-mm-ddTHH:MM:SS');
endtime   = datestr(t2,'yyyy-mm-ddTHH:MM:SS');

r_time = datevec(Vdate(1,1));
starttime = datestr(r_time);

year = num2str(r_time(1));
month = sprintf('%02d',r_time(2));
day = sprintf('%02d',r_time(3));

hour   = sprintf('%02d',r_time(4));
minute = sprintf('%02d',r_time(5));
second = sprintf('%02.f',r_time(6));

intper_mean = subsetmean(diff(matfile.data.utime,1,2));        % mean integration time (of quartile range 2 and 3)
intper_mean_str = num2str(round(intper_mean,3,'significant'));

e_time = datevec(Vdate(2,end));
endtime = datestr(e_time);

[~,filename,~] = fileparts(matfile_vel); 
pulses = {'arc','beata','bella','cp','CP','folke','hilde','ipy','manda','steffe','taro','tau','othia','tyko','gup0','gup3'};
ants   = {'32m','42m','uhf','vhf','esa','esr','eis','kir','sod','tro','lyr'};
if ~exist('name_expr','var')
    name_expr='missing'; 
    for pp = 1:length(pulses)
        if contains(filename,pulses{pp})
            name_expr = pulses{pp};
        end
    end
end
if ~exist('name_ant','var')
    name_ant='missing';
    for pp = 1:length(ants)
        if contains(filename,ants{pp})
            name_ant = ants{pp};
        end
    end
end
if ~exist('name_sig','var'), name_sig='missing'; end

datafolder = ['EISCAT_' year '-' month '-' day '_' name_expr '_V' intper_mean_str '@' name_ant];
storepath = fullfile(datapath,datafolder);

while exist(storepath)
    endnum = str2double(datafolder(end));
    if isnan(endnum)
        datafolder = [datafolder(1:end) num2str(1)];
    else
        datafolder = [datafolder(1:end-1) num2str(endnum+1)];            
    end
    storepath = fullfile(datapath,datafolder);
end

% Hdf5File = sprintf('%s%s',datafolder,'.hdf5');
% MatFile = sprintf('%s%s',datafolder,'.mat');
Hdf5File = [datafolder '.hdf5'];
MatFile  = [datafolder '.mat'];
hdffilename = fullfile(storepath,Hdf5File);
matfilename = fullfile(storepath,MatFile);
EISCAThdf5file = {hdffilename};
Storepath = {storepath};
GuisdapParFile = fullfile(path_GUP,'matfiles','Guisdap_Parameters.xlsx'); % path to the .xlsx file

if exist(hdffilename)==2, delete(hdffilename); end

[~,text] = xlsread(GuisdapParFile);
parameters_list = text(:,1);    % list that includes all Guisdap parameters and keep their positions from the excel arc

matfile.metadata.header= text(1,1:7)';

% TAI time (leapseconds)
if exist('leaps','var')
    if length(unique(leaps)) == 1
        if isfield(matfile.data,'par0d')
            ll0 = length(matfile.data.par0d);
        else	
            ll0 = 0;
        end
        matfile.data.par0d(ll0+1) = leaps(1);  
        a = find(strcmp('leaps',parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        matfile.metadata.par0d(:,ll0+1) = info';
    else
        ll1 = length(matfile.data.par1d(1,:));
        matfile.data.par1d(:,ll1+1:ll1+2) = leaps;
        a = find(strcmp('leaps1',parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        matfile.metadata.par1d(:,ll1+1) = info';
        a = find(strcmp('leaps2',parameters_list)==1);
        [~,~,info] = xlsread(GuisdapParFile,1,['A' num2str(a) ':E' num2str(a)]);
        info(6:7) = {num2str(xlsread(GuisdapParFile,1,['F' num2str(a)])) num2str(xlsread(GuisdapParFile,1,['G' num2str(a)]))};
        matfile.metadata.par1d(:,ll1+2) = info';
    end
end

nn = 0;
if exist('name_expr','var'); nn = nn + 1;
    infoname(1) = {'name_expr'};
    infoname(2) = {name_expr};
    a = find(strcmp('name_expr',parameters_list)==1);                
    [~,~,infodesc] = xlsread(GuisdapParFile,1,['B' num2str(a)]);
    infoname(3) = infodesc;
    matfile.metadata.names(:,nn) = infoname';
end
if exist('name_ant','var'); nn = nn + 1; 
    infoname(1) = {'name_ant'};
    infoname(2) = {name_ant};
    infoname(3) = {'name of receiving antenna, or alternative code name for multistatic analyses (tro: Tromso, kst: combination of Kiruna, Sodankyla, Tromso, esr: Svalbard, esa: combination of mainland and Svalbard'};
    matfile.metadata.names(:,nn) = infoname';
end
if exist('name_sig','var'); nn = nn + 1; 
    infoname(1) = {'name_sig'};
    infoname(2) = {name_sig};
    a = find(strcmp('name_sig',parameters_list)==1); 
    [~,~,infodesc] = xlsread(GuisdapParFile,1,['B' num2str(a)]);
    infoname(3) = infodesc;
    matfile.metadata.names(:,nn) = infoname';   
end
if exist('name_exps','var'); nn = nn + 1; 
    infoname(1) = {'name_exps'};
    for i = 1:length(name_exps(:,1))
        if i == 1
            nameexps = name_exps(1,:);
        else
            nameexps = [nameexps ' ' name_exps(i,:)];
        end
    end
    infoname(2) = {nameexps};
    a = find(strcmp('name_expr',parameters_list)==1);                
    [~,~,infodesc] = xlsread(GuisdapParFile,1,['B' num2str(a)]);
    infoname(3) = infodesc;
    matfile.metadata.names(:,nn) = infoname';
end
if exist('name_ants','var'); nn = nn + 1; 
    infoname(1) = {'name_ants'};
    for i = 1:length(name_ants(:,1))
        if i == 1
            nameants = name_ants(1,:);
        else
            nameants = [nameants ' ' name_ants(i,:)];
        end
    end
    infoname(2) = {nameants};
    a = find(strcmp('name_ant',parameters_list)==1); 
    [~,~,infodesc] = xlsread(GuisdapParFile,1,['B' num2str(a)]);
    infoname(3) = infodesc;
    matfile.metadata.names(:,nn) = infoname';
end

% Software
matfile.metadata.software.software_link = {software};
matfile.metadata.software.EISCAThdf5_ver = {hdf5ver};
if ~isempty(level2_link)
    matfile.metadata.software.level2_links = {level2_link};
end
if exist('name_strategies','var')
    for i = 1:length(name_strategies(:,1))
        if i == 1
            namestrategies = name_strategies(1,:);
        else
            namestrategies = [namestrategies ' ' name_strategies(i,:)];
        end
    end
    matfile.metadata.software.strategies = {namestrategies};
end

if isfield(matfile.metadata,'par0d')
    aa = find(cellfun('isempty',matfile.metadata.par0d(6,:)));    matfile.metadata.par0d(6,aa)= {'0'};
    aa = find(cellfun('isempty',matfile.metadata.par0d(7,:)));    matfile.metadata.par0d(7,aa)= {'0'};
end
if isfield(matfile.metadata,'par1d')
    aa = find(cellfun('isempty',matfile.metadata.par1d(6,:)));    matfile.metadata.par1d(6,aa)= {'0'};
    aa = find(cellfun('isempty',matfile.metadata.par1d(7,:)));    matfile.metadata.par1d(7,aa)= {'0'};
end
if isfield(matfile.metadata,'par2d')
    aa = find(cellfun('isempty',matfile.metadata.par2d(6,:)));    matfile.metadata.par2d(6,aa)= {'0'};
    aa = find(cellfun('isempty',matfile.metadata.par2d(7,:)));    matfile.metadata.par2d(7,aa)= {'0'};
end
  
symbols = ['a':'z' 'A':'Z' '0':'9'];
strLength = 10;
nums = randi(numel(symbols),[1 strLength]);
randstr = symbols(nums);
PID = ['doi://eiscat.se/3a/' year month day hour minute second '/' randstr];

matfile.metadata.schemes.DataCite.Identifier = {PID};
matfile.metadata.schemes.DataCite.Creator = {name_ant};
matfile.metadata.schemes.DataCite.Title = {['EISCAT_' filename]};
matfile.metadata.schemes.DataCite.Publisher = {'EISCAT Scientific Association'};
matfile.metadata.schemes.DataCite.ResourceType.Dataset = {'Level 4a Derived ionospheric data'};
matfile.metadata.schemes.DataCite.Date.Collected = {[starttime '/' endtime]};

% Find the smallest box (4 corners and mid-point) to enclose the data, 
% or 1 or unique 2 points that describes the data.
% If area of convhull (or distance between 2 points) < 10-4 deg^2, 
% define all points as one (average)
% imag = 1 to plot the data and the corresponding box if there is a box
im = 0;
gg_sp = [Vpos(:,2) Vpos(:,1)];

[plonlat,PointInPol] = polygonpoints([gg_sp(:,2) gg_sp(:,1)],im);
matfile.metadata.schemes.DataCite.GeoLocation.PolygonLon = plonlat(:,1);
matfile.metadata.schemes.DataCite.GeoLocation.PolygonLat = plonlat(:,2);
if ~isempty(PointInPol)
    matfile.metadata.schemes.DataCite.GeoLocation.PointInPolygonLon = PointInPol(1);
    matfile.metadata.schemes.DataCite.GeoLocation.PointInPolygonLat = PointInPol(2);
end

% Delete any empty fields from the structure
sFields = fieldnames(matfile);
for sf = sFields.' 
    tFields = fieldnames(matfile.(char(sf)));
    for tf = tFields.'
        if isempty(matfile.(char(sf)).(char(tf)))
            matfile.(char(sf)) = rmfield(matfile.(char(sf)),char(tf));
        end
    end
end

mkdir(storepath);
%save(matfilename,'matfile')

% Generate an HDF5-file from the MAT-file
chunklim = 10;
sFields = fieldnames(matfile);
for sf = sFields.'
    group1 = ['/' char(sf)];
    tFields = fieldnames(matfile.(char(sf)));
    for tf = tFields.'
        if strcmp('data',char(sf)) && (strcmp('par0d',char(tf)) || strcmp('par1d',char(tf)) || strcmp('par2d',char(tf)) || strcmp('par2d_pp',char(tf)) || strcmp('acf',char(tf)) || strcmp('ace',char(tf)) || strcmp('lag',char(tf)) || strcmp('freq',char(tf)) || strcmp('spec',char(tf)) || strcmp('om',char(tf)))
            npar  = length(matfile.data.(char(tf))(1,:));
            ndata = length(matfile.data.(char(tf))(:,1));
            if ge(ndata,chunklim) && ge(npar,chunklim), csize = [chunklim chunklim];
            elseif ge(ndata,chunklim), csize = [chunklim npar];
            elseif ge(npar,chunklim), csize = [ndata chunklim];
            else csize = [ndata npar]; end    
            h5create(hdffilename,['/' char(sf) '/' char(tf)],size([matfile.(char(sf)).(char(tf))]),'ChunkSize',csize,'Deflate',9,'Datatype','single');
            h5write(hdffilename,['/' char(sf) '/' char(tf)],single(matfile.(char(sf)).(char(tf))));
        elseif strcmp('metadata',char(sf)) 
            if isstruct(matfile.(char(sf)).(char(tf)))
                group2 = [group1 '/' char(tf)];
                uFields = fieldnames(matfile.(char(sf)).(char(tf)));
                for uf = uFields.'
                    if isstruct(matfile.(char(sf)).(char(tf)).(char(uf)))
                        group3 = [group2 '/' char(uf)];
                        vFields = fieldnames(matfile.(char(sf)).(char(tf)).(char(uf)));
                        for vf = vFields.'
                            if isstruct(matfile.(char(sf)).(char(tf)).(char(uf)).(char(vf)))
                                group4 = [group3 '/' char(vf)];
                                wFields = fieldnames(matfile.(char(sf)).(char(tf)).(char(uf)).(char(vf)));
                                for wf = wFields.'
                                    strdata = matfile.(char(sf)).(char(tf)).(char(uf)).(char(vf)).(char(wf));
                                    dsname = char(wf);
                                    strds2hdf5(hdffilename,group4,dsname,strdata)
                                end
                            else
                                strdata = matfile.(char(sf)).(char(tf)).(char(uf)).(char(vf));
                                dsname = char(vf);
                                strds2hdf5(hdffilename,group3,dsname,strdata)
                            end
                        end
                    else
                        strdata = matfile.(char(sf)).(char(tf)).(char(uf));
                        dsname = char(uf);
                        strds2hdf5(hdffilename,group2,dsname,strdata)
                    end
                end
            else
                strdata = matfile.(char(sf)).(char(tf));
                dsname = char(tf);
                strds2hdf5(hdffilename,group1,dsname,strdata)
            end
        else
            h5create(hdffilename,['/' char(sf) '/' char(tf)],size([matfile.(char(sf)).(char(tf))]));
            h5write(hdffilename,['/' char(sf) '/' char(tf)],matfile.(char(sf)).(char(tf)));
        end
    end   
end

matpath = fileparts(matfile_vel);
if addfigs
    image_filelist = [dir(fullfile(matpath,'*.png'));dir(fullfile(matpath,'*.pdf'))];
    npdf = 0;
    for ii = 1:length(image_filelist)
        figurefile = fullfile(matpath,image_filelist(ii).name);
        [~,filename,ext] = fileparts(figurefile);
        if strcmp(ext,'.png')
            image2hdf5(figurefile,hdffilename)
        elseif strcmp(ext,'.pdf')
            npdf = npdf + 1;
            pdf_forHDF5(npdf) = {[filename ext]};          
        end
    end
    if npdf>0
        strds2hdf5(hdffilename,'/figures','figure_links',pdf_forHDF5');
    end
end

if addnotes
    notesfiles = dir(fullfile(matpath,'notes*txt'));
    for nn = 1:length(notesfiles)
        notesfile = fullfile(matpath,notesfiles(nn).name);
        note2hdf5(notesfile,EISCAThdf5file{1},nn)
    end
end
