function [brukerdata]=readspendosy_GNAT(FilePath,brukerdata)
%*************************************************************************%
%[brukerdata]=readspendosy_GNAT(FilePath,brukerdata)
%reads the SPEN DOSY data and save them in brukerdata structure
%       -Input parameters
%              - FilePath - the path of the file
%              - brukerdata - uses the brukerdata structure which is saved
%                            via inputpspendosy_GNAT and brukerimport_SPEN
%  
%       -Outut parameters
%              - brukerdata - parameters (spectral width, dwell time,time
%                             domain size, big delta, gradient duration,no.  
%                             of loops, no. of interleaves, etc.)are saved.
%*************************************************************************%
brukerdata.acqu = readacqus([FilePath '/acqu']);
brukerdata.numscan = brukerdata.acqu.NS;

% Load the data
% complexdata_full = readfid_RM([pathexp num2str(expno) '/ser']);
if exist([FilePath '/acqu2s'], 'file') == 2
    brukerdata.DataDim=2;
elseif exist([FilePath '/acqus'], 'file') == 2
    brukerdata.DataDim=1;
end

if brukerdata.DataDim>=1
    brukerdata.acqus = readacqus([FilePath '/acqus']);
    brukerdata.acqus.PULPROG(1)=[];
    brukerdata.acqus.PULPROG(end)=[];
    brukerdata.pulprog=brukerdata.acqus.PULPROG;
    logtext=['pulse programme is: ' brukerdata.pulprog];
    disp(logtext)
    
    % % % % % Added by Rituraj % % % % %
    %     if strfind(parameters.pulprog,'dosyuf')~=0
    if exist([FilePath '/spnam1'], 'file') == 2
        brukerdata.shape=readacqus([FilePath '/spnam1']);
        brukerdata.sf=str2double(brukerdata.acqus.BF1);
        brukerdata.sfo1=str2double(brukerdata.acqus.SFO1);
        % to prepare data for choosing Geff in spen2GNAT according to 
        % linear or quadratic pulse
        A = textread([FilePath '/spnam1'],'%s','whitespace','\n');
        Index = find(contains(A,'XYPOINTS'));
        A = A (Index+1:end-1);
        A = string(A);
        amp = str2num(A((size(A,1)/2)));
        amp1 = str2num(A((size(A,1)/2+10)));
        amp2 = str2num(A((size(A,1)/2-10)));
        amp = amp(1); amp1 = amp1(1); amp2 = amp2(1);
        if (amp-amp1==0)|| (amp-amp2==0)
            brukerdata.islinearchirp = 1;
        else
            brukerdata.islinearchirp = 0;
        end
     % ending here
    else
        brukerdata.SW=str2double(brukerdata.acqus.SW);
    end
    % % % % % end % % % % %
end

if brukerdata.DataDim>=2
    brukerdata.acqu2s = readacqus([FilePath '/acqu2s']);
    if isfield(brukerdata.acqus, 'SPW') == 1
        brukerdata.SP1=brukerdata.acqus.SPW(2);
        brukerdata.PL1=brukerdata.acqus.PLW(2);
    else
        brukerdata.SP1=brukerdata.acqus.SP(2);
        brukerdata.PL1=brukerdata.acqus.PL(2);
    end
end

% chirp parameters
if exist([FilePath '/spnam1'], 'file') == 2
    if isfield(brukerdata.shape,'SHL_SW') == 1
        brukerdata.BW=brukerdata.shape.SHL_SW(1);
        brukerdata.Te=brukerdata.shape.SHAPE_LENGTH(1)*10^(-6);
    else
        shapeparam = textscan(brukerdata.shape.SHAPE_PARAMETERS, '%s');
        BWind = find(strcmp(shapeparam{1},'[Hz]'),1,'first');
        brukerdata.BW = str2double(cell2mat(shapeparam{1}(BWind+1)));
        tauind = find(strcmp(shapeparam{1},'[usec]'),1,'first');
        brukerdata.tau = str2double(cell2mat(shapeparam{1}(tauind+1)))*1e-6;
    end
    
    chirpfile_temp=brukerdata.shape.XYPOINTS;
    chirpfile_temp(1:9)=[];
    chirpfile_temp=str2num(chirpfile_temp);
    brukerdata.chirpfile=[chirpfile_temp(1:2:end); chirpfile_temp(2:2:end)]';
end
%Add by CJ to select automatically the axis of DOSY
Gdosy = [brukerdata.acqus.GPX(3) brukerdata.acqus.GPY(3) brukerdata.acqus.GPZ(3)];
if sum(Gdosy(:)) ==0 % Check that diffusion has been done on only one axis
    error('The acquisition has been done on a double axis. The calculated coefficient values would be wrong.')
elseif Gdosy(1) >0
    brukerdata.DOSYaxis = 'X';
elseif Gdosy(2) >0
    brukerdata.DOSYaxis = 'Y';
elseif Gdosy(3) >0
    brukerdata.DOSYaxis = 'Z';
end
[GmaxZ,GmaxY,GmaxX] = gradparlist(brukerdata.acqus.PROBHD);
switch brukerdata.DOSYaxis
    case 'X'
        Gacq = abs(brukerdata.acqus.GPX(3));
        brukerdata.Gdmax= GmaxX; %(T/m=10^-2*G/cm) maximum gradient amplitude
        brukerdata.GPZ0=brukerdata.acqus.GPX(1)/100; % strength for excitation gradient
    case 'Y'
        Gacq = abs(brukerdata.acqus.GPY(3));
        brukerdata.Gdmax= GmaxY; %(T/m=10^-2*G/cm) maximum gradient amplitude
        brukerdata.GPZ0=brukerdata.acqus.GPY(1)/100; % strength for excitation gradient
    case 'Z'
        Gacq = abs(brukerdata.acqus.GPZ(3));
        brukerdata.Gdmax= GmaxZ; %(T/m=10^-2*G/cm) maximum gradient amplitude
        brukerdata.GPZ0=brukerdata.acqus.GPZ(1)/100; % strength for excitation gradient
end
if brukerdata.DataDim>=1
    brukerdata.SW = brukerdata.acqus.SW_h; % spectral width
    brukerdata.NP = brukerdata.acqus.TD/2; % total number of points in one fid
    dw = 1/(brukerdata.SW); % dwell time
    brukerdata.dw=dw;
    brukerdata.byte_format = brukerdata.acqus.BYTORDA;
    brukerdata.byte_size = brukerdata.acqus.DTYPA;
    brukerdata.Gacq=Gacq; % acquisition gradient
    brukerdata.gamma1H = nuclist(brukerdata.acqus.NUC1); % rad/s.T
end
if brukerdata.DataDim>1
    % Some parameters are needed for processing and display
    Ni = brukerdata.acqu2s.TD;  % number of interleaves
    brukerdata.Ni=Ni;
    N2 = brukerdata.acqus.L(4); % number of loops in the acquisition
    brukerdata.N2=N2;
    brukerdata.B0 = brukerdata.acqus.BF1; % spectrometer frequency
    offsetconv=brukerdata.acqus.O1/brukerdata.B0; % offset
    D6=brukerdata.acqus.D(7); % acqusition gradient duration
    brukerdata.D6=D6;
    brukerdata.td=brukerdata.acqus.TD; % time domaine size
    brukerdata.offsetconv=offsetconv;
    brukerdata.D30=brukerdata.acqus.D(31); % diffusion time big DELTA (s)
    brukerdata.P1=brukerdata.acqus.P(2)*10^(-6);
    brukerdata.d11= brukerdata.acqus.D(12); % extra dephasing period (s)
    brukerdata.dscale=1e-10;% diffusion coefficient scale
    brukerdata.spoff=brukerdata.acqus.SPOFFS(2);
    brukerdata.Ta= ((brukerdata.NP*dw)/(2*N2));
end
end