% Script to build tables for Lz and Zave (cooling rates and average charge state)
source = 'web';
year   = [];

% If RAPTOR_path has been called before (when testing) skip the call
strpath = which('writestruct');
if isempty(strpath)
 run(fullfile(pwd,'../../..','RAPTOR_path.m')); % add RAPTOR path
end

% List of elements defined in chargemass
[A,Z,name] = chargemass;

missing_element = {};

% Start to compute D0/D+ equilibrium 
[Te1,Lz1,Zave1,Z2ave1,Zl,rep_H1,year_out,source_out,quality] = compute_lz('H',year,source);            
[Te2,Lz2,Zave2,Z2ave2,Zl,rep_H2,year_out,source_out,quality] = compute_lz('T',year,source);    
if all(Te1 == Te2)
    Te_D = Te1;
    rep_D = (rep_H1 + rep_H2) ./ 2;
nH0ne = rep_D(:,1);
else
     error('????????????????')
end
       
% Initialisation
for k=1:length(A)
    try 
      if name{k} == 'H'
	  name_loc = 'H';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('H',year,source);    
      elseif name{k} == 'D'
	  name_loc = 'H';
	  [Te1,Lz1,Zave1,Z2ave1,Zl,rep,year_out,source_out,quality] = compute_lz('H',year,source);            
	  [Te2,Lz2,Zave2,Z2ave2,Zl,rep,year_out,source_out,quality] = compute_lz('T',year,source);    
	  if all(Te1 == Te2)
	      Te = Te1;
	      Lz = 0.5 .* (Lz1 + Lz2);
	      Zave = 0.5 .* (Zave1 + Zave2);
	      Z2ave = 0.5 .* (Z2ave1 + Z2ave2);
	  else
		error('????????????????')
	  end
      elseif name{k} == 'T'
	  name_loc = 'H';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('T',year,source);    
      elseif strcmp(name{k},'Ar')
	  name_loc = 'Ar';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('Ar',[85,89],source);      
      elseif strcmp(name{k},'Ni')
	  name_loc = 'Ni';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('Ni',[85,89],source);      
      elseif strcmp(name{k},'He3')
	  name_loc = 'He';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('He',year,source,nH0ne);      
      elseif strcmp(name{k},'He4')
	  name_loc = 'He';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('He',year,source,nH0ne);      
      elseif strcmp(name{k},'He')
	  name_loc = 'He';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('He',year,source,nH0ne);      
      elseif strcmp(name{k},'Si')
	  name_loc = 'Si';
	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('Si',[96,96],source,nH0ne);      
      elseif strcmp(name{k},'W')
  	  name_loc = 'W';
  	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('W',50,source,nH0ne);    
      elseif strcmp(name{k},'Fe')
  	  name_loc = 'Fe';
  	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('Fe',[85,89],source,nH0ne);      
      elseif strcmp(name{k},'Cl')
  	  name_loc = 'Cl';
  	  [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz('Cl',89,source,nH0ne);    
      else
	  name_loc = name{k};
	  fprintf('Element: %s\n',name_loc);
          [Te,Lz,Zave,Z2ave,Zl,rep,year_out,source_out,quality] = compute_lz(name{k},year,source,nH0ne);
	  if isempty(Lz)
              disp('Problem with ADAS web data')
	  end
      end
	  
      if ~isempty(Lz)
	tabmat.(name{k}).A = A(k);
	tabmat.(name{k}).Z = Z(k);
	tabmat.(name{k}).data = cat(2,Te(:) ./ 1e3,Lz(:) ./ 1e6,Zave(:),Z2ave(:));
	tabmat.(name{k}).Zl = Zl;
	tabmat.(name{k}).rep = rep;
	tabmat.(name{k}).year = year_out;
	tabmat.(name{k}).source = source_out;
	tabmat.(name{k}).quality = quality(:);
	tabmat.(name{k}).Te = Te(:) ./ 1e3;
      else
	fprintf('No data for %s\n',name{k});
	missing_element{end+1} = name{k};
      end
      
    catch
	fprintf('Error for element: %s ->\n',name{k});
	disp(lasterr);
    end
    fclose('all');
end

% For compatibility
tabmat.He4 = tabmat.He;
  
save Lz_zave_HDT0.mat tabmat;

writestruct(tabmat,'Lz_zave_HDT0');
