%FBTWTCV  Write FBT specific data in TCV tree
% FBTWTCV(SHOT,TREE,BRANCH,L,LY) writes FBT specific data from LY in
% TCV SHOT TREE BRANCH. L contains ancillary data as obtained by FBTP. To
% be used with MEQWTCV. See also MEQWTCV,FBTP,FBTT.
%
% [ FreeBoundaryTokamak equilibrium ] Sw+ssPlasmaCenter EPF+Lausanne
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function fbtwtcv(shot,tree,branch,L,LY)

 %% 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)
 
 % MGAMS.EQUIL nodes do not have a version_num subnodes.
 
 % RZE/ZZE is the center of the plasma as defined by the boundary control points
 Rlim = [L.P.rlim1;L.P.rlia1]; % FBTPTCV converts analytic shape to discrete boundary points
 Zlim = [L.P.zlim1;L.P.zlia1]; % id.
 if L.P.idoublet
  % Temporary fix for doublets, FBTE-FORTRAN, uses upper half only
  Rlim(Zlim <= 0) = NaN;
  Zlim(Zlim <= 0) = NaN;
 end
 rze = (max(Rlim,[],1)+min(Rlim,[],1))/2;
 zze = (max(Zlim,[],1)+min(Zlim,[],1))/2;
 
 % SHVSEC is the flux from external currents at (rze, zze)
 ir = iround(L.G.rx,rze); iz = iround(L.G.zx,zze);
 ix = L.nzx*(ir-1)+iz; % Get nearest grid point
 shvsec = NaN(L.P.nruns,1);
 for irun = 1:L.P.nruns
   shvsec(irun) = L.G.Mxa(ix(irun),:)*LY.Ia(:,irun);
 end
 
 % NOTE: FBTE nodes have correct signs, FBT nodes are written for Ip>0
 if L.P.idoublet
  meqput('\FBT:RAX'        ,'pack(data(\MGAMS.EQUIL:R_AXIS)[0,*],$TRUE)','x',{'*'},'','')
  meqput('\FBT:ZAX'        ,'pack(data(\MGAMS.EQUIL:Z_AXIS)[0,*],$TRUE)','x',{'*'},'','')
  meqput('\FBT:RAX2'       ,'pack(data(\MGAMS.EQUIL:R_AXIS)[1,*],$TRUE)','x',{'*'},'','')
  meqput('\FBT:ZAX2'       ,'pack(data(\MGAMS.EQUIL:Z_AXIS)[1,*],$TRUE)','x',{'*'},'','')
 else
  meqput('\FBT:RAX'        ,'\MGAMS.EQUIL:R_AXIS'                       ,'x',{'*'},'','')
  meqput('\FBT:ZAX'        ,'\MGAMS.EQUIL:Z_AXIS'                       ,'x',{'*'},'','')
  meqput('\FBT:RAX2'       ,'0*\MGAMS.EQUIL:R_AXIS'                     ,'x',{'*'},'','')
  meqput('\FBT:ZAX2'       ,'0*\MGAMS.EQUIL:Z_AXIS'                     ,'x',{'*'},'','')
 end
 meqput('\FBT:SHVSEC'     ,L.P.iohfb*shvsec                             ,'f',{'*'},'','')
 meqput('\FBT:GCURR'      ,'\MGAMS.EQUIL:I_POL*\MGAMS:DATA.IOHFB'       ,'x',{'*'},'','')
 meqput('\FBT:BEE'        ,'\MGAMS.EQUIL:B_PROBE*\MGAMS:DATA.IOHFB'     ,'x',{'*'},'','')
 meqput('\FBT:FLUX'       ,'\MGAMS.EQUIL:PSI_LOOP*\MGAMS:DATA.IOHFB'    ,'x',{'*'},'','') % NB: FORTRAN version used Ff-FB
 meqput('\FBT:PLACAMP'    ,'\MGAMS.EQUIL:I_PL*\MGAMS:DATA.IOHFB'        ,'x',{'*'},'','')
 meqput('\FBTE:I_P'       ,'\MGAMS.EQUIL:I_PL'                          ,'x',{'*'},'','')
 meqput('\FBTE:L_I'       ,'\MGAMS.EQUIL:L_I_3*\MGAMS.EQUIL:R0/\FBT:RZE','x',{'*'},'','')
 if L.P.idoublet
  meqput('\FBTE:AREA'      ,'2*\MGAMS.EQUIL:AREA_EDGE'                   ,'x',{'*'},'','')
 else
  meqput('\FBTE:AREA'      ,'\MGAMS.EQUIL:AREA_EDGE'                     ,'x',{'*'},'','')
 end
 meqput('\FBTE:KAPPA_EDGE','\MGAMS.EQUIL:KAPPA_EDGE'                    ,'x',{'*'},'','')
 % Remaining FBT nodes
 meqput('\FBT:RZE'        ,rze                                          ,'f',{'*'},'','')
 meqput('\FBT:ZZE'        ,zze                                          ,'f',{'*'},'','')
 
end
