From 53df129841e1368c8643bd3b8f34b38d37af4001 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 1 Jun 2016 07:05:21 +0000
Subject: [PATCH] for generic gui as well
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5865 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/generic_fit.m | 158 +++++++++++++++++++++++++++++++++++++++++
crpptbx/line_average.m | 99 ++++++++++++++++++++++++++
2 files changed, 257 insertions(+)
create mode 100644 crpptbx/generic_fit.m
create mode 100644 crpptbx/line_average.m
diff --git a/crpptbx/generic_fit.m b/crpptbx/generic_fit.m
new file mode 100644
index 00000000..54a7f8c2
--- /dev/null
+++ b/crpptbx/generic_fit.m
@@ -0,0 +1,158 @@
+function [gen_signals_data] = generic_fit(shot,data_in,dataname,time_array_in,tension_xx_in,coeff_xx_in,in_option);
+%
+% [gen_signals_data,gen_signals_fit,filename_withfits] = generic_fit(shot,data_in,dataname,time_interval_in,tension_xx_in,coeff_xx_in,filename_withfits_in,in_option);
+%
+% Required inputs:
+%
+% shot: shot number
+% data_in{i}: array of structure from each "diagnostic" to combine, containing the following fields:
+% data_in{i}.data(:,it), data_in{i}.t, data_in{i}.error_bar(:,it)
+% and a subfield .x for the radial coordinate with:
+% data_in{i}.x(irho,it) or data_in{i}.x(irho) if the rho values do not change with time
+% and a subfield data_in{i}.xname which is used to describe if it is rhopol, rhotornorm, etc
+% dataname: name of data and of subfield to fill in in the gen_signals_data structure:
+%
+% It uses data_in{1}.t as the main time base thus as gen_signals_data.(dataname).t
+%
+% in_option: 0 (default) computes the fits from the provided data
+% 1 (combined data already performed thus uses data_in{1}.combined_data or data_in.combined_data to compute the fits
+%
+% output: gen_signals_data.(dataname).data, gen_signals_fit.(dataname).data
+%
+
+in_option_eff = 0;
+if exist('in_option') && ~isempty(in_option)
+ in_option_eff = in_option;
+end
+
+if in_option_eff == 0
+ % assumes receives just raw data in data_in{}
+ gen_signals_data.(dataname).raw_data = data_in;
+end
+if in_option_eff == 1
+ % assume one receives already gen_signals_data structure
+ if iscell(data_in)
+ gen_signals_data = data_in{1};
+ else
+ gen_signals_data = data_in;
+ end
+end
+
+if ~exist('coeff_xx_in') || length(coeff_xx_in)~=length(data_in)
+ coeff_xx_eff = ones(length(data_in),1);
+else
+ coeff_xx_eff = coeff_xx_in;
+end
+if exist('tension_xx_in') && ~isempty(tension_xx_in)
+ tension_xx_eff = tension_xx_in;
+else
+ tension_xx_eff=-3;
+end
+
+% combine the data to be able to do the fits
+if in_option_eff == 0
+
+ if isempty(data_in) || length(data_in)<1 || ~isfield(data_in{1},'t') || ~isfield(data_in{1},'data') || ~isfield(data_in{1},'x') || isempty(data_in{1}.data)
+ disp('no main data')
+ return
+ end
+
+ main_time_base = data_in{1}.t;
+
+ time_interval_in = [main_time_base(1) main_time_base(end)];
+ ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
+
+ gen_signals_data.(dataname).combined_data.t = main_time_base(ij);
+ gen_signals_data.(dataname).t = gen_signals_data.(dataname).combined_data.t;
+ 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
+ gen_signals_data.(dataname).combined_data.time_interval(i,1:2) = [t1(i) t2(i)];
+ end
+
+ for idata=1:length(data_in)
+ gen_signals_data.(dataname).combined_data.rholabel{idata} = data_in{idata}.xname;
+ end
+
+ for i=1:length(ij)
+ % for each time interval, collect the data
+ rhotor_data_tofit = [];
+ xxdata_tofit = [];
+ xxerr_tofit = [];
+ data_provenance = [];
+ for idata=1:length(data_in)
+ it_data{idata} = find(data_in{idata}.t>=t1(i) & data_in{idata}.t<t2(i));
+ % construct 1D array with data from data_in{idata}
+ if ~isempty(it_data{idata})
+ for it=1:length(it_data{idata})
+ jrho_data = find(data_in{idata}.x(:,it_data{idata}(it))>0 & data_in{idata}.x(:,it_data{idata}(it))<1.01);
+ if length(jrho_data)>0
+ rhotor_data_tofit(end+1:end+length(jrho_data)) = data_in{idata}.x(jrho_data,it_data{idata}(it));
+ xxdata_tofit(end+1:end+length(jrho_data)) = data_in{idata}.data(jrho_data,it_data{idata}(it));
+ xxerr_tofit(end+1:end+length(jrho_data)) = data_in{idata}.error_bar(jrho_data,it_data{idata}(it))./coeff_xx_eff(idata);
+ data_provenance(end+1:end+length(jrho_data)) = idata.*ones(size(jrho_data));
+ end
+ end
+ end
+ end
+
+ if isempty(xxdata_tofit)
+ disp('xxdata_tofit empty')
+ gen_signals_data.(dataname).combined_data.perDt{i}.rho = [];
+ gen_signals_data.(dataname).combined_data.perDt{i}.data = [];
+ gen_signals_data.(dataname).combined_data.perDt{i}.error_bar = [];
+ gen_signals_data.(dataname).combined_data.perDt{i}.provenance = [];
+ return
+ else
+ [rhoeff,irhoeff] = sort(rhotor_data_tofit);
+ xxdata_tofit_eff = xxdata_tofit(irhoeff);
+ xxerr_tofit_eff = xxerr_tofit(irhoeff);
+ data_provenance_eff = data_provenance(irhoeff);
+ gen_signals_data.(dataname).combined_data.perDt{i}.rho = rhoeff;
+ gen_signals_data.(dataname).combined_data.perDt{i}.data = xxdata_tofit_eff;
+ gen_signals_data.(dataname).combined_data.perDt{i}.error_bar = xxerr_tofit_eff;
+ gen_signals_data.(dataname).combined_data.perDt{i}.provenance = data_provenance_eff;
+ end
+ end
+end
+
+% do the fit
+doplot=0; % for debugging
+nb_points = 201;
+main_rho = linspace(0,1,nb_points);
+xxfit = [];
+dxxfitdrn = [];
+for i=1:length(gen_signals_data.(dataname).combined_data.perDt)
+ rhoeff = gen_signals_data.(dataname).combined_data.perDt{i}.rho;
+ xxdata_tofit_eff = gen_signals_data.(dataname).combined_data.perDt{i}.data;
+ xxerr_tofit_eff = gen_signals_data.(dataname).combined_data.perDt{i}.error_bar;
+ xeff = [0. rhoeff];
+ yeffxx = [xxdata_tofit_eff(1) xxdata_tofit_eff];
+ erreffxx = [100.*max(xxerr_tofit_eff) xxerr_tofit_eff];
+ [xxfit(:,i),dxxfitdrn(:,i)]=interpos(xeff,yeffxx,main_rho,tension_xx_eff,[1 0],[0 0],erreffxx);
+ if doplot
+ figure(11);clf
+ errorbar(rhoeff,xxdata_tofit_eff,xxerr_tofit_eff,'*');
+ hold all
+ plot(main_rho,xxfit(:,i))
+ pause(0.01)
+ end
+end
+
+gen_signals_data.(dataname).fit.rhofit = main_rho;
+gen_signals_data.(dataname).fit.data = xxfit;
+gen_signals_data.(dataname).fit.dydrho = dxxfitdrn;
+gen_signals_data.(dataname).shot = shot;
+gen_signals_data.(dataname).fit.tension = tension_xx_eff;
+gen_signals_data.(dataname).fit.coeff_datas = coeff_xx_eff;
+gen_signals_data.(dataname).fit.rholabel = gen_signals_data.(dataname).combined_data.rholabel{1};
+
+
diff --git a/crpptbx/line_average.m b/crpptbx/line_average.m
new file mode 100644
index 00000000..fb0d64c2
--- /dev/null
+++ b/crpptbx/line_average.m
@@ -0,0 +1,99 @@
+function [nel_out,s_out,nes_out,varargout] = line_average(shot,time,Rdiag,Zdiag,ne_in,rhopol_in,Rmesh,Zmesh,psi_RZmesh,psi_axis,psi_edge);
+%
+% compute s and ne(s) in s_out and nes_out; where s is distance from first point into the plasma along Rdiag,Zdiag chord
+% compute int[ne(s)ds]/int[ds] in nel_out
+%
+% can get Rmesh, Zmesh and psi_RZmesh from eqdsk at shot,time from gdat if not given, but for many times, should use equil output from gdat
+%
+
+nel_out = [];
+s_out = [];
+nes_out = [];
+
+time_eff = time;
+
+if ~exist('Rmesh') || isempty(Rmesh)
+ % do only 1st time
+ time_eff = time(1);
+ ab=gdat(shot,'eqdsk','time',time_eff);
+ eqdsk=ab.eqdsk;
+ if isempty(Rdiag)
+ % assume H-1 line
+ par=sf2par('DCN',shot,'H1LFS_R','DCNgeo');
+ Rlfs=par.value/1e3;
+ par=sf2par('DCN',shot,'H1LFS_z','DCNgeo');
+ Zlfs=par.value/1e3;
+ par=sf2par('DCN',shot,'H1HFS_R','DCNgeo');
+ Rhfs=par.value/1e3;
+ par=sf2par('DCN',shot,'H1HFS_z','DCNgeo');
+ Zhfs=par.value/1e3;
+ Rdiag=[Rlfs Rhfs];
+ Zdiag=[Zlfs Zhfs];
+ end
+ if length(Rdiag)==2
+ nbpoints=50;
+ Rout=linspace(Rdiag(1),Rdiag(2),nbpoints);
+ Zout=linspace(Zdiag(1),Zdiag(2),nbpoints);
+ elseif length(Rdiag)>10
+ Rout=Rdiag;
+ Zout=Zdiag;
+ else
+ disp(['unexpected nb of diag points: length(Rdiag) = ',length(Rdiag)]);
+ return
+ end
+ Rmesh = eqdsk.rmesh;
+ Zmesh = eqdsk.zmesh;
+ psi_RZmesh(:,:,1) = eqdsk.psirz;
+ psi_axis(1) = eqdsk.psiaxis;
+ psi_edge(1) = eqdsk.psiedge;
+end
+
+if prod(size(Rmesh)) == max(size(Rmesh))
+ Rmesh_eff = repmat(reshape(Rmesh,length(Rmesh),1),1,length(time_eff));
+ Zmesh_eff = repmat(reshape(Zmesh,length(Zmesh),1),1,length(time_eff));
+end
+
+for it=1:length(time_eff)
+ psi_at_routzout(:,it) = interpos2Dcartesian(Rmesh_eff(:,it),Zmesh(:,it),psi_RZmesh(:,:,it),Rout,Zout);
+ ij_in=find((psi_at_routzout(:,it)-psi_axis(it)).*(psi_at_routzout(:,it)-psi_edge(it))<=0);
+ if ~isempty(ij_in)
+ if ij_in(1)==1
+ Rstart=Rout(ij_in(1));
+ Zstart=Zout(ij_in(1));
+ else
+ Rstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Rout(ij_in(1)-1:ij_in(1))) + Rout(ij_in(1)-1);
+ Zstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Zout(ij_in(1)-1:ij_in(1))) + Zout(ij_in(1)-1);
+ end
+ if ij_in(end)==length(Rout)
+ Rend=Rout(ij_in(end));
+ Zend=Zout(ij_in(end));
+ else
+ Rend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Rout(ij_in(end):ij_in(end)+1)) + Rout(ij_in(end));
+ Zend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Zout(ij_in(end):ij_in(end)+1)) + Zout(ij_in(end));
+ end
+ s_eff = sqrt((Rout(ij_in)-Rstart).^2 + (Zout(ij_in)-Zstart).^2);
+ rhopol_eff=sqrt((psi_at_routzout(ij_in,it)-psi_axis(it))./(psi_edge(it)-psi_axis(it)));
+ if ij_in(1)==1
+ rhopol_eff = reshape(rhopol_eff,length(rhopol_eff),1);
+ s_eff = reshape(s_eff,length(s_eff),1);
+ else
+ rhopol_eff=[1 ;reshape(rhopol_eff,length(rhopol_eff),1); 1];
+ s_eff=[0 ;reshape(s_eff,length(s_eff),1); sqrt((Rend-Rstart).^2 + (Zend-Zstart).^2)];
+ end
+ [rho_in_eff,isort]=sort(rhopol_in(:,it));
+ if rho_in_eff(1)<1e-3
+ % assume 1st point is center
+ nes = interpos(rho_in_eff,ne_in(isort,it),rhopol_eff,-0.1,[1 0],[0 0]);
+ else
+ rho_in_eff = [0 ; rho_in_eff];
+ ne_in_eff = [ne_in(isort(1),it); ne_in(isort,it)];
+ sigma_in = ones(size(ne_in_eff));
+ sigma_in(1)=100;
+ nes = interpos(rho_in_eff,ne_in_eff,rhopol_eff,-1,[1 0],[0 0],sigma_in);
+ end
+ nes_out(:,it) = nes;
+ s_out(:,it) = s_eff;
+ [dum1,~,~,nel_out_prof] = interpos(s_eff,nes,-0.1);
+ nel_out(it) = nel_out_prof(end)./(s_eff(end)-s_eff(1));
+ end
+end
--
GitLab