diff --git a/matlab/AUG/aug_requests_mapping.m b/matlab/AUG/aug_requests_mapping.m index 43a991fbd99a2a8e3d054eec5134ee65c273898c..94e6f6755be28c2bd0a6c9eb1d5a7c40db97787c 100644 --- a/matlab/AUG/aug_requests_mapping.m +++ b/matlab/AUG/aug_requests_mapping.m @@ -276,8 +276,11 @@ switch lower(data_request) mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ... 'gdat_tmp=gdat_aug(shot,params_eff);params_eff.data_request=''n_greenwald'';' ... 'gdat_tmp2=gdat_aug(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);']; + 'if numel(gdat_tmp2.t)>1;tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... + 'else;tmp_data2=gdat_tmp2.data;end;' ... + 'gdat_tmp.ng=gdat_tmp2;gdat_tmp.ng.data_t_nel=tmp_data2;' ... + 'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);' ... + 'gdat_tmp.units='''';gdat_tmp.dimunits{1}=''s'';']; % $$$ case 'ni' % $$$ mapping.method = 'switchcase'; % especially since might have option fit, etc case 'pellet' diff --git a/matlab/CHDF/cdf2mat.m b/matlab/CHDF/cdf2mat.m index c9d7a14e693f33ce81ffa9f00d2d501cb6eb6fcc..f2c1155bd9ac11239d13c77fc09e08f60c2a946f 100644 --- a/matlab/CHDF/cdf2mat.m +++ b/matlab/CHDF/cdf2mat.m @@ -1,4 +1,5 @@ -function cdf2mat_out = cdf2mat(pfname) +function cdf2mat_out = cdf2mat(pfname,varargin) +% cdf2mat_out = cdf2mat(pfname,varargin) % % reads all variables and coordinates from netcdf % some trials with netcdf, using 50725c01.cdf as test @@ -7,6 +8,9 @@ function cdf2mat_out = cdf2mat(pfname) % % Uses matlab netcdf.open, ncinfo, netcdf.inqVarIDs, netcd.getVar, etc % +% varargin{1}: shot number to add (in particular when cdf file does not have it defined +% shot from file: allinfo=ncinfo(pfname);if ~isempty(allinfo.Attributes),top_attr_names = {allinfo.Attributes(:).Name};allinfo.Attributes(strmatch('shot',top_attr_names,'exact')).Value, else, disp('no Attributes'),end +% % if ~exist(pfname,'file') @@ -33,6 +37,7 @@ if length(allvarnames) ~= length(allvarids) end [varnames_sorted,~]=sort(allvarnames); +cdf2mat_out.allvarnames_sorted = varnames_sorted; % to find a variable: % strmatch('GFUN',allvarnames,'exact') @@ -56,6 +61,9 @@ for i=1:length(coordnames_sorted) for ij=1:length(fields_variables_to_copy) matcdf.coords(i).(fields_variables_to_copy{ij}) = allinfo.Variables(matcdf.coords(i).index_allvarnames).(fields_variables_to_copy{ij}); end + % define defaults before fetching the values in Attributes (if provided) + matcdf.coords(i).units = ''; + matcdf.coords(i).long_name = [matcdf.coords(i).name ' empty']; for jj=1:numel(allinfo.Variables(matcdf.coords(i).index_allvarnames).Attributes) if strcmp(lower(allinfo.Variables(matcdf.coords(i).index_allvarnames).Attributes(jj).Name),'units') matcdf.coords(i).units = strtrim(allinfo.Variables(matcdf.coords(i).index_allvarnames).Attributes(jj).Value); @@ -80,6 +88,7 @@ for i=1:length(varnames_sorted) matcdf.vars(i).name = varnames_sorted{i}; matcdf.vars(i).index_allvarnames = strmatch(matcdf.vars(i).name,allvarnames,'exact'); matcdf.vars(i).varid = allvarids(matcdf.vars(i).index_allvarnames); + % if i==strmatch('ZEFFC',varnames_sorted,'exact'), keyboard; end if ~isempty(matcdf.vars(i).index_allvarnames) if strcmp(allinfo.Variables(matcdf.vars(i).index_allvarnames).Datatype,'single') matcdf.vars(i).data = netcdf.getVar(funnetcdf,matcdf.vars(i).varid,'double'); @@ -89,6 +98,9 @@ for i=1:length(varnames_sorted) for ij=1:length(fields_variables_to_copy) matcdf.vars(i).(fields_variables_to_copy{ij}) = allinfo.Variables(matcdf.vars(i).index_allvarnames).(fields_variables_to_copy{ij}); end + % define defaults before fetching the values in Attributes (if provided) + matcdf.vars(i).units = ''; + matcdf.vars(i).long_name = [matcdf.vars(i).name ' empty']; for jj=1:numel(allinfo.Variables(matcdf.vars(i).index_allvarnames).Attributes) if strcmp(lower(allinfo.Variables(matcdf.vars(i).index_allvarnames).Attributes(jj).Name),'units') matcdf.vars(i).units = strtrim(allinfo.Variables(matcdf.vars(i).index_allvarnames).Attributes(jj).Value); @@ -126,9 +138,17 @@ end if ~isempty(allinfo.Attributes) top_attr_names = {allinfo.Attributes(:).Name}; ij = strmatch('shot',top_attr_names,'exact'); +else + ij = []; +end +if ~isempty(ij) cdf2mat_out.shot = allinfo.Attributes(ij).Value; else - cdf2mat_out.shot = NaN; + if nargin>1 && isnumeric(varargin{1}) + cdf2mat_out.shot = varargin{1}; + else + cdf2mat_out.shot = NaN; + end end cdf2mat_out.fname = pfname; [a1,a2,a3]=fileparts(pfname); diff --git a/matlab/D3D/d3d_requests_mapping.m b/matlab/D3D/d3d_requests_mapping.m index 84879d6ffc32a7347a6dabd35dd18df55c10334a..a40639a32184a880c4c2f6ec16157a9cfd1ad119 100644 --- a/matlab/D3D/d3d_requests_mapping.m +++ b/matlab/D3D/d3d_requests_mapping.m @@ -224,6 +224,11 @@ switch lower(data_request) mapping.label = 'NEBAR\_R0'; mapping.method = 'signal'; mapping.expression = [{'efit01'},{'\NEBAR_R0'}]; + case 'nel_v1' + mapping.timedim = 1; + mapping.label = 'NEBAR\_V1'; + mapping.method = 'signal'; + mapping.expression = [{'efit01'},{'\NEBAR_V1'}]; case {'ne', 'ne_rho'} mapping.timedim = 2; mapping.label = 'ne'; @@ -254,7 +259,9 @@ switch lower(data_request) mapping.expression = ['params_eff = gdat_data.gdat_params;params_eff.data_request=''nel'';' ... 'gdat_tmp=gdat_d3d(shot,params_eff);params_eff.data_request=''n_greenwald'';' ... 'gdat_tmp2=gdat_d3d(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);' ... + 'if numel(gdat_tmp2.t)>1;tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... + 'else;tmp_data2=gdat_tmp2.data;end;' ... + 'gdat_tmp.ng=gdat_tmp2;gdat_tmp.ng.data_t_nel=tmp_data2;' ... 'gdat_tmp.data = gdat_tmp.data*1e6./(tmp_data2+1e-5);' ... 'gdat_tmp.units='''';gdat_tmp.dimunits{1}=''s'';']; case {'f2b', 'pcf2b'} diff --git a/matlab/D3D/gdat_d3d.m b/matlab/D3D/gdat_d3d.m index c3b5c90df0453e60fb371d841279af6b0afd3fe1..dd7414586065ab67e4cd6f2ce9d7f9cefeb96c1f 100644 --- a/matlab/D3D/gdat_d3d.m +++ b/matlab/D3D/gdat_d3d.m @@ -49,6 +49,7 @@ function [gdat_data,gdat_params,error_status,varargout] = gdat_d3d(shot,data_req % a5=gdat(32827,'ip'); % standard call % a6=gdat(32827,'ip','Opt1',123,'Doplot',1,'opt2','Abc'); % standard call with a few options (note all lowercase in output) % +% From local: mdsconnect('atlas.gat.com'); mdsvalue('1+3') % From remote: % ssh -p 2039 -C -l sautero -L 8002:atlas.gat.com:8000 cybele.gat.com % And in another session in matlab: @@ -715,6 +716,12 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'equil'} + if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source) && ischar(gdat_data.gdat_params.source) && strcmp(lower(gdat_data.gdat_params.equil),'efit01') + % assume source provided instead of equil for choice of equil, since equil=default and source provided + warning('choice of equil should be provided in ''equil'' not in ''source'', moved entry') + gdat_data.gdat_params.equil = gdat_data.gdat_params.source; + gdat_data.gdat_params = rmfield(gdat_data.gdat_params,'source'); + end if ~isfield(gdat_data.gdat_params,'equil') || isempty(gdat_data.gdat_params.equil) || ~ischar(gdat_data.gdat_params.equil) gdat_data.gdat_params.equil = 'efit03'; else @@ -844,9 +851,9 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') if any(strcmp(fieldnames(a),'units')) gdat_data.(extra_source{i}).units=a.units; end - gdat_data.(extra_source{i}).r = mdsdata([main_source 'r']); - gdat_data.(extra_source{i}).z = mdsdata([main_source 'z']); - gdat_data.(extra_source{i}).error_bar = mdsdata([nodenameeff '_e'])'; + gdat_data.(extra_source{i}).r = mdsvalue([main_source 'r']); + gdat_data.(extra_source{i}).z = mdsvalue([main_source 'z']); + gdat_data.(extra_source{i}).error_bar = mdsvalue([nodenameeff '_e'])'; gdat_data.(extra_source{i}).data_fullpath=[data_request_eff 'from electrons ' nodenameeff]; gdat_data.(extra_source{i}).dim=[{gdat_data.x};{gdat_data.t}]; if strcmp(extra_source{i},'tangential') @@ -875,7 +882,10 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'neint_ts'} if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) || ~isstruct(gdat_data.gdat_params.source) - error('requires ''ne_rho'' gdat structure as input in ''source'' option'); + warning('requires ''ne_rho'' gdat structure as input in ''source'' option. Loading it...'); + params_eff = gdat_data.gdat_params; + params_eff.data_request = 'ne_rho'; + [gdat_data.gdat_params.source,params_nerho,error_status]=gdat_d3d(shot,params_eff); end ne_rho = gdat_data.gdat_params.source; % Computes R0 and V1 line integrated signal based on ne fit (so inside LCFS), provide also R(psi_ts) and Z(psi_ts) on these lines @@ -1049,11 +1059,16 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') return end % add rho coordinates - params_eff.data_request='equil'; - [equil,params_equil,error_status]=gdat_d3d(shot,params_eff); - if error_status>0 - if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end - return + if ~isfield(gdat_data.gdat_params,'equil') || isempty(gdat_data.gdat_params.equil) || ~isstruct(gdat_data.gdat_params.equil) + params_eff.data_request='equil'; + [equil,params_equil,error_status]=gdat_d3d(shot,params_eff); + if error_status>0 + if gdat_params.nverbose>=3; disp(['problems with ' params_eff.data_request]); end + return + end + else + params_equil.equil = gdat_data.gdat_params.equil.gdat_params.equil; + equil = gdat_data.gdat_params.equil; end gdat_data.gdat_params.equil = params_equil.equil; gdat_data.equil = equil; @@ -1071,24 +1086,26 @@ elseif strcmp(mapping_for_d3d.method,'switchcase') rhotornorm_out_source = NaN*ones(inb_chord_source,inb_time_source); rhovolnorm_out_source = NaN*ones(inb_chord_source,inb_time_source); % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)] - time_equil=[min(gdat_data.(extra_source{i}).t(1)-0.1,equil.t(1)-0.1); 0.5.*(equil.t(1:end-1)+equil.t(2:end)) ;max(equil.t(end)+0.1,gdat_data.(extra_source{i}).t(end)+0.1)]; - clear psi_out_core rhopolnorm_out_core rhotornorm_out_core rhovolnorm_out_core - for itequil=1:numel(time_equil)-1 - rr=equil.Rmesh(:,itequil); - zz=equil.Zmesh(:,itequil); - psirz_in = equil.psi2D(:,:,itequil); - it_core_inequil = find(gdat_data.(extra_source{i}).t>time_equil(itequil) & gdat_data.(extra_source{i}).t<=time_equil(itequil+1)); - if ~isempty(it_core_inequil) - rout_core=gdat_data.(extra_source{i}).r; - zout_core=gdat_data.(extra_source{i}).z; - psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core); - psi_out_core(:,it_core_inequil) = repmat(psi_at_routzout,1,numel(it_core_inequil)); - rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); - for it_cx=1:length(it_core_inequil) - rhotornorm_out_core(:,it_core_inequil(it_cx)) = ... - interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); - rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ... - interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); + if ~isempty(gdat_data.(extra_source{i}).t) + time_equil=[min(gdat_data.(extra_source{i}).t(1)-0.1,equil.t(1)-0.1); 0.5.*(equil.t(1:end-1)+equil.t(2:end)) ;max(equil.t(end)+0.1,gdat_data.(extra_source{i}).t(end)+0.1)]; + clear psi_out_core rhopolnorm_out_core rhotornorm_out_core rhovolnorm_out_core + for itequil=1:numel(time_equil)-1 + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_core_inequil = find(gdat_data.(extra_source{i}).t>time_equil(itequil) & gdat_data.(extra_source{i}).t<=time_equil(itequil+1)); + if ~isempty(it_core_inequil) + rout_core=gdat_data.(extra_source{i}).r; + zout_core=gdat_data.(extra_source{i}).z; + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core); + psi_out_core(:,it_core_inequil) = repmat(psi_at_routzout,1,numel(it_core_inequil)); + rhopolnorm_out_core(:,it_core_inequil) = sqrt((psi_out_core(:,it_core_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_core_inequil) + rhotornorm_out_core(:,it_core_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out_core(:,it_core_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_core(:,it_core_inequil(it_cx)),-3,[2 2],[0 1]); + end end end end diff --git a/matlab/IMAS/gdat_imas.m b/matlab/IMAS/gdat_imas.m index 88b72b2f647f2f8ca0bb8b09601439ac62e11e9d..d2b06fefcbdf8d2fcc5fc3b3289cc933cdc28ddd 100644 --- a/matlab/IMAS/gdat_imas.m +++ b/matlab/IMAS/gdat_imas.m @@ -633,7 +633,7 @@ end function [gdat_data] = get_thomson_raw_data(shot,data_request_eff,gdat_data,doedge,nverbose); % try - time=mdsdata('\results::thomson:times'); + time=mdsvalue('\results::thomson:times'); catch gdat_data.error_bar = []; if strcmp(data_request_eff(1:2),'ne') @@ -695,7 +695,7 @@ if strcmp(data_request_eff(1:2),'ne') disp('***********************************************************************') end end -z=mdsdata(['\diagz::thomson_set_up' edge_str_dot ':vertical_pos']); +z=mdsvalue(['\diagz::thomson_set_up' edge_str_dot ':vertical_pos']); gdat_data.dim=[{z};{time}]; gdat_data.dimunits=[{'Z [m]'} ; {'time [s]'}]; gdat_data.x=z; @@ -716,7 +716,7 @@ function [psiscatvol,psi_max] = psi_scatvol_zshift(gdat_data, zshift, system_str % Use psitbx t_th = gdat_data.t; -th_z = mdsdata(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}? +th_z = mdsvalue(['dim_of(\thomson' system_str_dot ':te,1)']); % TODO: Isn't this gdat_data.dim{1}? th_r = 0.9*ones(size(th_z)); ntt = numel(t_th); nch = numel(th_z); % number of chords diff --git a/matlab/TCV/gdat_tcv.m b/matlab/TCV/gdat_tcv.m index 43a7b34b78067a271255050f71e197dc3ec2a5fa..8acfc6f38c78d0457b3e330c4ed1c09b7b2be85e 100644 --- a/matlab/TCV/gdat_tcv.m +++ b/matlab/TCV/gdat_tcv.m @@ -866,6 +866,94 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.gdat_params.time_out = []; end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + case {'bfields'} + params_eff = gdat_data.gdat_params; + + try + if params_eff.liuqe == 1 + % load LIUQE data + [L,LY] = mds2meq(gdat_data.shot,'LIUQE.M',[],'ifield',true); % flag to allow field computation + LY = meqreprocess(L,LY); % compute fields etc, since not directly loaded from mds2meq + datapath = 'Outputs from [L,LY]=mds2meq(gdat_data.shot,''LIUQE.M'',[],''ifield'',true); LY = meqreprocess(L,LY);'; + Br = permute(LY.Brx,[2,1,3]); % permute for (Z,R)->(R,Z) + Bz = permute(LY.Bzx,[2,1,3]); + Btor = permute(LY.Btx,[2,1,3]); + dim_in = {L.rx,L.zx,LY.t}; + time_in = LY.t; + x_in = {L.rx,L.zx}; + else % compute B-fields from psi + ip_gdat = gdat(gdat_data.shot,'ip'); + [B,coord,time,psi,fsd] = tcv_mag_field (gdat_data.shot,ip_gdat.t,[],[],[],psitbx_str); + datapath ='Outputs from [B,coord,time,psi,fsd] = tcv_mag_field(shot,time,[],[],[],psitbx_str);'; + Br = B.R; + Bz = B.z; + Btor = B.phi; + dim_in = {coord.R,coord.z,time}; + time_in = time; + x_in = {coord.R,coord.z}; + end + + gdat_data.data_fullpath = datapath; + + % radial magnetic field + gdat_data.Br.data = Br; + gdat_data.Br.units = 'T'; + gdat_data.Br.dim = dim_in; + gdat_data.Br.dimunits = {'m','m','s'}; + gdat_data.Br.t = time_in; + gdat_data.Br.x = x_in; + gdat_data.Br.data_fullpath = datapath; + gdat_data.Br.label = 'Radial magnetic field map in (R,Z)'; + + % vertical magnetic field + gdat_data.Bz.data = Bz; + gdat_data.Bz.units = 'T'; + gdat_data.Bz.dim = dim_in; + gdat_data.Bz.dimunits = {'m','m','s'}; + gdat_data.Bz.t = time_in; + gdat_data.Bz.x = x_in; + gdat_data.Bz.data_fullpath = datapath; + gdat_data.Bz.label = 'Vertical magnetic field map in (R,Z)'; + + % toroidal magnetic field + gdat_data.Btor.data = Btor; + gdat_data.Btor.units = 'T'; + gdat_data.Btor.dim = dim_in; + gdat_data.Btor.dimunits = {'m','m','s'}; + gdat_data.Btor.t = time_in; + gdat_data.Btor.x = x_in; + gdat_data.Btor.data_fullpath = datapath; + gdat_data.Btor.label = 'Toroidal magnetic field map in (R,Z)'; + + % total magnetic field + gdat_data.Btot.data = sqrt(gdat_data.Btor.data.^2 + gdat_data.Bz.data.^2 + gdat_data.Br.data.^2); + gdat_data.Btot.units = 'T'; + gdat_data.Btot.dim = dim_in; + gdat_data.Btot.dimunits = {'m','m','s'}; + gdat_data.Btot.t = time_in; + gdat_data.Btot.x = x_in; + gdat_data.Btot.data_fullpath = datapath; + gdat_data.Btot.label = 'Total magnetic field map in (R,Z)'; + + gdat_data.data = gdat_data.Btot.data; + gdat_data.units = 'T'; + gdat_data.dim = dim_in; + gdat_data.dimunits = {'m','m','s'}; + gdat_data.t = time_in; + gdat_data.x = x_in; + gdat_data.label = 'Total magnetic field map in (R,Z)'; + + catch + + warning('Problem obtaining B-fields, check if requested nodes are filled.') + gdat_data.Btor = []; + gdat_data.Br = []; + gdat_data.Bz = []; + gdat_data.Btot = []; + + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% case {'betan', 'beta_tor_norm'} % 100*beta / |Ip[MA] * B0[T]| * a[m] @@ -2076,33 +2164,37 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') n3=tdi('abs(mhdmode("LFS",3,1))'); gdat_data.data_fullpath='abs(mhdmode("LFS",n,1));% for 1, 2, 3'; else + % handle source input if isfield(gdat_data.gdat_params,'source') && ~isempty(gdat_data.gdat_params.source) if ~any(contains(mhd_source_list,lower(gdat_data.gdat_params.source))) t1=sprintf('bad source for mhd: ''%s'', should be either ',gdat_data.gdat_params.source); t2=sprintf('''%s'', ',mhd_source_list{1:end-1}); t3=sprintf('or ''%s''',mhd_source_list{end}); - disp(sprintf('%s %s %s',t1,t2,t3)); + fprintf('%s %s %s\n',t1,t2,t3); return end else + % if no source input decide source based on average z_axis z_axis=gdat_tcv([],'z_axis'); z_axis_av = 0.; - if numel(z_axis.data > 10) && isnumeric(z_axis.data) - z_axis_av = nanmean(z_axis.data([round(numel(z_axis.data)/3):round(numel(z_axis.data)*0.85)])); + if numel(z_axis.data) > 10 && isnumeric(z_axis.data) + z_axis_av = nanmean(z_axis.data(round(numel(z_axis.data)/3):round(numel(z_axis.data)*0.85))); end if z_axis_av > 0.12 - gdat_data.gdat_params.source = '23'; + gdat_data.gdat_params.source = '23'; % probe array at z = +23cm elseif z_axis_av < -0.12 - gdat_data.gdat_params.source = '-23'; + gdat_data.gdat_params.source = '-23'; % probe array at z = -23cm else - gdat_data.gdat_params.source = '0'; + gdat_data.gdat_params.source = '0'; % probe array at z = 0cm end t1=sprintf('source set to ''%s'', can be ',gdat_data.gdat_params.source); t2=sprintf('''%s'', ',mhd_source_list{1:end-1}); t3=sprintf('or ''%s''',mhd_source_list{end}); - disp(sprintf('%s %s %s',t1,t2,t3)) + fprintf('%s %s %s',t1,t2,t3) end - if length(gdat_data.gdat_params.source)>=2 && strcmp(gdat_data.gdat_params.source(1:2),'23') + + % load data dependent on requested source + if numel(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; @@ -2135,14 +2227,14 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % 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=aaHFSz0_sect3; gdat_data.n1_HFS.data = gdat_data.n1_HFS.data - aaHFSz0_sect11.data; - gdat_data.n2_HFS=aaHFSz0_sect11; + gdat_data.n2_HFS=aaHFSz0_sect3; 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 - elseif length(gdat_data.gdat_params.source)>=3 && strcmp(gdat_data.gdat_params.source(1:3),'-23') + elseif numel(gdat_data.gdat_params.source)>=3 && strcmp(gdat_data.gdat_params.source(1:3),'-23') aaLFSzm23_sect3=tdi('\atlas::DT196_MHD_002:channel_003'); aaLFSzm23_sect3.data = - aaLFSzm23_sect3.data; % opposite polarity aaLFSzm23_sect11=tdi('\atlas::DT196_MHD_002:channel_011'); @@ -2153,7 +2245,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') % n3=n1; gdat_data.data_fullpath=['\atlas::DT196_MHD_002:channel_003 -+ \atlas::DT196_MHD_002:channel_011 for n=1,2,' ... 'LFS_sect_3/11, z=-23cm']; - elseif strcmp(lower(gdat_data.gdat_params.source(1:4)),'ltcc') + elseif strcmpi(gdat_data.gdat_params.source(1:4),'ltcc') switch lower(gdat_data.gdat_params.source) case {'ltcc', 'ltcc-dbpol'} % since other probes measure dbpol make "ltcc" default to ltcc-dbpol @@ -2166,8 +2258,8 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') aaLTCC=tdi('\atlas::dt132_ltcc_001:channel_003'); gdat_data.data_fullpath='ltcc-dbtor: \atlas::dt132_ltcc_001:channel_003'; otherwise - error(sprintf('case %s not known, choose between ltcc-dbpol(=ltcc), ltcc-dbrad or ltcc-dbtor', ... - lower(gdat_data.gdat_params.source))); + error('case %s not known, choose between ltcc-dbpol(=ltcc), ltcc-dbrad or ltcc-dbtor', ... + lower(gdat_data.gdat_params.source)); end n1 = aaLTCC; n2.data = []; @@ -2233,7 +2325,7 @@ elseif strcmp(mapping_for_tcv.method,'switchcase') gdat_data.units = 'particles/m^3'; nel.dim{1} = neft.fit.t'; nel.data = ones(size(nel.dim{1})); - gdat_data.r_center_chord = 0.903; + gdat_data.r_center_chord = 0.903; % mdsvalue('\diagz::fir_array:radii[5]') (6th chord is central) [ratio,fir_times,z_lineint]=fir_ratio_from_fits(neft.fit.data,neft.fit.x,psitbx01,nel,gdat_data.r_center_chord,gdat_data.t',1,1e3); gdat_data.data = reshape(1./ratio./z_lineint,numel(ratio),1); gdat_data.data_fullpath = ['line-averaged density from fit using ' neft.fit.data_fullpath]; diff --git a/matlab/TCV/tcv_requests_mapping.m b/matlab/TCV/tcv_requests_mapping.m index ce8a01441a5d585848dde445379a870f69c3b106..cbae8c287f7114b0b333059ef5fa974a874b3cb8 100644 --- a/matlab/TCV/tcv_requests_mapping.m +++ b/matlab/TCV/tcv_requests_mapping.m @@ -69,6 +69,11 @@ switch lower(data_request) mapping.label = 'B_0'; mapping.method = 'switchcase'; mapping.expression = ''; + case 'bfields' + mapping.timedim = 1; + mapping.label = 'bfields'; + mapping.method = 'switchcase'; + mapping.expression = ''; case {'beta', 'beta_tor'} mapping.timedim = 1; mapping.label = '\beta'; @@ -251,11 +256,21 @@ switch lower(data_request) mapping.label = 'line integrated el. density'; mapping.method = 'tdi'; mapping.expression = '\results::fir:lin_int_dens'; + case 'neints' + mapping.timedim = 1; + mapping.label = 'Array line integrated el. density'; + mapping.method = 'tdi'; + mapping.expression = '\tcv_shot::top.results.fir.fir_array:lin_int_dens'; case 'nel' mapping.timedim = 1; mapping.label = 'line-averaged el. density'; mapping.method = 'switchcase'; mapping.expression = '\results::fir:n_average'; + case 'nels' + mapping.timedim = 1; + mapping.label = 'Array line-averaged el. density'; + mapping.method = 'tdi'; + mapping.expression = '\tcv_shot::top.results.fir.fir_array:n_average'; case 'ne_rho' mapping.timedim = 2; mapping.label = 'ne'; @@ -297,6 +312,7 @@ switch lower(data_request) 'gdat_tmp2=gdat_tcv(shot,params_eff);ij=find(gdat_tmp2.data==0);gdat_tmp2.data(ij)=NaN;' ... 'if numel(gdat_tmp2.t)>1;tmp_data2=interp1(gdat_tmp2.t,gdat_tmp2.data,gdat_tmp.t,[],NaN);' ... 'else;tmp_data2=gdat_tmp2.data;end;' ... + 'gdat_tmp.ng=gdat_tmp2;gdat_tmp.ng.data_t_nel=tmp_data2;' ... 'gdat_tmp.data = gdat_tmp.data./(tmp_data2+1e-5);' ... 'gdat_tmp.units='''';gdat_tmp.dimunits{1}=''s'';']; case {'ec_data', 'aux', 'h_cd', 'nbi_data', 'ic_data', 'lh_data', 'ohm_data', 'bs_data'} @@ -442,6 +458,11 @@ switch lower(data_request) mapping.expression = '\results::r_axis'; mapping.expression = '\tcv_shot::top.results.equil_1.results:r_axis'; mapping.expression = 'tcv_eq("r_axis","LIUQE.M")'; + case {'r_xpts', 'rxpts', 'r_xpt'} + mapping.timedim = 2; + mapping.label = 'R\_Xpts'; + mapping.method = 'tdiliuqe'; + mapping.expression = 'tcv_eq("r_xpts","LIUQE.M")'; case {'rtc','scd'} if strcmp('scd',lower(data_request)); error('********* should call with request****** ''rtc'''); end mapping.timedim = 1; @@ -574,6 +595,11 @@ switch lower(data_request) mapping.expression = '\results::z_axis'; mapping.expression = '\tcv_shot::top.results.equil_1.results:z_axis'; mapping.expression = 'tcv_eq("z_axis","LIUQE.M")'; + case {'z_xpts', 'zxpts', 'z_xpt'} + mapping.timedim = 2; + mapping.label = 'Z\_Xpts'; + mapping.method = 'tdiliuqe'; + mapping.expression = 'tcv_eq("z_xpts","LIUQE.M")'; % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % extra TCV cases (not necessarily in official data_request name list) diff --git a/matlab/TCV_IMAS/tcv_get_ids_ec_launchers.m b/matlab/TCV_IMAS/tcv_get_ids_ec_launchers.m index e6f3b040e8ce238e03ed89f3baaeb51aebaa838a..56a9a5008965f260f6eb21fbc8ac27d5f3bc67aa 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_ec_launchers.m +++ b/matlab/TCV_IMAS/tcv_get_ids_ec_launchers.m @@ -76,9 +76,6 @@ for iant=1:nb_launchers % non time-dependent quantities, take 1st ok values ids_ec_launchers.beam{iant}.frequency.time = [pow.ec.t(1) pow.ec.t(end)]; ids_ec_launchers.beam{iant}.frequency.data =[launch_params{iant}{it_ok{iant}(1)}.freq launch_params{iant}{it_ok{iant}(end)}.freq]; - % ids_ec_launchers.beam{iant}.mode.time = [pow.ec.t(1) pow.ec.t(end)]; - % ids_ec_launchers.beam{iant}.mode.data = [-1 -1]; % at this stage assume X mode always, to change when available - % ids_ec_launchers.beam{iant}.time = [pow.ec.t(1) pow.ec.t(end)]; ids_ec_launchers.beam{iant}.mode = -1; % at this stage assume X mode always, to change when available for i=1:length(it_ok{iant}) r0 = sqrt(launch_params{iant}{it_ok{iant}(i)}.x0.^2 + launch_params{iant}{it_ok{iant}(i)}.y0.^2) / 100.; % in [m] @@ -87,15 +84,17 @@ for iant=1:nb_launchers ids_ec_launchers.beam{iant}.launching_position.z(i) = launch_params{iant}{it_ok{iant}(i)}.z0/100.; ids_ec_launchers.beam{iant}.launching_position.phi(i) = atan2(launch_params{iant}{it_ok{iant}(i)}.y0/100,r0); ids_ec_launchers.beam{iant}.time(i) = time_launch; - kz = cos(launch_params{iant}{it_ok{iant}(i)}.theta_toray*pi/180.); - kmr = -sin(launch_params{iant}{it_ok{iant}(i)}.theta_toray*pi/180.).*cos(launch_params{iant}{it_ok{iant}(i)}.phi_toray*pi/180.); - kphi = sin(launch_params{iant}{it_ok{iant}(i)}.theta_toray*pi/180.).*sin(launch_params{iant}{it_ok{iant}(i)}.phi_toray*pi/180.); %*sigma_Rphiz (=+1 for TCV cocos=17) - if (kz==0 && kmr==0) - ids_ec_launchers.beam{iant}.steering_angle_pol(i) = 0.; - else - ids_ec_launchers.beam{iant}.steering_angle_pol(i) = atan2(-kz,kmr); - end - ids_ec_launchers.beam{iant}.steering_angle_tor(i) = asin(kphi); + + % define angles based on TORAY angles, as discussed with F. Poli Oct 2024 + % Reminder: theta_toray, zero on z-axis, pointing inwards with 90deg + % phi_toray, 0deg=outwards radial, 180deg=pointing + % towards core, +: clockwise, -: counter-clockwise + % IDS angles have same convention as TORBEAM and GRAY + pol_angle_ids = launch_params{iant}{it_ok{iant}(i)}.theta_toray-90; + tor_angle_ids = launch_params{iant}{it_ok{iant}(i)}.phi_toray-180; + + ids_ec_launchers.beam{iant}.steering_angle_pol(i) = pol_angle_ids*pi/180; + ids_ec_launchers.beam{iant}.steering_angle_tor(i) = tor_angle_ids*pi/180; ids_ec_launchers.beam{iant}.spot.size(1,i) = 0.023; ids_ec_launchers.beam{iant}.spot.size(2,i) = 0.012; ids_ec_launchers.beam{iant}.spot.angle(i) = 0.0; diff --git a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m index b0b22e2d593c954f9091642d07135f68c40c2588..102ff032e284243b7eea96afa8d26d828f8f69c9 100644 --- a/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m +++ b/matlab/TCV_IMAS/tcv_get_ids_equilibrium.m @@ -5,7 +5,7 @@ function [ids_equilibrium,ids_equilibrium_description,varargout] = ... % tcv_get_ids_equilibrium(shot,ids_equil_empty,gdat_params,varargin); % % -% gdat_params: gdat_data.gdat_params to get all params passed from original call, +% gdat_params: gdat_data.gdat_params to get all params passed from original call, % in particular error_bar and cocos_out options % @@ -527,19 +527,31 @@ profiles_2d.r.data = repmat(repmat(profiles_2d.psi.dim{1},1,numel(profiles_2d.ps profiles_2d_desc.r = 'from dim{1} of ''psi'' repeated'; profiles_2d.z.data = repmat(repmat(profiles_2d.psi.dim{2}',numel(profiles_2d.psi.dim{1}),1),1,1,numel(profiles_2d.psi.dim{3})); profiles_2d_desc.z = 'from dim{2} of ''psi'' repeated'; +params_eff.data_request = 'bfields'; +bfields_gdat = gdat(params_equilibrium.shot,params_eff); +if ~isempty(bfields_gdat.Br) + profiles_2d.b_field_r = bfields_gdat.Br; + profiles_2d_desc.b_field_r = bfields_gdat.Br.label; + profiles_2d.b_field_z= bfields_gdat.Bz; + profiles_2d_desc.b_field_z = bfields_gdat.Bz.label; + profiles_2d.b_field_tor= bfields_gdat.Btor; + profiles_2d_desc.b_field_tor = bfields_gdat.Btor.label; +else + warning('B-fields could not be load from gdat(shot,''bfields''). B-fields only available for LIUQE.M, trial 1.') +end % theta = gdat(params_equilibrium.shot,'theta','machine',gdat_params.machine); profiles_2d_fieldnames = fieldnames(profiles_2d); special_fields = {'grid', 'grid_type'}; % fields needing non-automatic treatments for it=1:ntime - for i=1:numel(profiles_2d_fieldnames) - if ~any(strcmp(profiles_2d_fieldnames{i},special_fields)) - if ~isstruct(ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i})) - if ~ischar(profiles_2d.(profiles_2d_fieldnames{i}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i}).data) ... - && size(profiles_2d.(profiles_2d_fieldnames{i}).data,3)>=it - ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i}) = ... - profiles_2d.(profiles_2d_fieldnames{i}).data(:,:,it); + for i_name=1:numel(profiles_2d_fieldnames) + if ~any(strcmp(profiles_2d_fieldnames{i_name},special_fields)) + if ~isstruct(ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i_name})) + if ~ischar(profiles_2d.(profiles_2d_fieldnames{i_name}).data) && ~isempty(profiles_2d.(profiles_2d_fieldnames{i_name}).data) ... + && size(profiles_2d.(profiles_2d_fieldnames{i_name}).data,3)>=it + ids_equilibrium.time_slice{it}.profiles_2d{1}.(profiles_2d_fieldnames{i_name}) = ... + profiles_2d.(profiles_2d_fieldnames{i_name}).data(:,:,it); end else special_fields{end+1} = profiles_2d_fieldnames{i}; @@ -547,7 +559,6 @@ for it=1:ntime end end end - % special cases for it=1:ntime ids_equilibrium.time_slice{it}.profiles_2d{1}.grid_type.name = profiles_2d.grid_type.name; @@ -596,7 +607,6 @@ end % $$$ ids_equilibrium.time_slice{it}.profiles_2d{2}.psi_error_upper(:,:) = 12.*ones(ldim1,ldim2); % $$$ ids_equilibrium.time_slice{it}.profiles_2d{2}.psi_error_lower(:,:) = 10.*ones(ldim1,ldim2); % $$$ end - % cocos automatic transform if ~isempty(which('ids_generic_cocos_nodes_transformation_symbolic')) [ids_equilibrium,cocoscoeff]=ids_generic_cocos_nodes_transformation_symbolic(ids_equilibrium,'equilibrium',gdat_params.cocos_in, ... diff --git a/matlab/TCV_IMAS/tcv_ids_pf_active_definition.m b/matlab/TCV_IMAS/tcv_ids_pf_active_definition.m index 413340424beb1618110ad6bee605be8d466516d6..a53efd5a0f68fe05366377208bd8e9ea85e81030 100644 --- a/matlab/TCV_IMAS/tcv_ids_pf_active_definition.m +++ b/matlab/TCV_IMAS/tcv_ids_pf_active_definition.m @@ -88,11 +88,10 @@ mds_paths = {... '\magnetics::ipol[*,"F_006"]';...% Circuit 16 '\magnetics::ipol[*,"F_007"]';...% Circuit 17 '\magnetics::ipol[*,"F_008"]';... % Circuit 18 - 'is_in("G_001",dim_of(\magnetics::ipol,1)) ? \magnetics::ipol[*,"G_001"] : make_signal(zero(shape(data(\magnetics::ipol))[0],1.0),*,dim_of(\magnetics::ipol,0))';... % G coils % Circuit 19 + '"G_001" is_in \magnetics::ipol:dim ? \magnetics::ipol[*,"G_001"] : make_signal(zero(shape(data(\magnetics::ipol))[0],1.0),*,dim_of(\magnetics::ipol,0))';... % G coils % Circuit 19 '\magnetics::iphi';... % Connection between tf coils % Circuit 20 }; - % Combined structure combined_structure = struct(); combined_structure.coil_names = coil_names; @@ -127,7 +126,7 @@ circuit_struct.connection_matrix = struct([]); coil_column_index = 2*circuit_struct.ntotpowersupplies ; for ii=1:circuit_struct.ntotcircuits circuit_connection_matrix = zeros(circuit_struct.nnodespercircuit(ii), 2*circuit_struct.ntotelements); - + % Put power supply connection power_supply_index = ii; if circuit_struct.power_supply_current_signs{power_supply_index} == 1 @@ -137,27 +136,27 @@ for ii=1:circuit_struct.ntotcircuits circuit_connection_matrix(1,2*(power_supply_index-1)+1) = 1; circuit_connection_matrix(circuit_struct.nnodespercircuit(ii),2*(power_supply_index-1)+2) = 1; end - + % Put coil connection for jj=1:circuit_struct.ncoilpercircuit(ii) if circuit_struct.coil_current_signs{ii}(jj) == 1 circuit_connection_matrix(jj, coil_column_index + 2*(jj-1) + 1 ) = 1; circuit_connection_matrix(jj + 1, coil_column_index + 2*(jj-1) + 2 ) = 1; - + elseif circuit_struct.coil_current_signs{ii}(jj) == -1 circuit_connection_matrix(jj, coil_column_index + 2*(jj-1) + 2 ) = 1; circuit_connection_matrix(jj + 1, coil_column_index + 2*(jj-1) + 1 ) = 1; end end - + coil_column_index = coil_column_index + 2*circuit_struct.ncoilpercircuit(ii); circuit_struct.connection_matrix{ii} = circuit_connection_matrix; - + % Plot all the connaction matrices as a check if doplot plot_connection_matrix(circuit_connection_matrix, circuit_struct.power_supply_names, circuit_struct.coil_names); end - + end %% Plot connection matrix @@ -184,7 +183,7 @@ for ii=1:numel(psnames) end for ii=1:numel(cnames) for jj=1:numel(cnames{ii}) - + index = index +1; xlab{index} = [cnames{ii}{jj} 'in']; index = index +1; diff --git a/matlab/cell2str.m b/matlab/cell2str.m new file mode 100644 index 0000000000000000000000000000000000000000..464684efdf629c390bac48e1f3e6df871d274ead --- /dev/null +++ b/matlab/cell2str.m @@ -0,0 +1,50 @@ +function str = cell2str(C,n) +% function str = cell2str(C,n) +% Convert cell structure to string such that C = eval(str). +% +% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. + +assert(iscell(C),'C must be a cell'); +assert(ndims(C)<=2,'Cell arrays with more than 2 dimensions are not supported'); + +[nrow,ncol] = size(C); + +str = '{'; +for ii=1:nrow + for jj=1:ncol + celldata = C{ii,jj}; + if isstruct(celldata) + structstrcell = struct2str(celldata,n); % returns cell array of strings + datastr = ''; + for kk = 1:numel(structstrcell) + datastr = [datastr,structstrcell{kk},'\n']; %#ok<AGROW> + end + datastr = sprintf(datastr); + % contatenate them to get a single string again + elseif iscell(celldata) + datastr = cell2str(celldata,n); % recursive call + else + datastr = field2str(celldata,n); + end + str = [str,datastr]; %#ok<AGROW> + if jj == ncol && ii~=nrow + str = [str,';']; %#ok<AGROW> + elseif jj~=ncol + str = [str,',']; %#ok<AGROW> + end + end +end +str = [str,'}']; + +end + +function datastr = field2str(data,n) + +if isnumeric(data) || islogical(data) + datastr = mat2str(data,n,'class'); +elseif ischar(data) + datastr = ['''',data,'''']; +else + error('%s datatype not supported:',class(data)) +end +end