function [brukerdata]=brukerimport(DataType, inpath)
%   [brukerdata]=brukerimport(DataType, inpath)
%   Imports PFG-NMR data in Bruker format
%               DataType specifies wheter it is raw or Processed data.
%               DataType=0 is raw data, 1 processed, 2 is processed
%               data for an array, and 3 is raw data for a array
%
%   Useage: Point to the fid/ser file that contains the raw data. The
%           imported data will be returned in the structure brukerdata.
%
%   Example:
%   See also: DOSYToolbox, 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 GNAT.
%   Copyright 2017  <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

%% Changes
% MN      26 June    2017   Handling import and PSYCHEiDOSY experiment,
%                           given correct name for sequence
% MN      17 July    2017   Handling import and PSYCHEiDOSY experiment for
%                           both 2D and 3D versions
% GDP     15 August  2017   Importing 3D corrected for most of our experiments
% MN      23 October 2017   Added errordlg for import error
% MN      09 January 2018   Fixed and error calculating brukerdata.at for
%                           processed data
% GDP     25 October 2018   Importing 4D experiments
% MN      24 April   2019   Added support for bruker convection compensated
%                           pulse sequence dstebpgp
% GDP     26 April   2019   Added support for pure array import (DataType=3)
% MN      11 July    2019   Partial fix for too short pulse sequence names
% MJS     01 March   2021   Fix of import of pulprog (413)
% GDP     08 Feb     2021   Added DEPT/INEPT options for BPP and Oneshot
% GDP     02 March   2021   Added import for GPNAM6
% MN      12 May     2021   Added ST calculations according to gradient
%                           shape (GPNAM) - still needs doing for higher order experiments
% GDP     13 May     2021   Modifications for higher order experiments
% GDP     19 May     2021   House keeping on the code
% MN 29 March correcting chunkSize to accomodate for Neo. Thanks to Goran
% Carlstrom (Lund) for this fix 1712
%MN       17 April   2023   Sorted scaling problems by remembering to
%                           multuply FID data with the NC parameter
%MN       1  November2023   Sorted imprt array of vd and vclist
%MN       25 April   2024   Fixed import of processed data so that the DSP
%                           is taken into account and so it works for
%                           Reference Deconvolution
%MN       16 May   2024     Fixed import of 2D processed data so that the DSP
%                           is taken into account and so it works for
%                           Reference Deconvolution


%% Initialisation
%Check the OS so the paths gets written the right way
if isunix==1
    OSslash='/';
else
    OSslash='\';
end
% FilePath=[FilePath OSslash];
%a=b
if nargin==2
    if DataType==0 %raw data
        %inpath should point to the directory containing the ser/fid file
        FilePath=[inpath OSslash];
        ProcPathRoot=[FilePath 'pdata' OSslash '1' OSslash];
    elseif DataType==1 %processed data
        %inpath should point to the directory containing the 1r/2rr file
        ProcPathRoot=[inpath OSslash];
        StrInd=strfind(inpath,'pdata');
        if isempty(StrInd)
            warndlg('Please select the directory containing the processed data (e.g 1r, 2rr)','Import Bruker Processed Data')
            uiwait(gcf)
        end
        FilePath=inpath(1:(StrInd-1));
    elseif DataType==2 %processed data for a large array
        FilePath=[inpath OSslash];
        ProcPathRoot=[FilePath 'pdata' OSslash '1' OSslash];
    elseif DataType==3 %raw data for a large array (e.g. Vladmir's data)
        %inpath should point to the directory containing the ser/fid file
        FilePath=[inpath OSslash];
        ProcPathRoot=[FilePath 'pdata' OSslash '1' OSslash];
    else
        error('Unknown DataType')
    end
else
    error('brukerimport takes exactly 2 arguments')
end
%% Main
if FilePath
    %% set some default parameters
    brukerdata.sp1=[];
    brukerdata.ni=[];
    brukerdata.logfile={'brukerimport'};
    %% Check data dimension
    if exist([FilePath OSslash 'acqu5s'], 'file') == 2
        error('Data of order >4 is not supported')
    elseif exist([FilePath OSslash 'acqu4s'], 'file') == 2
        DataDim=4;
    elseif exist([FilePath OSslash 'acqu3s'], 'file') == 2
        DataDim=3;
    elseif exist([FilePath OSslash 'acqu2s'], 'file') == 2
        DataDim=2;
    elseif exist([FilePath OSslash 'acqus'], 'file') == 2
        DataDim=1;
    else
        errordlg( {'Can not find the acqus file','It is normally found in the numbered data directory'},'File error')
        uiwait;
        brukerdata.acqus_error=1;
        return
    end
    %% Read in acquisition parameters
    if DataDim>=1
        brukerdata.acqus = GetParameters([FilePath  'acqus']);
        brukerdata.acqus.PULPROG(1)=[];
        brukerdata.acqus.PULPROG(end)=[];
        brukerdata.pulprog=brukerdata.acqus.PULPROG;
        logtext=['pulse programme is: ' brukerdata.pulprog];
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        brukerdata.sw=str2double(brukerdata.acqus.SW);
    end
    
    if DataDim>=2
        brukerdata.acqu2s = GetParameters([FilePath  'acqu2s']);
    end
    
    if DataDim>=3
        brukerdata.acqu3s = GetParameters([FilePath  'acqu3s']);
    end
    
    if DataDim>=4
        brukerdata.acqu4s = GetParameters([FilePath  'acqu4s']);
    end
    %% Read Processing parameters
    %get processing parameters, if there are any
    if DataDim>=1
        Procpath=[ProcPathRoot 'procs'];
        if exist(Procpath, 'file') == 2
            brukerdata.procs = GetParameters(Procpath);
            brukerdata.sp=str2double(brukerdata.procs.OFFSET)-brukerdata.sw;
            %brukerdata.sp=1
            brukerdata.lp=str2double(brukerdata.procs.PHC1);
            brukerdata.rp=str2double(brukerdata.procs.PHC0)*-1 - brukerdata.lp;
        else
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='estimating processing parameters from acquisition parameters';
            brukerdata.logfile{end +1}=logtext;
            brukerdata.sf=str2double(brukerdata.acqus.BF1);
            brukerdata.sfo1=str2double(brukerdata.acqus.SFO1);
            brukerdata.sp=(brukerdata.sfo1/brukerdata.sf-1)*1e6+0.5*brukerdata.sw*(brukerdata.sfo1/brukerdata.sf)-brukerdata.sw;
            brukerdata.lp=0;
            brukerdata.rp=0;
            brukerdata.procs=[];
        end
    end
    
    if DataDim>=2
        Procpath=[ProcPathRoot 'proc2s'];
        if exist(Procpath, 'file') == 2
            brukerdata.proc2s = GetParameters(Procpath);
            
        else
            
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='No processing paramaters for the second dimension';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.proc2s=[];
        end
    end
    
    if DataDim>=3
        Procpath=[ProcPathRoot 'proc3s'];
        if exist(Procpath, 'file') == 2
            brukerdata.proc3s = GetParameters(Procpath);
            
        else
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='No processing paramaters for the third dimension';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.proc3s=[];
        end
    end
    
    if DataDim>=4
        Procpath=[ProcPathRoot 'proc4s'];
        if exist(Procpath, 'file') == 2
            brukerdata.proc4s = GetParameters(Procpath);
            
        else
            logtext=['cannot open file: ' Procpath];
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            logtext='No processing paramaters for the forth dimension';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.proc4s=[];
        end
    end
    %% Sort out parameters in DOSY Toolbox format
    if DataDim==1
        logtext='1D data';
        brukerdata.arraydim=1;
        brukerdata.ngrad=1;
        fidpath=[FilePath 'fid'];
        brukerdata.acqu2s=[];
        brukerdata.acqu3s=[];
    elseif DataDim==2
        logtext='2D data';
        brukerdata.acqu3s=[];
        brukerdata.arraydim=str2num(brukerdata.acqu2s.TD); %#ok<*ST2NM>
        fidpath=[FilePath 'ser'];
        brukerdata.sw1=str2double(brukerdata.acqu2s.SW);
    elseif DataDim==3
        logtext='3D data';
        brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
        fidpath=[FilePath 'ser'];
        brukerdata.sw1=str2double(brukerdata.acqu3s.SW);
    elseif DataDim==4
        logtext='4D data';
        brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD).*str2num(brukerdata.acqu4s.TD);
        fidpath=[FilePath 'ser'];
        brukerdata.sw1=str2double(brukerdata.acqu4s.SW);
    else
        error('Unknown data dimension')
    end
    disp(logtext)
    brukerdata.logfile{end +1}=logtext;
    %Get the parameters I need for now (in units I like)
    brukerdata.filename=FilePath(1:end-1);
    brukerdata.np=str2double(brukerdata.acqus.TD)/2;
    brukerdata.sfrq=str2double(brukerdata.acqus.SFO1);
    brukerdata.at=brukerdata.np/(brukerdata.sw*brukerdata.sfrq);
    
    if DataType==1
        if isempty(brukerdata.procs)
            error('cannot import processed data if the file "procs" is missing')
        end
        brukerdata.SI=round(str2double(brukerdata.procs.SI));
        brukerdata.at=0.5*brukerdata.at/(brukerdata.np/brukerdata.SI);
    elseif DataType==2
        if isempty(brukerdata.procs)
            error('cannot import processed data if the file "procs" is missing')
        end
        brukerdata.SI=round(str2double(brukerdata.procs.SI));
        brukerdata.at=0.5*brukerdata.at/(brukerdata.np/brukerdata.SI);
    
    end
    
    
    switch brukerdata.acqus.NUC1
        case '<1H>'
            brukerdata.gamma=267524618.573;
        case '<2H>'
            brukerdata.gamma=41066000;
        case '<10B>'
            brukerdata.gamma=28747000;
        case '<11B>'
            brukerdata.gamma=85847000;
        case '<13C>'
            brukerdata.gamma=67283000;
        case '<14N>'
            brukerdata.gamma=19338000;
        case '<15N>'
            brukerdata.gamma=-27126000;
        case '<17O>'
            brukerdata.gamma=-36281000;
        case '<19F>'
            brukerdata.gamma=251815000;
        case '<23Na>'
            brukerdata.gamma=70808000;
        case '<27Al>'
            brukerdata.gamma=69764000;
        case '<29Si>'
            brukerdata.gamma=-53190000;
        case '<31P>'
            brukerdata.gamma=108394000;
            
        otherwise
            logtext='unknown nucleus - defaulting to proton';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.gamma=267524618.573;
    end
    
    
    %% Import lists
    
    % Checking for a GNATimport list
    if exist([FilePath 'GNATimport'],'file') == 2
        % will use GNATimport file to determine how to import the data
        logtext='Using GNATimport file to determine import parameters';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        fileid=fopen([FilePath  'GNATimport'],'rt'); %#ok<NASGU>
    else
        %Using default import options
        logtext=['cannot open file: ' [FilePath  'GNATimport']];
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        logtext='Using default import (normal for 1D, DOSY and Relaxation data';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        if DataDim==1
            logtext='1D experiment';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            brukerdata.ngrad=1;
            brukerdata.arraydim=1;

            if exist([FilePath 'difflist'],'file') == 2 || exist([FilePath 'difflist.txt'],'file') == 2 %Diffusion experiment                
                %Lets get the proper value of TD for 1D arrays
                if exist([FilePath 'difflist.txt'],'file') == 2
                    G=fopen([FilePath 'difflist.txt'],'rt');
                elseif exist([FilePath 'difflist'],'file') == 2
                    G=fopen([FilePath 'difflist'],'rt');
                end                
                TD=1;
                while ~feof(G)
                   textline=fgetl(G);
                   if isempty(textline) || length(textline)==0 %#ok<ISMT>
                      continue;
                   end 
                   A(TD)=str2double(textline); %#ok<AGROW>
                   TD=TD+1;
                end 
                TD=TD-1; %#ok<NASGU>
                              
                logtext='Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
%                 brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
%                 brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                brukerdata.Gzlvl=A*0.01; %conversion to T/m
                brukerdata.ngrad=length(brukerdata.Gzlvl);
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
            elseif exist([FilePath 'vdlist'],'file') == 2 %relaxation experiment
                 %Lets get the proper value of TD for 1D arrays

                 disp('hi')
                 G=fopen([FilePath 'vdlist'],'rt');
                TD=1;
                while ~feof(G)
                    textline=fgetl(G);
                    if isempty(textline) || length(textline)==0 %#ok<ISMT>
                        continue;
                    end 
                   A(TD)=str2double(textline); %#ok<AGROW>
                    TD=TD+1;
                end 
                TD=TD-1;
                logtext='T1 relaxation experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD);
                logtext='importing vdlist';
                disp(logtext)
                %                 brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
                %                 brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                brukerdata.Gzlvl=A; %conversion to T/m
                brukerdata.ngrad=length(brukerdata.Gzlvl);
                logtext=['Relaxation pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                               
              elseif exist([FilePath 'vclist'],'file') == 2 %relaxation experiment
                 %Lets get the proper value of TD for 1D arrays
                 G=fopen([FilePath 'vclist'],'rt');
                TD=1;
                while ~feof(G)
                 textline=fgetl(G);

                if isempty(textline) || length(textline)==0 %#ok<ISMT>
                  continue;
                end 
                   A(TD)=str2double(textline); %#ok<AGROW>
                TD=TD+1;
                end 
                TD=TD-1;
                
               
                logtext='T2 relaxation experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                 brukerdata.vdlist = ImportList([FilePath 'vclist'], TD);
                logtext='importing vclist';
                disp(logtext)
%                 brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
%                 brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                brukerdata.Gzlvl=A; %conversion to T/m
                brukerdata.ngrad=length(brukerdata.Gzlvl);
                logtext=['Relaxation pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        elseif DataDim == 2 %likely a difusion or relaxation experiment
            logtext='2D experiment: trying diffusion and relaxation type';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            brukerdata.ngrad=1;
            TD= str2num(brukerdata.acqu2s.TD);
            if exist([FilePath 'nuslist'],'file') == 2 %NUS file
                logtext='NUS acquired data (nuslist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.nuslist=ImportList([FilePath 'nuslist'],TD);
                brukerdata.nuslist(isnan(brukerdata.nuslist))=0;
                brukerdata.nuslist=brukerdata.nuslist(1:1+(nnz(brukerdata.nuslist)));
                brukerdata.nuslist=brukerdata.nuslist';
            end
            if exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment       
                logtext='Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
                brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                brukerdata.ngrad=length(brukerdata.Gzlvl);
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                % Check for gradient shape GDP 08iii2021
                %MN 12 MAy 2021 - trying to make this more logic
                
                % set the default as 'SINE'
                brukerdata.gpnam='SINE';
                if isfield(brukerdata.acqus,'GPNAM6')
                    if contains(brukerdata.acqus.GPNAM6,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(brukerdata.acqus.GPNAM6,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(brukerdata.acqus.GPNAM6,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                else
                    logtext='Reading gradient shape from the acqus';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    [tmp]=GetGPNAM([FilePath  'acqus']);
                    if contains(tmp,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(tmp,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(tmp,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                    
                end
                if isempty(brukerdata.gpnam) || strcmp(brukerdata.gpnam,'<>')
                    logtext='DOSY gradient shape missing. I will set it to default';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.gpnam='SINE';
                end
                logtext=(['DOSY gradient shape is: ' brukerdata.gpnam]);
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient, which is done automatically with the TopSpin "dosy" setup macro';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                %check for pulseprogram
                %MN 12 May 2021 Using the correct definitions for calculation of dosyconstant and tau
                
                if strncmpi('Doneshot',brukerdata.pulprog,8)
                    logtext='Doneshot sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.alpha=brukerdata.acqus.CNST(15);
                    %MN use the actual values ranther than CNST18
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(5-3*brukerdata.alpha.^2)*(brukerdata.delta/16)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-brukerdata.delta*((3367*pi.^2-207)*brukerdata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                            (brukerdata.tau/2)*(1-brukerdata.alpha.^2));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(2-brukerdata.alpha.^2)*(brukerdata.delta/6)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    end
                    
                    
                elseif contains(brukerdata.pulprog,'psycheidosy','IgnoreCase',true) || contains(brukerdata.pulprog,'psyche_idosy','IgnoreCase',true)
                    logtext='psyche-idosy sequence; setting paramaters to monopolar';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/4);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(4860*pi.^2)));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    end
                    
                elseif strfind(brukerdata.pulprog,'oneshot')
                    logtext='oneshot sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    if  ~strcmp('<1H>',brukerdata.acqus.NUC1)
                        if strfind(brukerdata.pulprog,'dept')  %If encoding is at 1H
                            if strfind(brukerdata.pulprog,'oneshot') < strfind(brukerdata.pulprog,'dept')
                                logtext='Using DEPT before acquisition - I will assume diffusion was encoded in 1H';
                                disp(logtext)
                                brukerdata.logfile{end +1}=logtext;
                                brukerdata.gamma=267524618.573;
                            else
                                %encoding on X nuclei
                            end
                        elseif strfind(brukerdata.pulprog,'inept')
                            if strfind(brukerdata.pulprog,'oneshot') < strfind(brukerdata.pulprog,'inept')
                                logtext='Using INEPT before acquisition - I will assume diffusion was encoded in 1H';
                                disp(logtext)
                                brukerdata.logfile{end +1}=logtext;
                                brukerdata.gamma=267524618.573;
                            else
                                %encoding on X nuclei
                            end
                        else
                            % X Oneshot
                        end
                    end
                    
                    brukerdata.alpha=brukerdata.acqus.CNST(15);
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(5-3*brukerdata.alpha.^2)*(brukerdata.delta/16)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-brukerdata.delta*((3367*pi.^2-207)*brukerdata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                            (brukerdata.tau/2)*(1-brukerdata.alpha.^2));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(2-brukerdata.alpha.^2)*(brukerdata.delta/6)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    end
                    
                elseif contains(brukerdata.pulprog,'dstebpgp','IgnoreCase',true)
                    logtext='double stimulated echo for convection compensation';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21)- 2*brukerdata.acqus.D(17) - 2*brukerdata.acqus.P(31)*1e-6 ;
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(2430*pi.^2)));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-2*brukerdata.delta/3);
                    end
                    
                elseif contains(brukerdata.pulprog,'stebp','IgnoreCase',true) || contains(brukerdata.pulprog,'ledbp','IgnoreCase',true)
                    %bipolar sequence
                    logtext='bipolar sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                    end
                    
                elseif contains(brukerdata.pulprog,'ste','IgnoreCase',true) || contains(brukerdata.pulprog,'led','IgnoreCase',true)
                    logtext='monopolar sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/4);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(4860*pi.^2)));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    end
                    
                elseif strfind(brukerdata.pulprog,'proj')
                    logtext='PROJECTED sequence';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    logtext=('Attenuation function calculated in the pulse sequence');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                elseif strfind(brukerdata.pulprog,'zg')
                    logtext='Looks like a simulation, I will just import the values';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA);
                elseif strfind(brukerdata.pulprog,'sel_idiff')
                    logtext='Bart`s iDOSY for selective experiments';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    logtext=('Attenuation function calculated in the pulse sequence');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;

                elseif (numel(strfind(brukerdata.pulprog,'dstebpgp')))>0 && (numel(brukerdata.pulprog)>7) %&& (strfind(brukerdata.pulprog,'dstebpgp',8))&&  (numel(brukerdata.pulprog)>7)
                    logtext='double stimulated echo for convection compensation';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21)-brukerdata.acqus.D(17)/2;
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(31)*1e-6;
                    brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*(brukerdata.DELTA-2*brukerdata.delta/3-2*brukerdata.tau/2);
                elseif strfind(brukerdata.pulprog,'bpp')
                    %bipolar sequence
                    logtext='bipolar sequence';
                    disp(logtext)
                    
                    if  ~strcmp('<1H>',brukerdata.acqus.NUC1)
                        if strfind(brukerdata.pulprog,'dept')  %If encoding is at 1H
                            if strfind(brukerdata.pulprog,'bpp') < strfind(brukerdata.pulprog,'dept')
                                logtext='Using DEPT before acquisition - I will assume diffusion was encoded in 1H';
                                disp(logtext)
                                brukerdata.logfile{end +1}=logtext;
                                brukerdata.gamma=267524618.573;
                            else
                                %encoding on X nuclei
                            end
                        elseif strfind(brukerdata.pulprog,'inept')
                            if strfind(brukerdata.pulprog,'bpp') < strfind(brukerdata.pulprog,'inept')
                                logtext='Using INEPT before acquisition - I will assume diffusion was encoded in 1H';
                                disp(logtext)
                                brukerdata.logfile{end +1}=logtext;
                                brukerdata.gamma=267524618.573;
                            else
                                %encoding on X nuclei
                            end
                        else
                            
                        end
                    end
                    
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21)-brukerdata.acqus.D(17)/2;
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(31)*1e-6;
                    
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                    end
                else
                    if brukerdata.acqus.CNST(19)~=0
                        logtext='Probably an UoM pulse sequence; I will import the constants';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.alpha=brukerdata.acqus.CNST(15);
                        brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                        brukerdata.tau=brukerdata.acqus.CNST(18);
                        brukerdata.DELTA=brukerdata.acqus.D(21);
                        brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    else
                        logtext='unknown sequence; setting parameters to bipolar';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.DELTA=brukerdata.acqus.D(21);
                        brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse;
                        logtext=('Gradient shape will be used for calculation of the attenuation function');
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        if strcmpi(brukerdata.gpnam,"SINE")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                        elseif strcmpi(brukerdata.gpnam,"SMSQ")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                        elseif strcmpi(brukerdata.gpnam,"RECT")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                        end
                    end
                end
            elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                logtext='T2 experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.vclist = ImportList([FilePath 'vclist'], TD);
                if exist([FilePath 'vdlist'],'file') == 2
                    brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD);
                    logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                else
                    logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                    uiwait(gcf)
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                end
            elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                logtext='T1 experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD);
                logtext='using vdlist directly for delay values';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            else
                logtext='No recognizable files with increment values exist: importing empty arrays';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
            %else %higher order experiments, GDP 17 Jan 2018
        elseif DataDim == 3 && DataType~=3 %higher order experiments, for PARAFAC
            logtext='3D Data. Probably for PARAFAC. I will try to sort experimental parameters'; disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
            brukerdata.ngrad=str2num(brukerdata.acqu2s.TD);
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.dosyconstant='non existing';
            if exist([FilePath 'nuslist'],'file') == 2 %NUS file
                logtext='NUS acquired data (nuslist file exists)';
                nusmax=1024; %All the other dimensions must be smaller than this
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.nuslist=ImportDoubleList([FilePath 'nuslist'],nusmax,DataDim);
                brukerdata.nuslist=brukerdata.nuslist(1:(find(isnan(brukerdata.nuslist(:,1))))-1,:);
                % GDP Only works for 100% NUS
                TD2=max(brukerdata.nuslist(:,1))+1;
                TD3=max(brukerdata.nuslist(:,2))+1;
                if exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                    logtext='3D Diffusion experiment (difflist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    % Check for gradient shape
                    % set the default as 'SINE'
                    brukerdata.gpnam='SINE';
                    if isfield(brukerdata.acqus,'GPNAM6')
                        if contains(brukerdata.acqus.GPNAM6,"SINE",'IgnoreCase',true)
                            brukerdata.gpnam='SINE';
                        elseif contains(brukerdata.acqus.GPNAM6,"SMSQ",'IgnoreCase',true)
                            brukerdata.gpnam='SMSQ';
                        elseif contains(brukerdata.acqus.GPNAM6,"RECT",'IgnoreCase',true)
                            brukerdata.gpnam='RECT';
                        end
                    else
                        logtext='Reading gradient shape from the acqus';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        [tmp]=GetGPNAM([FilePath  'acqus']);
                        if contains(tmp,"SINE",'IgnoreCase',true)
                            brukerdata.gpnam='SINE';
                        elseif contains(tmp,"SMSQ",'IgnoreCase',true)
                            brukerdata.gpnam='SMSQ';
                        elseif contains(tmp,"RECT",'IgnoreCase',true)
                            brukerdata.gpnam='RECT';
                        end
                    end
                    if isempty(brukerdata.gpnam) || strcmp(brukerdata.gpnam,'<>')
                        logtext='DOSY gradient shape missing. I will set it to default';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.gpnam='SINE';
                    end
                    logtext=(['DOSY gradient shape is: ' brukerdata.gpnam]);
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient, which is done automatically with the TopSpin "dosy" setup macro';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    if  contains(brukerdata.pulprog,'psycheidosy','IgnoreCase',true) || contains(brukerdata.pulprog,'psyche_idosy','IgnoreCase',true)
                        logtext='psyche-idosy sequence; setting paramaters to monopolar';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                        Gzlvl3(isnan(Gzlvl3))=0;
                        Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                        Gzlvl3(isnan(Gzlvl3))=0;
                        if nnz(Gzlvl3)>nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl3;
                        elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl2;
                        elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                            if (length(Gzlvl3)-nnz(Gzlvl3))==0
                                brukerdata.Gzlvl=Gzlvl3;
                            else
                                brukerdata.Gzlvl=Gzlvl2;
                            end
                        end
                        brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                        brukerdata.ngrad=length(brukerdata.Gzlvl);
                        brukerdata.DELTA=brukerdata.acqus.D(21);
                        brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                        
                        %Set the ST values according to the gradient shape
                        logtext=('Gradient shape will be used for calculation of the attenuation function');
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        if strcmpi(brukerdata.gpnam,"SINE")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/4);
                        elseif strcmpi(brukerdata.gpnam,"SMSQ")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(4860*pi.^2)));
                        elseif strcmpi(brukerdata.gpnam,"RECT")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                        end
                        brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
                        
                    elseif strfind(brukerdata.pulprog,'oneshot')
                        logtext='oneshot sequence';
                        disp(logtext)
                        Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                        Gzlvl3(isnan(Gzlvl3))=0;
                        Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                        Gzlvl3(isnan(Gzlvl3))=0;
                        if nnz(Gzlvl3)>nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl3;
                        elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl2;
                        elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                            if (length(Gzlvl3)-nnz(Gzlvl3))==0
                                brukerdata.Gzlvl=Gzlvl3;
                            else
                                brukerdata.Gzlvl=Gzlvl2;
                            end
                        end
                        brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                        brukerdata.ngrad=length(brukerdata.Gzlvl);
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.alpha=brukerdata.acqus.CNST(15);
                        brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                        brukerdata.DELTA=brukerdata.acqus.D(21);
                        brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                        
                        %Set the ST values according to the gradient shape
                        logtext=('Gradient shape will be used for calculation of the attenuation function');
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        if strcmpi(brukerdata.gpnam,"SINE")
                            brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                                (brukerdata.DELTA-(5-3*brukerdata.alpha.^2)*(brukerdata.delta/16)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                        elseif strcmpi(brukerdata.gpnam,"SMSQ")
                            brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                                (brukerdata.DELTA-brukerdata.delta*((3367*pi.^2-207)*brukerdata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                                (brukerdata.tau/2)*(1-brukerdata.alpha.^2));
                        elseif strcmpi(brukerdata.gpnam,"RECT")
                            brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                                (brukerdata.DELTA-(2-brukerdata.alpha.^2)*(brukerdata.delta/6)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                        end
                    else
                        logtext='unknown sequence; setting parameters to bipolar';
                        disp(logtext)
                        Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                        Gzlvl3(isnan(Gzlvl3))=0;
                        Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                        Gzlvl2(isnan(Gzlvl2))=0;
                        if nnz(Gzlvl3)>nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl3;
                        elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                            brukerdata.Gzlvl=Gzlvl2;
                        elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                            if (length(Gzlvl3)-nnz(Gzlvl3))==0
                                brukerdata.Gzlvl=Gzlvl3;
                            else
                                brukerdata.Gzlvl=Gzlvl2;
                            end
                        end
                        brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                        brukerdata.ngrad=length(brukerdata.Gzlvl);
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.DELTA=brukerdata.acqus.D(21);
                        brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse
                        brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                        logtext=('Gradient shape will be used for calculation of the attenuation function');
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        if strcmpi(brukerdata.gpnam,"SINE")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                        elseif strcmpi(brukerdata.gpnam,"SMSQ")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                        elseif strcmpi(brukerdata.gpnam,"RECT")
                            brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                        end
                    end
                    if exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                        logtext='T2 experiment (vclist file exists)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        TD2= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vclist=ImportList([FilePath 'vclist'], TD2);
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(3) + brukerdata.acqus.P(3)*1e-6);
                        if exist([FilePath 'vdlist'],'file') == 2
                            brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                            logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                            disp(logtext)
                            brukerdata.logfile{end +1}=logtext;
                        else
                            logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                            disp(logtext)
                            brukerdata.logfile{end +1}=logtext;
                            brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                            warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                            uiwait(gcf)
                        end
                    elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                        logtext='T1 experiment (vdlist file exists)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        TD2= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                    else
                        logtext='only diffusion list found, I will sort like this';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                    end
                elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                    logtext='T2 experiment (vclist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    vclist3=ImportList([FilePath 'vclist'],TD3);
                    vclist3(isnan(vclist3))=0;
                    vclist2=ImportList([FilePath 'vclist'],TD2);
                    vclist2(isnan(vclist2))=0;
                    if nnz(vclist3)>nnz(vclist2)
                        brukerdata.vclist=vclist3;
                    elseif nnz(vclist3)<nnz(vclist2)
                        brukerdata.vclist=vclist2;
                    elseif nnz(vclist3)==nnz(vclist2)
                        if (length(vclist3)-nnz(vclist3))==0
                            brukerdata.vclist=vclist3;
                        else
                            brukerdata.vclist=vclist2;
                        end
                    end
                    brukerdata.ngrad=length(brukerdata.vclist);
                    if exist([FilePath 'vdlist'],'file') == 2
                        TD2= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD2);
                        logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    else
                        
                        logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                        uiwait(gcf)
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                        brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                    end
                elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                    logtext='T1 experiment (vdlist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    vdlist3=ImportList([FilePath 'vdlist'],TD3);
                    vdlist3(isnan(vdlist3))=0;
                    vdlist2=ImportList([FilePath 'vdlist'],TD2);
                    vdlist2(isnan(vdlist2))=0;
                    if nnz(vdlist3)>nnz(vdlist2)
                        brukerdata.vdlist=vdlist3;
                    elseif nnz(vdlist3)<nnz(vdlist2)
                        brukerdata.vdlist=vdlist2;
                    elseif nnz(vdlist3)==nnz(vdlist2)
                        if (length(vdlist3)-nnz(vdlist3))==0
                            brukerdata.vdlist=vdlist3;
                        else
                            brukerdata.vdlist=vdlist2;
                        end
                    end
                    brukerdata.ngrad=length(brukerdata.vdlist);
                    logtext='using vdlist directly for delay values';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                else
                    logtext='Any list was found';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                end
            elseif exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                logtext='3D Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                % Check for gradient shape
                % set the default as 'SINE'
                brukerdata.gpnam='SINE';
                if isfield(brukerdata.acqus,'GPNAM6')
                    if contains(brukerdata.acqus.GPNAM6,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(brukerdata.acqus.GPNAM6,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(brukerdata.acqus.GPNAM6,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                else
                    logtext='Reading gradient shape from the acqus';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    [tmp]=GetGPNAM([FilePath  'acqus']);
                    if contains(tmp,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(tmp,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(tmp,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                end
                if isempty(brukerdata.gpnam) || strcmp(brukerdata.gpnam,'<>')
                    logtext='DOSY gradient shape missing. I will set it to default';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.gpnam='SINE';
                end
                logtext=(['DOSY gradient shape is: ' brukerdata.gpnam]);
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient, which is done automatically with the TopSpin "dosy" setup macro';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                if  contains(brukerdata.pulprog,'psycheidosy','IgnoreCase',true) || contains(brukerdata.pulprog,'psyche_idosy','IgnoreCase',true)
                    logtext='psyche-idosy sequence; setting paramaters to monopolar';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD= str2num(brukerdata.acqu2s.TD);
                    brukerdata.Gzlvl=ImportList([FilePath 'difflist'],TD);
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/4);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(4860*pi.^2)));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    end
                    brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
                    
                elseif strfind(brukerdata.pulprog,'oneshot')
                    logtext='oneshot 3D sequence';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.alpha=brukerdata.acqus.CNST(15);
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(5-3*brukerdata.alpha.^2)*(brukerdata.delta/16)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-brukerdata.delta*((3367*pi.^2-207)*brukerdata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                            (brukerdata.tau/2)*(1-brukerdata.alpha.^2));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(2-brukerdata.alpha.^2)*(brukerdata.delta/6)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    end
                else
                    logtext='unknown sequence; setting parameters to bipolar';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                    end
                    
                end
                if exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                    logtext='T2 experiment (vclist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vclist=ImportList([FilePath 'vclist'], TD2);
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(3) + brukerdata.acqus.P(3)*1e-6);
                    if exist([FilePath 'vdlist'],'file') == 2
                        brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                        logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                    else
                        logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                        warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                        uiwait(gcf)
                    end
                elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                    logtext='T1 experiment (vdlist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                else
                    logtext='only diffusion list found, I will sort like this';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                end
            elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                logtext='T2 experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vclist3=ImportList([FilePath 'vclist'],TD3);
                vclist3(isnan(vclist3))=0;
                vclist2=ImportList([FilePath 'vclist'],TD2);
                vclist2(isnan(vclist2))=0;
                if nnz(vclist3)>nnz(vclist2)
                    brukerdata.vclist=vclist3;
                elseif nnz(vclist3)<nnz(vclist2)
                    brukerdata.vclist=vclist2;
                elseif nnz(vclist3)==nnz(vclist2)
                    if (length(vclist3)-nnz(vclist3))==0
                        brukerdata.vclist=vclist3;
                    else
                        brukerdata.vclist=vclist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vclist);
                if exist([FilePath 'vdlist'],'file') == 2
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD2);
                    logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                else
                    
                    logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                    uiwait(gcf)
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                end
            elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                logtext='T1 experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vdlist3=ImportList([FilePath 'vdlist'],TD3);
                vdlist3(isnan(vdlist3))=0;
                vdlist2=ImportList([FilePath 'vdlist'],TD2);
                vdlist2(isnan(vdlist2))=0;
                if nnz(vdlist3)>nnz(vdlist2)
                    brukerdata.vdlist=vdlist3;
                elseif nnz(vdlist3)<nnz(vdlist2)
                    brukerdata.vdlist=vdlist2;
                elseif nnz(vdlist3)==nnz(vdlist2)
                    if (length(vdlist3)-nnz(vdlist3))==0
                        brukerdata.vdlist=vdlist3;
                    else
                        brukerdata.vdlist=vdlist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vdlist);
                logtext='using vdlist directly for delay values';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            else
                logtext='Any list was found';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        elseif DataDim ~= 1 && DataType==3 % large array
            logtext='Array data for pure acquisition mode'; disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            if DataDim == 2
                brukerdata.arraydim=str2num(brukerdata.acqu2s.TD);
            elseif DataDim ==3
                brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
            else
                %not supported
                return
            end
            brukerdata.ngrad=brukerdata.arraydim;
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.tau='non existing';
            brukerdata.dosyconstant='non existing';
        else %higher order experiments, GDP 05 Nov 2018 - Probably for PARAFAC 
            % MJS editting as gets the Number right but associated with
            % wrong FIDs
            logtext='higher order experiment (>3D). Probably PARAFAC. I will try to sort experimental parameters'; disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD).*str2num(brukerdata.acqu4s.TD);
            brukerdata.ngrad=str2num(brukerdata.acqu2s.TD);
            brukerdata.ngrad2=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
            brukerdata.Gzlvl='non existing';
            brukerdata.DELTA='non existing';
            brukerdata.delta='non existing';
            brukerdata.tau='non existing';
            brukerdata.dosyconstant='non existing';
            if exist([FilePath 'nuslist'],'file') == 2 %NUS file
                logtext='NUS acquired data (nuslist file exists)';
                nusmax=1024; %All the other dimensions must be smaller than this
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                brukerdata.nuslist=ImportDoubleList([FilePath 'nuslist'],nusmax,DataDim);
                brukerdata.nuslist=brukerdata.nuslist(1:(find(isnan(brukerdata.nuslist(:,1))))-1,:,:);
                % GDP Only works for 100% NUS
                TD2=max(brukerdata.nuslist(:,1))+1;
                TD3=max(brukerdata.nuslist(:,2))+1;
                TD4=max(brukerdata.nuslist(:,3))+1; %#ok<NASGU>
                if exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                    logtext='difflist file exists';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    % Check for gradient shape
                    % set the default as 'SINE'
                    brukerdata.gpnam='SINE';
                    if isfield(brukerdata.acqus,'GPNAM6')
                        if contains(brukerdata.acqus.GPNAM6,"SINE",'IgnoreCase',true)
                            brukerdata.gpnam='SINE';
                        elseif contains(brukerdata.acqus.GPNAM6,"SMSQ",'IgnoreCase',true)
                            brukerdata.gpnam='SMSQ';
                        elseif contains(brukerdata.acqus.GPNAM6,"RECT",'IgnoreCase',true)
                            brukerdata.gpnam='RECT';
                        end
                    else
                        logtext='Reading gradient shape from the acqus';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        [tmp]=GetGPNAM([FilePath  'acqus']);
                        if contains(tmp,"SINE",'IgnoreCase',true)
                            brukerdata.gpnam='SINE';
                        elseif contains(tmp,"SMSQ",'IgnoreCase',true)
                            brukerdata.gpnam='SMSQ';
                        elseif contains(tmp,"RECT",'IgnoreCase',true)
                            brukerdata.gpnam='RECT';
                        end
                    end
                    if isempty(brukerdata.gpnam) || strcmp(brukerdata.gpnam,'<>')
                        logtext='DOSY gradient shape missing. I will set it to default';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.gpnam='SINE';
                    end
                    logtext=(['DOSY gradient shape is: ' brukerdata.gpnam]);
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient, which is done automatically with the TopSpin "dosy" setup macro';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    logtext=('Attenuation function calculated in the pulse sequence');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.dosyconstant=brukerdata.acqus.CNST(19).*brukerdata.gamma.^2;
                    brukerdata.tau=brukerdata.acqus.CNST(18);
                    if exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                        logtext='vclist file exists';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        TDn= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vclist=ImportList([FilePath 'vclist'], TDn);
                        brukerdata.vclist=brukerdata.vclist(1:(find(isnan(brukerdata.vclist)))-1);
                        TDn=length(brukerdata.vclist);
                        brukerdata.ngrad2=TDn*brukerdata.ngrad;
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(3) + brukerdata.acqus.P(3)*1e-6);
                        if exist([FilePath 'vdlist'],'file') == 2
                            brukerdata.vdlist=ImportList([FilePath 'vdlist'], TDn);
                            logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                            disp(logtext)
                            brukerdata.logfile{end +1}=logtext;
                        else
                            logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                            disp(logtext)
                            brukerdata.logfile{end +1}=logtext;
                            brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                            warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                            uiwait(gcf)
                        end
                    elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                        logtext='vdlist file exists';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        TDn= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vdlist=ImportList([FilePath 'vdlist'], TDn);
                        brukerdata.vdlist=brukerdata.vdlist(1:(find(isnan(brukerdata.vdlist)))-1);
                        TDn=length(brukerdata.vdlist);
                        brukerdata.ngrad2=TDn*brukerdata.ngrad;
                    else
                        logtext='only diffusion list found, I will sort like this';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                    end
                elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                    logtext='T2 experiment (vclist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    vclist3=ImportList([FilePath 'vclist'],TD3);
                    vclist3(isnan(vclist3))=0;
                    vclist2=ImportList([FilePath 'vclist'],TD2);
                    vclist2(isnan(vclist2))=0;
                    if nnz(vclist3)>nnz(vclist2)
                        brukerdata.vclist=vclist3;
                    elseif nnz(vclist3)<nnz(vclist2)
                        brukerdata.vclist=vclist2;
                    elseif nnz(vclist3)==nnz(vclist2)
                        if (length(vclist3)-nnz(vclist3))==0
                            brukerdata.vclist=vclist3;
                        else
                            brukerdata.vclist=vclist2;
                        end
                    end
                    brukerdata.ngrad=length(brukerdata.vclist);
                    if exist([FilePath 'vdlist'],'file') == 2
                        TD2= brukerdata.arraydim/brukerdata.ngrad;
                        brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD2);
                        logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    else
                        
                        logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                        uiwait(gcf)
                        brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                        brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                    end
                elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                    logtext='T1 experiment (vdlist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    vdlist3=ImportList([FilePath 'vdlist'],TD3);
                    vdlist3(isnan(vdlist3))=0;
                    vdlist2=ImportList([FilePath 'vdlist'],TD2);
                    vdlist2(isnan(vdlist2))=0;
                    if nnz(vdlist3)>nnz(vdlist2)
                        brukerdata.vdlist=vdlist3;
                    elseif nnz(vdlist3)<nnz(vdlist2)
                        brukerdata.vdlist=vdlist2;
                    elseif nnz(vdlist3)==nnz(vdlist2)
                        if (length(vdlist3)-nnz(vdlist3))==0
                            brukerdata.vdlist=vdlist3;
                        else
                            brukerdata.vdlist=vdlist2;
                        end
                    end
                    brukerdata.ngrad=length(brukerdata.vdlist);
                    logtext='using vdlist directly for delay values';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                else
                    logtext='Any list was found';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                end
            elseif exist([FilePath 'difflist'],'file') == 2 %Diffusion experiment
                logtext='3D Diffusion experiment (difflist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                logtext=['DOSY pulse programme is: ' brukerdata.pulprog];
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                
                % Check for gradient shape
                % set the default as 'SINE'
                brukerdata.gpnam='SINE';
                if isfield(brukerdata.acqus,'GPNAM6')
                    if contains(brukerdata.acqus.GPNAM6,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(brukerdata.acqus.GPNAM6,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(brukerdata.acqus.GPNAM6,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                else
                    logtext='Reading gradient shape from the acqus';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    [tmp]=GetGPNAM([FilePath  'acqus']);
                    if contains(tmp,"SINE",'IgnoreCase',true)
                        brukerdata.gpnam='SINE';
                    elseif contains(tmp,"SMSQ",'IgnoreCase',true)
                        brukerdata.gpnam='SMSQ';
                    elseif contains(tmp,"RECT",'IgnoreCase',true)
                        brukerdata.gpnam='RECT';
                    end
                end
                if isempty(brukerdata.gpnam) || strcmp(brukerdata.gpnam,'<>')
                    logtext='DOSY gradient shape missing. I will set it to default';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.gpnam='SINE';
                end
                logtext=(['DOSY gradient shape is: ' brukerdata.gpnam]);
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                logtext='Warning! - GNAT assumes that the difflist is corrected for the shape of the gradient, which is done automatically with the TopSpin "dosy" setup macro';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                
                if  strfind(brukerdata.pulprog,'oneshot')
                    logtext='oneshot 3D sequence';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.alpha=brukerdata.acqus.CNST(15);
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(5-3*brukerdata.alpha.^2)*(brukerdata.delta/16)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-brukerdata.delta*((3367*pi.^2-207)*brukerdata.alpha.^2+(6344*pi.^2-207))/(19440*pi.^2) -...
                            (brukerdata.tau/2)*(1-brukerdata.alpha.^2));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2*brukerdata.delta.^2*...
                            (brukerdata.DELTA-(2-brukerdata.alpha.^2)*(brukerdata.delta/6)-(1-brukerdata.alpha.^2)*(brukerdata.tau/2));
                    end
                elseif  strfind(brukerdata.pulprog,'psyche_idosy') %GDP 13v2021 This needs testing
                    logtext='psyche_idosy sequence; setting paramaters to monopolar';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=brukerdata.acqus.P(31)*1e-6;
                    
                    %Set the ST values according to the gradient shape
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/4);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((1484*pi.^2-207)/(4860*pi.^2)));
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3);
                    end
                    brukerdata.arraydim=str2num(brukerdata.acqu2s.TD).*str2num(brukerdata.acqu3s.TD);
                else
                    logtext='unknown sequence; setting parameters to bipolar';
                    disp(logtext)
                    TD3= str2num(brukerdata.acqu3s.TD);
                    TD2= str2num(brukerdata.acqu2s.TD);
                    Gzlvl3=ImportList([FilePath 'difflist'],TD3);
                    Gzlvl3(isnan(Gzlvl3))=0;
                    Gzlvl2=ImportList([FilePath 'difflist'],TD2);
                    Gzlvl2(isnan(Gzlvl2))=0;
                    if nnz(Gzlvl3)>nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl3;
                    elseif nnz(Gzlvl3)<nnz(Gzlvl2)
                        brukerdata.Gzlvl=Gzlvl2;
                    elseif nnz(Gzlvl3)==nnz(Gzlvl2)
                        if (length(Gzlvl3)-nnz(Gzlvl3))==0
                            brukerdata.Gzlvl=Gzlvl3;
                        else
                            brukerdata.Gzlvl=Gzlvl2;
                        end
                    end
                    brukerdata.Gzlvl=brukerdata.Gzlvl*0.01; %conversion to T/m
                    brukerdata.ngrad=length(brukerdata.Gzlvl);
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.DELTA=brukerdata.acqus.D(21);
                    brukerdata.delta=2*brukerdata.acqus.P(31)*1e-6; %probably assuming bipolar pulse
                    brukerdata.tau=brukerdata.acqus.D(17)+brukerdata.acqus.P(3)*1e-6;
                    logtext=('Gradient shape will be used for calculation of the attenuation function');
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    if strcmpi(brukerdata.gpnam,"SINE")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*(5/16) - brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"SMSQ")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta*((6344*pi.^2-207)/(19440*pi.^2))- brukerdata.tau/2);
                    elseif strcmpi(brukerdata.gpnam,"RECT")
                        brukerdata.dosyconstant=brukerdata.gamma.^2.*brukerdata.delta.^2.*(brukerdata.DELTA-brukerdata.delta/3 - brukerdata.tau/2);
                    end
                end
                if exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                    logtext='T2 experiment (vclist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vclist=ImportList([FilePath 'vclist'], TD2);
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(3) + brukerdata.acqus.P(3)*1e-6);
                    if exist([FilePath 'vdlist'],'file') == 2
                        brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                        logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                    else
                        logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                        disp(logtext)
                        brukerdata.logfile{end +1}=logtext;
                        brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                        warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                        uiwait(gcf)
                    end
                elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                    logtext='T1 experiment (vdlist file exists)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist=ImportList([FilePath 'vdlist'], TD2);
                else
                    logtext='only diffusion list found, I will sort like this';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                end
            elseif exist([FilePath 'vclist'],'file') == 2 %T2 experiment
                logtext='T2 experiment (vclist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vclist3=ImportList([FilePath 'vclist'],TD3);
                vclist3(isnan(vclist3))=0;
                vclist2=ImportList([FilePath 'vclist'],TD2);
                vclist2(isnan(vclist2))=0;
                if nnz(vclist3)>nnz(vclist2)
                    brukerdata.vclist=vclist3;
                elseif nnz(vclist3)<nnz(vclist2)
                    brukerdata.vclist=vclist2;
                elseif nnz(vclist3)==nnz(vclist2)
                    if (length(vclist3)-nnz(vclist3))==0
                        brukerdata.vclist=vclist3;
                    else
                        brukerdata.vclist=vclist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vclist);
                if exist([FilePath 'vdlist'],'file') == 2
                    TD2= brukerdata.arraydim/brukerdata.ngrad;
                    brukerdata.vdlist = ImportList([FilePath 'vdlist'], TD2);
                    logtext='using vdlist directly for delay values, overriding vclist, but saving it for manual conversion';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                else
                    
                    logtext='Assuming vclist is for a T2 experiment. Converting to vdlist by vclist*(2*d20 + p2)';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    warndlg('CPMG total echo time is used by default. To change it go to "Edit > Settings"')
                    uiwait(gcf)
                    brukerdata.vc_constant=(2*brukerdata.acqus.D(21) + brukerdata.acqus.P(3)*1e-6);
                    brukerdata.vdlist=brukerdata.vclist*brukerdata.vc_constant;
                end
            elseif exist([FilePath 'vdlist'],'file') == 2 %T1 experiment
                logtext='T1 experiment (vdlist file exists)';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                TD3= str2num(brukerdata.acqu3s.TD);
                TD2= str2num(brukerdata.acqu2s.TD);
                vdlist3=ImportList([FilePath 'vdlist'],TD3);
                vdlist3(isnan(vdlist3))=0;
                vdlist2=ImportList([FilePath 'vdlist'],TD2);
                vdlist2(isnan(vdlist2))=0;
                if nnz(vdlist3)>nnz(vdlist2)
                    brukerdata.vdlist=vdlist3;
                elseif nnz(vdlist3)<nnz(vdlist2)
                    brukerdata.vdlist=vdlist2;
                elseif nnz(vdlist3)==nnz(vdlist2)
                    if (length(vdlist3)-nnz(vdlist3))==0
                        brukerdata.vdlist=vdlist3;
                    else
                        brukerdata.vdlist=vdlist2;
                    end
                end
                brukerdata.ngrad=length(brukerdata.vdlist);
                logtext='using vdlist directly for delay values';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            else
                logtext='No list was found';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        end
    end
    
    %% Read in the binary data
    %check the format of the data
    %The byte order will be big from SGI (MIPS processor) and little from
    %Linux/Windows (Intel), this is stored as the status parameters BYTORDA
    %(acquisition data) and BYTORDP (processed data).  Processed data is always
    %stored as scaled 32 bit integers, but acquisition data can in principle be
    %present as double precision (64 bit) data - this is indicated by the
    %status parameter DTYPA (0=integer, 2=double).  The double state is
    %automatically triggered if the raw data would overflow (but even in the
    %modern age of large effective ADC resolution this is probably a rare event!)

    switch str2double(brukerdata.acqus.BYTORDA)
        case 0 %little endian
            byte_format='l';
        case 1  %big endian
            byte_format='b';
        otherwise
            error('unknown data format (BYTORDA)')
    end
     switch str2double(brukerdata.acqus.DTYPA) %MN 29 March correcting chunkSize to accomodate for Neo
        case 0 %32 bit integer
            byte_size='int32';
            chunkSize = 256;
        case 1  %double
            byte_size='double';
            chunkSize = 128;
        case 2  %double
            byte_size='double'; %MN 13March2017 Not sure if 1 or 2 does any difference to Matlab
            chunkSize = 128;
        otherwise
            error('unknown data format (DTYPA)')
    end
    
    %Check that np is a correct multiple  - otherwise correct
    % In a ser file each 1D fid will start at a 1024 byte block boundary even
    % if its size is not a multiple of 1k (256 points)
    corrpoints=rem(brukerdata.np*2,chunkSize);
    if corrpoints>0
        corrpoints=chunkSize-corrpoints;
        brukerdata.np=brukerdata.np+corrpoints/2;
    end
    
    FID=zeros(brukerdata.np,brukerdata.arraydim);
    
    
    %if importing the processed spectrum
    if DataType==1 %Processed data
        byte_format_proc='l';
        %Seems like it should always be l
        byte_size='int32';
        if DataDim==1
            ProcFIDpath=[ProcPathRoot '1r'];
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                realfid=fread(fileid_ProcFIDreal,brukerdata.SI*brukerdata.arraydim,byte_size);

                ProcFIDpath=[ProcPathRoot '1i'];
                if exist(ProcFIDpath, 'file') == 2
                    fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                    imagfid=fread(fileid_ProcFIDreal,brukerdata.SI*brukerdata.arraydim,byte_size);
                    % realfid=realfid*2^str2num(brukerdata.procs.NC_proc);
                    totfid=complex(realfid,imagfid);
                    totfid=totfid*2^str2num(brukerdata.procs.NC_proc);

                    FID=ifft(ifftshift(totfid));
                
                    brukerdata.FID=FID(1:round(brukerdata.SI));
              
                else
                    error(['Processed data 1D data 1i does not exist in:' ProcFIDpath])

                end

                % FID=ifft(ifftshift(realfid));
                % FID=FID(1:round(brukerdata.SI/2));
                fclose(fileid_ProcFIDreal);
            else
                error(['Processed data 1D data 1r does not exist in:' ProcFIDpath])
            end
            %size(FID)
        elseif DataDim==2
            ProcFIDpath=[ProcPathRoot '2rr'];
            
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                %  ProcFIDpath
                SI1=str2num(brukerdata.procs.SI);
                SI2=str2num(brukerdata.proc2s.SI);
                brukerdata.np=SI1;
                XDIM1=str2num(brukerdata.procs.XDIM);
                XDIM2=str2num(brukerdata.proc2s.XDIM);
                
                realfid=fread(fileid_ProcFIDreal,brukerdata.np*2*SI2,byte_size);
                realfid=realfid*2^str2num(brukerdata.proc2s.NC_proc);
                
                NoSM = SI1*SI2/(XDIM1*XDIM2);
                NoSM2 = SI2/XDIM2;
                
                %recipe from Nils Nyberg - works most of the time. Should
                %make it work 100% of the time when I feel like it.
                realfid = reshape(...
                    permute(...
                    reshape(...
                    permute(...
                    reshape(realfid,XDIM1,XDIM2,NoSM),...
                    [2 1 3]),...
                    XDIM2,SI1,NoSM2),...
                    [2 1 3]),...
                    SI1,SI2)';
                
                realfid=realfid';
                
                %This saves the processed 2D data for use with e.g. DOSY3D
                brukerdata.X2DSPEC=realfid;
                

                ProcFIDpath=[ProcPathRoot '2ii'];
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);

                imagfid=fread(fileid_ProcFIDreal,brukerdata.np*2*SI2,byte_size);
                imagfid=imagfid*2^str2num(brukerdata.proc2s.NC_proc);

                NoSM = SI1*SI2/(XDIM1*XDIM2);
                NoSM2 = SI2/XDIM2;

                %recipe from Nils Nyberg - works most of the time. Should
                %make it work 100% of the time when I feel like it.
                imagfid = reshape(...
                    permute(...
                    reshape(...
                    permute(...
                    reshape(imagfid,XDIM1,XDIM2,NoSM),...
                    [2 1 3]),...
                    XDIM2,SI1,NoSM2),...
                    [2 1 3]),...
                    SI1,SI2)';
                
                imagfid=imagfid';


                totfid=complex(realfid,imagfid);

                FID=totfid(:,1:brukerdata.arraydim);

                
                %FID=realfid(:,1:brukerdata.arraydim);
                
                
                for p=1:brukerdata.arraydim
                    FID(:,p)=ifft(ifftshift(FID(:,p)));
                end
                %brukerdata.FID=FID(1:brukerdata.np/2,:);
                brukerdata.FID=FID;
            else
                error(['Processed data 2D data 2rr does not exist in:' ProcFIDpath])
            end
        else
            error('Import of processed data for more that two dimensions is not implemented')
        end
        
        %brukerdata.FID=FID(1:brukerdata.np/2,:);
        brukerdata.lp=0;
        brukerdata.rp=0;
    elseif DataType==2 %Processed data array
        byte_format_proc='l';
        %Seems like it should always be l
        byte_size='int32';
        if DataDim==1

            ProcFIDpath=[ProcPathRoot '1r'];
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                realfid=fread(fileid_ProcFIDreal,brukerdata.SI*brukerdata.arraydim,byte_size);

                ProcFIDpath=[ProcPathRoot '1i'];
                if exist(ProcFIDpath, 'file') == 2
                    fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                    imagfid=fread(fileid_ProcFIDreal,brukerdata.SI*brukerdata.arraydim,byte_size);
                    % realfid=realfid*2^str2num(brukerdata.procs.NC_proc);
                    totfid=complex(realfid,imagfid);
                    totfid=totfid*2^str2num(brukerdata.procs.NC_proc);

                    FID=ifft(ifftshift(totfid));
size(FID)
                    brukerdata.FID=FID(1:round(brukerdata.SI));

                else
                    error(['Processed data 1D data 1i does not exist in:' ProcFIDpath])

                end


                fclose(fileid_ProcFIDreal);
            else
                error(['Processed data 1D data 1r does not exist in:' ProcFIDpath])
            end

        elseif DataDim==2
            ProcFIDpath=[ProcPathRoot '2rr'];
            
            if exist(ProcFIDpath, 'file') == 2
                fileid_ProcFIDreal=fopen(ProcFIDpath ,'r',byte_format_proc);
                %  ProcFIDpath
                SI1=str2num(brukerdata.procs.SI);
                SI2=str2num(brukerdata.proc2s.SI);
                brukerdata.np=SI1;
                XDIM1=str2num(brukerdata.procs.XDIM);
                XDIM2=str2num(brukerdata.proc2s.XDIM);
                
                realfid=fread(fileid_ProcFIDreal,brukerdata.np*2*SI2,byte_size);
                realfid=realfid*2^str2num(brukerdata.proc2s.NC_proc);
                
                NoSM = SI1*SI2/(XDIM1*XDIM2);
                NoSM2 = SI2/XDIM2;
                
                %recipe from Nils Nyberg - works most of the time. Should
                %make it work 100% of the time when I feel like it.
                realfid = reshape(...
                    permute(...
                    reshape(...
                    permute(...
                    reshape(realfid,XDIM1,XDIM2,NoSM),...
                    [2 1 3]),...
                    XDIM2,SI1,NoSM2),...
                    [2 1 3]),...
                    SI1,SI2)';
                
                realfid=realfid';
                
                %This saves the processed 2D data for use with e.g. DOSY3D
                brukerdata.X2DSPEC=realfid;
                
                
                FID=realfid(:,1:brukerdata.arraydim);
                
                
                for p=1:brukerdata.arraydim
                    FID(:,p)=ifft(ifftshift(FID(:,p)));
                end
            else
                error(['Processed data 2D data 2rr does not exist in:' ProcFIDpath])
            end
        else
            error('Import of processed data for more that two dimensions is not implemented')
        end
        
        %brukerdata.FID=FID(1:brukerdata.np/2,:);
        brukerdata.lp=0;
        brukerdata.rp=0;

    else %Raw FID data
        
        fileid_FID=fopen(fidpath,'r',byte_format);
        
        impfid=fread(fileid_FID,brukerdata.np*2*brukerdata.arraydim,byte_size);
        compfid=complex(impfid(1:2:end),-impfid(2:2:end));
        compfid=compfid*2^str2num(brukerdata.acqus.NC); %MN 17 April 2023
        [m,~]=size(compfid);
        
        if m<brukerdata.np
            logtext='Actual points in FID and number of points resported (TD) are different. This can happen after the use of the "convdta" macro. I will use the actual number and continue import';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            brukerdata.np=m;
            FID=zeros(brukerdata.np,brukerdata.arraydim);
        end
        
        for k=1:brukerdata.arraydim
            FID(:,k)=compfid((k-1)*brukerdata.np+1:k*brukerdata.np);
        end
        
        brukerdata.np=brukerdata.np-corrpoints/2;
        brukerdata.FID=FID((1:brukerdata.np),:);
        fclose(fileid_FID);
    end


    
    %% Sort out the digital filtering
    % This is needed for converion of digitally filtered data - se longer
    % explanation at the end of the file
    BrukDigital=[
        2       44.750       46.000       46.311;
        3       33.500       36.500       36.530;
        4       66.625       48.000       47.870;
        6       59.083       50.167       50.229;
        8       68.563       53.250       53.289;
        12       60.375       69.500       69.551;
        16       69.531       72.250       71.600;
        24       61.021       70.167       70.184;
        32       70.016       72.750       72.138;
        48       61.344       70.500       70.528;
        64       70.258       73.000       72.348;
        96       61.505       70.667       70.700;
        128       70.379       72.500       72.524;
        192       61.586       71.333            0;
        256       70.439       72.250            0;
        384       61.626       71.667            0;
        512       70.470       72.125            0;
        768       61.647       71.833            0;
        1024       70.485       72.063            0;
        1536       61.657       71.917            0;
        2048       70.492       72.031            0];
    
    % first check if GRPDLY exists and use that if so
    decim=1;
    dspfvs=1;
    if isfield(brukerdata.acqus,'GRPDLY') && str2double(brukerdata.acqus.GRPDLY)~=-1
        brukerdata.digshift=str2double(brukerdata.acqus.GRPDLY);
        %disp('Hi')
    else
        %brukerdata.acqus.DECIM=0;
        
        if isfield(brukerdata.acqus,'DECIM')
            decim=str2double(brukerdata.acqus.DECIM);
            switch decim
                case 2
                    decimrow=1;
                case 3
                    decimrow=2;
                case 4
                    decimrow=3;
                case 6
                    decimrow=4;
                case 8
                    decimrow=5;
                case 12
                    decimrow=6;
                case 16
                    decimrow=7;
                case 24
                    decimrow=8;
                case 32
                    decimrow=9;
                case 48
                    decimrow=10;
                case 64
                    decimrow=11;
                case 96
                    decimrow=12;
                case 128
                    decimrow=13;
                case 192
                    decimrow=14;
                case 256
                    decimrow=15;
                case 384
                    decimrow=16;
                case 512
                    decimrow=17;
                case 768
                    decimrow=18;
                case 1024
                    decimrow=19;
                case 1536
                    decimrow=20;
                case 2048
                    decimrow=21;
                otherwise
                    logtext='unknown value of DECIM parameter in acqus - cannot set compensation for digital filtering';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    decim=0;
                    decimrow=Inf;
            end
        else
            logtext='no DECIM parameter in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            decim=0;
            decimrow=Inf;
        end
        if isfield(brukerdata.acqus,'DSPFVS')
            dspfvs=str2double(brukerdata.acqus.DSPFVS);
            switch dspfvs
                case 10
                    dspfvsrow=2;
                case 11
                    dspfvsrow=3;
                case 12
                    dspfvsrow=4;
                otherwise
                    logtext='unknown value of DSPVFS parameter in acqus - cannot set compensation for digital filtering';
                    disp(logtext)
                    brukerdata.logfile{end +1}=logtext;
                    dspfvs=0;
                    dspfvsrow=0;
            end
        else
            logtext='no DSPFVS parameter in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            dspfvs=0;
        end
        if (decimrow>14) && (dspfvsrow==4)
            logtext='unknown combination of DSPVFS and DECIM parameters in acqus - cannot set compensation for digital filtering';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
            dspfvs=0;
            decim=0;
        end
        
    end
    if (decim==0) || (dspfvs==0)
        %No digital filtering
        logtext='Parameters for digital filtering unknown - assumed to be data without digital filtering';
        disp(logtext)
        brukerdata.logfile{end +1}=logtext;
        brukerdata.digshift=0;
    elseif (decim==1) && (dspfvs==1)
        %digital filtering set by GRPDLY
        %do nothing
    else
        brukerdata.digshift=BrukDigital(decimrow,dspfvsrow);
    end
    %brukerdata.digshift=0
    
    if DataType==1
        %brukerdata.digshift=0;
       %shift the DSP point to the front
       brukerdata.digshift=round(brukerdata.digshift);
       brukerdata.FID= circshift(brukerdata.FID,+brukerdata.digshift);
       brukerdata.np=brukerdata.SI;
       brukerdata.FID=brukerdata.FID(1:brukerdata.np,:);
       brukerdata.at=brukerdata.np/(brukerdata.sw*brukerdata.sfrq);
    elseif DataType==2
        brukerdata.digshift=0;
       brukerdata.FID= circshift(brukerdata.FID,+brukerdata.digshift);
       brukerdata.np=brukerdata.SI;
       brukerdata.FID=brukerdata.FID(1:brukerdata.np,:);
       brukerdata.at=brukerdata.np/(brukerdata.sw*brukerdata.sfrq);
    end    
    
    % if any(size(dir([FilePath '/*.m' ]),1))==1
    % 
    %     brukerdata.digshift=0;
    % end
    
    brukerdata.digshift=round(brukerdata.digshift);
    brukerdata.np=length(brukerdata.FID(:,1));
    %% Check the acquisition order for 3D Data
    if DataDim==3 && DataType~=3
        brukerdata.procs;
        brukerdata.acqus;
        brukerdata.acqus.CNST(1);
        brukerdata.acqus.CNST(2);
        brukerdata.acqus.CNST(3);
        brukerdata.acqus.CNST(4);
        brukerdata.acqus.CNST(5);
        if isfield(brukerdata.acqus,'AQSEQ')
            if str2double(brukerdata.acqus.AQSEQ)==0
                %all is fine - this is the structure we expect
                logtext='Normal aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            elseif str2double(brukerdata.acqus.AQSEQ)==1
                % reversed order - lets fix it
                logtext='Reversed aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
                orgfid=brukerdata.FID;
                
                % brukerdata.FID=brukerdata.FID.*0;
                
                array2=brukerdata.arraydim/brukerdata.ngrad;
                
                for k=1:array2
                    for p=1:brukerdata.ngrad
                        % brukerdata.FID(:,k+array2*(p-1))=orgfid(:,(k-1)*brukerdata.ngrad +p);
                        brukerdata.FID(:,(k-1)*brukerdata.ngrad +p)=orgfid(:,k+array2*(p-1));
                    end
                end
            else
                %we're not sure
                logtext='Unknown aquisition order';
                disp(logtext)
                brukerdata.logfile{end +1}=logtext;
            end
        else
            logtext='no AQSEQ parameter. Cannot determine aquisition mode. Sticking with default';
            disp(logtext)
            brukerdata.logfile{end +1}=logtext;
        end
    else
        %not a 3D set
    end
    
else
    brukerdata=[];
end
end
%% Utility Functions
function OutParams=GetParameters(FilePath)
%acquspath=[FilePath '/acqus'];
%Reading in all the acquisition or processing parameters from the acquXs or
% procXs files
fileid=fopen(FilePath,'rt');
if fileid==-1
    error(['cannot open file: ' acquspath])
end
kk=1;
parmflag=0;
while 1
    if parmflag==1
        parmflag=0;
    else
        parmline=fgetl(fileid);
    end
    if parmline==-1;  break;  end
    
    if length(parmline)<2
        logtext='probably a cut of line - I will continue anyway';
        disp(logtext)
    elseif strcmp(parmline(1:2),'##')
        if strcmp(parmline(3),'$')
            parmend=find(parmline=='=');
            if strcmp(parmline(parmend+2),'(')
                parmname=parmline(4:(parmend(1))-1);
                parmlist=[];
                while 1
                    parmline=fgetl(fileid);
                    if strcmp(parmline(1),'#') || strcmp(parmline(1),'<') || strcmp(parmline(1),'>') || strcmp(parmline(1),'$')
                        parmflag=1;
                        break
                    else
                        parmcell=textscan(parmline,'%q',1024);
                        tadam=parmcell{1};
                        tmplist=zeros(length(tadam),1);
                        for m=1:length(tadam)
                            tmplist(m)=str2num(cell2mat(tadam(m)));
                        end
                        parmlist=[parmlist; tmplist];                   %#ok<AGROW>
                    end
                end
                procpar.(parmname)=parmlist;
            else
                procpar.(parmline(4:(parmend(1))-1))=parmline((parmend+2):end);
            end
        else
            parmend=find(parmline=='=');
            procpar.(parmline(3:(parmend(1)-1)))=parmline((parmend+2):end);
        end
    elseif strcmp(parmline(1:2),'$$')
        str=['acqinfo' num2str(kk)];
        procpar.(str)=parmline(4:end);
        kk=kk+1;
    else
        %do nothing
    end
end
OutParams=procpar;
fclose(fileid);

end
function GPNAM6=GetGPNAM(FilePath)
%Reading gradient shape from the acqus
fileid=fopen(FilePath,'rt');
if fileid==-1
    error(['cannot open file: ' FilePath])
end

parmcell=textscan(fileid,'%s','Delimiter',{'##$GPNAM= (0..31)','##$GPX= (0..31)'});
tadam=parmcell{1};
INDEX=find(cellfun('isempty',tadam'));

if numel(INDEX)~=2
    disp('This is probably an older version of TopSpin')
    GPNAM6=[];
    return
end

tadam=tadam(INDEX(1)+1:INDEX(2)-1);
newStr = cell2mat(join(tadam));
newStr=split(newStr,' ',1);
GPNAM6=cell2mat(newStr(7));

fclose(fileid);

end
function Vals = ImportList(filepath, TD)
% import gradient values from difflist
fileid_gzlvl=fopen(filepath,'rt');
if (fileid_gzlvl~=-1)
    Vals=zeros(1,TD);
    for k=1:TD
        parmline=fgetl(fileid_gzlvl);
        Vals(k)=str2double(parmline);
        if parmline==-1
            logtext='Warning! -number of  values in list does not correspond to the number ponits in TD';
            disp(logtext)
            break;
        end
    end
    
    fclose(fileid_gzlvl);
else
    logtext='no list file; cannot read invalues';
    disp(logtext)
    %brukerdata.logfile{end +1}=logtext;
end
end
function Vals = ImportDoubleList(filepath, TD, DataDim)
fileid_gzlvl=fopen(filepath,'rt');
if (fileid_gzlvl~=-1)
    Vals=zeros(TD,(DataDim-1));
    for k=1:TD
        parmline=fgetl(fileid_gzlvl);
        if parmline==-1
            logtext='Warning! -number of  values in list does not correspond to the number points in TD';
            disp(logtext)
            Vals(k,:)=NaN;
            break;
        else
            tmp=str2num(parmline);
            Vals(k,:)=tmp;
        end
    end
    
    fclose(fileid_gzlvl);
else
    logtext='no list file; cannot read invalues';
    disp(logtext)
    %brukerdata.logfile{end +1}=logtext;
end
end
%% INFORMATION ABOUT BRUKER DIGITAL FILTER
% For older versions of the Bruker hardware:
%
% A nice piece was found on the internet on how to calculate the number of points
% semi-automatically. Note that currently matNMR doesn't allow for the necessary
% negative-time apodization.
%
%
%    W. M. Westler and F.  Abildgaard
%    July 16, 1996
%
%    The introduction of digital signal processing by Bruker in their DMX
%    consoles also introduced an unusual feature associated with the data. The
%    stored FID no longer starts at its maximum followed by a decay, but is
%    prepended with an increasing signal that starts from zero at the
%    first data point and rises to a maximum after several tens of data points.
%    On transferring this data to a non-Bruker processing program such as FELIX,
%    which is used at NMRFAM, the Fourier transform leads to an unusable spectrum
%    filled with wiggles. Processing the same data with Bruker's Uxnmr
%    program yields a correct spectrum. Bruker has been rather reluctant
%    to describe what tricks are implemented during their processing protocol.
%
%    They suggest the data should be first transformed in Uxnmr and then inverse
%    transformed, along with a GENFID command, before shipping the data to another
%    program. Bruker now supplies a piece of software to convert the digitally
%    filtered data to the equivalent analog form.
%    We find it unfortunate that the vendor has decided to complicate
%    the simple task of Fourier transformation. We find that the procedure
%    suggested by Bruker is cumbersome, and more so, we are irritated since
%    we are forced to use data that has been treated with an unknown procedure.
%    Since we do not know any details of Bruker's digital filtration procedure
%    or the "magic" conversion routine that is used in Uxnmr, we have been forced
%    into observation and speculation. We have found a very simple, empirical
%    procedure that leads to spectra processed in FELIX that are identical,
%    within the noise limits, to spectra processed with Uxnmr. We deposit
%    this information here in the hope that it can be of some
%    use to the general community.
%    The application of a nonrecursive (or recursive) digital filter to time
%    domain data is accomplished by performing a weighted running average of
%    nearby data points. A problem is encountered at the beginning of
%    the data where, due to causality, there are no prior values. The
%    weighted average of the first few points, therefore, must include data
%    from "negative" time. One naive procedure, probably appropriate to NMR
%    data, is to supply values for negative time points is to pad the data with
%    zeros. Adding zeros (or any other data values) to the beginning of
%    the FID, however, shifts the beginning of the time domain data (FID) to
%    a later positive time. It is well known that a shift in the time
%    domain data is equivalent to the application of a frequency-dependent,
%    linear phase shift in the frequency domain. The 1st order phase shift
%    corresponding to a time shift of a single complex dwell point is 360 degrees
%    across the spectral width. The typical number of prepended points
%    found in DMX digitally filtered data is about 60 data points (see below),
%
%    the corresponding 1st order phase correction is ~21,000 degrees.
%    This large linear phase correction can be applied to the transformed data
%    to obtain a normal spectrum. Another, equivalent approach is to time
%    shift the data back to its original position. This results in the need
%    of only a small linear phase shift on the transformed data.
%    There is a question as what to do with the data preceding the actual
%    FID. The prepended data can be simply eliminated with the addition
%    of an equal number of zeros at the end of the FID (left shift). This
%    procedure, however, introduces "frowns" (some have a preference to refer
%    to these as "smiles") at the edge of the spectrum. If the sweep
%    width is fairly wide this does not generally cause a problem. The
%    (proper) alternative is to retain this data by applying a circular left
%    shift of the data, moving the first 60 or so points (see recommendations
%    below) to the end of the FID. This is identical to a Fourier transformation
%    followed by the large linear phase correction mentioned above. The
%    resulting FID is periodic with the last of the data rising to meet the
%    first data point (in the next period). Fourier transform of this
%    data results in an approximately phased spectrum. Further linear
%    phase corrections of up to 180 degrees are necessary. A zero fill applied
%    after a circular shift of the data will cause a discontinuity and thus
%    introduce sinc wiggles on the peaks. The usual correction for DC
%    offset and apodization of the data, if not done correctly, also results
%    in the frowns at the edges of the spectrum.
%
%    In our previous document on Bruker digital filters, we presented deduced
%    rules for calculating the appropriate number of points to be circular left
%    shifted. However, since then, newer versions of hardware (DQD) and software
%    has introduced a new set of values. Depending on the firmware versions
%    (DSPFVS) and the decimation rate (DECIM), the following lookup table will
%    give the circular shift values needed to correct the DMX data. The values
%    of DECIM and DSPFVS can be found in the acqus file in the directory containing
%    the data.
%
%     DECIM           DSPFVS 10       DSPFVS 11      DSPFVS 12
%
%       2              44.7500         46.0000        46.311
%       3              33.5000         36.5000        36.530
%       4              66.6250         48.0000        47.870
%       6              59.0833         50.1667        50.229
%       8              68.5625         53.2500        53.289
%      12              60.3750         69.5000        69.551
%      16              69.5313         72.2500        71.600
%      24              61.0208         70.1667        70.184
%      32              70.0156         72.7500        72.138
%      48              61.3438         70.5000        70.528
%      64              70.2578         73.0000        72.348
%      96              61.5052         70.6667        70.700
%     128              70.3789         72.5000        72.524
%     192              61.5859         71.3333
%     256              70.4395         72.2500
%     384              61.6263         71.6667
%     512              70.4697         72.1250
%     768              61.6465         71.8333
%    1024              70.4849         72.0625
%    1536              61.6566         71.9167
%    2048              70.4924         72.0313
%
%
%    The number of points obtained from the table are usually not integers.  The appropriate procedure is to circular shift (see protocol for details) by the integer obtained from truncation of the obtained value and then the residual 1st order phase shift that needs to be applied can be obtained by multiplying the decimal portion of the calculated number of points by 360.
%
%    For example,
%
%    If DECIM = 32, and DSPFVS = 10,
%    then #points 70.0156
%
%    The circular shift performed on the data should be 70 complex points and the linear
%    phase correction after Fourier transformation is approximately 0.0156*360= 5.62 degrees.
%
%    Protocol:
%
%       1. Circular shift (rotate) the appropriate number of points in the data indicated by
%       the  DECIM parameter. (see above formulae).
%
%       2. After the circular shift, resize the data to the original size minus
%       the number of shifted points. This will leave only the part of the
%       data that looks like an FID. Baseline correction (BC) and/or apodization
%       (EM etc.) should be applied only on this data, otherwise "In come the frowns."
%
%       Since the first part of the data (the points that are shifted) represents
%       negative time, a correct apodization would also multiply the shifted points
%       by a negative time apodization. The data size is now returned to
%       its original size to reincorporate the shifted points. There may
%       still be a discontinuity between the FID portion and the shifted points
%       if thelast point of the FID portion is not at zero. This will cause
%       sinc wiggles in the peaks.
%
%       3. Applying a zero fill to this data will lead to a discontinuity in the data
%       between the rising portion of the shifted points and the zero padding.
%       To circumvent this problem, the shifted points are returned (by circular
%       shift) to the front of the data, the data is zero filled, and then the
%       first points are shifted again to the end of the zero filled data.
%
%       4) The data can now be Fourier transformed and the residual calculated
%       1st order phase correction can be applied.
%
%
%
% For newer versions of the Bruker hardware:
%
%     For firmware versions of 20 <= DSPFVS <= 23 the GRPDLY parameter directly shows the number of
%     points that need to be shifted.
%
%     Thanks for Bruker for supplying this information!
%
%
%
%26-09-'00
