function [time_psi,time_fast] = liuqe_times(shot,psi_type,fast_type,psi_extra,fast_extra)
%LIUQE_TIMES  Determine time array for LIUQE anasrv call
%
% [TIME_PSI,TIME_FAST] = LIUQE_TIMES(SHOT[,PSI_TYPE,FAST_TYPE,PSI_EXTRA,FAST_EXTRA])
%
% LIUQE_TIMES defines the time arrays for LIUQE reconstruction according 
% to the different TYPEs. TIME_PSI is the list of times for which the PSI 
% map will be stored and the PSITBX will be run. TIME_FAST contains all 
% the times for which LIUQE will be run and is a superset of TIME_PSI.
% PSI_EXTRA and FAST_EXTRA are 2 arrays containing extra times to be added
% to TIME_PSI and TIME_FAST respectively. Note elements of both PSI_EXTRA
% and FAST_EXTRA will be added to TIME_PSI.
%
% PSI_TYPE can currently be either:
%    - refined: equally spaced points between start and end of shot with 1/1000 s separation
%    - thomson: use \results::thomson:times if not empty, otherwise 
%               equally spaced points between start and end of shot with 1/60 s separation
%    - mixed: union of fast and thomson [default]
%
% FAST_TYPE can currently be either:
%    - none: TIME_FAST will be equal to TIME_PSI [default]
%    - magnetics: all data points from the magnetics in the [0,3s] window
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(nargin > 0,'First argument (shot) is not specified');
if nargin < 5, fast_extra =      [];end
if nargin < 4, psi_extra  =      [];end
if nargin < 3, fast_type  =  'none';end
if nargin < 2, psi_type   = 'mixed';end
validatestring(lower(psi_type) ,{'refined','thomson','mixed'},mfilename,'psi_type');
validatestring(lower(fast_type),{'none','magnetics'}         ,mfilename,'fast_type');

has_refined = true;
has_thomson = true;
switch lower(psi_type(1))
 case 'r'
  has_thomson = false;
 case 't'
  has_refined = false;
end

% Open shot
meqmdsopen(shot,'tcv_shot');

% Select interval
tmin = 0;
[td,st] = mdsvalue('time_disrupted()');
assert(isodd(st),'Could not retrieve time of disruption for shot %05d',shot);
switch td
 case 512
  warning('TCV shot %05d was aborted and has no useful data, LIUQE time selection will be empty',shot);
  tmax = -Inf;
 case 256
  warning('TCV shot %05d was a PCS/Stray shot, LIUQE time selection will be empty',shot);
  tmax = -Inf;
 case 128
  [tf,st] = mdsvalue('fbte_time("FBTE")');
  assert(isodd(st),'Could not retrieve FBTE times for shot %05d',shot);
  tmax = tf(end) + 0.050; % Add 50ms after planned end of shot
 otherwise
  tmax = td + 0.050; % Add 50ms after time of disruption
end

standard_times = (floor(tmin*60):ceil(tmax*60))/60;

refined_times = [];
if has_refined, refined_times = (floor(tmin*1000):ceil(tmax*1000))/1000;end

thomson_times = [];
if has_thomson
 % Get list of Thomson times ...
 [thomson_times,status] = mdsvalue('\results::thomson:times');
 if ~isodd(status), thomson_times = [];end
 % ... within interval of interest
 mask = thomson_times>=tmin & thomson_times<=tmax;
 thomson_times = thomson_times(mask);
 % Convert to row vector
 thomson_times = thomson_times(:).';
end

% Construct full time array
time_psi = unique([thomson_times,refined_times,psi_extra]);

% If time_psi is empty, replace with the standard times (can only happen
% for type_psi='thomson' and \results::thomson:times is empty)
if isempty(time_psi), time_psi = standard_times;end

% Round them to the nearest magnetics time
mag_times = mdsvalue('\magnetics::bpol_avg:dim0').';
k = iround(mag_times,time_psi);
k(diff(k) == 0) = []; % Remove duplicates (vector is already sorted).
time_psi = mag_times(k);

switch lower(fast_type(1))
 case 'n'
  % Merge time_psi and fast_extra and round to the nearest magnetics time
  time_fast = unique([time_psi,fast_extra]);
  k = iround(mag_times,time_fast);
  k(diff(k) == 0) = []; % Remove duplicates (vector is already sorted).
  time_fast = mag_times(k);
 case 'm'
  % No need to add extras, already all available times
  mask = mag_times >= tmin & mag_times <= tmax;
  time_fast = mag_times(mask);
end

end

