function [L,LY] = mds2meq(shot,source,server,varargin)
% function [L,LY] = mds2meq(shot,source,server,varargin)
% Load data from mds into meq structure
% source: source of equilibrium results. 'FBTE', 'LIUQE.M' load equilibrium results from 
%   resp. FBTE and LIUQE.M nodes. 'MAG' loads currents and voltages from magnetics data directly
%   'RTLIUQE' loads equilibria from RTC mems for shots when this was available (use server='scd')
% server: optional non-default server (empty gives default)
% varargin: optional additional arguments passed to liuqe() or fbt() calls used
%     to generate L
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

assert(any(contains(upper(source),{'FBTE','LIUQE.M','MAG','RTLIUQE'})),...
  'source %s not supported',source);

if nargin==2, server = []; end

switch source
  case 'RTLIUQE'
    tree = 'RTC';
  case 'FBTE'
    tree = 'PCS';
  otherwise
    tree = 'TCV_SHOT';
end
meqmdsopen(shot,tree,[],server);

LY.shot = shot;

% Data
switch upper(source)
  case 'LIUQE.M'
    % Determine version of code used for stored results
    [lver,st] = mdsvalue('\results::top.equil_1.results:i_pl:version_num');
    assert(isodd(st),'Could not retrieve LIUQE version information for this run');
    % Prior to version 4 deployment, the legacy LIUQE tree was used for parameters
    if lver < 4, psrc = 'LIUQE';
    else,        psrc = 'LIUQE.M';
    end
    % L structure
    L = liuqe(shot,[],'psrc',psrc,varargin{:});
    meqmdsopen(shot,tree,[],server);
    % time & mask
    LY.t = meqy('TIME_PSI',source); % time to microsecond precision
    fast = meqy('TIME_FAST',source);
    imask = iround(fast,LY.t);

    LY.rB = meqy('R_SURF',source,imask);
    LY.zB = meqy('Z_SURF',source,imask);
    LY.lB = ~isnan(LY.rB)&~isnan(LY.zB);
    LY.FB = meqy('PSI_SURF',source,imask);
    LY.dr2FA = meqy('DR2_PSI_AXIS',source,imask);
    LY.dz2FA = meqy('DZ2_PSI_AXIS',source,imask);
    LY.drzFA = meqy('DRZ_PSI_AXIS',source,imask);
    LY.lX    = meqy('HAS_XPTS',source,imask);
    
    % profiles
    LY.PpQg  = meqy( 'PPRIME_BFCT',source,imask);
    LY.TTpQg = meqy('TTPRIME_BFCT',source,imask);
    
    LY.PpQ  = meqy( 'PPRIME_RHO',source,imask);
    LY.TTpQ = meqy('TTPRIME_RHO',source,imask);
    LY.PQ   = meqy('P_RHO',source,imask);
    LY.TQ   = meqy('RBTOR_RHO',source,imask);
    
    % currents
    LY.Ia = meqy('I_POL',source,imask);
    LY.Iu = meqy('I_VESSEL',source,imask);
    % Fitted signals
    LY.Bm = meqy('B_PROBE',source,imask);
    LY.Ff = meqy('PSI_LOOP',source,imask);
    LY.Ft = meqy('TOR_FLUX_DML',source,imask);
    
    % Volume integrals
    LY.rBt = meqy('RBTOR_VAC',source,imask);
    LY.Ip  = meqy('I_PL',source,imask);
    LY.Ft0 = meqy('TOR_FLUX_VAC',source,imask);
    % LY.Ft  = meqy('TOR_FLUX_PL',source,mask);
    LY.Wk = meqy('W_MHD',source,imask);
    LY.Wt0 = meqy('W_TOR_VAC',source,imask);
    LY.Wt = meqy('W_TOR_PL',source,imask);
    LY.Wp = meqy('W_POL_PL',source,imask);
    LY.WN = meqy('W_NORM',source,imask);
    LY.bp = meqy('BETA_POL',source,imask);
    LY.bt = meqy('BETA_TOR',source,imask);
    LY.qA = meqy('Q_AXIS',source,imask);
    LY.li = meqy('L_I_3',source,imask);
    LY.mu = meqy('DIAMAG',source,imask);
    LY.bpli2 = meqy('LAMBDA',source,imask);
    LY.Vp = meqy('VOL_EDGE',source,imask);
    % Algorithm outputs
    LY.dz = meqy('DELTA_Z',source,imask);
  case 'FBTE'
    % find FBTE version
    meqmdsopen(shot,'PCS',[],'tcvdata')
    FBTEversion = mdsvalue('get_fbte_version()');
    if strcmp('FBTE',FBTEversion)
      npq = mdsvalue('\pcs::mgams.data:npsi');
      pq = sqrt((1:npq)/npq); % specific FBT fortran rho grid
      % Scale to SI units for p' and tt' profiles
      PpScale = 1/2/pi;
      TTpScale = 2e-7;
    elseif strcmp('FBTE.M',FBTEversion) % fbte matlab has rho grid node
      pq = mdsvalue('\pcs::mgams.equil:rho'); npq=numel(pq)-1;
      % Scale to SI units for p' and tt' profiles
      PpScale = 1;
      TTpScale = 1;
    else
      error('unknown FBTE version of shot is not open?')
    end
    % L structure
    L = fbt('TCV',shot,[],'pq',pq,'npq',npq,varargin{:});
    meqmdsopen(shot,tree,[],server) % previous operation leaves BASE/shot=-1 open

    % time & mask
    [t,status] = mdsvalue('dim_of(tcv_eq("I_P",$1))',source); % there must be a better way..
    assert(isodd(status),'Could not retrieve time for source %s',source);
    LY.t = reshape(t,1,numel(t));
    nt = numel(LY.t);
    imask= 1:nt;

    % flux
    LY.FB = zeros(1,nt);
    
    % profiles (includes scale correction for FBTE-FORTRAN)
    LY.PpQ   = meqy( 'PPRIME_PSI',source,imask)*PpScale;
    LY.TTpQ  = meqy('TTPRIME_PSI',source,imask)*TTpScale;
    
    % integrals
    LY.Ip = meqy('I_P',source,imask);
    LY.Ft = meqy('TOR_FLUX',source,imask);
    LY.bp = meqy('BETA_POL',source,imask);
    LY.qA = meqy('Q_ZERO',source,imask);
    LY.bt = meqy('BETA_TOR',source,imask);
    LY.li = meqy('L_I',source,imask);
    LY.rBt = 0.88*meqy('BZERO',source,imask); % don't know where else to get this
    LY.Vp = meqy('VOLUME',source,imask);
    
    % currents - value stored is always assuming positive Ip so need to correct
    Ipol = meqy('I_POL',source,imask);
    LY.Ia = Ipol(1:16,:); % don't know if there's a better way
    LY.Iu = zeros(0,nt);
  case 'MAG'
    % time
    t = mdsvalue('\MAGNETICS::IPHI:DIM0')';
    
    % L structure plus LIUQE measurements
    [L,LY] = liuqe(shot,t,'selu','n','ivesm',0,varargin{:});
    nt = numel(LY.t);
    
    % flux
    LY.FB = zeros(0,nt);
    
    % profiles (correct for TT' definition in nodes and 2pi error(?)
    LY.PpQ   = zeros(11,nt);
    LY.TTpQ  = zeros(11,nt);
    
    % integrals
    LY.Ip = zeros(1,nt);
    LY.Ft = zeros(1,nt);
    LY.bp = zeros(1,nt);
    LY.qA = zeros(0,nt);
    LY.bt = zeros(1,nt);
    LY.li = nan(1,nt);
    LY.Vp = zeros(1,nt);
    LY.Iv = zeros(L.G.nv,nt);
    
    meqmdsopen(shot,tree,[],server);
    for ii=1:8
      Uref_E(:,ii) = mdsvalue(sprintf('\\top.power.measurements.red.udref.e_%03d',ii)); %#ok<*AGROW>
      Uref_F(:,ii) = mdsvalue(sprintf('\\top.power.measurements.red.udref.f_%03d',ii));
    end
    for ii=1:2
      U_OH(:,ii) = mdsvalue(sprintf('\\top.power.measurements.red.udref.oh_%03d',ii));
    end
    V_G(:,1) = mdsvalue('\POWER::FPS:UFOS');
    
    t_ef(:,1) = mdsvalue('dim_of(\top.power.measurements.red.udref.e_001)');
    t_oh(:,1) = mdsvalue('dim_of(\top.power.measurements.red.udref.oh_001)');
    t_V_G(:,1) = mdsvalue('dim_of(\power::fps:UFOS)');
    
    e_scale = mdsvalue('\hybrid::output_scales_e_u');
    f_scale = mdsvalue('\hybrid::output_scales_f_u');
    oh_scale =mdsvalue('\hybrid::output_scales_oh_u');
    
    LY.Va = [e_scale*interp1(t_ef,Uref_E,LY.t,'linear',0)';
      f_scale*interp1(t_ef,Uref_F,LY.t,'linear',0)';
      interp1(t_V_G,V_G,LY.t,'linear',0);
      oh_scale*interp1(t_oh,U_OH,LY.t,'linear',0)'];
  
  case 'RTLIUQE'
    if shot<72453
      L = liuqe(shot,[],'bfct',@bf3pmex,'bfp',false,...
        'idml',1, ...
        'itert',1,... % Picard iterations per time step
        'iterq',1,... % we want post-processing
        'slx',1,... % Simulink version
        'npq',8,... % radial grid excluding LCFS and axis
        'pq',[],...
        'iterh',0,...
        varargin{:});
    else
      L = liuqe(shot,[],...
        'psrc','RTLIUQE',... % Use parameters stored in \PCS::TOP.RTLIUQE_1.PARAMETERS
        'slx',1,...          % Enforce Simulink version
        'bslvdprec',true,... % Use double precision for (Ie,ag) LSQ fit
        varargin{:});
    end
    
    assert(~isempty(which('get_scd_mems')),...
      'liuxrtc:noscd','did not find get_scd_mems.m - need to add SCD tools to the path')
    LYts = getfield(get_scd_mems(shot,'LY','verbosity',0),'LY');
    assert(~isempty(LYts),'no LY found in mems for rtc shot %d',shot);
    LY = ts2meq(L,LYts);
    % customizations for missing fields
    assert(~isfield(LY,'qA'),'RTLIUQE now giving qA directly?')
    if isfield(LY,'iqQ') && L.pQ(1)==0
      LY.qA = 1./LY.iqQ(1,:);
    end
end

switch source
  case {'FBTE','LIUQE.M'}
    %% Common quantities
    
    % Flux
    FxB   = meqy('PSI',source); % TIME_PSI only
    LY.Fx = permute(FxB+reshape(LY.FB,1,1,[]),[2 1 3]);
    LY.rA = meqy('R_AXIS',source,imask);
    LY.zA = meqy('Z_AXIS',source,imask);
    FAB   = meqy('PSI_AXIS',source,imask);
    LY.FA = FAB + LY.FB;
    LY.nA = int32(sum(~isnan(LY.rA),1));
    
    % Source terms
    
    Jx = meqy('J_TOR',source); % TIME_PSI only
    LY.Iy    = permute(Jx(2:end-1,2:end-1,:),[2,1,3])*L.dsx;
    
    LY.Iv = L.G.Tvu*LY.Iu;
    
  case 'MAG'
    % Flux
    LY.Fx = zeros(L.nzx,L.nrx,nt);
    LY.Iy = zeros(L.nzy,L.nry,nt);
    LY.rA = zeros(0,nt);
    LY.zA = zeros(0,nt);
    LY.FB = zeros(0,nt);
    LY.FA = zeros(0,nt);
    L.pQ = linspace(0,1,11);
    
  otherwise
    % nothing
end
  
end

function Y = meqy(node,source,imask)
[Y,stat] = mdsdata('tcv_eq($1,$2)',node,source);
assert(isodd(stat) && ~isempty(Y),'LIUQE:MDS:Cannot read result %s for source %s',node,source)
applymask = nargin>2;
if applymask
  % Apply mask on the time dimension (the last one)
  ndim = find(size(Y) == numel(imask),1,'last');
  if ~isequal(imask,1)
    assert(~isempty(ndim),'no dimension size matches time dimension %d',numel(imask))
    S = substruct('()',[repmat({':'},1,ndim-1),{imask(:)}]);
    Y = subsref(Y,S);
  end
end

% reorient (nt x 1) vectors to be row vectors (meq convention)
if (~applymask || numel(imask)>1) && iscolumn(Y)
  Y = Y';
end
end
