function [complexdata,parameters,inputp,shapepulse]=readpsycheparms_RM(inputp)
% This program has been created by Ludmilla Guduff, Jean-Nicolas Dumez,
% Corentin Jacquemmoz and Achille Marchand
% CEISAM, CNRS, Universit� de Nantes, France
% The authors are not responsible for the results obtained with this program
% This program can be used, modified and shared as long as the authors are
% credited and the program is share under the same conditions
% Please do not delete this header from the file.
%
%Chapter I : read the data and give a structure parameters with all that can be usefull
%Chapter II : load the data for it to be used under the name complexdata 
%-----------------------
    
%Read the data (Fid, ser) from an UfspenDOSY

% Chapter I : to load the parameters


acquparams = readacqus([inputp.pathexp num2str(inputp.diffexpno) '/acqus']); % structure for acquisition parameter ?? why was it working without de C:/ before?
acquparams2 = readacqus([inputp.pathexp num2str(inputp.diffexpno) '/acqu2s']); % structure of ??? if 2D
shapepulse= readacqus([inputp.pathexp num2str(inputp.diffexpno) '/spnam1']); % Structure of shapepulse

% shapepulse parameters for now only chirp

if isfield(shapepulse,'SHL_SW') == 1
    parameters.BW=shapepulse.SHL_SW(1);
    parameters.tau=shapepulse.SHAPE_LENGTH(1)*10^(-6);
else 
    shapeparam = textscan(shapepulse.SHAPE_PARAMETERS, '%s');
    BWind = find(strcmp(shapeparam{1},'[Hz]'),1,'first');
    parameters.BW = str2double(cell2mat(shapeparam{1}(BWind+1)));
    tauind = find(strcmp(shapeparam{1},'[usec]'),1,'first');
    parameters.tau = str2double(cell2mat(shapeparam{1}(tauind+1)))*1e-6;
end

%do not work good for QS yet
 shapepulse_temp=shapepulse.XYPOINTS; 
 shapepulse_temp(1:9)=[]; % ??? 
 shapepulse_temp=str2num(shapepulse_temp);
 parameters.chirpfile=[shapepulse_temp(1:2:end); shapepulse_temp(2:2:end)]';

Gdosy = [acquparams.GPX(1) acquparams.GPY(1) acquparams.GPZ(1)];

if length(Gdosy(Gdosy~=0)) == 1 % Check that acquisition has been done on only one axis
    if Gdosy(1) >0
        parameters.DOSYaxis = 'X';
    elseif Gdosy(2) >0
        parameters.DOSYaxis = 'Y';
    elseif Gdosy(3) >0
        parameters.DOSYaxis = 'Z';
    end
else
    error('The acquisition has been done on a double axis. This option is not implemented yet.')
end

%Chapter II : to merge parameters into the omnistructure of parameters

%F = acquparams2.TD;  % number of interleaves
parameters.Ni=1; % number of interleaves TD2
parameters.N2 = acquparams.L(4); % number of loops in the acquisition
parameters.NPS = acquparams2.TD;
parameters.SW = acquparams.SW_h;% spectral width
parameters.NP = acquparams.TD/2; % total number of points in one fid
parameters.dw = 1/(parameters.SW); % dwell time
parameters.B0 = acquparams.BF1; % spectrometer frequency
parameters.offsetconv=acquparams.O1/parameters.B0; % offset
parameters.D6=acquparams.D(7); % acqusition gradient duration
parameters.td=acquparams.TD; % time domaine size
% offsetconv=acquparams.O1/parameters.B0; % offset
byte_format = acquparams.BYTORDA;
byte_size = acquparams.DTYPA;
parameters.Ta= ((parameters.NP*parameters.dw)/(2*parameters.N2)); %Acquisiton Time
% parameters.offsetconv=offsetconv; % ?
parameters.Gacq=abs(acquparams.GPZ(16)); % acquisition gradient
parameters.gamma1H = nuclist(acquparams.NUC1); % rad/s.T 
parameters.D30=acquparams.D(21)-2*acquparams.P(12)*1e-6; % diffusion time big DELTA (s)
parameters.P1=acquparams.P(2)*10^(-6); % hard pulse duration in us
parameters.dChunk = 2*acquparams.IN(1);
parameters.npChunk = floor(parameters.dChunk / 2 / parameters.Ta);
[GmaxZ,GmaxY,GmaxX] = gradparlist(acquparams.PROBHD);
parameters.Gacq = abs(acquparams.GPZ(3));
% parameters.Gdmax= GmaxZ; %(T/m=10^-2*G/cm) maximum gradient amplitude
% parameters.GPZ0=acquparams.GPZ(11)/100; % strength for excitation gradient
switch parameters.DOSYaxis % find witch gradient to include in parameters
    case 'X'    
        Gacq = abs(acquparams.GPX(3));
        parameters.Gdmax= GmaxX; %(T/m=10^-2*G/cm) maximum gradient amplitude
        parameters.GPZ0=acquparams.GPX(1)/100; % strength for excitation gradient
    case 'Y'
        Gacq = abs(acquparams.GPY(3));
        parameters.Gdmax= GmaxY; %(T/m=10^-2*G/cm) maximum gradient amplitude
        parameters.GPZ0=acquparams.GPY(1)/100; % strength for excitation gradient     
    case 'Z'
        Gacq = abs(acquparams.GPZ(16));
        parameters.Gdmax= GmaxZ; %(T/m=10^-2*G/cm) maximum gradient amplitude
        parameters.GPZ0=acquparams.GPZ(1)/100; % strength for excitation gradient
end


% parameters.offsetconv=offsetconv; % ?
parameters.Gacq=Gacq; % acquisition gradient
parameters.P1=acquparams.P(2)*10^(-6); % hard pulse duration in us
parameters.d11= acquparams.P(12)*1e-6; % extra dephasing period (s)
parameters.dscale=1e-10;% diffusion coefficient scale
parameters.spoff=acquparams.SPOFFS(2); % shapulse offset

if isfield(acquparams, 'SPW') == 1
    parameters.SPW1=acquparams.SPW(2); % watt hard pulse
    parameters.PLW1=acquparams.PLW(2); % watt shape pulse
    parameters.gammaB1=1/(parameters.P1*4)*(parameters.SPW1/parameters.PLW1)^(1/2);% pulse strength
else
    parameters.SP1=acquparams.SP(2);% ?
    parameters.PL1=acquparams.PL(2); % ? 
    parameters.gammaB1=1/(parameters.P1*4)*10^(-(parameters.SP1-parameters.PL1)/20);% ??? pulse strength in dB ???? 
end

if acquparams.GRPDLY>1
    
    
end

% Chapter III : read and load so they may grow 

    complexdata_full = readfid([inputp.pathexp num2str(inputp.diffexpno) '/ser'],byte_format,byte_size); %read the ser and the byte format/size
    %complexdata=complexdata_full((inputp.numscan-1)*parameters.td/2+1:(inputp.numscan)*parameters.td/2); % take only the "inputp.numscan" in the data and throw the rest
    
        %complexdata_full = complexdata ;

% Select the region of fitting based on the diffusion axis
%is it still necessary ? 

switch parameters.DOSYaxis
    case 'Z'
        parameters.fitLengthNeg = inputp.fitLengthLongN;
        parameters.fitLengthPos = inputp.fitLengthLongP;
    case {'X', 'Y'}
        parameters.fitLengthNeg = inputp.fitLengthTransN;
        parameters.fitLengthPos = inputp.fitLengthTransP;
end
% Remove the group delay
    nD=acquparams.GRPDLY;
    chop=floor(nD) ;
    complexdata = circshift(complexdata_full,-chop);
    parameters.p= (nD-chop);
    
end