From e6f104bc63b8aa5ff0f66323991ff4c2606ecceb Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Sun, 27 Oct 2019 19:02:48 +0100
Subject: [PATCH] corrected time_out with grids 1d and use interpos linear cst
for extrapolating
---
matlab/TCV_IMAS/tcv_get_ids_equilibrium.m | 3 +-
matlab/gdat2time_out.m | 201 +++++++++++++++++-----
matlab/get_grids_1d.m | 8 +-
3 files changed, 171 insertions(+), 41 deletions(-)
diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
index 195391d1..fdc6b9f8 100644
--- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
+++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m
@@ -29,7 +29,8 @@ params_eff.data_request='b0'; params_eff.source='liuqe'; % to get liuqe time arr
bb = gdat(params_equilibrium.shot,params_eff);
params_eff = rmfield(params_eff,'source'); % to get original magnetics data
vacuum_toroidal_field.b0=gdat(params_equilibrium.shot,params_eff);
-vacuum_toroidal_field.b0.data = interp1(vacuum_toroidal_field.b0.t,vacuum_toroidal_field.b0.data,bb.t);
+ij_ok = [isfinite(vacuum_toroidal_field.b0.data)];
+vacuum_toroidal_field.b0.data = interpos(21,vacuum_toroidal_field.b0.t(ij_ok),vacuum_toroidal_field.b0.data(ij_ok),bb.t);
vacuum_toroidal_field.b0.t = bb.t;
vacuum_toroidal_field.b0.dim = {vacuum_toroidal_field.b0.t};
vacuum_toroidal_field_desc.b0 = ['''b0'',''timing source'',''liuqe=' num2str(gdat_params.liuqe) ''''];
diff --git a/matlab/gdat2time_out.m b/matlab/gdat2time_out.m
index 16e4f442..d8a50238 100644
--- a/matlab/gdat2time_out.m
+++ b/matlab/gdat2time_out.m
@@ -50,6 +50,7 @@ end
try
time_out = sort(gdat_data_in.gdat_params.time_out);
+ time_out = reshape(time_out,1,numel(time_out)); % usual shape if input as [t1 t2] or [t1,t2];
catch ME
getReport(ME)
error(['expect gdat_data_in.gdat_params.time_out'])
@@ -57,6 +58,19 @@ end
gdat_data_out = gdat_data_in;
+% do this special part first otherwise ref time may have changed
+if any(option==[21,22])
+ % add extra fields then correct
+ ab_tmp_rho = gdat_data_out.grids_1d;
+ ab_tmp_rho.data = gdat_data_out.grids_1d.rhotornorm; ab_tmp_rho.t = gdat_data_in.t;
+ ab_tmp_rho.gdat_params = gdat_data_out.gdat_params; ab_tmp_rho.mapping_for = gdat_data_out.mapping_for;
+ extra_fields = {'psi','rhotornorm','rhovolnorm','rhopolnorm_equil','rhotornorm_equil', ...
+ 'rhovolnorm_equil','rhotor_edge','volume_edge','b0'};
+ if option==22; extra_fields{end+1} = 'rhopolnorm'; end
+ [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields);
+ gdat_data_out.grids_1d = rmfield(ab_tmp_rho,{'data','t','gdat_params','mapping_for'});
+end
+
if any(option == [1, 21, 22])
if numel(time_out) == 2
% time interval provided, take all existing time points within this interval
@@ -148,7 +162,7 @@ if any(option == [1, 21, 22])
% default time_out_interp='linear', will implement {'interpos','tension','option',..} later when needed
% time interval provided, take all existing time points within this interval
ij = iround_os(gdat_data_out.t,time_out);
- gdat_data_out.t = time_out;
+ gdat_data_out.t = time_out'; % usual shape n X 1 for time array in gdat, at least when 1D array
if isfield(gdat_data_out,'dim'); gdat_data_out.dim{dim_time} = gdat_data_out.t; end
nbdims = numel(size(gdat_data_out.data));
nbdims_eff = nbdims;
@@ -157,10 +171,25 @@ if any(option == [1, 21, 22])
end
switch nbdims_eff
case 1
- gdat_data_out.data = interp1(gdat_data_in.t,gdat_data_in.data,gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data));
+ if numel(ij_ok) > 1
+ gdat_data_out.data = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data = gdat_data_in.data(ij_ok) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data = NaN .* ones(size(gdat_data_out.t));
+ end
if iothers_as_data
for j=1:length(others_as_data_eff)
- eval(['gdat_data_out.' others_as_data_eff{j} ' = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} ',gdat_data_out.t);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ij_ok),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = gdat_data_in.' others_as_data_eff{j} '(ij_ok) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
case 2
@@ -168,65 +197,146 @@ if any(option == [1, 21, 22])
if dim_time ==1
gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2));
for ix=1:size(gdat_data_out.data,2)
- gdat_data_out.data([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix),gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data(:,ix)));
+ if numel(ij_ok) > 1
+ gdat_data_out.data([1:len_t],ix) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok,ix),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data([1:len_t],ix) = gdat_data_in.data(ij_ok,ix) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data([1:len_t],ix) = NaN .* ones(size(gdat_data_out.t));
+ end
end
if iothers_as_data
for j=1:length(others_as_data_eff)
- eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,size(gdat_data_out.' others_as_data_eff{j} ',2));']);
- end
- for j=1:length(others_as_data_eff)
- ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
- for ix=1:ix_len
- eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix),gdat_data_out.t);']);
+ eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
+ if ~iscell(atmp)
+ ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,ix_len);']);
+ for ix=1:ix_len
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(:,ix)));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ij_ok,ix),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = gdat_data_in.' ...
+ others_as_data_eff{j} '(ij_ok,ix) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix) = NaN .* ones(size(gdat_data_out.t));']);
+ end
+ end
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '{1} = NaN*ones(len_t,1);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '{1}));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '{1} = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '{1}(ij_ok),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '{1} = gdat_data_in.' ...
+ others_as_data_eff{j} '{1}(ij_ok) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '{1} = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
end
else
gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t);
for ix=1:size(gdat_data_out.data,1)
- gdat_data_out.data(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:),gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data(ix,:)));
+ if numel(ij_ok) > 1
+ gdat_data_out.data(ix,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,ij_ok),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data(ix,[1:len_t]) = gdat_data_in.data(ix,ij_ok) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data(ix,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));
+ end
end
if iothers_as_data
for j=1:length(others_as_data_eff)
- eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(size(gdat_data_out.' others_as_data_eff{j} ',1),len_t);']);
- end
- for j=1:length(others_as_data_eff)
- ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
- for ix=1:ix_len
- eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:),gdat_data_out.t);']);
+ eval(['atmp = gdat_data_out.' others_as_data_eff{j} ';']);
+ if ~iscell(atmp)
+ ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',1);']);
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(ix_len,len_t);']);
+ for ix=1:ix_len
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,:)));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,ij_ok),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,ij_ok) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));']);
+ end
+ end
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '{2} = NaN*ones(1,len_t);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '{2}));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '{2} = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '{2}(ij_ok),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '{2} = gdat_data_in.' ...
+ others_as_data_eff{j} '{2}(ij_ok) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '{2} = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
end
end
case 3
+ % do not assume can have a cell case until need to deal with it, since not sure what dim is left
len_t = numel(ij);
if dim_time == 1
gdat_data_out.data = NaN*ones(len_t,size(gdat_data_out.data,2),size(gdat_data_out.data,3));
for ix=1:size(gdat_data_out.data,2)
for jy=1:size(gdat_data_out.data,3)
- gdat_data_out.data([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.data(:,ix,jy),gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data(:,ix,jy)));
+ if numel(ij_ok) > 1
+ gdat_data_out.data([1:len_t],ix,jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ij_ok,ix,jy),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data([1:len_t],ix,jy) = gdat_data_in.data(ij_ok,ix,jy) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data([1:len_t],ix,jy) = NaN .* ones(size(gdat_data_out.t));
+ end
end
end
if iothers_as_data
- for j=1:length(others_as_data_eff)
- eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,size(gdat_data_out.' others_as_data_eff{j} ...
- ',2),size(gdat_data_out.' others_as_data_eff{j} ',3));']);
- end
for j=1:length(others_as_data_eff)
ix_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']);
+ eval(['gdat_data_out.' others_as_data_eff{j} ' = NaN*ones(len_t,ix_len,jy_len);']);
for ix=1:ix_len
for jy=1:jy_len
- eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(:,ix,jy),gdat_data_out.t);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(:,ix,jy)));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ij_ok,ix,jy),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = gdat_data_in.' ...
+ others_as_data_eff{j} '(ij_ok,ix,jy) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '([1:len_t],ix,jy) = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
end
+ for j=1:length(others_as_data_eff)
+ end
end
elseif dim_time == 2
gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),len_t,size(gdat_data_out.data,3));
for ix=1:size(gdat_data_out.data,1)
for jy=1:size(gdat_data_out.data,3)
- gdat_data_out.data(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.data(ix,:,jy),gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data(ix,:,jy)));
+ if numel(ij_ok) > 1
+ gdat_data_out.data(ix,[1:len_t],jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,ij_ok,jy),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data(ix,[1:len_t],jy) = gdat_data_in.data(ix,ij_ok,jy) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data(ix,[1:len_t],jy) = NaN .* ones(size(gdat_data_out.t));
+ end
end
end
if iothers_as_data
@@ -239,7 +349,16 @@ if any(option == [1, 21, 22])
jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',3);']);
for ix=1:ix_len
for jy=1:jy_len
- eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,:,jy),gdat_data_out.t);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,:,jy)));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,ij_ok,jy),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,ij_ok,jy) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,[1:len_t],jy) = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
end
@@ -248,7 +367,14 @@ if any(option == [1, 21, 22])
gdat_data_out.data = NaN*ones(size(gdat_data_out.data,1),size(gdat_data_out.data,2),len_t);
for ix=1:size(gdat_data_out.data,1)
for jy=1:size(gdat_data_out.data,2)
- gdat_data_out.data(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.data(ix,jy,:),gdat_data_out.t);
+ ij_ok = find(isfinite(gdat_data_in.data(ix,jy,:)));
+ if numel(ij_ok) > 1
+ gdat_data_out.data(ix,jy,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.data(ix,jy,ij_ok),gdat_data_out.t);
+ elseif numel(ij_ok) == 1
+ gdat_data_out.data(ix,jy,[1:len_t]) = gdat_data_in.data(ix,jy,ij_ok) .* ones(size(gdat_data_out.t));
+ else
+ gdat_data_out.data(ix,jy,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));
+ end
end
end
if iothers_as_data
@@ -261,7 +387,16 @@ if any(option == [1, 21, 22])
jy_len = eval(['size(gdat_data_out.' others_as_data_eff{j} ',2);']);
for ix=1:ix_len
for jy=1:jy_len
- eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = interp1(gdat_data_in.t,gdat_data_in.' others_as_data_eff{j} '(ix,jy,:),gdat_data_out.t);']);
+ ij_ok = eval(['find(isfinite(gdat_data_in.' others_as_data_eff{j} '(ix,jy,:)));']);
+ if numel(ij_ok) > 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = interpos(21,gdat_data_in.t(ij_ok),gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,jy,ij_ok),gdat_data_out.t);']);
+ elseif numel(ij_ok) == 1
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = gdat_data_in.' ...
+ others_as_data_eff{j} '(ix,jy,ij_ok) .* ones(size(gdat_data_out.t));']);
+ else
+ eval(['gdat_data_out.' others_as_data_eff{j} '(ix,jy,[1:len_t]) = NaN .* ones(size(gdat_data_out.t));']);
+ end
end
end
end
@@ -274,15 +409,3 @@ if any(option == [1, 21, 22])
else
error(['option = ' num2str(option) ' not yet implemented'])
end
-
-if any(option==[21,22])
- % add extra fields then correct
- ab_tmp_rho = gdat_data_out.grids_1d;
- ab_tmp_rho.data = gdat_data_out.grids_1d.rhotornorm; ab_tmp_rho.t = gdat_data_out.t;
- ab_tmp_rho.gdat_params = gdat_data_out.gdat_params; ab_tmp_rho.mapping_for = gdat_data_out.mapping_for;
- extra_fields = {'psi','rhotornorm','rhovolnorm','rhopolnorm_equil','rhotornorm_equil', ...
- 'rhovolnorm_equil','rhotor_edge','volume_edge','b0'};
- if option==22; extra_fields{end+1} = 'rhopolnorm'; end
- [ab_tmp_rho] = gdat2time_out(ab_tmp_rho,1,extra_fields);
- gdat_data_out.grids_1d = rmfield(ab_tmp_rho,{'data','t','gdat_params','mapping_for'});
-end
diff --git a/matlab/get_grids_1d.m b/matlab/get_grids_1d.m
index f89a3b0f..97931339 100644
--- a/matlab/get_grids_1d.m
+++ b/matlab/get_grids_1d.m
@@ -27,7 +27,13 @@ if (nopt == 0) || isempty(gdat_data.x) || isempty(gdat_data.t) || isempty(gdat_d
gdat_data.grids_1d.b0 = [];
return
end
-params_eff = gdat_data.gdat_params;params_eff.doplot=0;
+% time_out extraction is done at end within gdat, otherwise may loose information
+if isfield(gdat_data.gdat_params,'time_out')
+ params_eff = rmfield(gdat_data.gdat_params,'time_out');
+else
+ params_eff = gdat_data.gdat_params;
+end
+params_eff.doplot=0;
params_eff.data_request='rhotor_norm';
rhotor_norm = gdat(gdat_data.shot,params_eff);
ndim_x_rhotor = length(find(size(rhotor_norm.x)>1));
--
GitLab