function [ec,nb,lh,ic] = get_HCD_from_IDS(core_sources,ec,nb,ic,lh)


% ----> Read H&CD profiles from input core_sources IDS

nsource = numel(core_sources.source);
nrho    = numel(core_sources.source{1}.profiles_1d{1}.grid.rho_tor_norm);
n_ec_antenna = 0;
n_ic_antenna = 0;
n_nb_antenna = 0;
n_lh_antenna = 0;
ec.active = 0;
ic.active = 0;
nb.active = 0;
lh.active = 0;
ec_cs_idx = [];
ic_cs_idx = [];
nb_cs_idx = [];
lh_cs_idx = [];

core_source_power     = zeros(nsource,nrho);
core_source_power_el  = zeros(nsource,nrho);
core_source_power_ion = zeros(nsource,nrho);
for isource=1:nsource
  if numel(core_sources.source{isource}.profiles_1d{1}.electrons.energy) > 0
    core_source_power_el(isource,:) = core_sources.source{isource}.profiles_1d{1}.electrons.energy;
  end
  if numel(core_sources.source{isource}.profiles_1d{1}.total_ion_energy) > 0
    core_source_power_ion(isource,:) = core_sources.source{isource}.profiles_1d{1}.total_ion_energy;
  end
  core_source_power(isource,:) = core_source_power_el(isource,:) + core_source_power_ion(isource,:);
end
  
for isource=1:nsource
  if strcmp(core_sources.source{isource}.identifier.name,'ec')
    n_ec_antenna = n_ec_antenna + 1;
    ec.active = 1;
    ec_cs_idx = [ec_cs_idx isource]; % Indices in core_sources IDS
  end
  if strcmp(core_sources.source{isource}.identifier.name,'ic')
    n_ic_antenna = n_ic_antenna + 1;
    ic.active = 1;
    ic_cs_idx = [ic_cs_idx isource]; % Indices in core_sources IDS
  end
  if strcmp(core_sources.source{isource}.identifier.name,'nbi')
    n_nb_antenna = n_nb_antenna + 1;
    nb.active = 1;
    nb_cs_idx = [nb_cs_idx isource]; % Indices in core_sources IDS
  end
  if strcmp(core_sources.source{isource}.identifier.name,'lh')
    n_lh_antenna = n_lh_antenna + 1;
    lh.active = 1;
    lh_cs_idx = [lh_cs_idx isource]; % Indices in core_sources IDS
  end
end
core_source_power_el_tot = zeros(nsource,1);
core_source_power_ion_tot = zeros(nsource,1);
for isource=1:nsource
  core_source_power_el_tot(isource)  = trapz(core_sources.source{isource}.profiles_1d{1}.grid.volume,core_source_power_el(isource,:));
  core_source_power_ion_tot(isource) = trapz(core_sources.source{isource}.profiles_1d{1}.grid.volume,core_source_power_ion(isource,:));
  core_source_power_total(isource)   = core_source_power_el_tot(isource) + core_source_power_ion_tot(isource);
end

%% ECH
ec.cd_eff = ones(1,n_ec_antenna);
idx_ec_min = 2;
idx_ec_max = n_ec_antenna+1;
ec.uindices = [idx_ec_min:idx_ec_max];
ec.rdep = [];
ec.wdep = [];
for ii=1:numel(ec_cs_idx)
  isource = ec_cs_idx(ii);
  prof = core_sources.source{isource}.profiles_1d{1}.electrons.energy;
  rhotornorm = core_sources.source{isource}.profiles_1d{1}.grid.rho_tor_norm;
  if sum(prof) > 0
    [ec.rdep(ii),ec.wdep(ii)] = fit_gaussian(rhotornorm,prof);
    ec.rdep(find(isnan(ec.rdep)))=0;
    ec.wdep(find(isnan(ec.wdep)))=0.1;
  else
    ec.rdep(ii) = 0;
    ec.wdep(ii) = 1;
  end
end

%% ICRH
idx_ic_min = n_ec_antenna+2;
idx_ic_max = n_ec_antenna+2+n_ic_antenna-1;
ic.uindices = [idx_ic_min:idx_ic_max];
for ii=1:numel(ic_cs_idx)
  isource = ic_cs_idx(ii);
  prof = core_sources.source{isource}.profiles_1d{1}.electrons.energy;
  rhotornorm = core_sources.source{isource}.profiles_1d{1}.grid.rho_tor_norm;
  ic.rho  = rhotornorm;
  if core_source_power_total(isource) > 0
    ic.ich  = prof./core_source_power_total(isource);
  else
    ic.nbh  = zeros(size(prof));
  end
  ic.iccd = zeros(size(core_sources.source{isource}.profiles_1d{1}.j_parallel));
end
ic.eta_ic0 = 0.;
if mean(core_source_power_el_tot(ic_cs_idx)./core_source_power_ion_tot(ic_cs_idx))>0
  ic.frac_to_electrons = mean(core_source_power_el_tot(ic_cs_idx)./core_source_power_ion_tot(ic_cs_idx));
end
  
%% NBI
idx_nb_min = n_ec_antenna+2+n_ic_antenna;
idx_nb_max = n_ec_antenna+2+n_ic_antenna+n_nb_antenna-1;
nb.uindices = [idx_nb_min:idx_nb_max];
for ii=1:numel(nb_cs_idx)
  isource = nb_cs_idx(ii);
  prof = core_sources.source{isource}.profiles_1d{1}.electrons.energy;
  rhotornorm = core_sources.source{isource}.profiles_1d{1}.grid.rho_tor_norm;
  nb.rho  = rhotornorm;
  if core_source_power_total(isource) > 0
    nb.nbh  = prof./core_source_power_total(isource);
  else
    nb.nbh  = zeros(size(prof));
  end
  nb.nbcd = zeros(size(core_sources.source{isource}.profiles_1d{1}.j_parallel));
end
if  mean(core_source_power_el_tot(nb_cs_idx)./core_source_power_ion_tot(nb_cs_idx))>0
  nb.frac_to_electrons = mean(core_source_power_el_tot(nb_cs_idx)./core_source_power_ion_tot(nb_cs_idx));
end
  
%% LH
idx_lh_min = n_ec_antenna+2+n_ic_antenna+n_nb_antenna;
idx_lh_max = n_ec_antenna+2+n_ic_antenna+n_nb_antenna+n_lh_antenna-1;
lh.uindices = [idx_lh_min:idx_lh_max];

for ii=1:numel(lh_cs_idx) 
  isource = lh_cs_idx(ii);
  prof = core_sources.source{isource}.profiles_1d{1}.electrons.energy;
  rhotornorm = core_sources.source{isource}.profiles_1d{1}.grid.rho_tor_norm;
  lh.rho  = rhotornorm;
  if core_source_power_total(isource) > 0
    lh.lhh  = prof./core_source_power_total(isource);
  else
    lh.nbh  = zeros(size(prof));
  end
  lh.lhcd = zeros(size(core_sources.source{isource}.profiles_1d{1}.j_parallel));
end
if mean(core_source_power_el_tot(lh_cs_idx)./core_source_power_ion_tot(lh_cs_idx))>0
  lh.frac_to_electrons = mean(core_source_power_el_tot(lh_cs_idx)./core_source_power_ion_tot(lh_cs_idx));
end

return
