function [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(shot,times_in,varargin);
%
% [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(shot,times_in,varargin);
% 
% varargin{1}: torbeam_in (from a previous run, so avoids to reload all the data)
%              leave then times_in empty to rerun at all times in torbeam_in.inbeam.t_some_beam_on or select times points from that array
% varargin{2}: source for gdat profiles/equilibrium related: [] (default sources), 'IDA' (IDA/IDE sources), {'TRA','expname'} for TRANSP
% varargin{3}: source for torbeam angles: 0 (default) from conversion from ECS/ECN with aug_setval2tp_mex; 1 from TBP/git shotfile
% varargin{4}: nb of radial points for outputs in shotfile structure (default 401)
% varargin{5}: if >0: zeff default if no data in gdat(shot,'zeff') (2. is default)
%              if <0: zeff value (abs(varargin{5}) to impose for all times
%
% Examples:
%       [torbeam_out,torbeam_in] = torbeam_prepare_inputs_and_run(30594,linspace(0,10,101));
%       [torbeam_out2] = torbeam_prepare_inputs_and_run(30594,torbeam_in.inbeam.t_some_beam_on(1:30:end),torbeam_in); % to rerun a few cases
%

torbeam_in = [];
torbeam_out = [];

if nargin==0 || isempty(shot)
  error('requires shot number')
  return
end
times_for_torbeam_eff = [];
if nargin>=2 && ~isempty(times_in)
  times_for_torbeam_eff = times_in;
end

doload = 1;
if nargin>=3 && ~isempty(varargin{1})
  torbeam_in = varargin{1};
  doload = 0;
end
torbeam_in.prepare_run_options.doload = doload;

source_gdat_profiles = [];
source_gdat_profiles_exp_name = [];
if nargin>=4 && ~isempty(varargin{2})
  if ~iscell( varargin{2})
    source_gdat_profiles = varargin{2};
  else
    if length(varargin{2}) >= 2
      source_gdat_profiles = varargin{2}{1};
      source_gdat_profiles_exp_name = varargin{2}{2};
    else
      error('expects cell(2) if for source and exp_name if a cell is provided for profiles');
    end
  end
end
torbeam_in.prepare_run_options.source_gdat_profiles = source_gdat_profiles;
torbeam_in.prepare_run_options.source_gdat_profiles_exp_name = source_gdat_profiles_exp_name;

tbangles_source = 0;
if nargin>=5 && ~isempty(varargin{3})
  tbangles_source = varargin{3};
end
torbeam_in.prepare_run_options.tbangles_source = tbangles_source;

nb_points_shotfile = 401;
if nargin>=6 && ~isempty(varargin{4})
  nb_points_shotfile = varargin{4};
end
torbeam_in.prepare_run_options.nb_points_shotfile = nb_points_shotfile;

zeff_default = 2.0;
if nargin>=7 && ~isempty(varargin{5})
  zeff_default = varargin{5};
end
torbeam_in.prepare_run_options.zeff_default = zeff_default;

if doload
  tic
  if isempty(source_gdat_profiles_exp_name)
    torbeam_in.nete = gdat(shot,'nete_rho','machine','aug','source',source_gdat_profiles);
    torbeam_in.zeff = gdat(shot,'zeff','machine','aug');
  else
    torbeam_in.nete=gdat(shot,'nete_rho','machine','aug','source',source_gdat_profiles,'source_exp_name',source_gdat_profiles_exp_name);
    torbeam_in.zeff = gdat(shot,'zeff','machine','aug','source_exp_name',source_gdat_profiles_exp_name);
  end
  % qrho_and_grid=gdat(shot,'q_rho','machine','aug');
else
  if ~isfield(torbeam_in,'nete') || isempty(torbeam_in.nete)
    error('expects nete as field of torbeam_in (varargin{1}) from gdat(shot,''nete_rho'')')
  end
  if ~isfield(torbeam_in,'zeff') || isempty(torbeam_in.zeff)
    error('expects zeff as field of torbeam_in (varargin{1}) from gdat(shot,''zeff'')')
  end
end

% base all times on Te data (most important? although ne as well), note take Thomson by default to avoid cut-off
if doload && (~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff)),
  times_for_torbeam_eff = torbeam_in.nete.fit.t(torbeam_in.nete.fit.t>0);
end
% load beam releated data including theta,phi in torbeam reference
if doload, 
  [torbeam_in.inbeam] = read_ech_py(shot,times_for_torbeam_eff,[],tbangles_source);
  times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
else
   if ~isfield(torbeam_in,'inbeam') || isempty(torbeam_in.inbeam)
    error('expects inbeam as field of torbeam_in (varargin{1}) from read_ech_py(shot,times_for_torbeam)')
  end
end
if isempty(times_for_torbeam_eff)
  disp(char(10))
  warning(['no times with beams on, return' char(10)])
  return
end
% if torbeam_in provided as input (with its own time input sequence) and times_in provided as well: need to adapt times_in to match available inputs
% (otherwise would take closest eqdsk)
if doload==0
  if ~isempty(times_for_torbeam_eff)
    times_for_torbeam_eff = intersect(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
  else
    times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
  end
end
if doload,
  torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0,'source',source_gdat_profiles,'doplot',0);
  torbeam_in.dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
  unix(['mkdir -p ' torbeam_in.dir_torbeam]);
else
  if ~isfield(torbeam_in,'eqd_all') || isempty(torbeam_in.eqd_all)
    error('expects eqd_all as field of torbeam_in (varargin{1}) from gdat(shot,''eqdsk'',''time'',torbeam_in.inbeam.t_some_beam_on,''write'',0,''doplot'',0);')
  end
end
if length(torbeam_in.eqd_all.eqdsk) < 1
  error('expects non empty torbeam_in.eqd_all.eqdsk')
end
r0exp = torbeam_in.eqd_all.eqdsk{1}.r0;

if doload; disp(['time for reading data: ' num2str(toc)]); end

if doload,
  nb_points_prof = 46;
  rhofit = linspace(0,1,nb_points_prof)';
  
  tb_par_fields=fieldnames(torbeam_in.inbeam.tb_par);

  itime_te = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
  itime_ne = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
  it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
  cocos_out = 17;
  cocos_in = 17;
  if strcmp(upper(source_gdat_profiles),'TRA')
    cocos_in = 17;
  end
  tic
  for it=1:length(torbeam_in.inbeam.t_some_beam_on)
    % Te
    it_te = itime_te(it);
    torbeam_in.fname{it}.te = fullfile(torbeam_in.dir_torbeam,['Te' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
    fid=fopen(torbeam_in.fname{it}.te,'w'); % te.dat file etc per time
    fprintf(fid,'%d\n',nb_points_prof);
    % tefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_te) , torbeam_in.nete.fit.te.data(:,it_te),rhofit,-1,[1 2],[0 0]);
    tefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_te) , torbeam_in.nete.fit.te.data(:,it_te),rhofit,-0.1,[1 0],[0 0]);
    aa=[rhofit(end:-1:1) , abs(tefit(end:-1:1))/1e3]; % in keV
    fprintf(fid,'%.7f %.7f\n',aa');
    fclose(fid);
    % ne
    it_ne = itime_ne(it);
    torbeam_in.fname{it}.ne = fullfile(torbeam_in.dir_torbeam,['ne' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
    fid=fopen(torbeam_in.fname{it}.ne,'w'); % ne.dat file etc per time
    fprintf(fid,'%d\n',nb_points_prof);
    nefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_ne) , torbeam_in.nete.fit.ne.data(:,it_ne),rhofit,-0.1,[1 0],[0 0]);
    aa=[rhofit(end:-1:1) , abs(nefit(end:-1:1))/1e19]; % in 1e19
    fprintf(fid,'%.6f %.6f\n',aa');
    fclose(fid);
    % topfile
    % eqd_with_b = read_eqdsk(torbeam_in.eqd_all.eqdsk{it},cocos_in,1);
    eqd_with_b2=eqdsk_cocos_transform(torbeam_in.eqd_all.eqdsk{it},[cocos_in cocos_out]);
    eqd_with_b1=eqdsk_transform(eqd_with_b2,65,129,-0.1);
    eqd_with_b = read_eqdsk(eqd_with_b1,cocos_out,1); % to compute BR, etc
    torbeam_in.fname{it}.topfile = fullfile(torbeam_in.dir_torbeam,['topfile' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it))]);
    fid=fopen(torbeam_in.fname{it}.topfile,'w');
    fprintf(fid,'%s %d %.6f\n',' Number of radial and vertical grid points in ',shot,torbeam_in.inbeam.t_some_beam_on(it));
    fprintf(fid,'%d %d\n',length(eqd_with_b.rmesh),length(eqd_with_b.zmesh));
    fprintf(fid,'%s\n','Inside and Outside radius and psi_sep');
    psi_shift = 0.; % or eqd_with_b.psiedge for 0 value at edge
    fprintf(fid,'%.6f %.6f %.6f\n',min(eqd_with_b.rmesh),max(eqd_with_b.rmesh),eqd_with_b.psiedge-psi_shift);
    fprintf(fid,'%s\n','Radial grid coordinates');
    nb_rmesh = length(eqd_with_b.rmesh);
    nb_zmesh = length(eqd_with_b.zmesh);
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.rmesh(1:nb_rmesh));
    if mod(nb_rmesh,8) ~= 0; fprintf(fid,'\n'); end
    fprintf(fid,'%s\n','Vertical grid coordinates');
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.zmesh(1:nb_zmesh));
    if mod(nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
    fprintf(fid,'%s\n','B_r values');
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.BR(1:nb_rmesh,1:nb_zmesh));
    if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
    fprintf(fid,'%s\n','B_t values');
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.Bphi(1:nb_rmesh,1:nb_zmesh));
    if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
    fprintf(fid,'%s\n','B_z values');
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.BZ(1:nb_rmesh,1:nb_zmesh));
    if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
% $$$   fprintf(fid,'%s\n','psi norm values');
% $$$   fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',(eqd_with_b.psi(1:nb_rmesh,1:nb_zmesh)-eqd_with_b.psiaxis) ...
% $$$           ./ (eqd_with_b.psiedge-eqd_with_b.psiaxis));
    fprintf(fid,'%s\n',['psi cocos_out=' num2str(cocos_out)]);
    fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.psi(1:nb_rmesh,1:nb_zmesh));
    if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
    fclose(fid);
    
    if isempty(torbeam_in.zeff.data) || torbeam_in.prepare_run_options.zeff_default < 0.
      zeff_eff = abs(torbeam_in.prepare_run_options.zeff_default);
    else
      zeff_eff = interpos(21,torbeam_in.zeff.t,torbeam_in.zeff.data,torbeam_in.inbeam.t_some_beam_on(it));
    end
    torbeam_in.inbeam.zeff_t(it_inb(it)) = zeff_eff;
    % extract beam info at given time
    for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it)))',
      torbeam_in.fname{it}.inbeam{ibeam} = fullfile(torbeam_in.dir_torbeam,['inbeam' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) 'g' num2str(ibeam) '.dat']);
      [aa,bb]=unix(['cp /afs/ipp-garching.mpg.de/home/o/osauter/TORBEAM/inbeam_default.dat ' torbeam_in.fname{it}.inbeam{ibeam}]);
      fid=fopen(torbeam_in.fname{it}.inbeam{ibeam},'a');
      fprintf(fid,' xf = %f,	     ! wave frequency om=2*pi*xf\n',torbeam_in.inbeam.freq(ibeam));
      fprintf(fid,' nmod = %d,          ! mode selection: O-mode (1), X-mode (-1)\n',torbeam_in.inbeam.rfmod(ibeam));
      fprintf(fid,' xtordeg = %f,	     ! geom. optics injection angle\n',torbeam_in.inbeam.phi_t(ibeam,it_inb(it)));
      fprintf(fid,' xpoldeg = %f,	     ! geom. optics injection angle\n',torbeam_in.inbeam.theta_t(ibeam,it_inb(it)));
      % may be position and curvature, but get defaults from tb_gui. Check with tb_gui how several launchers are given
      for i=1:length(tb_par_fields)
        fprintf(fid,[' ' tb_par_fields{i} ' = %f,\n'],torbeam_in.inbeam.tb_par.(tb_par_fields{i})(ibeam));
      end
      fprintf(fid,' xpw0 = %f,\n',torbeam_in.inbeam.power_t(ibeam,it_inb(it))/1.e6);
      fprintf(fid,' xzeff = %f,\n',torbeam_in.inbeam.zeff_t(it_inb(it)));
      fprintf(fid,'/\n');
      fclose(fid);
    end
  end
  disp(['time for creating all input files: ' num2str(toc)])
end
  
% run TORBEAM
current_dir = pwd;
eval(['cd ' torbeam_in.dir_torbeam])
%                   namelist,  ol proj 3 rays(R,z)*100, not used now,  tor proj 3 rays(R,Y),  vol()   dP/dV(rhopol) j(rhopol)
%tb_out_files_to_cp = {'echo.dat', 't1_LIB.dat', 't2_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat', 'cntpr.dat'};
tb_out_files_to_cp = {'t1_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat'};
tb_out_files_to_cp_varname = strrep(tb_out_files_to_cp,'.dat','');
torbeam_out.run_torbeam = false .* ones(size(torbeam_in.inbeam.power_t,1),length(times_for_torbeam_eff));
it_inb_beam_on = iround_os(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
torbeam_out.t = torbeam_in.inbeam.t_some_beam_on(it_inb_beam_on);
torbeam_out.shotfile.rhopol = linspace(0,1,nb_points_shotfile);
torbeam_out.shot = shot;
torbeam_out.shotfile.shot = shot;
torbeam_out.shotfile.t = torbeam_out.t;
tic
icount = 0;
ld_library_path=getenv('LD_LIBRARY_PATH'); % need /afs/@cell/common/soft/intel/ics2013/14.0/compiler/lib/intel64 (ok on toks but not toki now)
ij_intel = findstr(ld_library_path,'intel');
if isempty(ij_intel)
  warning('no intel in LD_LIBRARY_PATH, torbeam may not run properly, add manually /afs/@cell/common/soft/intel/ics2013/14.0/compiler/lib/intel64');
  ld_library_path = [ld_library_path ':/afs/@cell/common/soft/intel/ics2013/14.0/compiler/lib/intel64']; % $INTEL_HOME/compiler/lib/intel64
end
for it_eff=1:length(times_for_torbeam_eff)
  it = it_inb_beam_on(it_eff);
  unix(['cp ' torbeam_in.fname{it}.te ' Te.dat']);
  unix(['cp ' torbeam_in.fname{it}.ne ' ne.dat']);
  unix(['cp ' torbeam_in.fname{it}.topfile ' topfile']);
  disp(['   t=' num2str(torbeam_in.inbeam.t(it_inb(it_eff))) '   beam=' num2str(find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))','%d')]);
  for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))',
    unix(['cp ' torbeam_in.fname{it}.inbeam{ibeam} ' inbeam.dat']);
    % [a,b]=unix(['setenv LD_LIBRARY_PATH ' ld_library_path '; /afs/ipp-garching.mpg.de/home/o/osauter/TORBEAM/lib-OUT/a.out']);
    [a,b]=unix('/afs/ipp-garching.mpg.de/home/o/osauter/TORBEAM/lib-OUT/a.out');
    icount = icount + 1;
    imw = findstr('(MW)',b);
    if ~isempty(imw)
      torbeam_out.run_torbeam(ibeam,it_eff) = true;
      torbeam_out.pabs(ibeam,it_eff) = 1e6 * sscanf(b(imw+4:end),'%f',1);
      ima = findstr('(MA)',b);
      torbeam_out.icd(ibeam,it_eff) = 1e6 * sscanf(b(ima+4:end),'%f',1);
      imaxat = findstr('max. at',b);
      torbeam_out.pabs_rhopolmax(ibeam,it_eff) = sscanf(b(imaxat(1)+7:end),'%f',1);
      torbeam_out.icd_rhopolmax(ibeam,it_eff) = sscanf(b(imaxat(2)+7:end),'%f',1);
      for j=1:length(tb_out_files_to_cp)
        torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it_eff} = load(fullfile(torbeam_in.dir_torbeam,[tb_out_files_to_cp{j}]));
      end
      nnn = load('cntpr.dat');
      torbeam_out.cntpr.tbegabs(ibeam,it_eff) = nnn(2);
      torbeam_out.cntpr.dtabs(ibeam,it_eff) = nnn(1);
      torbeam_out.pdens_jcd_rho(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,1);
      torbeam_out.pdens(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,2) * 1e6;
      torbeam_out.jcd(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,3) * 1e6;
      % structure for shot file, add rho=0 point, use a generic nb_points_shotfile
      x2 = [0.; torbeam_out.pdens_jcd_rho(:,ibeam,it_eff)];
      pdens2 = [torbeam_out.pdens(1,ibeam,it_eff); torbeam_out.pdens(:,ibeam,it_eff)];
      jcd2 = [torbeam_out.jcd(1,ibeam,it_eff); torbeam_out.jcd(:,ibeam,it_eff)];
      sig2 = ones(size(x2)); sig2(1)=1e3;
      sigvol=ones(size(torbeam_out.volumes{ibeam,it_eff}(1:end,1))); sigvol(2) = 1e3; %sigvol(2:3) = 1e3;
      torbeam_out.shotfile.volume(:,it_eff) = ...
          interpos([0.;torbeam_out.volumes{ibeam,it_eff}(2:end,1)],[0.;torbeam_out.volumes{ibeam,it_eff}(2:end,2)].^.5,torbeam_out.shotfile.rhopol,-0.001,[2 0],[0 0],sigvol);
      torbeam_out.shotfile_vol_ok(it_eff) = 1;
      if min(diff(torbeam_out.shotfile.volume(:,it_eff))) <=0
        sigvol(2:3) = 1e3;
        torbeam_out.shotfile.volume(:,it_eff) = ...
            interpos([0.;torbeam_out.volumes{ibeam,it_eff}(2:end,1)],[0.;torbeam_out.volumes{ibeam,it_eff}(2:end,2)].^.5,torbeam_out.shotfile.rhopol,-0.001,[2 0],[0 0],sigvol);
        torbeam_out.shotfile_vol_ok(it_eff) = 2;
        if min(diff(torbeam_out.shotfile.volume(:,it_eff))) <=0
          disp(char(10))
          warning(['non-monotonic volume for it_eff = ' num2str(it_eff) ' (ibeam= ' num2str(ibeam) ')' char(10)])
          torbeam_out.shotfile_vol_ok(it_eff) = 0;
        end
      end
      torbeam_out.shotfile.volume(1,it_eff) = 0.;
      torbeam_out.shotfile.volume(:,it_eff) = torbeam_out.shotfile.volume(:,it_eff).^2;
      torbeam_out.shotfile.pdens(:,ibeam,it_eff) = interpos(x2,pdens2,torbeam_out.shotfile.rhopol,-0.01,[1 0],[0 0],sig2);
      torbeam_out.shotfile.jcd(:,ibeam,it_eff) = interpos(x2,jcd2,torbeam_out.shotfile.rhopol,-0.01,[1 0],[0 0],sig2);
      if torbeam_out.shotfile_vol_ok(it_eff) > 0
        [~,~,~,torbeam_out.shotfile.pabs_int(:,ibeam,it_eff)] = interpos(torbeam_out.shotfile.volume(:,it_eff),torbeam_out.shotfile.pdens(:,ibeam,it_eff));
        [~,~,~,torbeam_out.shotfile.icd_int(:,ibeam,it_eff)] = interpos(torbeam_out.shotfile.volume(:,it_eff),torbeam_out.shotfile.jcd(:,ibeam,it_eff));
        torbeam_out.shotfile.icd_int(:,ibeam,it_eff) = torbeam_out.shotfile.icd_int(:,ibeam,it_eff) ./2 ./ pi ./ r0exp / 1.01; % 1.01 on average while waiting to insert proper metric
      else
        torbeam_out.shotfile.pabs_int(:,ibeam,it_eff) = nan(size(torbeam_out.shotfile.volume(:,it_eff)));
        torbeam_out.shotfile.icd_int(:,ibeam,it_eff) = nan(size(torbeam_out.shotfile.volume(:,it_eff)));
      end
      torbeam_out.shotfile.rz_r_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1_LIB{ibeam,it_eff}(:,1:2:end-1)/100.;
      torbeam_out.shotfile.rz_z_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1_LIB{ibeam,it_eff}(:,2:2:end)/100.;
      torbeam_out.shotfile.ry_r_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1tor_LIB{ibeam,it_eff}(:,1:2:end-1)/100.;
      torbeam_out.shotfile.ry_y_ray(1:size(torbeam_out.t1_LIB{ibeam,it_eff},1),1:size(torbeam_out.t1_LIB{ibeam,it_eff},2)/2,ibeam,it_eff) = torbeam_out.t1tor_LIB{ibeam,it_eff}(:,2:2:end)/100.;
      torbeam_out.shotfile.index_abs_ray(1,ibeam,it_eff) = nnn(2) + 1;
      torbeam_out.shotfile.index_abs_ray(2,ibeam,it_eff) = nnn(2) + 1 + nnn(1);
      torbeam_out.shotfile.torbeam_was_run = torbeam_out.run_torbeam;
      torbeam_out.shotfile.pabs = torbeam_out.pabs;
      torbeam_out.shotfile.icd = torbeam_out.icd;
      torbeam_out.shotfile.pabs_rhopolmax = torbeam_out.pabs_rhopolmax;
      torbeam_out.shotfile.icd_rhopolmax = torbeam_out.icd_rhopolmax;
    else
      disp(['TORBEAM did not run properly, no MW result'])
    end
  end
end
aaa=toc;
disp(['time for running all TORBEAMs: ' num2str(aaa) ' for ' num2str(icount) ' runs so ' num2str(aaa/icount) 's per TORBEAM run'])
for ibeam=find(any(torbeam_out.run_torbeam,2))'
  figure;
  subplot(2,1,1)
  plot(squeeze(torbeam_out.pdens_jcd_rho(:,ibeam,:)),squeeze(torbeam_out.pdens(:,ibeam,:)))
  title(['AUG shot #' num2str(shot) ': for beam = ' num2str(ibeam)])
  ylabel('Pdens')
  xlabel('rhopol')
  subplot(2,1,2)
  plot(squeeze(torbeam_out.pdens_jcd_rho(:,ibeam,:)),squeeze(torbeam_out.jcd(:,ibeam,:)))
  ylabel('jcd')
  xlabel('rhopol')
end

figure
subplot(3,1,1)
plot(torbeam_out.t,torbeam_out.pabs/1e6,'-')
ylabel('Pabs [MW]')

title(['AUG shot #' num2str(shot)])

subplot(3,1,2)
plot(torbeam_out.t,torbeam_out.icd/1e6,'-')
ylabel('Icd [MA]')

subplot(3,1,3)
plot(torbeam_out.t,torbeam_out.pabs_rhopolmax,'-')
ylabel('rhopol(max Pdens)')
xlabel('t [s]')


eval(['cd ' current_dir])