function [torbeam] = read_ech_py(shot,varargin);
%
% [torbeam] = read_ech_py(shot,varargin);
%
% loads angles and powers from ECS/ECN for given shot
% transforms ECS/ECN into theta,phi of Torbeam
%
% varargin{1}: time array for Torbeam relevant data output
% varargin{2}: min_power_on to decide if the beam was on at a given time (default P>=2e3)
% varargin{3}: source for TB angles, etc: 0 (default, aug_setval2tp_mex function useing setmirrors library); 1 (TBP/xxx/git shotfile)
%

if nargin==0
  error('read_ech_py requires at least shot number as input');
else
  torbeam.shot = shot;
end

ecsystems=[4,4];
nbgy = sum(ecsystems);
igy_seq=[1:ecsystems(1);ecsystems(1)+1:ecsystems(1)+ecsystems(2)]';

load_option = 0,
if nargin >= 4 && ~isempty(varargin{3})
  load_option = varargin{3};
end
torbeam.read_ech_options.load_option = load_option;

% Load parameters independent of angles, like power etc:
%
tec_min = NaN;
tec_max = NaN;
ecsystems=[4,4];
nbgy = sum(ecsystems);
igy_seq=[1:ecsystems(1);ecsystems(1)+1:ecsystems(1)+ecsystems(2)]';
for isys=1:length(ecsystems)
  isyseff = isys;
  if isys==1 && shot > 33725; isyseff = 3; end
  for igy=1:ecsystems(isys)
    ps_source = ['P_sy' num2str(isyseff) '_g' num2str(igy)];
    if isyseff == 1
      aaon=gdat(shot,{'ECS',ps_source,'','param:gy_on'},'machine','aug');
    else
      aaon=gdat(shot,{'ECS',ps_source,'','param:gy_HV_on'},'machine','aug');
    end
    if isempty(aaon.data)
      warning(['problem with detecting if gy_on for isyseff = ' num2str(isyseff) ', igy = ',num2str(igy)]);
    end
    torbeam.is_gyro_on(igy_seq(igy,isys)) = false;
    if aaon.data; torbeam.is_gyro_on(igy_seq(igy,isys)) = true; end
    if torbeam.is_gyro_on(igy_seq(igy,isys))
      if isys==1
        aa=gdat(shot,{'ECS',['PG' num2str(igy)]},'machine','aug');
      else
        aa=gdat(shot,{'ECS',['PG' num2str(igy) 'N']},'machine','aug');
      end
      torbeam.pow{igy_seq(igy,isys)}.data = aa.data;
      torbeam.pow{igy_seq(igy,isys)}.t = aa.t;
      tec_min = min(min(torbeam.pow{igy_seq(igy,isys)}.t),tec_min);
      tec_max = max(max(torbeam.pow{igy_seq(igy,isys)}.t),tec_max);
    else
      torbeam.pow{igy_seq(igy,isys)}.data = [];
      torbeam.pow{igy_seq(igy,isys)}.t = [];
    end
  end
end

if ~isfinite(tec_min) || ~isfinite(tec_max)
  warning('no good time interval')
  tec_min
  tec_max
  return
end

tec_min = max(tec_min, 0.);
if tec_max <= tec_min
  warning('no good time interval')
  tec_min
  tec_max
  return
end

if nargin>=2 && ~isempty(varargin{1})
  time_for_torbeam = varargin{1};
else
  time_for_torbeam = linspace(tec_min,tec_max,round(tec_max-tec_min)*3e2); % aim at 3ms
end
torbeam.t = time_for_torbeam;
torbeam.theta_t = NaN(nbgy,length(time_for_torbeam));
torbeam.phi_t = NaN(nbgy,length(time_for_torbeam));
torbeam.err_status_t = NaN(nbgy,length(time_for_torbeam));

isys = 1;
for igy=1:ecsystems(isys)
  igy_1_8 = igy_seq(igy,isys);
  if torbeam.is_gyro_on(igy_seq(igy,isys))
      % linear interpolation for power to avoid over/under shoot
    torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
  end
end
isys = 2;
for igy=1:ecsystems(isys)
  igy_1_8 = igy_seq(igy,isys);
  if torbeam.is_gyro_on(igy_seq(igy,isys))
    torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
  end
end

if nargin>=3 && ~isempty(varargin{2})
  min_power_on = varargin{2};
else
  min_power_on = 2e3;
end
torbeam.read_ech_options.min_power_on = min_power_on;

torbeam.beam_on = false.*ones(size(torbeam.power_t));
torbeam.beam_on(torbeam.power_t >= min_power_on) = true;
torbeam.t_some_beam_on = torbeam.t(any(torbeam.beam_on,1));

switch load_option
 case 0
  torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
  torbeam.tb_par.('xxb') =  [238., 238., 231.1, 231.1, 236.4, 236.4, 236.4 ,236.4];
  torbeam.tb_par.('xzb') =  [0.,   0.,   0.,    0.,    32.,   32.,   -32.,  -32.];
  torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 85.4,   85.4,   85.4,   85.4]; % 85.0 in GT's py, 85.4 in tb_gui saved files
  torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 85.4,   85.4,   85.4,   85.4];
  torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97,  2.97,  2.55,  2.55,  2.55,  2.55];
  torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97,  2.97,  2.55,  2.55,  2.55,  2.55];
  % from Ringberg 2018 Martin Schubert
  torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 122.7, 122.7, 139.7, 139.7]; % 85.0 in GT's py, 85.4 in tb_gui saved files
  torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 122.7, 122.7, 139.7, 139.7];
  torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97,  2.97,  3.589, 3.589, 3.196, 3.196];
  torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97,  2.97,  3.589, 3.589, 3.196, 3.196];
  
% $$$   torbeam.tr2as = { ...
% $$$       'XECECH', 'RR_n'; 'ZECECH', 'ZZ_n'; ...
% $$$       'ECB_WIDTH_HOR', 'xwyyb'; 'ECB_WIDTH_VERT', 'xwzzb'; ...
% $$$       'ECB_CURV_HOR' , 'xryyb'; 'ECB_CURV_VERT' , 'xrzzb' };

  % following python script:
  % WARNING seems alpha, beta inverted in read_ech.py script, so here use alpha for poloidal related angle and beta for tor
  % Thus use alpha_pol and beta_tor to avoid confusion (inverted with respect to .py alpha, beta script)
  for isys=1:length(ecsystems)
    isyseff = isys;
    if isys==1 && shot > 33725; isyseff = 3; end
    for igy=1:ecsystems(isys)
      ps_source = ['P_sy' num2str(isyseff) '_g' num2str(igy)];
      aa=gdat(shot,{'ECS',ps_source,'','param:GPolPos'},'machine','aug');
      torbeam.alpha_pol(igy_seq(igy,isys),1) = aa.data;
      aa=gdat(shot,{'ECS',ps_source,'','param:GTorPos'},'machine','aug');
      torbeam.beta_tor(igy_seq(igy,isys),1) = aa.data;
      if load_option == 0
        aa=gdat(shot,{'ECS',ps_source,'','param:gyr_freq'},'machine','aug');
        if aa.data > 1000
          torbeam.freq(igy_seq(igy,isys),1) = aa.data;
        else
          torbeam.freq(igy_seq(igy,isys),1) = 1.398765E+11;
        end
      else
        % assume already loaded in above switch case
      end
    end
  end

  % Compute theta,phi for TORBEAM
  torbeam.theta = NaN(nbgy,1);
  torbeam.phi = NaN(nbgy,1);
  torbeam.err_status = -ones(nbgy,1);
  torbeam.isystunit = NaN(nbgy,1);
  torbeam.datum = NaN(nbgy,1);
  for isys=1:length(ecsystems)
    for igy=1:ecsystems(isys)
      igy_1_8 = igy_seq(igy,isys);
      if torbeam.is_gyro_on(igy_seq(igy,isys))
        [torbeam.theta(igy_seq(igy,isys),1),torbeam.phi(igy_seq(igy,isys),1),torbeam.err_status(igy_seq(igy,isys),1),torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
            aug_setval2tp_mex(torbeam.alpha_pol(igy_seq(igy,isys),1),torbeam.beta_tor(igy_seq(igy,isys),1),igy_1_8,shot);
      end
    end
  end

  % get time dependent data
  if shot < 23187
    % Before ECRH2 installed
    return
  end

  isys=length(ecsystems);
  isyseff = isys;
  diagname = 'ECN';
  % 1st system angles fixed so use above values
  isys = 1;
  for igy=1:ecsystems(isys)
    igy_1_8 = igy_seq(igy,isys);
    if torbeam.is_gyro_on(igy_seq(igy,isys))
      torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_seq(igy,isys));
      torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_seq(igy,isys));
      torbeam.err_status_t(igy_1_8,:) = torbeam.err_status(igy_seq(igy,isys),1);
    end
  end
  isys = 2;
  for igy=1:ecsystems(isys)
    torbeam.alpha_polN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'POL']},'machine','aug');
    torbeam.alpha_polN{igy}.data = 0.01 .* torbeam.alpha_polN{igy}.data;
    torbeam.alpha_pol_t = interpos(-21,torbeam.alpha_polN{igy}.t,torbeam.alpha_polN{igy}.data,torbeam.t,-1e2);
% $$$   torbeam.beta_torN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'TOR']},'machine','aug');
    igy_1_8 = igy_seq(igy,isys);
    if torbeam.is_gyro_on(igy_seq(igy,isys))
      [theta_t,phi_t,err_status_t,torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
          aug_setval2tp_mex(torbeam.alpha_pol_t,torbeam.beta_tor(igy_seq(igy,isys),1).*ones(size(torbeam.alpha_pol_t)),igy_1_8,shot);
      torbeam.theta_t(igy_1_8,:) = theta_t';
      torbeam.phi_t(igy_1_8,:) = phi_t';
      torbeam.err_status_t(igy_1_8,:) = err_status_t';
    end
  end

 case 1
  torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
  tmp=gdat(shot,{'TBP','TB_par','git','param:Rlaunch'},'machine','aug');
  torbeam.tb_par.('xxb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:Zlaunch'},'machine','aug');
  torbeam.tb_par.('xzb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:wid_hor'},'machine','aug');
  torbeam.tb_par.('xwyyb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:wid_ver'},'machine','aug');
  torbeam.tb_par.('xwzzb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:curv_hor'},'machine','aug');
  torbeam.tb_par.('xryyb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:curv_ver'},'machine','aug');
  torbeam.tb_par.('xrzzb') =  reshape(tmp.data,1,length(tmp.data)) * 100.;
  tmp=gdat(shot,{'TBP','TB_par','git','param:the_ps'},'machine','aug');
  torbeam.tb_par_sf.('the_ps') =  reshape(tmp.data,1,length(tmp.data));
  tmp=gdat(shot,{'TBP','TB_par','git','param:phi_ps'},'machine','aug');
  torbeam.tb_par_sf.('phi_ps') =  reshape(tmp.data,1,length(tmp.data));
  tmp=gdat(shot,{'TBP','TB_par','git','param:freq'},'machine','aug');
  torbeam.tb_par_sf.('freq') =  reshape(tmp.data,1,length(tmp.data));
  tmp=gdat(shot,{'TBP','the_ec2','git'},'machine','aug');
  torbeam.tbp.the_ec2 = tmp;
  tmp=gdat(shot,{'TBP','phi_ec2','git'},'machine','aug');
  torbeam.tbp.phi_ec2 = tmp;
  
  % fill in relevant output tables
  torbeam.freq = torbeam.tb_par_sf.freq';
  torbeam.theta = torbeam.tb_par_sf.the_ps';
  torbeam.phi = torbeam.tb_par_sf.phi_ps';
  isys = 1;
  igy_1_8 = igy_seq(1:ecsystems(isys),isys); % sys1 has fixed angles vs time
  torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_1_8)*ones(1,size(torbeam.theta_t,2));
  torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_1_8)*ones(1,size(torbeam.phi_t,2));
  isys = 2;
  igy_1_8 = igy_seq(1:ecsystems(isys),isys); % sys1 has fixed angles vs time
  for igy=1:size(igy_seq,1)
    igy_1_8 = igy_seq(igy,isys);
    torbeam.theta_t(igy_1_8,:) = interpos(-21,torbeam.tbp.the_ec2.t,torbeam.tbp.the_ec2.data(igy,:),torbeam.t,-1e2);
    torbeam.phi_t(igy_1_8,:) = interpos(-21,torbeam.tbp.phi_ec2.t,torbeam.tbp.phi_ec2.data(igy,:),torbeam.t,-1e2);
  end
 
 otherwise
  error(['option for Torbeam parameters and angles = ' num2str(load_option) ' not known nor implemented yet. Ask O. Sauter']);
end

torbeam = orderfields(torbeam);
