%%NMRMasterProcessing    1-18-2022      Micah Jeppesen
#Get Files
display('Using typical unix filename syntax (i.e. ./* or ../your/file/is/here/* or ???). This will only use directories.')
filepath = input('What is the path to your data? ',"s");
F.fileglob = glob(filepath);  %% Number of ? is the number of digits of the file name

%% We only use directories here
count = 1;
sizeoffileglob = numel(F.fileglob);
currentdir = pwd;
for i = [1:sizeoffileglob]
  if isfolder(F.fileglob{i})
	  cd(F.fileglob{i});
    if (stat('ser') || stat('fid'))
      continue;
    end
    cd(currentdir);
    F.dirs{count++} = F.fileglob{i};
  end
end
checkforexistence = exist('F.dirs')
if checkforexistence
  error('There is no NMR data here, seems you have the wrong filepath.')
end


#Extract Data
%% This could be loadnmr or loaddmx - default should be loadnmr
display('With "@loadnmr" and "@loaddmx" as your choices. Refer to the manual for these functions.');
loadingfunc = input('Which would you like to use? ', "s");
if loadingfunc ~= '@loadnmr' || loadingfunc ~= '@loaddmx'
  error('A loading function was not chosen correctly');
end
[F.data, F.parms, F.t] = feval(loadingfunc, F.dirs); %% Assumes F.dirs were all run with similar parameters

#Apodization, Zerofill and Fourier Transform

apodizing = input('Would you like to apodize?(y,n) ', "s")
if apodizing == 'y' || apodizing == 'Y'
  display('With "@expwindow", "@gausswindow", and "@sinewindow" being the options.(Refer to the manual)')
  functionhandleofapodization = input('Which function would you like to use? ',"s"); %% User choice on function handle
  F.data = apodize(F.data,F.parms,functionhandleofapodization);
end

zerofilling = input('Would you like to zerofill?(y,n) ', "s");
if zerofilling == 'y' || zerofilling == 'Y'
  numberofzerofills = input('How many zerofills would you like? '); %% User choice on number
  F.data = zerofill(F.data,F.parms,numberofzerofills);
end
[S.data, S.ppm] = nmrft(F.data, F.parms);%% finally do something

#Phase and extract real
%% phase correction method objective can be decided
display('Automatic phasing will be performed, you can choose an objective function. Your options are "@simplex_minimum", "@simplex_entropy", "@simplex_whiten", or "@simplex_integral". (Refer to manual)')
phaseobjective = input('Which objective do you choose? ', "s");
[S.data, S.phc0, S.phc1] = autophase(S.data, F.parms, phaseobjective);
[S.data, S.phc0, S.phc1] = autophase(S.data, F.parms, phaseobjective);
X.data = realnmr(S.data, F.parms); %% extract numerical data
X.ppm = S.ppm;

#Alignment and Referencing
%% an if could be used to decide if icosshift or cosshift here**
display('With "@coshift" and "@icoshift" as your options.')
alignmentoption = input('Which alignment method would you like to use? ',"s");
X.data = feval(alignmentoption, X.data, X.ppm);
X.data = feval(alignmentoption, X.data, X.ppm);
X.data = feval(alignmentoption, X.data, X.ppm);
plot(X.ppm, X.data);

%% an automatic referencing method could be used instead of this (no automatic method currently available)
adjust_zero = 1
zero_adjustment = 0
  while adjust_zero

  	X.ppm = refadj(X.ppm, zero_adjustment, 0.0);

  	plot(X.ppm, X.data);

  	adjust_input = input('Do you need to zero?(y/n)  ',"s");
  	if adjust_input == "Y" || adjust_input == "y"
          	adjust_zero = 1;
  		      zero_adjustment = input('How much Adjustment?  ');
    else
          	adjust_zero = 0;
    end
  end

%% this can be selected as the manual removal method - otherwise ROI methods could be used?**
  remove_var = 1;
  while remove_var

  	remove_input = input('Do you need to remove data?(y/n)  ',"s");
  	if remove_input == "Y" || remove_input == "y"
          	remove_var = 1;
  		      kmin_input = input('What is the minumum of the range?  ');
            kmax_input = input('What is the maximum of the range?  ');
            idxmin = findnearest(X.ppm, kmax_input);
            idxmax = findnearest(X.ppm, kmax_input);
            X.rmvar = [idxmax:idxmin];
            [X.data, X.ppm] = rmvar(X.data, X.ppm, X.rmvar);
            plot(X.ppm, X.data);
    else
          	remove_var = 0;
    end
  end

%% binning can be done a number of ways we default to adaptive binning***
display('There are multiple ways of binning, this walkthrough uses binadapt.')
[X.bindata, X.binppm] = binadapt(X.data, X.ppm, F.parms);
%% normalization has multiple options these need to be considered ***
X.bindata = histmatch(X.bindata);

display('To this point, the important data is extracted and binned as "X.bindata" and can be analyzed');
