%PSIWTCV  Write PSITBX specific data in TCV tree
% PSIWTCV(SHOT,TREE,BRANCH,L,LY,IT_PSI,DEP) writes the
% results of the PSITBX analysis of MEQ data 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 psiwtcv(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);
 
 nthg = 129;
 if isfield(L.P,'nthg' ), nthg = L.P.nthg ;end
 
 psi = meq2psi(L,LY);
 fsd = psitbxp2p(subs_time(psi,it_psi),'FS',psitbxgrid('Flux','Grid',{L.pQ,linspace(-pi,pi,nthg),NaN}));
 rzs = psitbxg2g(fsd.grid,'C',fsd);
 fsg = psitbxfsg(fsd);
 khfs = iround(fsd.grid.x{2},        0);
 klfs = iround(fsd.grid.x{2},       pi);
 k95  = iround(fsd.grid.x{1},sqrt(.95));
 X.R_RHO        = 'SET_RANGE(SHAPE(RHO,0),SHAPE(THETA,0),SHAPE(TIME_PSI,0),REPLICATE(R_AXES[0,*],0,SHAPE(RHO,0)*SHAPE(THETA,0)))-SPREAD(SPREAD(COS(THETA),0,SHAPE(RHO,0)),2,SHAPE(TIME_PSI,0))*A_RHO';
 X.Z_RHO        = 'SET_RANGE(SHAPE(RHO,0),SHAPE(THETA,0),SHAPE(TIME_PSI,0),REPLICATE(Z_AXES[0,*],0,SHAPE(RHO,0)*SHAPE(THETA,0)))+SPREAD(SPREAD(SIN(THETA),0,SHAPE(RHO,0)),2,SHAPE(TIME_PSI,0))*A_RHO';
 X.A_EDGE       = squeeze(fsd.x   (end,   :,:));
 X.R_EDGE       = squeeze(rzs.x{1}(end,   :,:));
 X.Z_EDGE       = squeeze(rzs.x{2}(end,   :,:));
 X.R_IN_MID     = squeeze(rzs.x{1}(:  ,khfs,:));
 X.R_OUT_MID    = squeeze(rzs.x{1}(:  ,klfs,:));
 X.A_MINOR_MID  = (X.R_OUT_MID-X.R_IN_MID)*0.5;
 X.R_GEOM_MID   = (X.R_OUT_MID+X.R_IN_MID)*0.5;
 X.KAPPA_AXIS   = sqrt(fsd.dr2mag./fsd.dz2mag);
 X.KAPPA        = double((fsg.ztop - fsg.zbot) ./ (fsg.rout - fsg.rin)); X.KAPPA(1,:) = X.KAPPA_AXIS;
 X.KAPPA_95     = X.KAPPA(k95,:);
 X.KAPPA_EDGE   = X.KAPPA(end,:);
 X.DELTA_TOP    = double((fsg.rin + fsg.rout - 2*fsg.rtop) ./ (fsg.rout - fsg.rin)); X.DELTA_TOP(1,:) = 0;
 X.DELTA_BOT    = double((fsg.rin + fsg.rout - 2*fsg.rbot) ./ (fsg.rout - fsg.rin)); X.DELTA_BOT(1,:) = 0;
 X.DELTA        = (X.DELTA_TOP + X.DELTA_BOT)*0.5;
 X.DELTA_95_TOP = X.DELTA_TOP(k95,:);
 X.DELTA_95_BOT = X.DELTA_BOT(k95,:);
 X.DELTA_95     = X.DELTA    (k95,:);
 X.DELTA_ED_TOP = X.DELTA_TOP(end,:);
 X.DELTA_ED_BOT = X.DELTA_BOT(end,:);
 X.DELTA_EDGE   = X.DELTA    (end,:);
 % Temporary fix for doublets, treat only the upper half
 if L.P.idoublet
  X.Q           = fsg.ugpsiiri.x .* reshape(LY.TQ(:,1,it_psi),L.nQ,numel(it_psi))/(2*pi);
 else
  X.Q           = fsg.ugpsiiri.x .* LY.TQ(:,it_psi)/(2*pi);
 end
 %
 X.Q_AXIS       = X.Q(  1,:);
 X.Q_95         = X.Q(k95,:);
 X.Q_EDGE       = X.Q(end,:);
 X.TOR_FLUX_TOT = cumtrapz(2*mean(diff(fsd.grid.x{1}))*fsd.grid.x{1}(:)*fsd.psimag.*X.Q,1); X.TOR_FLUX_TOT(end,:) = LY.Ft0(it_psi) + LY.Ft(it_psi);

 %% Grids
 put('RHO'      ,fsd.grid.x{1},'f',{'*'},''   ,'rho grid, 0 on axis, 1 on LCFS'                 )
 put('THETA'    ,fsd.grid.x{2},'f',{'*'},'Rad','theta grid, 0 on HFS, CW [Rad]'                 )
 put('TIME_FAST',LY.t         ,'f',{'*'},'s'  ,'Time base for scalars, and some 1D profiles [s]')
 put('TIME_PSI' ,fsd.t        ,'f',{'*'},'s'  ,'Time base for 2D functions [s]'                 )
 
 %% Flux surface geometry
 put('R_AXES'      ,fsd.rmag      ,'f',{'*'          ,'TIME_PSI'},'m'     ,'r position of magnetic axes, sorted by absolute flux, 0 on LCFS (*,t) [m]'                                                              )
 put('Z_AXES'      ,fsd.zmag      ,'f',{'*'          ,'TIME_PSI'},'m'     ,'z position of magnetic axes, sorted by absolute flux, 0 on LCFS (*,t) [m]'                                                              )
 put('PSI_AXES'    ,fsd.psimag    ,'f',{'*'          ,'TIME_PSI'},'Wb'    ,'Poloidal flux on magnetic axes, sorted by absolute flux, 0 on LCFS, for 2*pi radian (*,t)[Wb]'                                          )
 put('DR2_PSI_AXES',fsd.dr2mag    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to r of poloidal flux on magnetic axes, sorted by largest absolute flux, for 2*pi radian (*,t)[Wb/m^2]'              )
 put('DZ2_PSI_AXES',fsd.dz2mag    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to z of poloidal flux on magnetic axes, sorted by largest absolute flux, for 2*pi radian (*,t)[Wb/m^2]'              )
 put('DRZ_PSI_AXES',fsd.drzmag    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to r and z of poloidal flux on magnetic axes, sorted by largest absolute flux, for 2*pi radian (*,t)[Wb/m^2]'        )
 put('R_XPTS'      ,fsd.rxpt      ,'f',{'*'          ,'TIME_PSI'},'m'     ,'r position of X points, sorted by absolute flux difference from LCFS (*,t) [m]'                                                         )
 put('Z_XPTS'      ,fsd.zxpt      ,'f',{'*'          ,'TIME_PSI'},'m'     ,'z position of X points, sorted by absolute flux difference from LCFS (*,t) [m]'                                                         )
 put('PSI_XPTS'    ,fsd.psixpt    ,'f',{'*'          ,'TIME_PSI'},'Wb'    ,'Poloidal flux on X points, sorted by absolute difference from LCFS, 0 on LCFS, for 2*pi radian (*,t)[Wb]'                               )
 put('DR2_PSI_XPTS',fsd.dr2xpt    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to r of poloidal flux on X points, sorted by absolute flux difference from LCFS, for 2*pi radian (*,t)[Wb/m^2]'      )
 put('DZ2_PSI_XPTS',fsd.dz2xpt    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to z of poloidal flux on magnetic axes, sorted by absolute flux difference from LCFS, for 2*pi radian (*,t)[Wb/m^2]' )
 put('DRZ_PSI_XPTS',fsd.drzxpt    ,'f',{'*'          ,'TIME_PSI'},'Wb/m^2','2nd derivative wrt to r and z of poloidal flux on X points, sorted by absolute flux difference from LCFS, for 2*pi radian (*,t)[Wb/m^2]')
 put('A_RHO'       ,fsd.x         ,'f',{'RHO','THETA','TIME_PSI'},'m'     ,'Distance to magnetic axis of flux surfaces (rho,theta,t) [m]'         )
 put('R_RHO'       ,X.R_RHO       ,'x',{'RHO','THETA','TIME_PSI'},'m'     ,'r position of flux surfaces (rho,theta,t) [m]'                        )
 put('Z_RHO'       ,X.Z_RHO       ,'x',{'RHO','THETA','TIME_PSI'},'m'     ,'z position of flux surfaces (rho,theta,t) [m]'                        )
 put('A_EDGE'      ,X.A_EDGE      ,'f',{      'THETA','TIME_PSI'},'m'     ,'Distance to magnetic axis of the LCFS (theta,t) [m]'                  )
 put('R_EDGE'      ,X.R_EDGE      ,'f',{      'THETA','TIME_PSI'},'m'     ,'r position of the LCFS (theta,t) [m]'                                 )
 put('Z_EDGE'      ,X.Z_EDGE      ,'f',{      'THETA','TIME_PSI'},'m'     ,'z position of the LCFS (theta,t) [m]'                                 )
 put('R_IN_MID'    ,X.R_IN_MID    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the flux surfaces on inner midplane (rho,t) [m]'        )
 put('R_OUT_MID'   ,X.R_OUT_MID   ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the flux surfaces on outer midplane (rho,t) [m]'        )
 put('A_MINOR_MID' ,X.A_MINOR_MID ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'Minor radius on midplane (rho,t) [m]'                                 )
 put('R_GEOM_MID'  ,X.R_GEOM_MID  ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'Geometrical major radius on midplane (rho,t) [m]'                     )
 put('R_IN'        ,fsg.rin.x     ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the inner most point of flux surfaces (rho,t) [m]'      )
 put('Z_IN'        ,fsg.zin.x     ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'z position of the inner most point of flux surfaces (rho,t) [m]'      )
 put('R_OUT'       ,fsg.rout.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the outer most point of flux surfaces (rho,t) [m]'      )
 put('Z_OUT'       ,fsg.zout.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'z position of the outer most point of flux surfaces (rho,t) [m]'      )
 put('R_TOP'       ,fsg.rtop.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the top most point of flux surfaces (rho,t) [m]'        )
 put('Z_TOP'       ,fsg.ztop.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'z position of the top most point of flux surfaces (rho,t) [m]'        )
 put('R_BOT'       ,fsg.rbot.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'r position of the bottom most point of flux surfaces (rho,t) [m]'     )
 put('Z_BOT'       ,fsg.zbot.x    ,'f',{'RHO'        ,'TIME_PSI'},'m'     ,'z position of the bottom most point of flux surfaces (rho,t) [m]'     )
 put('KAPPA'       ,X.KAPPA       ,'f',{'RHO'        ,'TIME_PSI'},''      ,'Elongation of flux surfaces (rho,t)'                                  )
 put('KAPPA_AXIS'  ,X.KAPPA_AXIS  ,'f',{              'TIME_PSI'},''      ,'Elongation on magnetic axis (t)'                                      )
 put('KAPPA_95'    ,X.KAPPA_95    ,'f',{              'TIME_PSI'},''      ,'Elongation of flux surface with rho**2=0.95 (t)'                      )
 put('KAPPA_EDGE'  ,X.KAPPA_EDGE  ,'f',{              'TIME_PSI'},''      ,'Elongation of LCFS (t)'                                               )
 put('DELTA_TOP'   ,X.DELTA_TOP   ,'f',{'RHO'        ,'TIME_PSI'},''      ,'Top triangularity of flux surfaces (rho,t)'                           )
 put('DELTA_BOT'   ,X.DELTA_BOT   ,'f',{'RHO'        ,'TIME_PSI'},''      ,'Bottom triangularity of flux surfaces (rho,t)'                        )
 put('DELTA'       ,X.DELTA       ,'f',{'RHO'        ,'TIME_PSI'},''      ,'Top/bottom average triangularity of flux surfaces (rho,t)'            )
 put('DELTA_95_TOP',X.DELTA_95_TOP,'f',{              'TIME_PSI'},''      ,'Top triangularity of flux surface with rho**2=0.95 (t)'               )
 put('DELTA_95_BOT',X.DELTA_95_BOT,'f',{              'TIME_PSI'},''      ,'Bottom triangularity of flux surface with rho**2=0.95 (t)'            )
 put('DELTA_95'    ,X.DELTA_95    ,'f',{              'TIME_PSI'},''      ,'Top/bottom average triangularity of flux surface with rho**2=0.95 (t)')
 put('DELTA_ED_TOP',X.DELTA_ED_TOP,'f',{              'TIME_PSI'},''      ,'Top triangularity of LCFS (t)'                                        )
 put('DELTA_ED_BOT',X.DELTA_ED_BOT,'f',{              'TIME_PSI'},''      ,'Bottom triangularity of LCFS (t)'                                     )
 put('DELTA_EDGE'  ,X.DELTA_EDGE  ,'f',{              'TIME_PSI'},''      ,'Top/bottom average triangularity of LCFS (t)'                         )
 
 %% Flux surface integrals
 put('AREA'        ,fsg.area.x       ,'f',{'RHO','TIME_PSI'},'m^2'       ,'Cross section of flux surfaces (rho,t) [m^2]'                     )
 put('AREA_EDGE'   ,fsg.area.x(end,:),'f',{      'TIME_PSI'},'m^2'       ,'Cross section of LCFS (t) [m^2]'                                  )
 put('VOL'         ,fsg.vol.x        ,'f',{'RHO','TIME_PSI'},'m^3'       ,'Volume of flux surfaces (rho,t) [m^3]'                            )
 put('VOL_EDGE'    ,fsg.vol.x(end,:) ,'f',{      'TIME_PSI'},'m^3'       ,'Volume of LCFS (t) [m^3]'                                         )
 put('SURF'        ,fsg.surf.x       ,'f',{'RHO','TIME_PSI'},'m^2'       ,'Surface of flux surfaces (rho,t) [m^2]'                           )
 put('DRHO_AREA'   ,fsg.darea.x      ,'f',{'RHO','TIME_PSI'},'m^2'       ,'Differential cross section of flux surfaces by drho (rho,t) [m^2]')
 put('DRHO_VOL'    ,fsg.dvol.x       ,'f',{'RHO','TIME_PSI'},'m^3'       ,'Differential volume of flux surfaces by drho (rho,t) [m^3]'       )
 put('GPSI0_R1_FSA',fsg.ngpsi0r1.x   ,'f',{'RHO','TIME_PSI'},'m^-1'      ,'<1/R> flux surface average (rho,t) [m^-1]'                        )
 put('GPSI0_R2_FSA',fsg.ngpsi0r2.x   ,'f',{'RHO','TIME_PSI'},'m^-2'      ,'<1/R^2> flux surface average (rho,t) [m^-2]'                      )
 put('GPSI2_R0_FSA',fsg.ngpsi2r0.x   ,'f',{'RHO','TIME_PSI'},'Wb^2/m^2'  ,'<|grad psi|^2> flux surface average (rho,t) [Wb^2/m^2]'           )
 put('GPSI2_R2_FSA',fsg.ngpsi2r2.x   ,'f',{'RHO','TIME_PSI'},'Wb^2/m^4'  ,'<|grad psi|^2/R^2> flux surface average (rho,t) [Wb^2/m^4]'       )
 put('Q'           ,X.Q              ,'f',{'RHO','TIME_PSI'},''          ,'Safety factor of flux surfaces (rho,t)'                           )
 put('Q_EDGE'      ,X.Q_EDGE         ,'f',{      'TIME_PSI'},''          ,'Safety factor of LCFS (t)'                                        )
 put('Q_95'        ,X.Q_95           ,'f',{      'TIME_PSI'},''          ,'Safety factor of flux surface with rho**2=0.95 (t)'               )
 put('Q_AXIS'      ,X.Q_AXIS         ,'f',{      'TIME_PSI'},''          ,'Safety factor of axis (t)'                                        )
 put('TOR_FLUX_TOT',X.TOR_FLUX_TOT   ,'f',{'RHO','TIME_PSI'},'Wb'        ,'Total toroidal flux in flux surfaces (rho,t) [Wb]'                )

 %% Unused nodes
 % put('DPSI_AREA'   ,,'f',{'RHO','TIME_PSI'},'m^2/Wb'       ,'Differential cross section of flux surfaces by dpsi (rho,t) [m^2/Wb]')
 % put('DPSI_VOL'    ,,'f',{'RHO','TIME_PSI'},'m^3/Wb'       ,'Differential volume of flux surfaces by dpsi (rho,t) [m^3/Wb]'       )
 
end
