diff --git a/crpptbx/JET/gdat_jet.m b/crpptbx/JET/gdat_jet.m index b3caeb64aadb8e7f5da469e4563353634b0a47a7..68b3ff4d01f32a8b088a42afc7595e0286686dec 100644 --- a/crpptbx/JET/gdat_jet.m +++ b/crpptbx/JET/gdat_jet.m @@ -387,129 +387,6 @@ elseif strcmp(mapping_for_jet.method,'switchcase') % First the request names valid for "all" machines: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'a_minor','rgeom'} - % compute average minor or major radius (on z=zaxis normally) - nodenameeff=['\results::r_max_psi' substr_liuqe]; - rmaxpsi=tdi(nodenameeff); - ijnan = find(isnan(rmaxpsi.data)); - if isempty(rmaxpsi.data) || isempty(rmaxpsi.dim) || ischar(rmaxpsi.data) || ... - ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rmaxpsi.data)) ) - if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end - if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end - return - end - nodenameeff2=['\results::r_min_psi' substr_liuqe]; - rminpsi=tdi(nodenameeff2); - ijnan = find(isnan(rminpsi.data)); - if isempty(rminpsi.data) || isempty(rminpsi.dim) || ... - ( ~isempty(ijnan) && prod(size(ijnan))==prod(size(rminpsi.data)) ) - if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff2 ' for data_request= ' data_request_eff]); end - if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end - return - end - ij=find(rmaxpsi.data<0.5 | rmaxpsi.data>1.2); - if ~isempty(ij); rmaxpsi.data(ij)=NaN; end - ij=find(rminpsi.data<0.5 | rminpsi.data>1.2); - if ~isempty(ij); rminpsi.data(ij)=NaN; end - if strcmp(data_request_eff,'a_minor') - gdat_data.data=0.5.*(rmaxpsi.data(end,:) - rminpsi.data(end,:)); - gdat_data.data_fullpath=[nodenameeff ' - ' nodenameeff2 ' /2']; - elseif strcmp(data_request_eff,'rgeom') - gdat_data.data=0.5.*(rmaxpsi.data(end,:) + rminpsi.data(end,:)); - gdat_data.data_fullpath=[nodenameeff ' + ' nodenameeff2 ' /2']; - else - if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end - return - end - gdat_data.dim = rmaxpsi.dim(2); - gdat_data.t = gdat_data.dim{1}; - if any(strcmp(fieldnames(rmaxpsi),'units')) - gdat_data.units = rmaxpsi.units; - end - gdat_data.dimunits = rmaxpsi.dimunits(2); - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'zgeom'} - % compute average minor or major radius (on z=zaxis normally) - nodenameeff=['\results::z_contour' substr_liuqe]; - zcontour=tdi(nodenameeff); - if isempty(zcontour.data) || isempty(zcontour.dim) % || ischar(zcontour.data) (to add?) - if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end - if (gdat_params.nverbose>=3); disp(['rerun LIUQE?']); end - return - end - if strcmp(data_request_eff,'zgeom') - gdat_data.data=0.5.*(max(zcontour.data,[],1) + min(zcontour.data,[],1)); - gdat_data.data_fullpath=['(max+min)/2 of ' nodenameeff]; - gdat_data.dim{1} = zcontour.dim{2}; - gdat_data.dimunits{1} = zcontour.dimunits{2}; - else - if (gdat_params.nverbose>=1); warning(['should not be in this case with data_request_eff = ' data_request_eff]); end - return - end - gdat_data.t = gdat_data.dim{mapping_for_jet.gdat_timedim}; - if any(strcmp(fieldnames(zcontour),'units')) - gdat_data.units = zcontour.units; - end - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'b0'} - % B0 at R0=0.88 - R0EXP=0.88; - if liuqe_version==-1 - nodenameeff = 'jet_eq("BZERO","FBTE")'; - tracetdi=tdi(nodenameeff); - gdat_data.data = tracetdi.data; - else - nodenameeff=['\magnetics::iphi']; - tracetdi=tdi(nodenameeff); - gdat_data.data=192.E-07 * 0.996 *tracetdi.data/R0EXP; - end - if isempty(tracetdi.data) || isempty(tracetdi.dim) % || ischar(tracetdi.data) (to add?) - if (gdat_params.nverbose>=1); warning(['problems loading data for ' nodenameeff ' for data_request= ' data_request_eff]); end - return - end - gdat_data.data_fullpath=[nodenameeff]; - gdat_data.dim = tracetdi.dim; - gdat_data.t = gdat_data.dim{1}; - if any(strcmp(fieldnames(tracetdi),'units')) - gdat_data.units = tracetdi.units; - end - gdat_data.dimunits = tracetdi.dimunits; - gdat_data.request_description = ['vacuum magnetic field at R0=' num2str(R0EXP) 'm; COCOS=17']; - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - case {'betan'} - % 100*beta / |Ip[MA] * B0[T]| * a[m] - % get B0 from gdat_jet, without re-opening the shot and using the same parameters except data_request - % easily done thanks to structure call for options - params_eff = gdat_data.gdat_params; - params_eff.data_request='b0'; - b0=gdat_jet([],params_eff); % note: no need to set .doplot=0 since gdat_jet does not call gdat_plot in any case - params_eff.data_request='ip'; - ip=gdat_jet([],params_eff); - params_eff.data_request='beta'; - beta=gdat_jet([],params_eff); - params_eff.data_request='a_minor'; - a_minor=gdat_jet([],params_eff); - % use beta as time base - if isempty(b0.data) || isempty(b0.dim) || isempty(ip.data) || isempty(ip.dim) || isempty(a_minor.data) || isempty(a_minor.dim) || isempty(beta.data) || isempty(beta.dim) - if (gdat_params.nverbose>=1); warning(['problems loading data for data_request= ' data_request_eff]); end - return - end - gdat_data.dim = beta.dim; - gdat_data.t = beta.dim{1}; - gdat_data.data = beta.data; - ij=find(~isnan(ip.data)); - ip_t = interp1(ip.dim{1}(ij),ip.data(ij),gdat_data.t); - ij=find(~isnan(b0.data)); - b0_t = interp1(b0.dim{1}(ij),b0.data(ij),gdat_data.t); - ij=find(~isnan(a_minor.data)); - a_minor_t = interp1(a_minor.dim{1}(ij),a_minor.data(ij),gdat_data.t); - gdat_data.data = 100.*beta.data ./ abs(ip_t).*1.e6 .* abs(b0_t) .* a_minor_t; - gdat_data.data_fullpath='100*beta/ip*1e6*b0*a_minor, each from gdat_jet'; - gdat_data.units = ''; - gdat_data.dimunits{1} = 's'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'cxrs'} @@ -747,11 +624,18 @@ elseif strcmp(mapping_for_jet.method,'switchcase') end if strcmp(lower(gdat_data.gdat_params.source),'hrt2') nodenameeff_rho = []; % profiles already on rho + gdat_data.gdat_params.fit = 0; % no need, chain2 is already a fit on rhopol end else gdat_data.gdat_params.source = 'hrts'; end - params_eff = gdat_data.gdat_params; + if ~isfield(gdat_data.gdat_params,'fit_ne_edge') || isempty(gdat_data.gdat_params.fit_ne_edge) || ~isnumeric(gdat_data.gdat_params.fit_ne_edge) + gdat_data.gdat_params.fit_ne_edge = 1e18; % default edge ne value (-1 to let free) + end + if ~isfield(gdat_data.gdat_params,'fit_te_edge') || isempty(gdat_data.gdat_params.fit_te_edge) || ~isnumeric(gdat_data.gdat_params.fit_te_edge) + gdat_data.gdat_params.fit_te_edge = 50; % default edge te value (-1 to let free) + end + 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}; @@ -795,9 +679,15 @@ elseif strcmp(mapping_for_jet.method,'switchcase') gdat_data.(data_request_eff).rhopol = gdat_data.x; gdat_data.dimunits{1} = 'rhopol'; end -keyboard - gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose); - + if length(gdat_data.x) == numel(gdat_data.rhopol) + % gdat_data.x is already rhopol and 1D + gdat_data = get_grids_1d(gdat_data,1,1,gdat_params.nverbose); + else + aa=gdat_data; aa.x = aa.rhopol; + gdat_data2 = get_grids_1d(aa,2,1,gdat_params.nverbose); + gdat_data.grids_1d = gdat_data2.grids_1d; + clear aa gdat_data2 + end % defaults for fits, so user always gets std structure gdat_data.fit.rhotornorm = []; % same for both ne and te gdat_data.fit.rhopolnorm = []; @@ -810,7 +700,7 @@ keyboard gdat_data.fit.raw.te.rhotornorm = []; gdat_data.fit.raw.ne.data = []; gdat_data.fit.raw.ne.rhotornorm = []; - fit_tension_default = -0.1; + fit_tension_default = -1; if isfield(gdat_data.gdat_params,'fit_tension') fit_tension = gdat_data.gdat_params.fit_tension; else @@ -853,17 +743,17 @@ keyboard 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]); + % make rhopol->rhotor transformation for each time since equilibrium might have changed + irhook=find(gdat_data.grids_1d.rhopolnorm(:,it)>0 & gdat_data.grids_1d.rhopolnorm(:,it)<1); % no need for ~isnan + [rhoeff isort]=sort(gdat_data.grids_1d.rhopolnorm(irhook,it)); + gdat_data.fit.rhotornorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.grids_1d.rhotornorm(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); + idatate = find(gdat_data.te.data(:,it)>1 & gdat_data.grids_1d.rhopolnorm(:,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.rhopolnorm(idatate,it) = gdat_data.grids_1d.rhopolnorm(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,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(idatate,it)); rhoeff = [0; rhoeff]; teeff = gdat_data.te.data(idatate(irhoeff),it); te_err_eff = gdat_data.te.error_bar(idatate(irhoeff),it); @@ -873,16 +763,22 @@ keyboard % 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); + if gdat_data.gdat_params.fit_te_edge < 0 + [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhopolnorm(:,it)] = interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 0],[0 ... + 0],te_err_eff); + else + [gdat_data.fit.te.data(:,it), gdat_data.fit.te.drhopolnorm(:,it)] = ... + interpos(rhoeff,teeff,rhopolfit,fit_tension.te,[1 2],[0 gdat_data.gdat_params.fit_te_edge],te_err_eff); + end 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); + idatane = find(gdat_data.ne.data(:,it)>1 & gdat_data.grids_1d.rhopolnorm(:,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.rhopolnorm(idatane,it) = gdat_data.grids_1d.rhopolnorm(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,irhoeff] = sort(gdat_data.grids_1d.rhopolnorm(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); @@ -892,7 +788,15 @@ keyboard % 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); + if gdat_data.gdat_params.fit_ne_edge < 0 + [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhopolnorm(:,it)] = interpos(rhoeff,neeff,rhopolfit,fit_tension.ne,[1 0],[0 ... + 0],ne_err_eff); + else + ij_in_1 = find(rhoeff<1); + [gdat_data.fit.ne.data(:,it), gdat_data.fit.ne.drhopolnorm(:,it)] = ... + interpos([rhoeff(ij_in_1); 1.],[neeff(ij_in_1); gdat_data.gdat_params.fit_ne_edge], ... + rhopolfit,fit_tension.ne,[1 2],[0 gdat_data.gdat_params.fit_ne_edge],[ne_err_eff(ij_in_1);ne_err_eff(1)]); + end end end end @@ -1364,6 +1268,7 @@ keyboard rad_bulk2=gdat_jet(shot,params_eff); gdat_data.rad_bulk_u = rad_bulk2.data; gdat_data.rad_bulk_avg = 0.5.*(rad_bulk1.data+rad_bulk2.data); + gdat_data.rad_bulk_t = rad_bulk1.t; end gdat_data.t = rad.t; gdat_data.dim{1} = gdat_data.t; @@ -1396,6 +1301,7 @@ keyboard %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'q_rho'} % q profile on psi from liuqe + nodenameeff=[{'ppf'},{'EFTM'},{'Q'}]; nodenameeff=[{'ppf'},{'EFIT'},{'Q'}]; ppftype = nodenameeff{1}; tracename = [nodenameeff{2} '/' nodenameeff{3}]; diff --git a/crpptbx/JET/get_grids_1d.m b/crpptbx/JET/get_grids_1d.m index 783436477bdd8c5e14195cec6e35dabfa5c5eebe..497eb6c51db33606e8584715d7ce1666d9a1cfe2 100644 --- a/crpptbx/JET/get_grids_1d.m +++ b/crpptbx/JET/get_grids_1d.m @@ -27,7 +27,7 @@ end params_eff = gdat_data.gdat_params;params_eff.doplot=0; params_eff.data_request='rhotor_norm'; -rhotor_norm = gdat(gdat_data.shot,params_eff) +rhotor_norm = gdat(gdat_data.shot,params_eff); params_eff.data_request='rhovol'; rhovol = gdat(gdat_data.shot,params_eff); params_eff.data_request='psi_axis'; diff --git a/crpptbx/JET/jet_help_parameters.m b/crpptbx/JET/jet_help_parameters.m index ee75d00d39b8332c95f07cd56136d6c09579f164..0635755a3be8c0195bbc7384644d3f411683e4d5 100644 --- a/crpptbx/JET/jet_help_parameters.m +++ b/crpptbx/JET/jet_help_parameters.m @@ -27,15 +27,17 @@ help_struct_all = struct(... % $$$ help_struct_all.time_interval = ['if provided sets a specific time interval [tstart tend].' ... % $$$ char(10) 'cxrs: (time_interval can have several nbs) take data and average over time interval(s) only, plots from CXRS_get_profiles are then provided' ... % $$$ ' as well']; -% $$$ help_struct_all.fit_tension = ['smoothing value used in interpos fitting routine, -30 means ''30 times default value'', thus -1 often a' ... -% $$$ ' good value' char(10) ... -% $$$ 'cxrs: if numeric, default for all cases, if structure, default for non given fields']; +help_struct_all.fit_tension = ['smoothing value used in interpos fitting routine, -30 means ''30 times default value'', thus -1 often a' ... + ' good value' char(10) ... + ]; help_struct_all.time = 'time(s) value(s) if relevant, for example eqdsk is provided by default only for time=1.0s'; help_struct_all.zshift = 'vertical shift of equilibrium, either for eqdsk (1 or -99 to shift to zaxis=0) or for mapping measurements on to rho surfaces [m]'; help_struct_all.cocos = ['cocos value desired in output, uses eqdsk_cocos_transform. Note should use latter if a specific Ip and/or B0 sign' ... 'is wanted. See O. Sauter et al Comput. Phys. Commun. 184 (2013) 293']; % $$$ help_struct_all.edge = '0 (default), 1 to get edge Thomson values'; -% $$$ help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type'; +help_struct_all.fit = '0, no fit profiles, 1 (default) if fit profiles desired as well, relevant for _rho profiles. See also fit_type'; +help_struct_all.fit_ne_edge = '-1, no edge value imposed, >0, value at boundary for all times'; +help_struct_all.fit_te_edge = '-1, no edge value imposed, >0, value at boundary for all times'; % $$$ help_struct_all.fit_type = 'provenance of fitted profiles ''conf'' (default) from conf nodes, ''avg'' or ''local'' for resp. proffit: nodes'; help_struct_all.source = 'sxr: ''H'' (default), camera name ''T'', ''V'' ; for Vloop: ''ppf''(default),''jpf'''; help_struct_all.camera = ['[] (default, all), [i1 i2 ...] chord nbs ([1 3 5] if only chords 1, 3 and 5 are desired), ''central'' uses H10']; diff --git a/crpptbx/JET/jet_requests_mapping.m b/crpptbx/JET/jet_requests_mapping.m index 3b579e5ece405661fa3c3044bae5a59d9003f845..f76d7674fa910e4597fc01aaa58ff44c93466ea7 100644 --- a/crpptbx/JET/jet_requests_mapping.m +++ b/crpptbx/JET/jet_requests_mapping.m @@ -126,6 +126,11 @@ switch lower(data_request) mapping.label = 'l_i_3'; mapping.method = 'signal'; mapping.expression = [{'ppf'},{'EFIT'},{'LI3M'}]; + case {'locked', 'locked_mode'} + mapping.timedim = 1; + mapping.label = 'locked mode'; + mapping.method = 'signal'; + mapping.expression = [{'jpf'},{'DA'},{'C2-LOCA'}]; case 'mhd' mapping.timedim = 1; mapping.label = 'n=1 and n=2 amplitude';