From 485172340721ce3e1aad3486e1bede4e3b0b91c4 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Fri, 12 Apr 2019 08:06:13 +0000
Subject: [PATCH] start adding icd sources

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11737 d63d8f72-b253-0410-a779-e742ad2e26cf
---
 crpptbx/TCV/gdat_tcv.m             | 85 ++++++++++++++++++++++++++++--
 crpptbx/TCV/tcv_requests_mapping.m |  7 ++-
 crpptbx/get_grids_1d.m             | 30 ++++++-----
 3 files changed, 106 insertions(+), 16 deletions(-)

diff --git a/crpptbx/TCV/gdat_tcv.m b/crpptbx/TCV/gdat_tcv.m
index c5bdef32..97177bd4 100644
--- a/crpptbx/TCV/gdat_tcv.m
+++ b/crpptbx/TCV/gdat_tcv.m
@@ -1138,10 +1138,10 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-   case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data'}
+   case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data','ohm_data', 'bs_data'}
     % note: same time array for all at main.data level, then individual at .ec, .nbi levels
     % At this stage fill just eccd, later add nbi
-    sources_avail = {'ec','nbi'}; % can be set in parameter source
+    sources_avail = {'ec','nbi','ohm','bs'}; % can be set in parameter source
     % create empty structures for all, so in return one always have same substructres
     for i=1:length(sources_avail)
       gdat_data.(sources_avail{i}).data = [];
@@ -1158,7 +1158,16 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     end
 
     if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source)
-      gdat_data.gdat_params.source = sources_avail;
+      switch data_request_eff
+       case 'ec_data'
+        gdat_data.gdat_params.source = {'ec'};
+       case 'ohm_data'
+        gdat_data.gdat_params.source = {'ohm'};
+       case 'bs_data'
+        gdat_data.gdat_params.source = {'bs'};
+       otherwise
+        gdat_data.gdat_params.source = sources_avail;
+      end
     elseif ~iscell(gdat_data.gdat_params.source)
       if ischar(gdat_data.gdat_params.source)
 	gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source);
@@ -1187,6 +1196,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
     % create structure for icd sources from params and complete with defaults
     source_icd.ec = 'toray';
     source_icd.nbi = '';
+    source_icd.ohm = 'ibs';
+    source_icd.bs = 'ibs';
     for i=1:length(gdat_data.gdat_params.source)
       if ~isfield(gdat_data.gdat_params,['source_' gdat_data.gdat_params.source{i}])
         gdat_data.gdat_params.(['source_' gdat_data.gdat_params.source{i}]) = source_icd.(gdat_data.gdat_params.source{i});
@@ -1388,6 +1399,74 @@ elseif strcmp(mapping_for_tcv.method,'switchcase')
       end
     end
     %
+    if any(strmatch('ohm',gdat_data.gdat_params.source))
+      data_fullpath = '';
+      ohm_help = '';
+      if strcmp(lower(source_icd.ohm),'ibs')
+        ohm_data.cd_tot = gdat_tcv(shot,'\results::ibs:iohm');
+        ohm_data.cd_tot.data = ohm_data.cd_tot.data';
+        ohm_data.cd_tot.units = strrep(ohm_data.cd_tot.units,'kA','A');
+        ohm_data.cd_dens = gdat_tcv(shot,'\results::ibs:johmav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
+        ohm_data.cd_dens.units = strrep(ohm_data.cd_dens.units,'kA','A');
+        abc=get_grids_1d(ohm_data.cd_dens,1,1);
+        ohm_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
+        ohm_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
+        ohm_help = 'from IBS ohm related nodes, from Iohm=Ip-Icd-Ibs assuming stationary state';
+      else
+        disp(['source_icd.ohm = ' source_icd.ohm ' not yet implemented, ask O. Sauter'])
+        for subfields={'cd_tot','cd_dens'};
+          for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
+            ohm_data.(subfields{1}).(subsubfields{1}) = [];
+          end
+        end
+      end
+      gdat_data.ohm.ohm_data = ohm_data;
+      data_fullpath = ['from ibs:ohmic related nodes, all results in .ohm_data, subfield=' field_for_main_data ...
+                       'in ohm.data, .x, .t, .dim, .dimunits, .label, .units'];
+      for subfields={'data','x','t','units','dim','dimunits','label'}
+        gdat_data.ohm.(subfields{1}) = gdat_data.ohm.ohm_data.(field_for_main_data).(subfields{1});
+      end
+      if isempty(gdat_data.t); gdat_data.t = gdat_data.ohm.t; end
+      gdat_data.ohm.data_fullpath = data_fullpath;
+      gdat_data.ohm.help = ohm_help;
+      gdat_data.data(end+1,:) = interpos(gdat_data.ohm.t,gdat_data.ohm.data,gdat_data.t,-0.1);
+      gdat_data.label{end+1}=gdat_data.ohm.label;
+    end
+    %
+    if any(strmatch('bs',gdat_data.gdat_params.source))
+      data_fullpath = '';
+      bs_help = '';
+      if strcmp(lower(source_icd.bs),'ibs')
+        bs_data.cd_tot = gdat_tcv(shot,'\results::ibs:ibs');
+        bs_data.cd_tot.data = bs_data.cd_tot.data';
+        bs_data.cd_tot.units = strrep(bs_data.cd_tot.units,'kA','A');
+        bs_data.cd_dens = gdat_tcv(shot,'\results::ibs:jbsav'); % =jtild=<j.B> / R0 / <Bphi/R>, with Bphi=F(psi)/R
+        bs_data.cd_dens.units = strrep(bs_data.cd_dens.units,'kA','A');
+        abc=get_grids_1d(bs_data.cd_dens,1,1);
+        bs_data.cd_dens.rhotor_norm = abc.grids_1d.rhotornorm;
+        bs_data.cd_dens.rhopol_norm = abc.grids_1d.rhopolnorm;
+        bs_help = 'from IBS bs related nodes, from Ibs=Ip-Icd-Ibs assuming stationary state';
+      else
+        disp(['source_icd.bs = ' source_icd.bs ' not yet implemented, ask O. Sauter'])
+        for subfields={'cd_tot','cd_dens'};
+          for subsubfields={'data','x','t','units','dim','dimunits','label','rhotor_norm','rhopol_norm'}
+            bs_data.(subfields{1}).(subsubfields{1}) = [];
+          end
+        end
+      end
+      gdat_data.bs.bs_data = bs_data;
+      data_fullpath = ['from ibs:bsic related nodes, all results in .bs_data, subfield=' field_for_main_data ...
+                       'in bs.data, .x, .t, .dim, .dimunits, .label, .units'];
+      for subfields={'data','x','t','units','dim','dimunits','label'}
+        gdat_data.bs.(subfields{1}) = gdat_data.bs.bs_data.(field_for_main_data).(subfields{1});
+      end
+      if isempty(gdat_data.t); gdat_data.t = gdat_data.bs.t; end
+      gdat_data.bs.data_fullpath = data_fullpath;
+      gdat_data.bs.help = bs_help;
+      gdat_data.data(end+1,:) = interpos(gdat_data.bs.t,gdat_data.bs.data,gdat_data.t,-0.1);
+      gdat_data.label{end+1}=gdat_data.bs.label;
+    end
+    %
     % add all to last index of .data(:,i)
     gdat_data.data(end+1,:) = nansum(gdat_data.data,1);
     gdat_data.x = [1:size(gdat_data.data,1)];
diff --git a/crpptbx/TCV/tcv_requests_mapping.m b/crpptbx/TCV/tcv_requests_mapping.m
index 1f989879..124ffeef 100644
--- a/crpptbx/TCV/tcv_requests_mapping.m
+++ b/crpptbx/TCV/tcv_requests_mapping.m
@@ -48,6 +48,11 @@ switch lower(data_request)
   mapping.label = 'area';
   mapping.method = 'tdiliuqe';
   mapping.expression = 'tcv_eq(''''area'''',''''LIUQE.M'''')';
+ case 'area_rho'
+  mapping.timedim = 2;
+  mapping.label = 'area_rho';
+  mapping.method = 'tdiliuqe';
+  mapping.expression = 'tcv_eq(''''area_rho'''',''''LIUQE.M'''')';
  case 'area_edge'
   mapping.timedim = 1;
   mapping.label = 'area\_lcfs';
@@ -229,7 +234,7 @@ switch lower(data_request)
 		    'gdat_tmp2=gdat_tcv(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ...
 		    'tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ...
 		    'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);'];
- case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data'}
+ case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data', 'ohm_data', 'bs_data'}
   mapping.timedim = 1;
   mapping.label = 'various Pdens, Icd, jcd';
   mapping.method = 'switchcase';
diff --git a/crpptbx/get_grids_1d.m b/crpptbx/get_grids_1d.m
index 12da7bd6..96c59a1f 100644
--- a/crpptbx/get_grids_1d.m
+++ b/crpptbx/get_grids_1d.m
@@ -1,5 +1,7 @@
 function [gdat_data] = get_grids_1d(gdat_data,nbdim_x,nopt,nverbose);
 %
+% [gdat_data] = get_grids_1d(gdat_data,nbdim_x,nopt,nverbose);
+%
 % add various rhos in grids_1d
 %
 % assume gdat_data.x=rhopol, gdat_data.t, and gdat_data(x,t)
@@ -37,7 +39,8 @@ params_eff.data_request='psi_axis';
 psi_axis = gdat(gdat_data.shot,params_eff);
 params_eff.data_request='psi_edge';
 psi_edge = gdat(gdat_data.shot,params_eff);
-psi0_edge = interpos(63,psi_axis.t,psi_axis.data-psi_edge.data,gdat_data.t,-0.01);
+ij=~isnan(psi_axis.t);
+psi0_edge = interpos(63,psi_axis.t(ij),psi_axis.data(ij)-psi_edge.data(ij),gdat_data.t,-0.01);
 if (nbdim_x == 1)
   gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2*reshape(psi0_edge,1,length(psi0_edge));
 elseif (nbdim_x == 2)
@@ -53,7 +56,6 @@ if (isempty(rhotor_norm.x) ||isempty(rhotor_norm.t) || isempty(rhotor_norm.data)
       || (isempty(rhovol.x) ||isempty(rhovol.t) || isempty(rhovol.data))
   return
 end
-
 it_rt = iround_os(rhotor_norm.t,gdat_data.t);
 it_vol = iround_os(rhovol.t,gdat_data.t);
 for it=1:length(gdat_data.t)
@@ -66,36 +68,40 @@ for it=1:length(gdat_data.t)
     ii=find(isfinite(gdat_data.grids_1d.rhopolnorm(:,it)));
   end
   if (nbdim_x == 1)
-    if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
+% $$$     if length(ii)==length(gdat_data.grids_1d.rhopolnorm)
+    if length(ii) >0
+      nb_ii = length(ii);
       if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
         if ndim_x_rhotor==1
-          gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
+          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
         else
-          gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm);
+          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii));
         end
       end
       if ~isempty(rhovol.x) && ~isempty(rhovol.data)
 	if ndim_x_rhovol==1
-	  gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
+	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
 	else
-	  gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm);
+	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii));
 	end
       end
     end
   else
-    if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
+% $$$     if length(ii)==size(gdat_data.grids_1d.rhopolnorm,1)
+    if length(ii) >0
       if ~isempty(rhotor_norm.x) && ~isempty(rhotor_norm.data)
+        nb_ii = length(ii);
         if ndim_x_rhotor==1
-          gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
         else
-          gdat_data.grids_1d.rhotornorm(:,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+          gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(-13,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
         end
       end
       if ~isempty(rhovol.x) && ~isempty(rhovol.data)
 	if ndim_x_rhovol==1
-	  gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
 	else
-	  gdat_data.grids_1d.rhovolnorm(:,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(:,it));
+	  gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(-13,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it));
 	end
       end
     end
-- 
GitLab