function [expeq_main_filename,tout,normalization_infos]=RAPTORout2expeq(out,t,filename,varargin)
%
% function [EXPEQname,tout] = RAPTORout2expeq(out,t,filename,varargin);
%
% Writes some profile quantities to ASCII for later use with CHEASE
% out: RAPTOR_out structure
% t: time
% filename: path of ASCII file to be written
%
% varargin{1}: chease_data structure for relevant time t (which contains eqdsk)
%
% At this stage requires varargin{1}, otherwise plasma boundary in particular and R0, B0 are not defined
%
%%
debug=0;

[~,it] = min(abs(t-out.time));
tout = out.time(it);

eqdsk_in=[];
if nargin>=4 && ~isempty(varargin{1})
  chease_data = varargin{1};
  if isfield(chease_data,'eqdsk')
    eqdsk_in = chease_data.eqdsk;
  end
  R0 = eqdsk_in.r0;
  B0 = eqdsk_in.b0;
else
  return
  R0 = 0.88;
  B0 = 1.43;
end

rhotor = out.rho;
rhotor_gauss = out.rhogauss;
rhopol = out.rhopol(:,it);
q = out.q(:,it);
psi = out.psi(:,it)-out.psi(end,it);
% jdotB = out.jdotB(:,it);
% Ipar = out.Ipar(:,it);
p = out.p(:,it);
F = out.F(:,it);
g3 = out.g3(:,it);
Ip = out.Ip(:,it);
jpar = out.jpar(:,it);
kappa = out.kappa(:,it);
%delta = out.delta(:,it);
volume = out.Volume(:,it);
pprime = out.pprime(:,it);
ttprime = out.ttprime(:,it);

%%
if debug==1
  fid = fopen(filename,'w');
  fprintf(fid,'%s\n',...
          '%rhotor rhopol Ip q psi-psib p RBphi <1/R^2> jpar kappa volume ');
  fprintf(fid,...
          '%.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e %.8e \n',...
          [rhotor rhopol Ip q psi p F g3 jpar kappa volume]'); % R0 B0
  fclose(fid);
  
  fprintf('\nwrote data for t=%03.3f to %s\n',tout,filename)
  %%
  
  fid = fopen([filename '_gauss'],'w');
  fprintf(fid,'%s\n',...
          '%rhotor_gauss pprime ttprime');
  fprintf(fid,...
          '%.8e %.8e %.8e \n',...
          [rhotor_gauss pprime ttprime]');
  fclose(fid);
  
  fprintf('\nwrote data for t=%03.3f to %s\n',tout,[filename '_gauss']);
end
%%

% make expeq
mu0=4e-7.*pi;
r_phys2chease = 1./R0;
b_phys2chease = 1./B0;
p_phys2chease = mu0./B0.^2;
pprime_phys2chease = mu0.*R0.^2./B0 .*2 .*pi; % Raptor cocos=12, thus as chease but extra 2pi in psi
ttprime_phys2chease = 1./B0 .*2 .*pi; % Raptor cocos=12, thus as chease but extra 2pi in psi
Ip_phys2chease = mu0./R0./B0;
j_phys2chease = mu0.*R0./B0;
psi_phys2chease = 1./B0./R0.^2./2./pi; % Raptor cocos=12, thus as chease but extra 2pi in psi
%p_phys2chease = ;

normalization_infos.r0 = R0;
normalization_infos.r0_mean_r = eqdsk_in.syscode_formulas.R0;
normalization_infos.b0 = B0;
normalization_infos.p_phys2chease = p_phys2chease;
normalization_infos.pprime_phys2chease = pprime_phys2chease;
normalization_infos.ttprime_phys2chease = ttprime_phys2chease;
normalization_infos.ip_phys2chease = Ip_phys2chease;
normalization_infos.j_phys2chease = j_phys2chease;
normalization_infos.psi_phys2chease = psi_phys2chease;
normalization_infos.ip = Ip(end);
normalization_infos.q0_raptor = q(1);
normalization_infos.qa_raptor = q(end);
normalization_infos.q0_eqdsk = eqdsk_in.q(1);
normalization_infos.qa_eqdsk = eqdsk_in.q(end);
normalization_infos.li3_raptor = out.li3(it);
normalization_infos.betan_raptor = out.betaN(it);
normalization_infos.betap_raptor = out.betapol(it);
normalization_infos.a_minor = eqdsk_in.syscode_formulas.a_minor;
normalization_infos.kappa = eqdsk_in.syscode_formulas.kappa;
normalization_infos.delta = eqdsk_in.syscode_formulas.delta;

if ~isempty(eqdsk_in)
  expeq.zgeom=eqdsk_in.zaxis .* r_phys2chease;
  expeq.raxis=eqdsk_in.raxis .* r_phys2chease;
  expeq.pedge = p(end) .* p_phys2chease;
  expeq.epsilon = out.epsilon(end,it);
  expeq.RZ_psi = [eqdsk_in.rplas  eqdsk_in.zplas] .* r_phys2chease;
  expeq.nsttp = 1;
  expeq.nrhotype = 0; % no rhotor for nsttp=1 yet
  expeq.rho = linspace(0,1,101);
  expeq.Pprime = interpos(rhotor_gauss,pprime.*pprime_phys2chease,rhotor,-1);
  expeq.Pprime = interpos(rhopol,expeq.Pprime,expeq.rho,-0.1);
  expeq.TTprime = interpos(rhotor_gauss,ttprime.*ttprime_phys2chease,rhotor,-1);
  expeq.TTprime = interpos(rhopol,expeq.TTprime,expeq.rho,-0.1);
  expeq.extralines = {eqdsk_in.fnamefull};
  expeq.extralines{end+1,1} = ['eqdsk: q0 = ' num2str(eqdsk_in.q(1)) ' ; qedge = ' num2str(eqdsk_in.q(end))];
  expeq.extralines{end+1,1} = ['eqdsk: q95 = ' num2str(eqdsk_in.syscode_formulas.q95)];
  expeq.extralines{end+1,1} = ['raptor li3 = ' num2str(out.li3(it))];
  expeq.extralines{end+1,1} = ['raptor betan = ' num2str(out.betaN(it))];
  expeq.extralines{end+1,1} = ['raptor betapol = ' num2str(out.betapol(it))];
else
  return
end

% adapt kappa if needed
do_adapt_kappa=1;
if do_adapt_kappa
    % re-scale z for new kappa from zgeom
    zgeom = (max(expeq.RZ_psi(:,2))+min(expeq.RZ_psi(:,2)))./2.;
    kappa_ratio = kappa(end)./normalization_infos.kappa;
    znew = zgeom + (expeq.RZ_psi(:,2)-zgeom).* kappa_ratio;
    expeq.RZ_psi(:,2) = znew;    
end

if debug==1; write_expeq(expeq,['EXPEQ_nsttp1_' filename]); end

expeq.nsttp = 4;
expeq.TTprime = [];
expeq.jdotb_over_b0 = interpos(rhopol,jpar.*j_phys2chease,expeq.rho,-1);
if debug==1; write_expeq(expeq,['EXPEQ_nsttp4_' filename]); end

expeq.nsttp = 4;
expeq.nppfun = 8;
expeq.Pprime = [];
expeq.pressure = interpos(rhopol,p.*p_phys2chease,expeq.rho,-1);
expeq.jdotb_over_b0 = interpos(rhopol,jpar.*j_phys2chease,expeq.rho,-1);
if debug==1; write_expeq(expeq,['EXPEQ_nsttp4_p_' filename]); end

expeq.nrhotype=1;
expeq.pressure = interpos(rhotor,p.*p_phys2chease,expeq.rho,-1);
expeq.jdotb_over_b0 = interpos(rhotor,jpar.*j_phys2chease,expeq.rho,-1);
expeq_main_filename = [filename '_EXPEQ_nsttp4_p_rhotor'];
write_expeq(expeq,expeq_main_filename);

% % option of using q instead of j_parr (uncomment lines below)
% expeq.nsttp = 5;
% expeq.q =             interpos(rhopol,q,expeq.rho,-1);
% expeq.jdotb_over_b0 = [];
% expeq_main_filename = [filename '_EXPEQ_nsttp5_p_rhotor'];
% write_expeq(expeq,expeq_main_filename);

