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) % % nbdim_x = 1: same x(1:length(x)) (rhopol) for all times % = 2: x(:,t), rhopol depends on time like Thomson projection % % nopt = 0: do not fill in, just make the empty structure % = 1: do compute various fields % % compute psi, rhotor_edge, rhotornorm, volume_edge and rhovol and b0 used for rhotor % % use gdat calls to psi_axis, psi_edge, rhotor, etc with same basic parameters as data.gdat_params % gdat_data.grids_1d.rhopolnorm = gdat_data.x; if (nopt == 0) || isempty(gdat_data.x) || isempty(gdat_data.t) || isempty(gdat_data.data) || ischar(gdat_data.data) gdat_data.grids_1d.rhotornorm = []; gdat_data.grids_1d.rhovolnorm = []; gdat_data.grids_1d.psi = []; gdat_data.grids_1d.rhotor_edge = []; gdat_data.grids_1d.volume_edge = []; gdat_data.grids_1d.b0 = []; return 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); ndim_x_rhotor = length(find(size(rhotor_norm.x)>1)); params_eff.data_request='rhovol'; rhovol = gdat(gdat_data.shot,params_eff); ndim_x_rhovol = length(find(size(rhotor_norm.x)>1)); % check that rhotor and rhovol come from same equil psi if sum(abs(size(rhotor_norm.x)-size(rhovol.x)))~=0 || sum(sum(abs(rhotor_norm.x-rhovol.x)))>1e-10 ... || sum(abs(rhotor_norm.t-rhovol.t))>1e-10; error(['get_grids_1d: x and t arrays for rhotor_norm and rhovol are different, although should refer to same equilibrium' char(10) ... 'size(rhotor_norm.x) = ' num2str(size(rhotor_norm.x)) '; size(rhovol.x) = ' num2str(size(rhovol.x))]); end 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); 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) gdat_data.grids_1d.psi = gdat_data.grids_1d.rhopolnorm.^2.*repmat(reshape(psi0_edge,1,length(psi0_edge)),size(gdat_data.grids_1d.rhopolnorm,1),1); else if nverbose>=0; disp(['option: nbdim_x = ' numstr(nbdim_x) ' not implemented, check with O. Sauter']); end return end gdat_data.grids_1d.rhotornorm = NaN*ones(size(gdat_data.data)); gdat_data.grids_1d.rhovolnorm = NaN*ones(size(gdat_data.data)); 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); tens0=0.; for it=1:length(gdat_data.t) % do an interpolation on closest point to avoid 2D interp it_rt_eff = it_rt(it); it_vol_eff = it_vol(it); if ndim_x_rhotor==1 gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x; % to have same sizes for all 3 reference rhos from equil else gdat_data.grids_1d.rhopolnorm_equil(:,it) = rhotor_norm.x(:,it_rt_eff); end gdat_data.grids_1d.rhotornorm_equil(:,it) = rhotor_norm.data(:,it_rt_eff); % can assume same .x for rhotor and rhovol since called with same params thus equil (checked with sum before) gdat_data.grids_1d.rhovolnorm_equil(:,it) = rhovol.data(:,it_vol_eff); % if (nbdim_x == 1) ii=find(isfinite(gdat_data.grids_1d.rhopolnorm)); else ii=find(isfinite(gdat_data.grids_1d.rhopolnorm(:,it))); end if (nbdim_x == 1) 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(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0); else gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0); end end if ~isempty(rhovol.x) && ~isempty(rhovol.data) if ndim_x_rhovol==1 gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0); else gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii),tens0); end end end else % $$$ 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(1:nb_ii,it)=interpos(23,rhotor_norm.x,rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0); else gdat_data.grids_1d.rhotornorm(1:nb_ii,it)=interpos(23,rhotor_norm.x(:,it_rt_eff),rhotor_norm.data(:,it_rt_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0); end end if ~isempty(rhovol.x) && ~isempty(rhovol.data) if ndim_x_rhovol==1 gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x,rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0); else gdat_data.grids_1d.rhovolnorm(1:nb_ii,it)=interpos(23,rhovol.x(:,it_vol_eff),rhovol.data(:,it_vol_eff),gdat_data.grids_1d.rhopolnorm(ii,it),tens0); end end end end end gdat_data.grids_1d.rhotor_edge=interpos(-63,rhotor_norm.t',rhotor_norm.rhotor_edge,gdat_data.t',-0.01); gdat_data.grids_1d.volume_edge=interpos(-63,rhovol.t',rhovol.volume_edge,gdat_data.t',-0.01); gdat_data.grids_1d.b0=interpos(-63,rhotor_norm.t',rhotor_norm.b0,gdat_data.t',-0.01);