% store_results: Internal routine to print and plot analysis results after each gate
% GUISDAP v.1.60 96-05-27 Copyright Asko Huuskonen and Markku Lehtinen
%
% See also: half_prof
% function store_results(aa,meas,var,result,alpha,chi2,status,kd2,p_coeffg,f_womega,p_om,pldfvv);
  function store_results(aa,meas,var,result,alpha,chi2,status,kd2,p_coeffg,f_womega,p_om,pldfvv);

 
global a_priori a_priorierror ad_lpg p_RECloc
global ch_el di_fit g_ind
global lpg_lag di_figures
global r_range r_status r_param r_dp r_error r_res r_apriori r_apriorierror r_ind

global sdate0 val_alt result_path d_time sv_acf
if val_alt, global plot_alt0, end

global height_arr1 r_lp_ind r_ac_ind
global sig_r_arr1 err_r_arr1 sig_i_arr1 
global err_i_arr1 fitted_r_arr1 fitted_i_arr1
global sig_r_arr2 err_r_arr2 sig_i_arr2 height_arr2
global err_i_arr2 fitted_r_arr2 fitted_i_arr2 sig_r_arr0
global status1 status2


% Scale residual and Xfer results to physical units
comp=aa(6);
aa=aa(1:5);
len=length(meas)-length(aa);
chi2=chi2/len;
err=covm2vec(alpha); 
res=scaled_to_real(result);
er=err;er(1:5)=scaled_to_real(err(1:5));
height=range_to_height(r_range(r_ind),ch_el(1));


% Print results to the console
if rem(r_ind,50)==1,
fprintf(' alt   Ne/1e11      Ti     Te/Ti      coll     velocity [O+]/Ne resid status\n')
end
fprintf('%5.1f',height);   % alt
fprintf(' %4.2f:%4.2f',res(1)/1e11,er(1)/1e11);   % Ne/1e11
fprintf(' %4.0f:%3.0f',res(2),er(2));   % Ti
fprintf(' %4.2f:%4.2f',res(3),er(3));   % Te/Ti
fprintf(' %4.1f:%4.1f',res(4)/1e3,er(4)/1e3);   % coll
fprintf(' %6.1f:%4.1f',res(5),er(5));   % velocity
fprintf(' %5.2f '     ,comp);
fprintf(' %4.2f  ',chi2);
if status==0, str='OK'; 
elseif status==1 str='Max iter';
elseif status==2 str='No fit done';
end
fprintf(str);
fprintf('\n')


% Store results to result variables
r_param(r_ind,:)=res;
r_dp(r_ind,:)=comp;
r_error(r_ind,:)=er; 
r_res(r_ind,:)=[chi2,sqrt(2/len)];
r_status(r_ind,:)=status;
r_apriori(r_ind,1:6)=scaled_to_real(a_priori(g_ind,:));
r_apriorierror(r_ind,1:6)=scaled_to_real(a_priorierror(g_ind,:));

p_di_fig3=di_figures(3);
if exist('plot_alt0') & di_figures(3),
 if round(height)~=plot_alt0,
  di_figures(3)=0;
 end
end


if height>=sv_acf(1) & height<=sv_acf(2),
%------------------------------------------------

 theo=dirthe([result,comp],p_coeffg,f_womega,kd2,p_om,pldfvv);
 res_err=err(1:5)';
 sig_err=sqrt(var);
 indr=1:len/2;
 indi=len/2+1:len;
 indp=len+(1:5);
 sig_r=meas(indr);err_r=sig_err(indr);
 sig_i=meas(indi);err_i=sig_err(indi);
 sig_p=meas(indp);err_p=sig_err(indp);
 fitted_r=theo(indr);
 fitted_i=theo(indi);
 fitted_p=theo(indp);
 indi=indi-length(indi);
 indp=1:5;

 if ~finite(er(1))|abs(er(1))>abs(res(1)) status=status+2^8; end
 if ~finite(er(2))|abs(er(2))>abs(res(2)) status=status+2^7; end
 if ~finite(er(3))|abs(res(2).*er(3)+res(3).*er(2))>abs(res(2).*res(3)) status=status+2^6; end
 if ~finite(er(5))|abs(er(5))>10000 status=status+2^5; end

 if isempty(sig_r_arr0) sig_r_arr0=sig_r; end

 if length(sig_r)~=length(sig_r_arr0), % for LP
  if isempty(height_arr2) r_lp_ind=1; end
  height_arr2(r_lp_ind)=height;
  status2(r_lp_ind)=status;
  sig_r_arr2(:,r_lp_ind)=sig_r;
  sig_i_arr2(:,r_lp_ind)=sig_i;
  err_r_arr2(:,r_lp_ind)=err_r;
  err_i_arr2(:,r_lp_ind)=err_i;
  fitted_r_arr2(:,r_lp_ind)=fitted_r;
  fitted_i_arr2(:,r_lp_ind)=fitted_i;
  r_lp_ind=r_lp_ind+1;
 end
 if length(sig_r)==length(sig_r_arr0), 
  if isempty(height_arr1) r_ac_ind=1; end
  height_arr1(r_ac_ind)=height;
  status1(r_ac_ind)=status;
  sig_r_arr1(:,r_ac_ind)=sig_r;
  sig_i_arr1(:,r_ac_ind)=sig_i;
  err_r_arr1(:,r_ac_ind)=err_r;
  err_i_arr1(:,r_ac_ind)=err_i;
  fitted_r_arr1(:,r_ac_ind)=fitted_r;
  fitted_i_arr1(:,r_ac_ind)=fitted_i;
  r_ac_ind=r_ac_ind+1;
 end

%------------------------------------------------
end %for (if height>=sv_acf(1) & height<=sv_acf(2))


if di_figures(3),
  theo=dirthe([result,comp],p_coeffg,f_womega,kd2,p_om,pldfvv);
  res_err=err(1:5)';
  sig_err=sqrt(var);
  indr=1:len/2;  
  indi=len/2+1:len;
  indp=len+(1:5);
  sig_r=meas(indr);err_r=sig_err(indr);
  sig_i=meas(indi);err_i=sig_err(indi);
  sig_p=meas(indp);err_p=sig_err(indp);
  fitted_r=theo(indr);
  fitted_i=theo(indi);
  fitted_p=theo(indp);
  indi=indi-length(indi);
  indp=1:5;
 
  figure(di_figures(3));clf
 
  set(gca,'Position',[.1 .1 .7 .8],'NextPlot','replace');
  plot(0,-1,'b.',indr,sig_r,'ro',indr,fitted_r,'g-'),hold on
  plot(indi+.2,sig_i,'bo',indi+.2,fitted_i,'g-');
  plot([indr;indr],[sig_r-err_r,sig_r+err_r]','r-')
  plot([indi+.2;indi+.2],[sig_i-err_i,sig_i+err_i]','b-')
  set(get(gca,'Children'),'MarkerSize',4)

  title([' Data  (o) and fit results (solid line) ' ....
         int2str(tosecs(d_time(2,:))) ' at ' int2str(height) ' km ']);
%  title(' Data  (o) and fit results (solid line)');
  ylabel('Power [K]'); xlabel(' # of data point'); 
 
  axes('Position',[.85 .1 .12 .8],'NextPlot','replace');
  plot(indp-0.15,sig_p,'ro',indp+0.15,fitted_p,'go'),hold on
  plot([indp;indp]-0.15,[sig_p-err_p,sig_p+err_p]','r:',...
         [indp;indp]+0.15,[fitted_p-res_err,fitted_p+res_err]','g-')
  set(get(gca,'Children'),'MarkerSize',4)
  axislim=[0 6 -1 max([5,ceil(fitted_p')])];
  if all(~isnan(axislim))
    axis(axislim);
  end
  set(gca,'Xtick',1:5,'XTickLabel',['N';'T';'r';'c';'v']);
  title('parameters')
  drawnow 

% if val_alt, 
 if exist('plot_alt0'), 
  ps_name=[result_path '/acf' sdate0];
  eval(['print -dpsc2 -append ' ps_name]),
  fitted_r,fitted_i,sig_r,sig_i,
 end
end


di_figures(3)=p_di_fig3;
