CEZ_CMZ_fit.m

function [filename_withfits,filename_data,cez_cmz_data,cez_cmz_fit] = CEZ_CMZ_fit(shot,time_interval_in,tension_ti_in,tension_vrot_in,coeff_ti_cmz_in,coeff_vrot_cmz_in);
%
% [filename_withfits,filename_data,cez_cmz_data,cez_cmz_fit] = CEZ_CMZ_fit(shot,time_interval_in,tension_ti_in,tension_vrot_in,coeff_ti_cmz_in,coeff_vrot_cmz_in);
%

filename_data = [num2str(shot) '_CEZ_CMZ.mat'];
filename_withfits = [num2str(shot) '_CEZ_CMZ_forgui.mat'];
if exist(filename_data,'file')
  load(filename_data)
else
  cxrs=gdat(shot,'cxrs_rho','doplot',1);
  cmz=gdat(shot,'cxrs_rho','doplot',1,'source','cmz');
end
cez_cmz_data = [];
cez_cmz_fit = [];

if isempty(cxrs.data) || isempty(cxrs.t) || ischar(cxrs.data)
  disp('no main cxrs')
  return
end
if ~exist(filename_data,'file')
  eval(['save ' filename_data ' cxrs cmz'])
end

main_time_base = cxrs.t;
if ~exist('time_interval_in') || isempty(time_interval_in)
  time_interval_in = [cxrs.t(1) cxrs.t(end)];
  ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
else
  ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
  if length(ij)==1
    disp('not yet ready')
    return
  end
end

cez_cmz_data.time = main_time_base(ij);
for i=1:length(ij)
  if i==1
    t1(i) = main_time_base(ij(i)) - 0.5*diff(main_time_base(ij([i:i+1])));
  else
    t1(i) = mean(main_time_base(ij([i-1:i])));
  end
  if i==length(ij)
    t2(i) = main_time_base(ij(i)) + 0.5*diff(main_time_base(ij([i-1:i])));
  else
    t2(i) = mean(main_time_base(ij([i:i+1])));
  end
  cez_cmz_data.time_interval(i,1:2) = [t1(i) t2(i)];
end

rhotorfit=linspace(0,1,200);
ticoeff_err_cxrs=1.;
if exist('coeff_ti_cmz_in') && ~isempty(coeff_ti_cmz_in)
  ticoeff_err_cmz=coeff_ti_cmz_in; % (if weight is different)
else
  ticoeff_err_cmz=1.0; % (if weight is different)
end
vrotcoeff_err_cxrs=1.;
if exist('coeff_vrot_cmz_in') && ~isempty(coeff_vrot_cmz_in)
  vrotcoeff_err_cmz=coeff_vrot_cmz_in; % (if weight is different)
else
  vrotcoeff_err_cmz=0.1; % (if weight is different)
end
if exist('tension_ti_in') && ~isempty(tension_ti_in)
  tension_ti_eff = tension_ti_in;
else
  tension_ti_eff=-3;
end
if exist('tension_vrot_in') && ~isempty(tension_vrot_in)
  tension_vrot_eff = tension_vrot_in;
else
  tension_vrot_eff=-3;
end

tifit = [];
dtifitdrn = [];

doplot=0;
for i=1:length(ij)
  it_cxrs = find(cxrs.t>=t1(i) & cxrs.t<t2(i));
  it_cmz = find(cmz.t>=t1(i) & cmz.t<t2(i));
  % construct 1D array with data from both cxrs, cmz
  rhotor_data_tofit = []; 
  tidata_tofit = [];
  vrotdata_tofit = [];
  tierr_tofit = [];
  vroterr_tofit = [];
  if ~isempty(it_cxrs)
    for it=1:length(it_cxrs)
      idata = find(cxrs.ti.data(:,it_cxrs(it))>0 & cxrs.rhotornorm(:,it_cxrs(it))<1.01);
      if length(idata)>0
	rhotor_data_tofit(end+1:end+length(idata)) = cxrs.rhotornorm(idata,it_cxrs(it));
	tidata_tofit(end+1:end+length(idata)) = cxrs.ti.data(idata,it_cxrs(it));
	vrotdata_tofit(end+1:end+length(idata)) = cxrs.vrot.data(idata,it_cxrs(it));
	tierr_tofit(end+1:end+length(idata)) = cxrs.ti.error_bar(idata,it_cxrs(it))./ticoeff_err_cxrs;
	vroterr_tofit(end+1:end+length(idata)) = cxrs.vrot.error_bar(idata,it_cxrs(it))./vrotcoeff_err_cxrs;
      end
    end
  end
  if ~isempty(it_cmz)
    for it=1:length(it_cmz)
      idata = find(cmz.ti.data(:,it_cmz(it))>0 & cmz.rhotornorm(:,it_cmz(it))<1.01);
      if length(idata)>0
	rhotor_data_tofit(end+1:end+length(idata)) = cmz.rhotornorm(idata,it_cmz(it));
	tidata_tofit(end+1:end+length(idata)) = cmz.ti.data(idata,it_cmz(it));
	vrotdata_tofit(end+1:end+length(idata)) = cmz.vrot.data(idata,it_cmz(it));
	tierr_tofit(end+1:end+length(idata)) = cmz.ti.error_bar(idata,it_cmz(it))./ticoeff_err_cmz;
	vroterr_tofit(end+1:end+length(idata)) = cmz.vrot.error_bar(idata,it_cmz(it))./vrotcoeff_err_cmz;
      end
    end
    rhotor_data_tofit_cmz = cmz.rhotornorm(:,it_cmz);
    tidata_tofit_cmz = cmz.ti.data(:,it_cmz);
    vrotdata_tofit_cmz = cmz.vrot.data(:,it_cmz);
    tierr_tofit_cmz = cmz.ti.error_bar(:,it_cmz);
    vroterr_tofit_cmz = cmz.vrot.error_bar(:,it_cmz);
  end
  [rhoeff,irhoeff] = sort(rhotor_data_tofit);
  tidata_tofit_eff = tidata_tofit(irhoeff);
  if isempty(tidata_tofit_eff)
    disp('tidata_tofit_eff empty')
    return
  else
    vrotdata_tofit_eff = vrotdata_tofit(irhoeff);
    tierr_tofit_eff = tierr_tofit(irhoeff);
    vroterr_tofit_eff = vroterr_tofit(irhoeff);
    cez_cmz_data.ti{i}.rhotornorm = rhoeff;
    cez_cmz_data.ti{i}.data = tidata_tofit_eff;
    cez_cmz_data.ti{i}.err = tierr_tofit_eff;
    cez_cmz_data.vrot{i}.rhotornorm = rhoeff;
    cez_cmz_data.vrot{i}.data = vrotdata_tofit_eff;
    cez_cmz_data.vrot{i}.err = vroterr_tofit_eff;
    xeff = [0. rhoeff];
    yeffti = [tidata_tofit_eff(1) tidata_tofit_eff];
    yeffvrot = [vrotdata_tofit_eff(1) vrotdata_tofit_eff];
    erreffti = [100.*max(tierr_tofit_eff) tierr_tofit_eff];
    erreffvrot = [100.*max(vroterr_tofit_eff) vroterr_tofit_eff];
    [tifit(:,i),dtifitdrn(:,i)]=interpos(xeff,yeffti,rhotorfit,tension_ti_eff,[1 0],[0 0],erreffti);
    [vrotfit(:,i),dvrotfitdrn(:,i)]=interpos(xeff,yeffvrot,rhotorfit,tension_vrot_eff,[1 0],[0 0],erreffvrot);
    if doplot
      figure(11);clf
      errorbar(rhoeff,tidata_tofit_eff,tierr_tofit_eff,'*');
      hold all
      plot(rhotorfit,tifit(:,i))
      figure(12);clf
      errorbar(rhoeff,vrotdata_tofit_eff,vroterr_tofit_eff,'*');
      hold all
      plot(rhotorfit,vrotfit(:,i))
      pause(0.01)
    end
  end
end

cez_cmz_fit.rhotornorm = rhotorfit;
cez_cmz_fit.ti = tifit;
cez_cmz_fit.dtidrhotornorm = dtifitdrn;
cez_cmz_fit.vrot = vrotfit;
cez_cmz_fit.dvrotdrhotornorm = dvrotfitdrn;
cez_cmz_fit.shot = shot;
cez_cmz_fit.tension_ti = tension_ti_eff;
cez_cmz_fit.tension_vrot = tension_vrot_eff;
cez_cmz_fit.ticoeff_err_cxrs = ticoeff_err_cxrs;
cez_cmz_fit.ticoeff_err_cmz = ticoeff_err_cmz;
cez_cmz_fit.vrotcoeff_err_cxrs = vrotcoeff_err_cxrs;
cez_cmz_fit.vrotcoeff_err_cmz = vrotcoeff_err_cmz;
cez_cmz_data.shot = shot;

eval(['save ' filename_withfits ' cez_cmz_fit cez_cmz_data'])