From 41e3662cc791e1d8ee0fd3d868c55fa97213d49c Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Tue, 17 Nov 2015 12:26:12 +0000
Subject: [PATCH] add several extra data_requests
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5275 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/aug_requests_mapping.m | 26 ++-
crpptbx/AUG/gdat_aug.m | 273 +++++++++++++++++++++++------
crpptbx/AUG/geteqdskAUG.m | 40 +++--
crpptbx/gdat_plot.m | 4 +-
4 files changed, 268 insertions(+), 75 deletions(-)
diff --git a/crpptbx/AUG/aug_requests_mapping.m b/crpptbx/AUG/aug_requests_mapping.m
index bbb93eb0..393c495f 100644
--- a/crpptbx/AUG/aug_requests_mapping.m
+++ b/crpptbx/AUG/aug_requests_mapping.m
@@ -175,8 +175,17 @@ switch lower(data_request)
case 'ni'
mapping.method = 'switchcase'; % especially since might have option fit, etc
case 'powers'
+ mapping.timedim = 1;
mapping.label = 'various powers';
mapping.method = 'switchcase';
+ case 'psi_axis'
+ mapping.timedim = 1;
+ mapping.method = 'switchcase'; % there is psi_axis-psi_edge in FPG but otherwise complicated to get from equil, thus needs swticth case
+ mapping.label ='psi_\axis' ;
+ case 'psi_edge'
+ mapping.timedim = 1;
+ mapping.method = 'switchcase'; % is set to zero, so not in tree nodes
+ mapping.label = 'psi\_edge';
case 'q0'
mapping.timedim = 1;
mapping.label = 'q_0';
@@ -219,9 +228,18 @@ switch lower(data_request)
mapping.timedim = 1;
mapping.method = 'signal';
mapping.expression = [{'FPG'},{'Raus'},{'AUGD'}];
+ case 'rhotor'
+ mapping.timedim = 2;
+ mapping.method = 'switchcase';
+ mapping.label = 'rhotor\_norm';
+ case 'rhotor_edge'
+ mapping.timedim = 1;
+ mapping.method = 'switchcase';
+ mapping.label = 'rhotor\_edge';
case 'rhovol'
+ mapping.timedim = 2;
mapping.label = 'rhovol\_norm';
- mapping.method = 'switchcase'; % from conf if exist otherwise computes it
+ mapping.method = 'switchcase';
case 'rmag'
mapping.label = 'R\_magaxis';
mapping.timedim = 1;
@@ -255,8 +273,12 @@ switch lower(data_request)
case 'volume'
mapping.label = 'Volume';
mapping.timedim = 1;
- mapping.method = 'signal';
+ mapping.method = 'switchcase';
mapping.expression = [{'FPG'},{'Vol'},{'AUGD'}];
+ case 'volume_rho'
+ mapping.timedim = 2;
+ mapping.label = 'volume\_norm';
+ mapping.method = 'switchcase';
case 'zeff'
mapping.label = 'zeff from cxrs';
mapping.timedim = 1;
diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m
index 33cf62c0..c382962a 100644
--- a/crpptbx/AUG/gdat_aug.m
+++ b/crpptbx/AUG/gdat_aug.m
@@ -699,31 +699,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'equil'}
- exp_name_eff = gdat_data.gdat_params.exp_name;
- gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
- DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
- M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
- M_Rmesh = M_Rmesh_par.value + 1; % nb of points
- N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
- N_Zmesh = N_Zmesh_par.value + 1; % nb of points
- NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
- NTIME = NTIME_par.value; % nb of points
- Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff);
- % since June, nb of time points in EQ results is not consistent with NTIME and time
- % It seems the first NTIME points are correct, so use this explicitely
- NTIME_Lpf = length(Lpf_par.value);
- if (NTIME < NTIME_Lpf)
- if gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
- elseif (NTIME > NTIME_Lpf)
- if gdat_params.nverbose >= 1
- disp('ERROR: nb of times points LARGER then equil results')
- disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
- end
- return
- end
- Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
- Lpf_SOL = fix(Lpf_tot/100000);
- Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
+ % get equil params and time array in gdat_data.t
+ [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff);
ndimrho = size(qpsi.data,2);
@@ -879,9 +856,6 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.Rcoils=equil_Rcoil.value;
[equil_Zcoil,e]=rdaAUG_eff(shot,DIAG,'Zcl',exp_name_eff);
gdat_data.Zcoils=equil_Zcoil.value;
- % get time values
- [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff);
- gdat_data.t = equil_time.value(1:NTIME);
%
[IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',exp_name_eff);
gdat_data.Ip = IpiPSI.value(1:NTIME);
@@ -1093,32 +1067,99 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- case {'q_rho'}
- exp_name_eff = gdat_data.gdat_params.exp_name;
- gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
- DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
- M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
- M_Rmesh = M_Rmesh_par.value + 1; % nb of points
- N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
- N_Zmesh = N_Zmesh_par.value + 1; % nb of points
- NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
- NTIME = NTIME_par.value; % nb of points
- Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff);
- % since June, nb of time points in EQ results is not consistent with NTIME and time
- % It seems the first NTIME points are correct, so use this explicitely
- NTIME_Lpf = length(Lpf_par.value);
- if (NTIME < NTIME_Lpf)
- if gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
- elseif (NTIME > NTIME_Lpf)
- if gdat_params.nverbose >= 1
- disp('ERROR: nb of times points LARGER then equil results')
- disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
- end
+ case {'powers'}
+ % case total only
+ params_eff = gdat_data.gdat_params;
+ % ohmic, use its time-base
+ params_eff.data_request={'TOT','P_OH'};
+ try
+ gdat_data=gdat_aug(shot,params_eff);
+ catch
+ end
+ gdat_data.ohm.data = gdat_data.data;
+ gdat_data.ohm.dim = gdat_data.dim;
+ gdat_data.ohm.t = gdat_data.t;
+ gdat_data.ohm.label = gdat_data.label;
+ if isempty(gdat_data.data) || isempty(gdat_data.dim)
+ if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end
return
end
- Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
- Lpf_SOL = fix(Lpf_tot/100000);
- Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
+ mapping_for_aug.timedim = 1; mapping_for_aug.gdat_timedim = 1;
+ taus = -10;
+ gdat_data.data = interpos(gdat_data.t,gdat_data.data,5.*taus);
+ gdat_data.data = reshape(gdat_data.data,length(gdat_data.t),1);
+ gdat_data.x =[1];
+ gdat_data.label={'P_{ohmic}'};
+ % nbi
+ params_eff.data_request={'NIS','PNI'};
+ try
+ nbi=gdat_aug(shot,params_eff);
+ catch
+ end
+ if ~isempty(nbi.data) && ~isempty(nbi.dim)
+ gdat_data.nbi.data = nbi.data;
+ gdat_data.nbi.dim = nbi.dim;
+ gdat_data.nbi.t = nbi.t;
+ gdat_data.nbi.label = nbi.label;
+ gdat_data.data(:,end+1) = interpos(gdat_data.nbi.t,gdat_data.nbi.data,gdat_data.t,taus);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{nbi}';
+ else
+ gdat_data.nbi.data = [];
+ gdat_data.nbi.dim = [];
+ gdat_data.nbi.t = [];
+ gdat_data.nbi.label = [];
+ end
+ % ic
+ params_eff.data_request={'ICP','PICRN'};
+ try
+ ic=gdat_aug(shot,params_eff);
+ catch
+ end
+ if ~isempty(ic.data) && ~isempty(ic.dim)
+ gdat_data.ic.data = ic.data;
+ gdat_data.ic.dim = ic.dim;
+ gdat_data.ic.t = ic.t;
+ gdat_data.ic.label = ic.label;
+ gdat_data.data(:,end+1) = interpos(gdat_data.ic.t,gdat_data.ic.data,gdat_data.t,taus);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{ic}';
+ else
+ gdat_data.ic.data = [];
+ gdat_data.ic.dim = [];
+ gdat_data.ic.t = [];
+ gdat_data.ic.label = [];
+ end
+ % ec
+ params_eff.data_request={'ECS','PECRH'};
+ try
+ ec=gdat_aug(shot,params_eff);
+ catch
+ end
+ if ~isempty(ec.data) && ~isempty(ec.dim)
+ gdat_data.ec.data = ec.data;
+ gdat_data.ec.dim = ec.dim;
+ gdat_data.ec.t = ec.t;
+ gdat_data.ec.label = ec.label;
+ gdat_data.data(:,end+1) = interpos(gdat_data.ec.t,gdat_data.ec.data,gdat_data.t,taus);
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ gdat_data.label{end+1}='P_{ec}';
+ else
+ gdat_data.ec.data = [];
+ gdat_data.ec.dim = [];
+ gdat_data.ec.t = [];
+ gdat_data.ec.label = [];
+ end
+ % add tot power
+ gdat_data.data(:,end+1) = interpos(gdat_data.t,nansum(gdat_data.data,2),100*taus);
+ gdat_data.label{end+1}='P_{tot}';
+ gdat_data.x(end+1) =gdat_data.x(end)+1;
+ 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'}
+ [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
% since Lpf depends on time, need to load all first and then loop over time for easier mapping
[qpsi,e]=rdaAUG_eff(shot,DIAG,'Qpsi',exp_name_eff);
ndimrho = size(qpsi.data,2);
@@ -1142,6 +1183,15 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
max_qValue = 30.; % Note could just put a NaN on LCFS value since ill-defined when diverted
for it=1:NTIME
Lpf1 = Lpf1_t(it);
+ % Qpsi and similar data is on (time,radius) with radius being: LCFS..Lpf_points dummy LCFS..SOL part
+ % change it to (radial,time) and use only Lpf+1 points up to LCFS
+ ijok=find(qpsi.value(:,1)); % note: eqr fills in only odd points radially
+ % set NaNs to zeroes
+ if qpsi.value(ijok(1),1)<0
+ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue);
+ else
+ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue);
+ end
% get x values
psi_it=psi_tree.value(it,Lpf1:-1:1)';
gdat_data.psi_axis(it)= psi_it(1);
@@ -1191,8 +1241,6 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
gdat_data.x = gdat_data.rhopolnorm;
% get time values
- [equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff);
- gdat_data.t = equil_time.value(1:NTIME);
gdat_data.data = gdat_data.qvalue; % put q in data
gdat_data.units=[]; % not applicable
gdat_data.data_fullpath = [DIAG ' from expname: ' gdat_data.gdat_params.exp_name '; q(rhopolnorm,t) in .data(.x,.t)'];
@@ -1202,6 +1250,86 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
gdat_data.dimunits{1} = 'rhopolnorm';
gdat_data.dimunits{2} = 'time [s]';
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'psi_axis', 'psi_edge'}
+ if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
+ gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
+ gdat_params_eff.data_request = [{'FPG'},{'fax-bnd'},{'AUGD'}];
+ gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
+ gdat_data.gdat_params = gdat_params_prev;
+ gdat_data.label = 'psi\_axis-psi\_edge';
+ gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only psi\_axis-psi\_edge from FPG/fax-bnd'];
+ else
+ [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
+ % since Lpf depends on time, need to load all first and then loop over time for easier mapping
+ [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff);
+ ndimrho = size(psi_tree.data,2);
+ if ndimrho==NTIME_Lpf
+ % data seems to be transposed
+ ndimrho = size(psi_tree.data,1);
+ itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
+ else
+ itotransposeback = 0;
+ end
+ psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
+ ijnan=find(isnan(psi_tree.value));
+ psi_tree.value(ijnan)=0;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.dimunits{1} = 'time [s]';
+ for it=1:NTIME
+ Lpf1 = Lpf1_t(it);
+ psi_it=psi_tree.value(it,Lpf1:-1:1)';
+ if strcmp(data_request_eff,'psi_axis')
+ gdat_data.data(it)= psi_it(1);
+ elseif strcmp(data_request_eff,'psi_edge')
+ gdat_data.data(it)= psi_it(end);
+ else
+ end
+ end
+ gdat_data.units = psi_tree.units;
+ gdat_data.data_fullpath = [DIAG '/PFL extract ' data_request_eff];
+ end
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ case {'rho_tor', 'rhotor_edge', 'rhovol', 'volume', 'volume_rho'}
+ if strcmp(upper(gdat_data.gdat_params.equil),'FPG')
+ gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev;
+ gdat_params_eff.data_request = [{'FPG'},{'fax-bnd'},{'AUGD'}];
+ gdat_data = gdat_aug(gdat_data.shot,gdat_params_eff);
+ gdat_data.gdat_params = gdat_params_prev;
+ gdat_data.label = 'psi\_axis-psi\_edge';
+ gdat_data.data_fullpath = [gdat_data.data_fullpath ' yields only psi\_axis-psi\_edge from FPG/fax-bnd'];
+ else
+ [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL] = get_EQ_params(gdat_data);
+ % since Lpf depends on time, need to load all first and then loop over time for easier mapping
+ [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff);
+ ndimrho = size(psi_tree.data,2);
+ if ndimrho==NTIME_Lpf
+ % data seems to be transposed
+ ndimrho = size(psi_tree.data,1);
+ itotransposeback = 1; % seems x,time inverted so transpose and exchange .x and .t
+ else
+ itotransposeback = 0;
+ end
+ psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback);
+ ijnan=find(isnan(psi_tree.value));
+ psi_tree.value(ijnan)=0;
+ gdat_data.dim{1} = gdat_data.t;
+ gdat_data.dimunits{1} = 'time [s]';
+ for it=1:NTIME
+ Lpf1 = Lpf1_t(it);
+ psi_it=psi_tree.value(it,Lpf1:-1:1)';
+ if strcmp(data_request_eff,'psi_axis')
+ gdat_data.data(it)= psi_it(1);
+ elseif strcmp(data_request_eff,'psi_edge')
+ gdat_data.data(it)= psi_it(end);
+ else
+ end
+ end
+ gdat_data.units = psi_tree.units;
+ gdat_data.data_fullpath = [DIAG '/PFL extract ' data_request_eff];
+ end
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case {'sxr'}
% sxr from B by default or else if 'source','else' is provided
@@ -1338,3 +1466,40 @@ error_status=0;
return
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% functions for portions called several times
+
+function [gdat_data,exp_name_eff,DIAG,NTIME_Lpf,NTIME,Lpf1_t,Lpf_SOL,M_Rmesh,N_Zmesh] = get_EQ_params(gdat_data);
+% get basic params to be able to read results in EQ-like shotfiles
+% M_Rmesh,N_Zmesh only needed for equil when 2D quantities are required
+%
+shot=gdat_data.shot;
+exp_name_eff = gdat_data.gdat_params.exp_name;
+gdat_data.gdat_params.equil = upper(gdat_data.gdat_params.equil);
+DIAG = gdat_data.gdat_params.equil; % 'EQI' by default at this stage, should become EQH?
+M_Rmesh_par = sf2par(DIAG,shot,'M','PARMV');
+M_Rmesh = M_Rmesh_par.value + 1; % nb of points
+N_Zmesh_par = sf2par(DIAG,shot,'N','PARMV');
+N_Zmesh = N_Zmesh_par.value + 1; % nb of points
+NTIME_par = sf2par(DIAG,shot,'NTIME','PARMV');
+NTIME = NTIME_par.value; % nb of points
+Lpf_par = rdaAUG_eff(shot,DIAG,'Lpf',exp_name_eff);
+% since June, nb of time points in EQ results is not consistent with NTIME and time
+% It seems the first NTIME points are correct, so use this explicitely
+NTIME_Lpf = length(Lpf_par.value);
+if (NTIME < NTIME_Lpf)
+ if gdat_data.gdat_params.nverbose>=3; disp('WARNING: nb of times points smaller then equil results, use first NTIME points'); end
+elseif (NTIME > NTIME_Lpf)
+ if gdat_data.gdat_params.nverbose >= 1
+ disp('ERROR: nb of times points LARGER then equil results')
+ disp('this is unexpected, so stop there and ask Olivier.Sauter@epfl.ch')
+ end
+ return
+end
+Lpf_tot = Lpf_par.value(1:NTIME); % nb of points: 100000*nb_SOL points + nb_core
+Lpf_SOL = fix(Lpf_tot/100000);
+Lpf1_t = mod(Lpf_tot,100000)+1; % nb of Lpf points
+
+[equil_time,e]=rdaAUG_eff(shot,DIAG,'time',exp_name_eff);
+gdat_data.t = equil_time.value(1:NTIME);
diff --git a/crpptbx/AUG/geteqdskAUG.m b/crpptbx/AUG/geteqdskAUG.m
index 7a661a18..ac5731f4 100644
--- a/crpptbx/AUG/geteqdskAUG.m
+++ b/crpptbx/AUG/geteqdskAUG.m
@@ -91,26 +91,32 @@ hold on
psiedge_eff = 0.001*eqdsk.psiaxis + 0.999*eqdsk.psiedge;
[hh1 hh2]=contour(eqdsk.rmesh,eqdsk.zmesh,eqdsk.psi',[psiedge_eff psiedge_eff],'k');
axis equal
-ij=1;
-ij_prev = 0;
-nbhh1(ij)=hh1(2,ij_prev+1);
-Rbnd{ij}=hh1(1,2:1+nbhh1(ij));
-Zbnd{ij}=hh1(2,2:1+nbhh1(ij));
-ij_prev = ij_prev+1+nbhh1(ij);
-while (ij_prev+1<size(hh1,2))
- % next
- ij=ij + 1;
+if ~isempty(hh1)
+ ij=1;
+ ij_prev = 0;
nbhh1(ij)=hh1(2,ij_prev+1);
- Rbnd{ij}=hh1(1,ij_prev+2:ij_prev+1+nbhh1(ij));
- Zbnd{ij}=hh1(2,ij_prev+2:ij_prev+1+nbhh1(ij));
+ Rbnd{ij}=hh1(1,2:1+nbhh1(ij));
+ Zbnd{ij}=hh1(2,2:1+nbhh1(ij));
ij_prev = ij_prev+1+nbhh1(ij);
+ while (ij_prev+1<size(hh1,2))
+ % next
+ ij=ij + 1;
+ nbhh1(ij)=hh1(2,ij_prev+1);
+ Rbnd{ij}=hh1(1,ij_prev+2:ij_prev+1+nbhh1(ij));
+ Zbnd{ij}=hh1(2,ij_prev+2:ij_prev+1+nbhh1(ij));
+ ij_prev = ij_prev+1+nbhh1(ij);
+ end
+ % assume LCFS with most points
+ [zzz irz]=max(nbhh1);
+ eqdsk.nbbound = nbhh1(irz);
+ eqdsk.rplas = Rbnd{irz}';
+ eqdsk.zplas = Zbnd{irz}';
+ plot(eqdsk.rplas,eqdsk.zplas,'k-')
+else
+ eqdsk.nbbound = [];
+ eqdsk.rplas = [];
+ eqdsk.zplas = [];
end
-% assume LCFS with most points
-[zzz irz]=max(nbhh1);
-eqdsk.nbbound = nbhh1(irz);
-eqdsk.rplas = Rbnd{irz}';
-eqdsk.zplas = Zbnd{irz}';
-plot(eqdsk.rplas,eqdsk.zplas,'k-')
[aget]=which('geteqdskAUG');
[path1,name2,ext3]=fileparts(aget);
diff --git a/crpptbx/gdat_plot.m b/crpptbx/gdat_plot.m
index ad09d006..147f1c46 100644
--- a/crpptbx/gdat_plot.m
+++ b/crpptbx/gdat_plot.m
@@ -89,7 +89,7 @@ if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty
imagesc(T+tmhdm(1),F/1e3,20*log10(abs(B)));axis xy;colormap jet;
ylabel('freq')
xlabel(gdat_data.dimunits{1})
- title(gdat_data.label{i})
+ title([upper(gdat_data.gdat_params.machine) '#' num2str(gdat_data.shot) ' ' gdat_data.label{i}])
mhd_sum_data = mhd_sum_data + gdat_data.data(:,i);
end
[B,F,T]=specgram(mhd_sum_data./size(gdat_data.data,2),nfft,1/mean(diff(gdat_data.t)),hanning(nfft),nfft/2);
@@ -97,7 +97,7 @@ if prod(isfield(gdat_data,{'data','t'})) && ~isempty(gdat_data.data) && ~isempty
imagesc(T+tmhdm(1),F/1e3,20*log10(abs(B)));axis xy;colormap jet;
ylabel('freq')
xlabel(gdat_data.dimunits{1})
- title(['sum of ' gdat_data.label{:}])
+ title([upper(gdat_data.gdat_params.machine) '#' num2str(gdat_data.shot) ' sum of ' gdat_data.label{:}])
end
else
disp('cannot plot gdat_data, has empty data or t field')
--
GitLab