From 2dd41eba541f64ae725b4283d861e58b72d83bad Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Wed, 5 Jul 2017 16:42:05 +0000
Subject: [PATCH] add mhd odd/even for JET

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@7777 d63d8f72-b253-0410-a779-e742ad2e26cf
---
 crpptbx/JET/gdat_jet.m             | 268 +++++++++++++++++------------
 crpptbx/JET/jet_requests_mapping.m |  33 ++--
 crpptbx/gdat_plot.m                |   6 +-
 3 files changed, 190 insertions(+), 117 deletions(-)

diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m
index 925ec592..2a0ca4e7 100644
--- a/crpptbx/JET/gdat_jet.m
+++ b/crpptbx/JET/gdat_jet.m
@@ -706,95 +706,52 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'mhd'}
-    % load n=1, 2 and 3 Bdot from magnetic measurements
-    if shot< 50926
-      n1=tdi('abs(mhdmode("LFS",1,1))');
-      n2=tdi('abs(mhdmode("LFS",2,1))');
-      n3=tdi('abs(mhdmode("LFS",3,1))');
-      gdat_data.data_fullpath='abs(mhdmode("LFS",n,1));% for 1, 2, 3';
-    else
-      if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
-        % gdat_data.gdat_params.source;
-      else
-        gdat_data.gdat_params.source = '23';
-      end      
-      if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23')
-        aaLFSz23_sect3=tdi('\atlas::DT196_MHD_001:channel_067');
-        aaLFSz23_sect11=tdi('\atlas::DT196_MHD_001:channel_075');
-        n1 = aaLFSz23_sect3;
-        n1.data = aaLFSz23_sect3.data - aaLFSz23_sect11.data;
-        n2 = aaLFSz23_sect3;
-        n2.data = aaLFSz23_sect3.data + aaLFSz23_sect11.data;
-        n3=n1;
-        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm';
-        if strcmp(gdat_data.gdat_params.source,'23full')
-          % HFS from sec 3 and 11
-          aaHFSz23_sect3=tdi('\atlas::DT196_MHD_002:channel_018');
-          aaHFSz23_sect11=tdi('\atlas::DT196_MHD_002:channel_022');
-          gdat_data.n1_HFS=aaHFSz23_sect3;
-          gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz23_sect11.data;
-          gdat_data.n2_HFS=aaHFSz23_sect3;
-          gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz23_sect11.data;
-          gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_067 -+ \atlas::DT196_MHD_001:channel_075 for n=1,2, LFS_sect_3/11, z=23cm; _HFS' ...
-                    ' same for sector HFS: MHD_002:channel_018-+MHD_002:channel_022'];
-        end
-      elseif strcmp(gdat_data.gdat_params.source(1),'0')
-        aaLFSz0_sect3=tdi('\atlas::DT196_MHD_001:channel_083');
-        aaLFSz0_sect11=tdi('\atlas::DT196_MHD_001:channel_091');
-        n1 = aaLFSz0_sect3;
-        n1.data = aaLFSz0_sect3.data - aaLFSz0_sect11.data;
-        n2 = aaLFSz0_sect3;
-        n2.data = aaLFSz0_sect3.data + aaLFSz0_sect11.data;
-        n3=n1;
-        gdat_data.data_fullpath='\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect_3/11, z=0cm';
-        if strcmp(gdat_data.gdat_params.source,'0full')
-          % sect 11 180deg from sec 3
-          aaHFSz0_sect3=tdi('\atlas::DT196_MHD_002:channel_034');
-          aaHFSz0_sect11=tdi('\atlas::DT196_MHD_002:channel_038');
-          gdat_data.n1_HFS=aaHFSz0_sect11;
-          gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz0_sect11.data;
-          gdat_data.n2_HFS=aaHFSz0_sect11;
-          gdat_data.n2_HFS.data = gdat_data.n2_HFS.data + aaHFSz0_sect11.data;
-          gdat_data.data_fullpath=['\atlas::DT196_MHD_001:channel_083 -+ \atlas::DT196_MHD_001:channel_091 for n=1,2, LFS_sect3/11, z=0cm; _HFS' ...
-                    ' same for HFS: MHD_002:channel_034-+MHD_002:channel_038'];
-          
-        end
-
-      
-      else
-        disp('should not be here in ''mhd'', ask O. Sauter')
-        return
-      end
+    params_eff = gdat_data.gdat_params;
+    if ~isfield(params_eff,'nfft') || isempty(params_eff.nfft)
+      params_eff.nfft = 4096;
     end
-    if ~isempty(n1.data)
-      gdat_data.data(:,1) = reshape(n1.data,length(n1.data),1);
-      if length(n2.data)==length(n1.data); gdat_data.data(:,2) = reshape(n2.data,length(n2.data),1); end
-      if length(n3.data)==length(n1.data); gdat_data.data(:,3) = reshape(n3.data,length(n3.data),1); end
-      gdat_data.dim{1} = n1.dim{1};
-      gdat_data.t = gdat_data.dim{1};
-      gdat_data.dim{2} = [1; 2; 3]; 
-      gdat_data.dimunits{1} = n1.dimunits{1};
-      gdat_data.dimunits{2} = 'n number';
-      gdat_data.units = 'T/s';
-      gdat_data.request_description = 'delta_Bdot from magnetic probes to get n=1, 2 and 3';
-      gdat_data.label = {'n=1','n=2','n=3'}; % can be used in legend(gdat_data.label)
+    gdat_data.gdat_params = params_eff;
+    params_eff.data_request={'JPF','DA','C1M-H305'};
+    gdat_tmp = gdat_jet(shot,params_eff);
+    gdat_data.sig1.data = reshape(gdat_tmp.data,length(gdat_tmp.data),1 );
+    gdat_data.t = gdat_tmp.t;
+    gdat_data.dim{1} = gdat_data.t;
+    gdat_data.dim{2} = [1 2];
+    gdat_data.dimunits = {'s', 'odd/even'};
+    gdat_data.x = gdat_data.dim{2};
+    gdat_data.sig1.data_request = params_eff.data_request;
+    % other coil 180deg apart toroidally
+    params_eff.data_request = {'JPF','DA','C1M-T009'};
+    gdat_tmp=gdat_jet(shot,params_eff);
+    gdat_data.sig2.data = reshape(gdat_tmp.data,length(gdat_tmp.data),1);
+    gdat_data.sig2.data_request=params_eff.data_request;
+    gdat_data.label={'n=odd','n=even'};
+    gdat_data.full_path='n=1 in data(:,1) and .sig1; .sig2';
+    gdat_data.units = 'T/s';
+    gdat_data.label = {'n=odd','n=even'}; % can be used in legend(gdat_data.label)
+    %
+    if ~isempty(gdat_data.sig1.data) && ~isempty(gdat_data.sig2.data)
+      gdat_data.data(:,1) = (gdat_data.sig1.data - gdat_data.sig2.data)./2;
+      gdat_data.data(:,2) = (gdat_data.sig1.data + gdat_data.sig2.data)./2;
     end
-    
+    keyboard
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ne','te'}
     nodenameeff_rho = 'rho';
-    params_eff = gdat_data.gdat_params;
-    if isfield(params_eff,'source') && ~isempty(params_eff.source)
-      if strcmp(lower(params_eff.source),'chain2')
-        params_eff.source = 'hrt2';
+    if ~isfield(gdat_data.gdat_params,'fit') || isempty(gdat_data.gdat_params.fit) || ~isnumeric(gdat_data.gdat_params.fit)
+      gdat_data.gdat_params.fit = 1; % default do fit
+    end
+    if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source)
+      if strcmp(lower(gdat_data.gdat_params.source),'chain2')
+        gdat_data.gdat_params.source = 'hrt2';
       end
-      if strcmp(lower(params_eff.source),'hrt2')
+      if strcmp(lower(gdat_data.gdat_params.source),'hrt2')
         nodenameeff_rho = []; % profiles already on rho
       end
     else
-      params_eff.source = 'hrts';
+      gdat_data.gdat_params.source = 'hrts';
     end
-    gdat_data.gdat_params = params_eff;
+   params_eff = gdat_data.gdat_params;
     params_eff.doplot = 0;
     % ne or Te from Thomson data on raw z mesh vs (z,t)
     params_eff.data_request = {'ppf',params_eff.source,data_request_eff};
@@ -838,6 +795,106 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
       gdat_data.(data_request_eff).rhopol = gdat_data.x;
       gdat_data.dimunits{1} = 'rhopol';
     end
+
+    % defaults for fits, so user always gets std structure
+    gdat_data.fit.rhotornorm = []; % same for both ne and te
+    gdat_data.fit.rhopolnorm = [];
+    gdat_data.fit.t = [];
+    gdat_data.fit.te.data = [];
+    gdat_data.fit.te.drhotornorm = [];
+    gdat_data.fit.ne.data = [];
+    gdat_data.fit.ne.drhotornorm = [];
+    gdat_data.fit.raw.te.data = [];
+    gdat_data.fit.raw.te.rhotornorm = [];
+    gdat_data.fit.raw.ne.data = [];
+    gdat_data.fit.raw.ne.rhotornorm = [];
+    fit_tension_default = -0.1;
+    if isfield(gdat_data.gdat_params,'fit_tension')
+      fit_tension = gdat_data.gdat_params.fit_tension;
+    else
+      fit_tension = fit_tension_default;
+    end
+    if ~isstruct(fit_tension)
+      fit_tension_eff.te = fit_tension;
+      fit_tension_eff.ne = fit_tension;
+      fit_tension = fit_tension_eff;
+    else
+      if ~isfield(fit_tension,'te'); fit_tension.te = fit_tension_default; end
+      if ~isfield(fit_tension,'ne '); fit_tension.ne = fit_tension_default; end
+    end
+    gdat_data.gdat_params.fit_tension = fit_tension;
+    if isfield(gdat_data.gdat_params,'fit_nb_rho_points')
+      fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points;
+    else
+      fit_nb_rho_points = 201;
+    end
+    gdat_data.gdat_params.fit_nb_rho_points = fit_nb_rho_points;
+    %
+    if gdat_data.gdat_params.fit==1
+      % add fits
+      gdat_data.fit.t = gdat_data.t;
+      rhopolfit = linspace(0,1,fit_nb_rho_points)';
+      gdat_data.fit.rhotornorm = NaN*ones(length(rhopolfit),length(gdat_data.fit.t));
+      gdat_data.fit.rhopolnorm = rhopolfit;
+      if any(strfind(data_request_eff,'te'))
+	gdat_data.fit.raw.te.data = NaN*ones(size(gdat_data.te.data));
+	gdat_data.fit.raw.te.error_bar = NaN*ones(size(gdat_data.te.data));
+	gdat_data.fit.raw.te.rhopolnorm = NaN*ones(size(gdat_data.te.data));
+	gdat_data.fit.te.data = gdat_data.fit.rhotornorm;
+	gdat_data.fit.te.drhotornorm = gdat_data.fit.rhotornorm;
+      end
+      if any(strfind(data_request_eff,'ne'))
+	gdat_data.fit.raw.ne.data = NaN*ones(size(gdat_data.ne.data));
+	gdat_data.fit.raw.ne.error_bar = NaN*ones(size(gdat_data.ne.data));
+	gdat_data.fit.raw.ne.rhopolnorm = NaN*ones(size(gdat_data.ne.data));
+	gdat_data.fit.ne.data =gdat_data.fit.rhotornorm;
+	gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhotornorm;
+      end
+      for it=1:length(gdat_data.t)
+	% make rhotor->rhopol transformation for each time since equilibrium might have changed
+	irhook=find(gdat_data.grids_1d.rhotornorm(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<1); % no need for ~isnan
+	[rhoeff isort]=sort(gdat_data.grids_1d.rhotornorm(irhook,it));
+	gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhopolnorm(irhook(isort),it); 1],rhopolfit,-0.1,[2 2],[0 1]);
+	if any(strfind(data_request_eff,'te'))
+	  idatate = find(gdat_data.te.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
+	  if length(idatate)>0
+	    gdat_data.fit.raw.te.rhotornorm(idatate,it) = gdat_data.grids_1d.rhotornorm(idatate,it);
+	    gdat_data.fit.raw.te.data(idatate,it) = gdat_data.te.data(idatate,it);
+	    gdat_data.fit.raw.te.error_bar(idatate,it) = gdat_data.te.error_bar(idatate,it);
+	    [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatate,it));
+	    rhoeff = [0; rhoeff];
+	    teeff = gdat_data.te.data(idatate(irhoeff),it);
+	    te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it);
+	    % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
+	    ij=find(teeff./te_err_eff>10.);
+	    if ~isempty(ij); te_err_eff(ij) = teeff(ij)./0.1; end
+	    %
+	    teeff =  [teeff(1); teeff];
+	    te_err_eff =  [1e4; te_err_eff];
+	    [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhotornorm(:,it)] = interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 0],[0 0],te_err_eff);
+	  end
+	end
+	if any(strfind(data_request_eff,'ne'))
+	  idatane = find(gdat_data.ne.data(:,it)>0 & gdat_data.grids_1d.rhotornorm(:,it)<=1.05);
+	  if length(idatane)>0
+	    gdat_data.fit.raw.ne.rhotornorm(idatane,it) = gdat_data.grids_1d.rhotornorm(idatane,it);
+	    gdat_data.fit.raw.ne.data(idatane,it) = gdat_data.ne.data(idatane,it);
+	    gdat_data.fit.raw.ne.error_bar(idatane,it) = gdat_data.ne.error_bar(idatane,it);
+	    [rhoeff,irhoeff] = sort(gdat_data.grids_1d.rhotornorm(idatane,it));
+	    rhoeff = [0; rhoeff];
+	    % they are some strange low error_bars, so remove these by max(Te/error_bar)<=10; and changing it to large error bar
+	    neeff = gdat_data.ne.data(idatane(irhoeff),it);
+	    ne_err_eff = gdat_data.ne.error_bar(idatane(irhoeff),it);
+	    ij=find(neeff./ne_err_eff>10.);
+	    if ~isempty(ij); ne_err_eff(ij) = neeff(ij)./0.1; end
+	    %
+	    neeff =  [neeff(1); neeff];
+	    ne_err_eff =  [1e21; ne_err_eff];
+	    [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhotornorm(:,it)] = interpos(rhoeff,neeff,rhopolfit,fit_tension.ne,[1 0],[0 0],ne_err_eff);
+	  end
+	end
+      end
+    end
     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'ni','ti'}
@@ -1280,33 +1337,6 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
     gdat_data.dim{2} = gdat_data.x; gdat_data.dimunits{2} = '';
     gdat_data.data_fullpath = 'tot powers from each sources, and total power in .data(:,end)';
 
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-   case {'q_rho'}
-    % q profile on psi from liuqe
-    nodenameeff=[{'ppf'},{'EFIT'},{'Q'}];
-    ppftype = nodenameeff{1};
-    tracename = [nodenameeff{2} '/' nodenameeff{3}];
-    [a,x,t,d,e]=rda_jet(shot,ppftype,tracename);
-    if e==0
-      if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with data_request_eff= ' data_request_eff]); end
-      return
-    end
-    aatmp.data = a; aatmp.t = t; aatmp.x = sqrt(x); % psi to rhopol
-    gdat_data.data = aatmp.data;
-    gdat_data.t = aatmp.t;
-    gdat_data.x = aatmp.x;
-    if isempty(gdat_data.data)
-      return
-    end
-    gdat_data.dim = [{gdat_data.x} {gdat_data.t}];
-    gdat_data.data_fullpath=[nodenameeff ' on rhopol'];
-    gdat_data.dimunits{1} = '';
-    gdat_data.dimunits{2} = 's';
-    gdat_data.units = '';
-    gdat_data.request_description = 'q(rhopol\_norm)';
-    % add grids_1d to have rhotor, etc
-    gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
-    
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'p_rad', 'prad'}
     params_eff = gdat_data.gdat_params;
@@ -1362,6 +1392,32 @@ elseif strcmp(mapping_for_jet.method,'switchcase')
     gdat_data.request_description = '0 since LIUQE construct psi to be zero at LCFS';
 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+   case {'q_rho'}
+    % q profile on psi from liuqe
+    nodenameeff=[{'ppf'},{'EFIT'},{'Q'}];
+    ppftype = nodenameeff{1};
+    tracename = [nodenameeff{2} '/' nodenameeff{3}];
+    [a,x,t,d,e]=rda_jet(shot,ppftype,tracename);
+    if e==0
+      if gdat_params.nverbose>=3; disp(['error after rda_jet in signal with data_request_eff= ' data_request_eff]); end
+      return
+    end
+    aatmp.data = a; aatmp.t = t; aatmp.x = sqrt(x); % psi to rhopol
+    gdat_data.data = aatmp.data;
+    gdat_data.t = aatmp.t;
+    gdat_data.x = aatmp.x;
+    if isempty(gdat_data.data)
+      return
+    end
+    gdat_data.dim = [{gdat_data.x} {gdat_data.t}];
+    gdat_data.data_fullpath=[nodenameeff ' on rhopol'];
+    gdat_data.dimunits{1} = '';
+    gdat_data.dimunits{2} = 's';
+    gdat_data.units = '';
+    gdat_data.request_description = 'q(rhopol\_norm)';
+    % add grids_1d to have rhotor, etc
+    gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose);
+    
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    case {'rhotor', 'rhotor_edge', 'rhotor_norm', 'rhovol', 'volume_rho'}
     % Phi(LCFS) = int(Bphi dSphi), can use Eq.(11) of "Tokamak coordinate conventions: COCOS" paper:
diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m
index 03bde4ee..dd010bf4 100644
--- a/crpptbx/JET/jet_requests_mapping.m
+++ b/crpptbx/JET/jet_requests_mapping.m
@@ -129,16 +129,29 @@ switch lower(data_request)
  case 'mhd'
   mapping.timedim = 1;
   mapping.label = 'n=1 and n=2 amplitude';
-  mapping.method = 'expression'; % {'JPF','DA','C1M-H306'} {'JPF','DA','C1M-T009'} (big) {'JPF','DA','C1M-MIT27'}
-  mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''JPF'',''DA'',''C1H-I301''}; ' ...
-                    'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ...
-		    'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ...
-		    'gdat_tmp.n_1.data = gdat_tmp.data;gdat_tmp.n_1.data_request=params_eff.data_request;' ...
-		    'params_eff.data_request={''JPF'',''DA'',''C1H-I302''};' ...
-		    'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data(:,2)=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ...
-		    'gdat_tmp.n_2.data = gdat_tmp2.data;gdat_tmp.n_2.data_request=params_eff.data_request;gdat_tmp.label={''n=1'',''n=2''};' ...
-		    'gdat_tmp.full_path=''n=1 in data and .n_1; .n_2'';' ...
-		    'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
+  mapping.timedim = 1;
+  mapping.method = 'switchcase';
+% $$$  mapping.method = 'expression'; % {'JPF','DA','C1M-H306'} {'JPF','DA','C1M-T009'} (big) {'JPF','DA','C1M-MIT27'}{'jpf'  'da'  'c1m-h304'} nfft=1024*3;
+% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''JPF'',''DA'',''C1H-I301''}; ' ...
+% $$$                     'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ...
+% $$$ 		    'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ...
+% $$$ 		    'gdat_tmp.n_1.data = gdat_tmp.data;gdat_tmp.n_1.data_request=params_eff.data_request;' ...
+% $$$ 		    'params_eff.data_request={''JPF'',''DA'',''C1H-I302''};' ...
+% $$$ 		    'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data(:,2)=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ...
+% $$$ 		    'gdat_tmp.n_2.data = gdat_tmp2.data;gdat_tmp.n_2.data_request=params_eff.data_request;gdat_tmp.label={''n=1'',''n=2''};' ...
+% $$$ 		    'gdat_tmp.full_path=''n=1 in data and .n_1; .n_2'';' ...
+% $$$ 		    'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;'];
+
+% $$$   mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request={''JPF'',''DA'',''C1M-T009''}; ' ...
+% $$$                     'gdat_tmp=gdat_jet(shot,params_eff);gdat_tmp.data=reshape(gdat_tmp.data,length(gdat_tmp.data),1 );' ...
+% $$$ 		    'gdat_tmp.dim{1}=gdat_tmp.t;gdat_tmp.dim{2}=[1 2];gdat_tmp.x=gdat_tmp.dim{2};' ...
+% $$$ 		    'gdat_tmp.n_1.data = gdat_tmp.data;gdat_tmp.n_1.data_request=params_eff.data_request;' ...
+% $$$ 		    'params_eff.data_request={''JPF'',''DA'',''C1M-H305''};' ...
+% $$$ 		    'gdat_tmp2=gdat_jet(shot,params_eff);gdat_tmp.data(:,2)=reshape(gdat_tmp2.data,length(gdat_tmp2.data),1);' ...
+% $$$ 		    'gdat_tmp.n_2.data = gdat_tmp2.data;gdat_tmp.n_2.data_request=params_eff.data_request;gdat_tmp.label={''n=1'',''n=2''};' ...
+% $$$ 		    'gdat_tmp.full_path=''n=1 in data and .n_1; .n_2'';' ...
+% $$$ 		    'gdat_tmp.gdat_request=''mhd'';gdat_tmp.gdat_params.data_request=gdat_tmp.gdat_request;' ...
+% $$$                     ''];
  case 'mhd_ampl'
   mapping.timedim = 1;
   mapping.label = 'n=1 and n=2 amplitude';
diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m
index 3377affc..edec2ea4 100644
--- a/crpptbx/gdat_plot.m
+++ b/crpptbx/gdat_plot.m
@@ -94,7 +94,11 @@ if all(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty(
     legend(gdat_data.label);
     % add spectrogram per signal
     mhd_sum_data = 0.;
-    nfft=1024;
+    if isfield(gdat_data.gdat_params,'nfft') && ~isempty(gdat_data.gdat_params.nfft)
+      nfft = gdat_data.gdat_params.nfft;
+    else
+      nfft=1024;
+    end
     tmhdm=mean(reshape(gdat_data.t(1:nfft*fix(length(gdat_data.t)/nfft)),nfft,fix(length(gdat_data.t)/nfft)));
     for i=1:size(gdat_data.data,2)
       [B,F,T]=specgram(gdat_data.data(:,i),nfft,1/mean(diff(gdat_data.t)),hanning(nfft),nfft/2);
-- 
GitLab