%MEQWTCV  Writes MEQ data in TCV tree
% MEQWTCV(SHOT,TREE,BRANCH,L,LY,IT_PSI,DEP) writes the
% data contained in L and LY in the BRANCH of TREE for SHOT.
% IT_PSI is the list of LY slices to be processed using PSITBX.
% DEP is the list of node paths to be written in NODES_USED
% subnodes (used in TCV to automatically rerun LIUQE)
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function meqwtcv(shot,tree,branch,L,LY,it_psi,dep)
 
 if nargin < 7, dep = {};end
 if nargin < 6, it_psi = 1:numel(LY.t);end
 
 %% Checks
 assert(isequal(numel(LY.t),numel(unique(LY.t))),'LY.t has some repeated values, this run cannot be stored in MDS');
 
 %%
 
 meqmdsopen(shot,tree,branch)
 
 v = meqver(); % Get meq version number
 
 put = @(node,x,f,d,u,h) meqput(node,x,f,d,u,h,v,dep);
 
 %% grids
 put('R'        ,L.G.rx      ,'f',{'*'},'m','r grid [m]'                                     )
 put('Z'        ,L.G.zx      ,'f',{'*'},'m','z grid [m]'                                     )
 put('RHO'      ,L.pQ        ,'f',{'*'},'' ,'rho grid, 0 on axis, 1 on LCFS'                 )
 put('COCOS'    ,17          ,'s',{'*'},'' ,'COCOS = 17'                                     )
 put('R0'       ,L.P.r0      ,'f',{'*'},'m','Nominal major radius [m]'                       )
 put('TIME_FAST',LY.t        ,'f',{'*'},'s','Time base for scalars, and some 1D profiles [s]')
 put('TIME_PSI' ,LY.t(it_psi),'f',{'*'},'s','Time base for 2D functions [s]'                 )
 
 %% Flux map
 FxB = permute(LY.Fx(:,:,it_psi),[2,1,3]);
 FAB = LY.FA;
 FxB_help = '0 at r=0, for 2*pi radian';
 if L.nD == 1
  FxB = FxB - reshape(LY.FB(it_psi),1,1,[]);
  FAB = FAB - LY.FB;
  FxB_help = '0 on LCFS, for 2*pi radian';
 end
 put('PSI'         ,FxB     ,'f',{'R' 'Z' 'TIME_PSI'},'Wb',['Poloidal flux (r,z,t)[Wb], ',FxB_help]                                                       )
 put('R_SURF'      ,LY.rB   ,'f',{'TIME_FAST'}       ,'m' , 'r position of point defining LCFS (t)[m]'                                                    )
 put('Z_SURF'      ,LY.zB   ,'f',{'TIME_FAST'}       ,'m' , 'z position of point defining LCFS (t)[m]'                                                    )
 put('PSI_SURF'    ,LY.FB   ,'f',{'TIME_FAST'}       ,'Wb', 'Poloidal flux on LCFS (t)[Wb], 0 at r=0, for 2*pi radian'                                    )
 put('R_AXIS'      ,LY.rA   ,'f',{'TIME_FAST'}       ,'m' , 'r position of magnetic axis (t)[m]'                                                          )
 put('Z_AXIS'      ,LY.zA   ,'f',{'TIME_FAST'}       ,'m' , 'z position of magnetic axis (t)[m]'                                                          )
 put('PSI_AXIS'    ,FAB     ,'f',{'TIME_FAST'}       ,'Wb',['Poloidal flux on magnetic axis (t)[Wb], ',FxB_help]                                          )
 put('DR2_PSI_AXIS',LY.dr2FA,'f',{'TIME_FAST'}       ,'Wb', '2nd derivative wrt to r of poloidal flux on magnetic axis (t)[Wb/m^2], for 2*pi radian'      )
 put('DZ2_PSI_AXIS',LY.dz2FA,'f',{'TIME_FAST'}       ,'Wb', '2nd derivative wrt to z of poloidal flux on magnetic axis (t)[Wb/m^2], for 2*pi radian'      )
 put('DRZ_PSI_AXIS',LY.drzFA,'f',{'TIME_FAST'}       ,'Wb', '2nd derivative wrt to r and z of poloidal flux on magnetic axis (t)[Wb/m^2], for 2*pi radian')
 put('HAS_XPTS'    ,LY.lX   ,'s',{'TIME_FAST'}       ,''  , 'Is diverted (t)'                                                                             )
 
 %% Source terms
 Jx = zeros(L.nrx,L.nzx,length(it_psi)); Jx(2:end-1,2:end-1,:) = permute(1/L.dsx*LY.Iy(:,:,it_psi),[2 1 3]);
 PA = LY.PQ(1,:)';
 % First domains are axis domains
 nA = size(LY.rA,1);
 JA = 2*pi*(LY.rA.*squeeze(LY.PpQ(1,1:nA,:))+squeeze(LY.TTpQ(1,1:nA,:))./mu0./LY.rA);
 put( 'PPRIME_BFCT',LY.PpQg ,'f',{'RHO' '*' 'TIME_FAST'},'Pa/Wb','p'' base functions multiplied by their coefficients (rho,*,t)[Pa/Wb]' )
 put('TTPRIME_BFCT',LY.TTpQg,'f',{'RHO' '*' 'TIME_FAST'},'T*m^4','TT'' base functions multiplied by their coefficients (rho,*,t)[T*m^4]')
 put( 'PPRIME_RHO' ,LY.PpQ  ,'f',{'RHO'     'TIME_FAST'},'Pa/Wb','p'' profile (rho,t)[Pa/Wb]'                                           )
 put('TTPRIME_RHO' ,LY.TTpQ ,'f',{'RHO'     'TIME_FAST'},'T*m^4','TT'' profile (rho,t)[T*m^4]'                                          )
 put('P_RHO'       ,LY.PQ   ,'f',{'RHO'     'TIME_FAST'},'Pa'   ,'Pressure profile (rho,t)[Pa], 0 on LCFS'                              )
 put('RBTOR_RHO'   ,LY.TQ   ,'f',{'RHO'     'TIME_FAST'},'m*T'  ,'T = r*Btor profile (rho,t)[m*T]'                                      )
 put('P_AXIS'      ,PA      ,'f',{          'TIME_FAST'},'Pa'   ,'Pressure on axis (t)[Pa]'                                             )
 put('J_TOR'       ,Jx      ,'f',{'R'   'Z' 'TIME_PSI '},'A/m^2','Toroidal plasma current density (r,z,t)[A/m^2]'                       )
 put('J_TOR_AXIS'  ,JA      ,'f',{          'TIME_FAST'},'A/m^2','Toroidal plasma current density on axis (t)[A/m^2]'                   )
 put('I_POL'       ,LY.Ia   ,'f',{      '*' 'TIME_FAST'},'A'    ,'Fitted poloidal field coil currents (*,t)[A]'                         )
 put('I_VESSEL'    ,LY.Iu   ,'f',{      '*' 'TIME_FAST'},'A'    ,'Fitted vessel segment currents (*,t)[A]'                              )
 
 %% Fitted signals
 put('B_PROBE'     ,LY.Bm   ,'f',{      '*' 'TIME_FAST'},'T'    ,'Fitted magnetic probe measurements (*,t)[T]'              )
 put('PSI_LOOP'    ,LY.Ff   ,'f',{      '*' 'TIME_FAST'},'Wb'   ,'Fitted flux loop poloidal flux (*,t)[Wb], for 2*pi radian')
 put('TOR_FLUX_DML',LY.Ft   ,'f',{          'TIME_FAST'},'Wb'   ,'Fitted toroidal diamagnetic flux (t)[Wb]'                 )
 
 %% Volume integrals
 put('RBTOR_VAC'   ,LY.rBt  ,'f',{          'TIME_FAST'},'m*T'  ,'Vacuum toroidal field * r (t)[m*T]'              )
 put('I_PL'        ,LY.Ip   ,'f',{          'TIME_FAST'},'A'    ,'Total plasma current (t)[A]'                     )
 put('VOL_PL'      ,LY.Vp   ,'f',{          'TIME_FAST'},'m^3'  ,'Volume delimited by LCFS (t)[m^3]'               )
 put('TOR_FLUX_VAC',LY.Ft0  ,'f',{          'TIME_FAST'},'Wb'   ,'Vacuum contribution to the toroidal flux (t)[Wb]')
 put('TOR_FLUX_PL' ,LY.Ft   ,'f',{          'TIME_FAST'},'Wb'   ,'Plasma contribution to the toroidal flux (t)[Wb]')
 put('W_MHD'       ,LY.Wk   ,'f',{          'TIME_FAST'},'J'    ,'Total plasma kinetic energy (t)[J]'              )
 put('W_TOR_VAC'   ,LY.Wt0  ,'f',{          'TIME_FAST'},'J'    ,'Vacuum toroidal magnetic energy (t)[J]'          )
 put('W_TOR_PL'    ,LY.Wt   ,'f',{          'TIME_FAST'},'J'    ,'Plasma toroidal magnetic energy (t)[J]'          )
 put('W_POL_PL'    ,LY.Wp   ,'f',{          'TIME_FAST'},'J'    ,'Plasma poloidal magnetic energy (t)[J]'          )
 put('W_NORM'      ,LY.WN   ,'f',{          'TIME_FAST'},'J'    ,'Normalisation energy (t)[J], mu0*r0*Ip^2/4'      )
 put('BETA_POL'    ,LY.bp   ,'f',{          'TIME_FAST'},''     ,'Beta poloidal (t)'                               )
 put('BETA_TOR'    ,LY.bt   ,'f',{          'TIME_FAST'},''     ,'Beta toroidal (t)'                               )
 put('L_I_3'       ,LY.li   ,'f',{          'TIME_FAST'},''     ,'Normalised internal inductance (t)'              )
 put('DIAMAG'      ,LY.mu   ,'f',{          'TIME_FAST'},''     ,'Normalised diamagnetism (t)'                     )
 put('LAMBDA'      ,LY.bpli2,'f',{          'TIME_FAST'},''     ,'Betap + li/2 (t)'                                )
 
 %% Iteration status
 put('RMS_ERROR'   ,LY.chi  ,'l',{          'TIME_FAST'},''     ,'Mean weighted error for fitter measurements'     )
 put('N_ITER'      ,LY.niter,'l',{          'TIME_FAST'},''     ,'Number of Picard iterations to reach convergence')
 
 %% Unused nodes
 % put('ERROR_CODE',,'l',{          'TIME_FAST'},''     ,'Some description of error')
 % put('P_FIT'     ,,'l',{          'TIME_FAST'},''     ,'Pressure measurements fit')
 % put('VLOOP_SURF',,'l',{          'TIME_FAST'},''     ,'Some time derivative of the LCFS flux')
 
end
