function [jasondata]=jasonimport(DataType,path,varargin)
%   [jasondata]=jasonimport()
%   Imports NMR data in JASON format
%
%   DataType specifies wheter it is raw or processed data.
%           DataType == 0 raw data
%                    == 1 processed, 
%                    == 2 raw data for an array
%                    == 3 processed data for an array
%   Path specifies wheter it is raw or processed data.
%           DataType == 0 single experiment
%                    == string raw data for an array, 
%
%   Checking optional inputs for jasonimport:
%          varargin{1} = expnr     = Checking the number of experiments 
%          varargin{2} = filename  = path of jjh5 file
%          varargin{3} = procpar   = general data inside the jjh5 file 
%
%   Point to the fid directory that contains the raw data.
%   The imported data will be returned in the structure variandata
%   containing the following members:
%               at: acquisition time (in seconds)
%               DELTA: diffusion time
%               delta: diffusion encoding time
%               dosyconstant: gamma.^2*delta^2*DELTAprime
%               FID: Free induction decays
%               filename: original file name and path
%               gamma: magnetogyric ratio of the nucleus
%               Gzlvl: gradient strengths
%               lp: first-order phase
%               ngrad: number of gradient levels
%               np: number of complex data points per gradient level
%               sfrq: spectrometer frequency (im MHz)
%               rp: zero-order phase
%               sp: start of spectrum (in ppm)
%               sp(n e.g 2,3,4): start of spectrum in n(nd) indirect dimension (in ppm)
%               sw: spectral width (in ppm)
%
%   See also: GNAT, dosy_mn, score_mn, decra_mn, mcr_mn, varianimport,
%             brukerimport, jeolimport, peakpick_mn, dosyplot_mn,
%             dosyresidual, dosyplot_gui, scoreplot_mn, decraplot_mn,
%             mcrplot_mn
%
%   This is a part of the GNAT
%   Copyright 2023  <Mathias Nilsson>

%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.

%   This program 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 General Public License along
%   with this program; if not, write to the Free Software Foundation, Inc.,
%   51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

%   Dr. Mathias Nilsson
%   School of Chemistry, University of Manchester,
%   Oxford Road, Manchester M13 9PL, UK
%   Telephone: +44 (0) 161 306 4465
%   Fax: +44 (0)161 275 4598
%   mathias.nilsson@manchester.ac.uk

%% Initialisation
% Extracting the information present in the jjh5 file
if nargin>3 % for arrays saved in one single file (same vendor)
    expnr = varargin{1};
    jasondata.filename = path;
    procpar = h5info(jasondata.filename); 
    clear varargin
elseif ~isnumeric(path) % for arrays saved in multiple file (same vendor)
    jasondata.filename=path;
    procpar = h5info(jasondata.filename);      
else % for single spectrum 
    [file, path]=uigetfile('*.jjh5','Choose the JASON header file (*.jjh5)');
    jasondata.filename=[path file];
    jasondata.filepath=path;
    if isnumeric(jasondata.filepath)
        return
    end
    procpar = h5info(jasondata.filename);    
end

% Filtering information useful for GNAT about the acquisition
% ('/JasonDocument/NMR/NMRData/n_experiment')
if exist('expnr','var') % for arrays in one file
    procpar = procpar.Groups.Groups(2).Groups.Groups(expnr);
else % for a single spectrum
    procpar = procpar.Groups.Groups(2).Groups.Groups;
    expnr=1;
end

% Check the type of data saved inside JASON file selected (At the moment there is no differentiation between thse methods)
if DataType==0 %raw data
    % Checking if inside JASON the imaginary component is saved 
    temp = h5info(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints']);
    sz = size(temp.Datasets,1);
    % Extracting the real part of the spectra
    real = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints/0']);
    if sz == 2 % real and imaginary data available
        % Extracting the imaginary part of the spectra
        img = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints/1']); 
    else % just the real data
        % To obtaing the imaginary part, we're going to use the Hilbert transformation 
        % see also [2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” 
        % IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.
        img = imag(hilbert(real));
    end   
    
elseif DataType==1 %processed data    
    % Checking if inside JASON the imaginary component is saved 
    temp = h5info(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints']);
    sz = size(temp.Datasets,1);
    % Extracting the real part of the spectra
    real = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints/0']);
    if sz == 2 % real and imaginary data available
        % Extracting the imaginary part of the spectra
        img = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints/1']); 
    else % just the real data
        % To obtaing the imaginary part, we're going to use the Hilbert transformation 
        % see also [2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” 
        % IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.
        img = imag(hilbert(real));
    end   

elseif DataType==2 %raw data for a large array
    % Checking if inside JASON the imaginary component is saved 
    temp = h5info(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints']);
    sz = size(temp.Datasets,1);
    % Extracting the real part of the spectra
    real = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints/0']);
    if sz == 2 % real and imaginary data available
        % Extracting the imaginary part of the spectra
        img = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/DataPoints/1']); 
    else % just the real data
        % To obtaing the imaginary part, we're going to use the Hilbert transformation 
        % see also [2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” 
        % IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.
        img = imag(hilbert(real));
    end   

elseif DataType==3 %processed data for a large array
    % Checking if inside JASON the imaginary component is saved 
    temp = h5info(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints']);
    sz = size(temp.Datasets,1);
    % Extracting the real part of the spectra
    real = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints/0']);
    if sz == 2 % real and imaginary data available
        % Extracting the imaginary part of the spectra
        img = h5read(jasondata.filename,['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/DataPoints/1']); 
    else % just the real data
        % To obtaing the imaginary part, we're going to use the Hilbert transformation 
        % see also [2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” 
        % IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.
        img = imag(hilbert(real));
    end 
else
    error('Unknown DataType')
end

% Extracting from .jjh5 the DOSY folder from specinfo
specinfo = procpar.Groups(6).Groups(2).Groups;  
p = find(matches({specinfo.Name}, ['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/SpecInfo/DosyParameters']));
if ~isempty(p) %The DosyParameters folder exists
    diffparameter = specinfo(p).Attributes;
else
    %I will exclusively extrate all the DOSY parameters from DosyParameters folder  
end

% Extracting lists (diffusion lists)
if sum(matches({procpar.Groups(6).Groups(2).Groups.Name},['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/SpecInfo/lists'])) == 1
    lists = procpar.Groups(6).Groups(2).Groups(4).Attributes.Value;
end

% set some default parameters
jasondata.procpar = procpar;
jasondata.rp=0;
jasondata.lp=0;
jasondata.sp=0;
jasondata.sp1=[];
jasondata.ni=[];
jasondata.logfile={'jasonimport'};

% Checking the data dimension
if isfield(procpar.Attributes,'Name')
    p = find(contains({procpar.Attributes.Name}, 'nDim'));
    DataDim = procpar.Attributes(p).Value;
else
    errordlg( {'Can not find the dimention information'},'File error')
    uiwait;
    jasondata.nDim_error=1;
    return
end

%% Read in acquisition parameters
specinfo = procpar.Groups(6).Groups(2).Attributes(1:end);
if DataDim>=1 % (I am only focusing on 1D parameters at the moment)
    p = find(matches({specinfo.Name}, {'PulseProgram','SW','SpectrometerFrequencies','appliedPh0','X'}));
    jasondata.pulprog = specinfo(p(1)).Value;
    logtext=['pulse programme is: ' jasondata.pulprog];
    disp(logtext)
    jasondata.logfile{end +1}=logtext;
    sw = specinfo(p(2)).Value;
    jasondata.sfrq = specinfo(p(3)).Value(1); %(in MHz)   
    jasondata.sw=sw(1)/jasondata.sfrq; %(in ppm)
    jasondata.rp=specinfo(p(4)).Value(1); 
    %---------------------------------------------------------------------%
    if DataType == 0 || DataType == 2
        p = find(matches({procpar.Groups(6).Attributes.Name}, 'Length'));
        jasondata.np=double(procpar.Groups(6).Attributes(p).Value(1));
    elseif DataType == 1 || DataType == 3
        p = find(matches({procpar.Attributes.Name}, 'Length'));
        jasondata.np=double(procpar.Attributes(p).Value(1));
    end
    p = find(matches({specinfo.Name}, {'SpectrumRef'}));
    jasondata.at=jasondata.np/(sw(1));
    jasondata.x_offset=(specinfo(p).Value(1))/jasondata.sfrq; % (in ppm)    
    jasondata.sp=jasondata.x_offset-jasondata.sw/2;
end
%% Read Processing parameters
% At the moment I am just extracting Raw data

%% Sort out parameters in GNAT format

if DataDim==1
    logtext='1D data';
    jasondata.arraydim=1;
    jasondata.ngrad=1;
    jasondata.Gzlvl='non existing';
    jasondata.DELTA='non existing';
    jasondata.delta='non existing';
    jasondata.dosyconstant='non existing';   
elseif DataDim==2
    logtext='2D data';
%     sw1 = specinfo(p(2)).Value(2); %(in Hz)
%     jasondata.sw1=sw1/jasondata.sfrq; %(in ppm)
elseif DataDim==3
    logtext='3D data';
%     sw1 = specinfo(p(2)).Value(3); %(in Hz)
%     jasondata.sw1=sw1/jasondata.sfrq; %(in ppm)
elseif DataDim==4
    logtext='4D data';
%     sw1 = specinfo(p(2)).Value(4); %(in Hz)
%     jasondata.sw1=sw1/jasondata.sfrq; %(in ppm)
else
    error('Unknown data dimension')
end
disp(logtext)
jasondata.logfile{end +1}=logtext;

% Extraction of the nucleis present in the Specinfo variable
p = find(matches({procpar.Groups(6).Groups(2).Groups.Name}, ['/JasonDocument/NMR/NMRData/',num2str(expnr-1),'/RawData/SpecInfo/Nucleides']));
nucleides = procpar.Groups(6).Groups(2).Groups(p).Groups; 
i=0;
while i<=8 
    i=i+1;
    temp=nucleides(i).Attributes(2).Value;
    if isempty(temp) == 1
        break
    else
        nucl(i)=temp;
    end
end

switch char(nucl(1))
    case 'H'
        jasondata.gamma=267524618.573;
    case 'C'
        jasondata.gamma=67283000;
    case 'N' % create condition for 14N
        jasondata.gamma=19338000;
%     case 'N' % create condition for 15N
%         jasondata.gamma=-27126000;
    case 'O'
        jasondata.gamma=-36281000;
    case 'F'
        jasondata.gamma=251815000;
    otherwise
        logtext='unknown nucleus - defaulting to proton';
        disp(logtext)
        jasondata.logfile{end +1}=logtext;
        jasondata.gamma=267524618.573;
end

%% Allocating DOSY parameters

% Peter, this is a temporary solution for variables not available in
% DosyParameters. (GRAD_SHAPE in JEOL and GPNAM in Bruker)
p = find(matches({specinfo.Name}, {'OrigFileFormat.str'}));
FileFormat = specinfo(p).Value;
% if isequal(FileFormat,'JCAMP-DX') || isequal(FileFormat,'JEOL Delta')
%     parameter = procpar.Groups(6).Groups(2).Groups(3).Groups(4).Attributes; 
% elseif isequal(FileFormat,'TopSpin')
%     parameter = procpar.Groups(6).Groups(2).Groups(2).Groups(2).Attributes;
% end
jasondata.FileFormat = FileFormat;

% For this part, the variables are extracted from DosyParameters folder
if exist([jasondata.filename 'GNATimport'],'file') == 2
    logtext='Using GNATimport file to determine import parameters';
    disp(logtext)
    jasondata.logfile{end +1}=logtext;
    fileid=fopen([jasondata.filename  'GNATimport'],'rt'); %#ok<NASGU>
else
    %Using default import options
    logtext=['cannot open file: ' [jasondata.filename  'GNATimport']];
    disp(logtext)
    jasondata.logfile{end +1}=logtext;
    logtext='Using default import (normal for 1D, DOSY and Relaxation data';
    disp(logtext)
    jasondata.logfile{end +1}=logtext;
    if DataDim==1
        % Do nothing
    elseif DataDim == 2 %likely a difusion or relaxation experiment
        logtext='2D experiment: trying diffusion and relaxation type';
        disp(logtext)
        jasondata.logfile{end +1}=logtext;
        % Default result for GNAT
        jasondata.Gzlvl='non existing';
        jasondata.DELTA='non existing';
        jasondata.delta='non existing';
        jasondata.dosyconstant='non existing';
        jasondata.ngrad=1;
        % Reading Diffusion experiment
        p = find(matches({specinfo.Name},'SpectrumTypeExt'));
        if specinfo(p).Value(2) == 6        
            logtext='Diffusion experiment';
            disp(logtext)
            jasondata.logfile{end +1}=logtext;
            jasondata.Gzlvl=(lists*0.001)'; % conversion from G/cm to T/m
            jasondata.ngrad=length(jasondata.Gzlvl);

            % set the default as 'SINE'
            jasondata.gpnam='SINE';
            % temporary solution for diffusion
            if sum(matches({parameter.Name}, 'GPNAM')) == 1 % Bruker 
                p = find(matches({parameter.Name}, 'GPNAM'));
                jasondata.gpnam=cell2str(parameter(p).Value(2));
            elseif sum(matches({parameter.Name}, 'GRAD_SHAPE')) == 1 % JEOL
                p = find(matches({parameter.Name}, 'GRAD_SHAPE'));
                jasondata.gpnam=parameter(p).Value;
            else
                logtext='DOSY gradient shape missing. I will set it to default (SINE)';
                disp(logtext)
                jasondata.logfile{end +1}=logtext;
            end
            logtext=(['DOSY gradient shape is: ' jasondata.gpnam]);
            disp(logtext)
            jasondata.logfile{end +1}=logtext;
            logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient';
            disp(logtext)
            jasondata.logfile{end +1}=logtext;
            
            % Extracting variables to the Stejskal-Tanner Equation
            p = find(matches({diffparameter.Name}, {'InterGradientDelay','DiffusionTime_nominal','GradientPulseDuration','OneShotAlpha'}));
            jasondata.tau=diffparameter(p(1)).Value(1);
            jasondata.DELTA=diffparameter(p(2)).Value(2);
            jasondata.delta=diffparameter(p(3)).Value(1);
            jasondata.alpha=diffparameter(p(4)).Value(1);

            % Set the ST equation values according to the gradient shape
            logtext=('Gradient shape will be used for calculation of the attenuation function');
            disp(logtext)
            jasondata.logfile{end +1}=logtext;
            if strcmpi(jasondata.gpnam,"SINE") || strcmpi(jasondata.gpnam,"<SINE>")
                jasondata.dosyconstant=jasondata.gamma.^2*jasondata.delta.^2*...
                    (jasondata.DELTA-(5-3*jasondata.alpha.^2)*(jasondata.delta/16)-(1-jasondata.alpha.^2)*(jasondata.tau/2));
            elseif strcmpi(jasondata.gpnam,"SMSQ") || strcmpi(jasondata.gpnam,"<SMSQ10.100>") || strcmpi(jasondata.gpnam,"SQUARE")
                jasondata.dosyconstant=jasondata.gamma.^2*jasondata.delta.^2*...
                    (jasondata.DELTA-jasondata.delta*((3367*pi.^2-207)*jasondata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                    (jasondata.tau/2)*(1-jasondata.alpha.^2));
            elseif strcmpi(jasondata.gpnam,"RECT") || strcmpi(jasondata.gpnam,"<RECT>")
                jasondata.dosyconstant=jasondata.gamma.^2*jasondata.delta.^2*...
                    (jasondata.DELTA-(2-jasondata.alpha.^2)*(jasondata.delta/6)-(1-jasondata.alpha.^2)*(jasondata.tau/2));
            end
        else
            logtext='No recognizable variable with increment values exist: importing empty arrays';
            disp(logtext)
            jasondata.logfile{end +1}=logtext;
        end
    end
end
        
%% Stracting the Spectra

compfid=complex(real,img);
% Checking the presence of the Group Delay
if sum(matches({specinfo.Name},'GroupDelay')) >= 1
    p = find(contains({specinfo.Name}, {'GroupDelay'}));
    jasondata.GroupDelay = specinfo(p(1)).Value;
    jasondata.lp = -2.0 * 3.14159265358979323846 * jasondata.GroupDelay;
    jasondata.lrfid = jasondata.GroupDelay;
    jasondata.rp=1.5157*-1 - jasondata.lp;
else
    jasondata.lrfid = 0;
end

jasondata.FID=compfid;

end