diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m index 73c36a27da6ff60bf0b887a35fc857710ef04a98..86f18fbe06e5bc068976d96b6c7ab1aa6ac66c36 100644 --- a/crpptbx/AUG/gdat_aug.m +++ b/crpptbx/AUG/gdat_aug.m @@ -217,7 +217,7 @@ data_request_eff = gdat_params.data_request; % in case was defined in pairs % default equil: if ~isfield(gdat_params,'equil'); gdat_params.equil = 'EQI'; end if isfield(gdat_params,'source') && (any(strcmp(lower(gdat_params.source),'ide')) || any(strcmp(lower(gdat_params.source),'idg')) ... - || any(strcmp(lower(gdat_params.source),'ida'))) + || any(strcmp(lower(gdat_params.source),'ida'))) gdat_params.equil = 'IDE'; end extra_arg_sf2sig = '[]'; @@ -310,7 +310,7 @@ if strcmp(mapping_for_aug.method,'signal') gdat_data.gdat_params.extra_arg_sf2sig = mapping_for_aug.expression{4}; end [aatmp,error_status]=rdaAUG_eff(shot,mapping_for_aug.expression{1},mapping_for_aug.expression{2},exp_location, ... - time_interval,gdat_data.gdat_params.extra_arg_sf2sig); + time_interval,gdat_data.gdat_params.extra_arg_sf2sig); if error_status~=0 if gdat_params.nverbose>=3; disp(['error after rdaAUG in signal with data_request_eff= ' data_request_eff]); end return @@ -329,9 +329,9 @@ if strcmp(mapping_for_aug.method,'signal') elseif length(size(aatmp.data))==2 % 2 true dimensions if length(aatmp.t) == size(aatmp.data,1) - mapping_for_aug.timedim = 1; + mapping_for_aug.timedim = 1; else - mapping_for_aug.timedim = 2; + mapping_for_aug.timedim = 2; end else % more than 2 dimensions, find the one with same length to define timedim @@ -353,9 +353,9 @@ if strcmp(mapping_for_aug.method,'signal') nbdims = length(size(gdat_data.data)); for i=1:nbdims if i~=mapping_for_aug.timedim - ieff=i; - if i>mapping_for_aug.timedim; ieff=i-1; end - gdat_data.dim{i} = gdat_data.x{ieff}; + ieff=i; + if i>mapping_for_aug.timedim; ieff=i-1; end + gdat_data.dim{i} = gdat_data.x{ieff}; end end end @@ -513,8 +513,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') a.name = 'vrot'; if isempty(a.data) || isempty(a.t) || error_status>0 if gdat_params.nverbose>=3; - a - disp(['with data_request= ' data_request_eff]) + a + disp(['with data_request= ' data_request_eff]) end end gdat_data.vrot.data = a.data; @@ -545,8 +545,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') params_equil.data_request = 'equil'; [equil,params_equil,error_status] = gdat_aug(shot,params_equil); if error_status>0 - if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end - return + if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end + return end gdat_data.gdat_params.equil = params_equil.equil; gdat_data.equil = equil; @@ -560,30 +560,30 @@ elseif strcmp(mapping_for_aug.method,'switchcase') time_equil=[min(gdat_data.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.t(end)+0.1)]; iok=find(~isnan(gdat_data.r(:,1))); for itequil=1:length(time_equil)-1 - rr=equil.Rmesh(:,itequil); - zz=equil.Zmesh(:,itequil); - psirz_in = equil.psi2D(:,:,itequil); - it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1)); - if ~isempty(it_cxrs_inequil) - if strcmp(upper(gdat_data.gdat_params.source),'CEZ') - rout=gdat_data.r(iok,it_cxrs_inequil); - zout=gdat_data.z(iok,it_cxrs_inequil); - else - rout=gdat_data.r(iok); - zout=gdat_data.z(iok); - end - psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); - if strcmp(upper(gdat_data.gdat_params.source),'CEZ') - psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil)); - else - psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]); - end - rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); - for it_cx=1:length(it_cxrs_inequil) - rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); - rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); - end - end + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_cxrs_inequil = find(gdat_data.t>time_equil(itequil) & gdat_data.t<=time_equil(itequil+1)); + if ~isempty(it_cxrs_inequil) + if strcmp(upper(gdat_data.gdat_params.source),'CEZ') + rout=gdat_data.r(iok,it_cxrs_inequil); + zout=gdat_data.z(iok,it_cxrs_inequil); + else + rout=gdat_data.r(iok); + zout=gdat_data.z(iok); + end + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); + if strcmp(upper(gdat_data.gdat_params.source),'CEZ') + psi_out(iok,it_cxrs_inequil) = reshape(psi_at_routzout,length(iok),length(it_cxrs_inequil)); + else + psi_out(iok,it_cxrs_inequil) = repmat(psi_at_routzout,[1,length(it_cxrs_inequil)]); + end + rhopolnorm_out(iok,it_cxrs_inequil) = sqrt((psi_out(iok,it_cxrs_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_cxrs_inequil) + rhotornorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out(iok,it_cxrs_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(iok,it_cxrs_inequil(it_cx)),-3,[2 2],[0 1]); + end + end end gdat_data.psi = psi_out; gdat_data.rhopolnorm = rhopolnorm_out; @@ -602,67 +602,67 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.fit.raw.vrot.data = []; fit_tension_default = -1; if isfield(gdat_data.gdat_params,'fit_tension') - fit_tension = gdat_data.gdat_params.fit_tension; + fit_tension = gdat_data.gdat_params.fit_tension; else - fit_tension = fit_tension_default; + fit_tension = fit_tension_default; end if ~isstruct(fit_tension) - fit_tension_eff.ti = fit_tension; - fit_tension_eff.vrot = fit_tension; - fit_tension = fit_tension_eff; + fit_tension_eff.ti = fit_tension; + fit_tension_eff.vrot = fit_tension; + fit_tension = fit_tension_eff; else - if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end - if ~isfield(fit_tension,'vrot '); fit_tension.vrot = fit_tension_default; end + if ~isfield(fit_tension,'ti'); fit_tension.ti = fit_tension_default; end + if ~isfield(fit_tension,'vrot '); fit_tension.vrot = 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; + fit_nb_rho_points = gdat_data.gdat_params.fit_nb_rho_points; else - fit_nb_rho_points = 201; + 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.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data)); - gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data)); - gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data)); - rhotornormfit = linspace(0,1,fit_nb_rho_points)'; - gdat_data.fit.rhotornorm = rhotornormfit; - gdat_data.fit.t = gdat_data.t; - for it=1:length(gdat_data.t) - % make rhotor->rhopol transformation for each time since equilibrium might have changed - irhook=find(gdat_data.rhotornorm(:,it)>0 & gdat_data.rhotornorm(:,it)<1); % no need for ~isnan - [rhoeff isort]=sort(gdat_data.rhotornorm(irhook,it)); - gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]); - idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01); - if length(idata)>0 - gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); - gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it); - gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it); - gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); - gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it); - gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it); - [rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it)); - rhoeff = [0; rhoeff]; - % they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar - tieff = gdat_data.ti.data(idata(irhoeff),it); - ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it); - ij=find(tieff./ti_err_eff>10.); - if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end - vroteff = gdat_data.vrot.data(idata(irhoeff),it); - vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it); - ij=find(vroteff./vrot_err_eff>10.); - if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end - % - tieff = [tieff(1); tieff]; - ti_err_eff = [1e4; ti_err_eff]; - vroteff = [vroteff(1); vroteff]; - vrot_err_eff = [1e5; vrot_err_eff]; - [gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff); - [gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff); - end - end + % add fits + gdat_data.fit.raw.rhotornorm = NaN*ones(size(gdat_data.ti.data)); + gdat_data.fit.raw.ti.data = NaN*ones(size(gdat_data.ti.data)); + gdat_data.fit.raw.vrot.data = NaN*ones(size(gdat_data.vrot.data)); + rhotornormfit = linspace(0,1,fit_nb_rho_points)'; + gdat_data.fit.rhotornorm = rhotornormfit; + gdat_data.fit.t = gdat_data.t; + for it=1:length(gdat_data.t) + % make rhotor->rhopol transformation for each time since equilibrium might have changed + irhook=find(gdat_data.rhotornorm(:,it)>0 & gdat_data.rhotornorm(:,it)<1); % no need for ~isnan + [rhoeff isort]=sort(gdat_data.rhotornorm(irhook,it)); + gdat_data.fit.rhopolnorm(:,it)=interpos([0; rhoeff; 1],[0; gdat_data.rhopolnorm(irhook(isort),it); 1],rhotornormfit,-0.1,[2 2],[0 1]); + idata = find(gdat_data.ti.data(:,it)>0 & gdat_data.rhotornorm(:,it)<1.01); + if length(idata)>0 + gdat_data.fit.ti.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); + gdat_data.fit.ti.raw.data(idata,it) = gdat_data.ti.data(idata,it); + gdat_data.fit.ti.raw.error_bar(idata,it) = gdat_data.ti.error_bar(idata,it); + gdat_data.fit.vrot.raw.rhotornorm(idata,it) = gdat_data.rhotornorm(idata,it); + gdat_data.fit.vrot.raw.data(idata,it) = gdat_data.vrot.data(idata,it); + gdat_data.fit.vrot.raw.error_bar(idata,it) = gdat_data.vrot.error_bar(idata,it); + [rhoeff,irhoeff] = sort(gdat_data.rhotornorm(idata,it)); + rhoeff = [0; rhoeff]; + % they are some strange low error_bars, so remove these by max(Ti/error_bar)<=10; and changing it to large error bar + tieff = gdat_data.ti.data(idata(irhoeff),it); + ti_err_eff = gdat_data.ti.error_bar(idata(irhoeff),it); + ij=find(tieff./ti_err_eff>10.); + if ~isempty(ij); ti_err_eff(ij) = tieff(ij)./0.1; end + vroteff = gdat_data.vrot.data(idata(irhoeff),it); + vrot_err_eff = gdat_data.vrot.error_bar(idata(irhoeff),it); + ij=find(vroteff./vrot_err_eff>10.); + if ~isempty(ij); vrot_err_eff(ij) = vroteff(ij)./0.1; end + % + tieff = [tieff(1); tieff]; + ti_err_eff = [1e4; ti_err_eff]; + vroteff = [vroteff(1); vroteff]; + vrot_err_eff = [1e5; vrot_err_eff]; + [gdat_data.fit.ti.data(:,it), gdat_data.fit.ti.drhotornorm(:,it)] = interpos(rhoeff,tieff,rhotornormfit,fit_tension.ti,[1 0],[0 0],ti_err_eff); + [gdat_data.fit.vrot.data(:,it), gdat_data.fit.vrot.drhotornorm(:,it)] = interpos(rhoeff,vroteff,rhotornormfit,fit_tension.vrot,[1 0],[0 0],vrot_err_eff); + end + end end end @@ -738,8 +738,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') if match_rz_to_time % interpolate R structure on ece data time array, to ease plot vs R for i=1:length(channels) - radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t); - zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t); + radius.data(i,:) = interp1(aR.t,aR.data(channels(i),:),gdat_data.t); + zheight.data(i,:) = interp1(aZ.t,aZ.data(channels(i),:),gdat_data.t); end radiuszheight.t=gdat_data.t; else @@ -759,8 +759,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') params_equil.data_request = 'equil'; [equil,params_equil,error_status] = gdat_aug(shot,params_equil); if error_status>0 - if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end - return + if gdat_params.nverbose>=3; disp(['problems with ' params_equil.data_request]); end + return end gdat_data.gdat_params.equil = params_equil.equil; gdat_data.equil = equil; @@ -774,29 +774,29 @@ elseif strcmp(mapping_for_aug.method,'switchcase') % constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)] time_equil=[min(gdat_data.rz.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.rz.t(end)+0.1)]; for itequil=1:length(time_equil)-1 - rr=equil.Rmesh(:,itequil); - zz=equil.Zmesh(:,itequil); - psirz_in = equil.psi2D(:,:,itequil); - it_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1)); - if ~isempty(it_ece_inequil) - r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil); - iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0); - if ~isempty(iok) - rout=r_it_ece_inequil(iok); - z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil); - zout=z_it_ece_inequil(iok); - psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); - [ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok); - for ij=1:length(iok) - psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij); - end - rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)))); - for it_cx=1:length(it_ece_inequil) - rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); - rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); - end - end - end + rr=equil.Rmesh(:,itequil); + zz=equil.Zmesh(:,itequil); + psirz_in = equil.psi2D(:,:,itequil); + it_ece_inequil = find(gdat_data.rz.t>time_equil(itequil) & gdat_data.rz.t<=time_equil(itequil+1)); + if ~isempty(it_ece_inequil) + r_it_ece_inequil = gdat_data.rz.r(:,it_ece_inequil); + iok = find(~isnan(r_it_ece_inequil) & r_it_ece_inequil>0); + if ~isempty(iok) + rout=r_it_ece_inequil(iok); + z_it_ece_inequil = gdat_data.rz.z(:,it_ece_inequil); + zout=z_it_ece_inequil(iok); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout,zout); + [ieff,jeff]=ind2sub(size(gdat_data.rz.r(:,it_ece_inequil)),iok); + for ij=1:length(iok) + psi_out(ieff(ij),it_ece_inequil(jeff(ij))) = psi_at_routzout(ij); + end + rhopolnorm_out(:,it_ece_inequil) = sqrt(abs((psi_out(:,it_ece_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil)))); + for it_cx=1:length(it_ece_inequil) + rhotornorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out(:,it_ece_inequil(it_cx)) = interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out(:,it_ece_inequil(it_cx)),-3,[2 2],[0 1]); + end + end + end end gdat_data.rhos.psi = psi_out; gdat_data.rhos.rhopolnorm = rhopolnorm_out; @@ -845,9 +845,9 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.gdat_params.cocos = cocos_out; end if equil.cocos ~= cocos_out - [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskAUG,[equil.cocos cocos_out]); + [eqdsk_cocosout, eqdsk_cocosout_IpB0pos,cocos_inout]=eqdsk_cocos_transform(eqdskAUG,[equil.cocos cocos_out]); else - eqdsk_cocosout = eqdskAUG; + eqdsk_cocosout = eqdskAUG; end % for several times, use array of structure for eqdsks, % cannot use it for psi(R,Z) in .data and .x since R, Z might be different at different times, @@ -859,10 +859,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.dim{1} = gdat_data.eqdsk{itime}.rmesh; gdat_data.dim{2} = gdat_data.eqdsk{itime}.zmesh; else - xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2)); - yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1); - aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ... - ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1); + xx=repmat(reshape(gdat_data.dim{1},length(gdat_data.dim{1}),1),1,size(gdat_data.eqdsk{itime}.psi,2)); + yy=repmat(reshape(gdat_data.dim{2},1,length(gdat_data.dim{2})),size(gdat_data.eqdsk{itime}.psi,1),1); + aa = interpos2Dcartesian(gdat_data.eqdsk{itime}.rmesh,gdat_data.eqdsk{itime}.zmesh ... + ,gdat_data.eqdsk{itime}.psi,xx,yy,-1,-1); gdat_data.data(:,:,itime) = aa; end else @@ -935,10 +935,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase') [B2ave,e]=rdaAUG_eff(shot,DIAG,'B2ave',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); B2ave=adapt_rda(B2ave,NTIME,ndimrho,itotransposeback); if strcmp(DIAG,'IDE') - FTRA.value=[]; + FTRA.value=[]; else - [FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); - FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback); + [FTRA,e]=rdaAUG_eff(shot,DIAG,'FTRA',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); + FTRA=adapt_rda(FTRA,NTIME,ndimrho,itotransposeback); end else Rinv.value=[]; R2inv.value=[]; Bave.value=[]; B2ave.value=[]; FTRA.value=[]; @@ -960,9 +960,9 @@ elseif strcmp(mapping_for_aug.method,'switchcase') ijok=find(abs(qpsi.value)>0); % note: eqr fills in only odd points radially % set NaNs to zeroes if qpsi.value(ijok(1))<0 - gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue); + 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); + gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue); end % get x values gdat_data.psi(:,it)=psi_tree.value(it,Lpf1:-1:1)'; @@ -970,24 +970,24 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.psi_lcfs(it)= gdat_data.psi(end,it); gdat_data.rhopolnorm(:,it) = sqrt(abs((gdat_data.psi(:,it)-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it)))); if strcmp(DIAG,'EQR'); - % q value has only a few values and from center to edge, assume they are from central rhopol values on - % But they are every other point starting from 3rd - ijk=find(gdat_data.qvalue(:,it)~=0); - if length(ijk)>2 - % now shots have non-zero axis values in eqr - rhoeff=gdat_data.rhopolnorm(ijk,it); - qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above - if ijk(1)>1 - rhoeff = [0.; rhoeff]; - qeff = [qeff(1) ;qeff]; - end - ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); - qfit = zeros(size(gdat_data.rhopolnorm(:,it))); - qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); - else - qfit = zeros(size(gdat_data.rhopolnorm(:,it))); - end - gdat_data.qvalue(:,it) = qfit; + % q value has only a few values and from center to edge, assume they are from central rhopol values on + % But they are every other point starting from 3rd + ijk=find(gdat_data.qvalue(:,it)~=0); + if length(ijk)>2 + % now shots have non-zero axis values in eqr + rhoeff=gdat_data.rhopolnorm(ijk,it); + qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above + if ijk(1)>1 + rhoeff = [0.; rhoeff]; + qeff = [qeff(1) ;qeff]; + end + ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); + else + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + end + gdat_data.qvalue(:,it) = qfit; end % get rhotor values gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)'; @@ -1004,40 +1004,40 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.pressure(:,it)=Pres.value(it,2*Lpf1-1:-2:1)'; gdat_data.dpressuredpsi(:,it)=Pres.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi if ~isempty(Jpol.value) - gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)'; - gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + gdat_data.jpol(:,it)=Jpol.value(it,2*Lpf1-1:-2:1)'; + gdat_data.djpolpsi(:,it)=Jpol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi else - gdat_data.jpol = []; - gdat_data.djpolpsi = []; + gdat_data.jpol = []; + gdat_data.djpolpsi = []; end gdat_data.ffprime(:,it) = FFP.value(it,Lpf1:-1:1)'; gdat_data.Xpoints.psi(1:LPFx.value(it)+1,it) = PFxx.value(it,1:LPFx.value(it)+1); gdat_data.Xpoints.Rvalue(1:LPFx.value(it)+1,it) = RPFx.value(it,1:LPFx.value(it)+1); gdat_data.Xpoints.Zvalue(1:LPFx.value(it)+1,it) = zPFx.value(it,1:LPFx.value(it)+1); if ~isempty(Rinv.value) - gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)'; + gdat_data.rinv(:,it) = Rinv.value(it,Lpf1:-1:1)'; else - gdat_data.rinv = []; + gdat_data.rinv = []; end if ~isempty(R2inv.value) - gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)'; + gdat_data.r2inv(:,it) = R2inv.value(it,Lpf1:-1:1)'; else - gdat_data.r2inv = []; + gdat_data.r2inv = []; end if ~isempty(Bave.value) - gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)'; + gdat_data.bave(:,it) = Bave.value(it,Lpf1:-1:1)'; else - gdat_data.bave = []; + gdat_data.bave = []; end if ~isempty(B2ave.value) - gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)'; + gdat_data.b2ave(:,it) = B2ave.value(it,Lpf1:-1:1)'; else - gdat_data.b2ave = []; + gdat_data.b2ave = []; end if ~isempty(FTRA.value) - gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)'; + gdat_data.ftra(:,it) = FTRA.value(it,Lpf1:-1:1)'; else - gdat_data.ftra = []; + gdat_data.ftra = []; end % end @@ -1051,8 +1051,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') if strcmp(DIAG,'IDE') [IpiPSI,e]=rdaAUG_eff(shot,'IDG','Itor',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); if length(IpiPSI.value)~=NTIME - disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter') - return + disp('Problem with nb time points of IDG/Itor with respect to IDE NTIME, check with O. Sauter') + return end else [IpiPSI,e]=rdaAUG_eff(shot,DIAG,'IpiPSI',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); @@ -1077,8 +1077,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') [a,error_status]=rdaAUG_eff(shot,'VTA',nodenameeff,exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); if isempty(a.data) || isempty(a.t) || error_status>0 if gdat_params.nverbose>=3; - a - disp(['with data_request= ' data_request_eff]) + a + disp(['with data_request= ' data_request_eff]) end return end @@ -1123,7 +1123,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase') [alow,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'elow_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); [aup,e]=rdaAUG_eff(shot,'VTA',[upper(data_request_eff(1)) 'eupp_e'],exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); gdat_data.(lower(node_child_nameeff_e)).error_bar = max(aup.value-gdat_data.(lower(node_child_nameeff_e)).data, ... - gdat_data.(lower(node_child_nameeff_e)).data-alow.value); + gdat_data.(lower(node_child_nameeff_e)).data-alow.value); gdat_data.error_bar(nb_core+1:nb_core+nb_edge,:) = gdat_data.(lower(node_child_nameeff_e)).error_bar(1:nb_edge,iaaa); else gdat_data.(lower(node_child_nameeff_e)).data = []; @@ -1183,32 +1183,32 @@ elseif strcmp(mapping_for_aug.method,'switchcase') psirz_in = equil.psi2D(:,:,itequil); it_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1)); if ~isempty(it_core_inequil) - rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil); - zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil); - psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core); - psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(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 + rout_core=gdat_data.(lower(node_child_nameeff)).r(:,it_core_inequil); + zout_core=gdat_data.(lower(node_child_nameeff)).z(:,it_core_inequil); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_core,zout_core); + psi_out_core(:,it_core_inequil) = reshape(psi_at_routzout,inb_chord_core,length(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 % edge it_edge_inequil = find(gdat_data.(lower(node_child_nameeff_e)).t>time_equil(itequil) & gdat_data.(lower(node_child_nameeff_e)).t<=time_equil(itequil+1)); if ~isempty(it_edge_inequil) - rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil); - zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil); - psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge); - psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil)); - rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); - for it_cx=1:length(it_edge_inequil) - rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ... - interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); - rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ... - interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); - end + rout_edge=gdat_data.(lower(node_child_nameeff_e)).r(:,it_edge_inequil); + zout_edge=gdat_data.(lower(node_child_nameeff_e)).z(:,it_edge_inequil); + psi_at_routzout = interpos2Dcartesian(rr,zz,psirz_in,rout_edge,zout_edge); + psi_out_edge(:,it_edge_inequil) = reshape(psi_at_routzout,inb_chord_edge,length(it_edge_inequil)); + rhopolnorm_out_edge(:,it_edge_inequil) = sqrt((psi_out_edge(:,it_edge_inequil)-equil.psi_axis(itequil))./(equil.psi_lcfs(itequil)-equil.psi_axis(itequil))); + for it_cx=1:length(it_edge_inequil) + rhotornorm_out_edge(:,it_edge_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhotornorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); + rhovolnorm_out_edge(:,it_edge_inequil(it_cx)) = ... + interpos(equil.rhopolnorm(:,itequil),equil.rhovolnorm(:,itequil),rhopolnorm_out_edge(:,it_edge_inequil(it_cx)),-3,[2 2],[0 1]); + end end end gdat_data.(lower(node_child_nameeff)).psi = psi_out_core; @@ -1311,62 +1311,62 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.fit.rhotornorm = rhotornormfit; gdat_data.fit.rhopolnorm = NaN*ones(length(rhotornormfit),length(gdat_data.fit.t)); 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.rhotornorm = NaN*ones(size(gdat_data.te.data)); - gdat_data.fit.te.data = gdat_data.fit.rhopolnorm; - gdat_data.fit.te.drhotornorm = gdat_data.fit.rhopolnorm; + 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.rhotornorm = NaN*ones(size(gdat_data.te.data)); + gdat_data.fit.te.data = gdat_data.fit.rhopolnorm; + gdat_data.fit.te.drhotornorm = gdat_data.fit.rhopolnorm; 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.rhotornorm = NaN*ones(size(gdat_data.ne.data)); - gdat_data.fit.ne.data =gdat_data.fit.rhopolnorm; - gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhopolnorm; + 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.rhotornorm = NaN*ones(size(gdat_data.ne.data)); + gdat_data.fit.ne.data =gdat_data.fit.rhopolnorm; + gdat_data.fit.ne.drhotornorm = gdat_data.fit.rhopolnorm; 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],rhotornormfit,-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,rhotornormfit,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,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff); - end - end + % 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],rhotornormfit,-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,rhotornormfit,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,rhotornormfit,fit_tension.ne,[1 0],[0 0],ne_err_eff); + end + end end end @@ -1402,122 +1402,122 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i} = gdat_data_i; - gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1); - gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; - if max(gdat_data_i.data) > 0. - labels{i} = ['EC_' num2str(i)]; - end + gdat_data.ec{i} = gdat_data_i; + gdat_data.data(:,i) = reshape(gdat_data_i.data,nb_timepoints,1); + gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; + if max(gdat_data_i.data) > 0. + labels{i} = ['EC_' num2str(i)]; + end end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy1_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']); % eval(['a = sf2par(''ECS'',shot,''gyr_freq'',''P_sy1_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - % gyr_freq not present (old shots for example) - a=[]; + % gyr_freq not present (old shots for example) + a=[]; end if isempty(a) else - gdat_data.ec{i}.freq = a; - gdat_data.freq_ec(i) = a.value; + gdat_data.ec{i}.freq = a; + gdat_data.freq_ec(i) = a.value; end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy1_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']); - % eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy1_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); + % eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy1_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - % GPolPos not present - a=[]; + % GPolPos not present + a=[]; end if isempty(a) else - gdat_data.ec{i}.polpos = a.value; - gdat_data.polpos_ec(i) = a.value; + gdat_data.ec{i}.polpos = a.value; + gdat_data.polpos_ec(i) = a.value; end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy1_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']); - % eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy1_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); + % eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy1_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - a=[]; + a=[]; end if isempty(a) else - gdat_data.ec{i}.torpos = a.value; - gdat_data.torpos_ec(i) = a.value; + gdat_data.ec{i}.torpos = a.value; + gdat_data.torpos_ec(i) = a.value; end % "new" ECRH2 gyrotrons: gyro 5 to 8 in pgyro params_eff.data_request={'ECS',['PG' num2str(i) 'N']}; gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i+4} = gdat_data_i; - gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1); - gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; - if max(gdat_data_i.data) > 0. - labels{i+4} = ['EC_' num2str(i+4)]; - end + gdat_data.ec{i+4} = gdat_data_i; + gdat_data.data(:,i+4) = reshape(gdat_data_i.data,nb_timepoints,1); + gdat_data.dim = [{gdat_data_i.t} {[1:9]}]; + if max(gdat_data_i.data) > 0. + labels{i+4} = ['EC_' num2str(i+4)]; + end end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'gyr_freq']); - % eval(['a = sf2par(''ECS'',shot,''gyr_freq'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); + % eval(['a = sf2par(''ECS'',shot,''gyr_freq'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - a=[]; + a=[]; end if isempty(a) else - gdat_data.ec{i+4}.freq = a; - gdat_data.freq_ec(i+4) = a.value; + gdat_data.ec{i+4}.freq = a; + gdat_data.freq_ec(i+4) = a.value; end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GPolPos']); - % eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); + % eval(['a = sf2par(''ECS'',shot,''GPolPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - a=[]; + a=[]; end if isempty(a) else - gdat_data.ec{i+4}.polpos = a.value; - gdat_data.polpos_ec(i+4) = a.value; + gdat_data.ec{i+4}.polpos = a.value; + gdat_data.polpos_ec(i+4) = a.value; end try [a,e]=rdaAUG_eff(shot,'ECS',['P_sy2_g' num2str(i)],gdat_data.gdat_params.exp_name,[], ... gdat_data.gdat_params.extra_arg_sf2sig,['param:' 'GTorPos']); - % eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); + % eval(['a = sf2par(''ECS'',shot,''GTorPos'',''P_sy2_g' num2str(i) '''' extra_arg_sf2sig_eff_string ');']); catch - a=[]; + a=[]; end if isempty(a) else - gdat_data.ec{i+4}.torpos = a.value; - gdat_data.torpos_ec(i+4) = a.value; + gdat_data.ec{i+4}.torpos = a.value; + gdat_data.torpos_ec(i+4) = a.value; end params_eff.data_request={'ECN',['G' num2str(i) 'POL']}; gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i+4}.gpol_ec = gdat_data_i; + gdat_data.ec{i+4}.gpol_ec = gdat_data_i; end params_eff.data_request={'ECN',['G' num2str(i) 'TOR']}; gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i+4}.gtor_ec = gdat_data_i; + gdat_data.ec{i+4}.gtor_ec = gdat_data_i; end params_eff.data_request={'ECN',['G' num2str(i) 'PO4']}; gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i+4}.gpo4_ec = gdat_data_i; + gdat_data.ec{i+4}.gpo4_ec = gdat_data_i; end params_eff.data_request={'ECN',['G' num2str(i) 'PO8']}; gdat_data_i=gdat_aug(shot,params_eff); if isempty(gdat_data_i.data) || isempty(gdat_data_i.dim) else - gdat_data.ec{i+4}.gpo8_ec = gdat_data_i; + gdat_data.ec{i+4}.gpo8_ec = gdat_data_i; end end if ~isempty(gdat_data.dim) @@ -1529,10 +1529,10 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.data_fullpath=['ECS/' 'PGi and PGiN, etc']; icount=0; for i=1:length(labels) - if ~isempty(labels{i}) - icount=icount+1; - gdat_data.label{icount} = labels{i}; - end + if ~isempty(labels{i}) + icount=icount+1; + gdat_data.label{icount} = labels{i}; + end end else gdat_data.freq_ech_units =[]'; @@ -1560,7 +1560,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.gdat_params.source = sources_avail; elseif ~iscell(gdat_data.gdat_params.source) if ischar(gdat_data.gdat_params.source) - gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source); + gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source); if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail))) if (gdat_params.nverbose>=1) warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]); @@ -1599,9 +1599,9 @@ elseif strcmp(mapping_for_aug.method,'switchcase') end if ~isempty(ohm.data) && ~isempty(ohm.dim) for i=1:length(fields_to_copy) - if isfield(ohm,fields_to_copy{i}) - gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i}); - end + if isfield(ohm,fields_to_copy{i}) + gdat_data.ohm.(fields_to_copy{i}) = ohm.(fields_to_copy{i}); + end end gdat_data.ohm.raw_data = gdat_data.ohm.data; else @@ -1627,19 +1627,19 @@ elseif strcmp(mapping_for_aug.method,'switchcase') params_eff.data_request={'ECS','PECRH'}; params_eff.data_request='pgyro'; try - ec=gdat_aug(shot,params_eff); + ec=gdat_aug(shot,params_eff); catch end if ~isempty(ec.data) && ~isempty(ec.dim) - for i=1:length(fields_to_copy) - % if has pgyro, use not_copy - if isfield(ec,fields_to_copy{i}) && ~any(strmatch(fields_to_not_copy,fields_to_copy{i})) - gdat_data.ec.(fields_to_copy{i}) = ec.(fields_to_copy{i}); - end - end - gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t); - gdat_data.x(end+1) =gdat_data.x(end)+1; - gdat_data.label{end+1}='P_{ec}'; + for i=1:length(fields_to_copy) + % if has pgyro, use not_copy + if isfield(ec,fields_to_copy{i}) && ~any(strmatch(fields_to_not_copy,fields_to_copy{i})) + gdat_data.ec.(fields_to_copy{i}) = ec.(fields_to_copy{i}); + end + end + gdat_data.data(:,end+1) = interpos(-21,gdat_data.ec.t,gdat_data.ec.data(:,end),gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{ec}'; end end % @@ -1647,21 +1647,21 @@ elseif strcmp(mapping_for_aug.method,'switchcase') % nbi params_eff.data_request={'NIS','PNIQ'}; try - nbi=gdat_aug(shot,params_eff); + nbi=gdat_aug(shot,params_eff); catch nbi.data = []; nbi.dim = []; end if ~isempty(nbi.data) && ~isempty(nbi.dim) - for i=1:length(fields_to_copy) - if isfield(nbi,fields_to_copy{i}) - gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i}); - end - end + for i=1:length(fields_to_copy) + if isfield(nbi,fields_to_copy{i}) + gdat_data.nbi.(fields_to_copy{i}) = nbi.(fields_to_copy{i}); + end + end % add to main with linear interpolation and 0 for extrapolated values - gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t); - gdat_data.x(end+1) =gdat_data.x(end)+1; - gdat_data.label{end+1}='P_{nbi}'; + gdat_data.data(:,end+1) = interpos(-21,gdat_data.nbi.t,gdat_data.nbi.data(:,end),gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{nbi}'; end end % @@ -1669,18 +1669,18 @@ elseif strcmp(mapping_for_aug.method,'switchcase') % ic params_eff.data_request={'ICP','PICRN'}; try - ic=gdat_aug(shot,params_eff); + ic=gdat_aug(shot,params_eff); catch end if ~isempty(ic.data) && ~isempty(ic.dim) - for i=1:length(fields_to_copy) - if isfield(ic,fields_to_copy{i}) - gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i}); - end - end - gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t); - gdat_data.x(end+1) =gdat_data.x(end)+1; - gdat_data.label{end+1}='P_{ic}'; + for i=1:length(fields_to_copy) + if isfield(ic,fields_to_copy{i}) + gdat_data.ic.(fields_to_copy{i}) = ic.(fields_to_copy{i}); + end + end + gdat_data.data(:,end+1) = interpos(-21,gdat_data.ic.t,gdat_data.ic.data,gdat_data.t); + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{ic}'; end end index_prad = []; @@ -1688,20 +1688,20 @@ elseif strcmp(mapping_for_aug.method,'switchcase') % radiated power params_eff.data_request='prad'; try - rad=gdat_aug(shot,params_eff); + rad=gdat_aug(shot,params_eff); catch end if ~isempty(rad.data) && ~isempty(rad.dim) - for i=1:length(fields_to_copy) - if isfield(rad,fields_to_copy{i}) - gdat_data.rad.(fields_to_copy{i}) = rad.(fields_to_copy{i}); - end - end + for i=1:length(fields_to_copy) + if isfield(rad,fields_to_copy{i}) + gdat_data.rad.(fields_to_copy{i}) = rad.(fields_to_copy{i}); + end + end % add to main with linear interpolation and 0 for extrapolated values - gdat_data.data(:,end+1) = interpos(-21,gdat_data.rad.t,gdat_data.rad.data,gdat_data.t); + gdat_data.data(:,end+1) = interpos(-21,gdat_data.rad.t,gdat_data.rad.data,gdat_data.t); index_prad = size(gdat_data.data,2); % to not add for ptot - gdat_data.x(end+1) =gdat_data.x(end)+1; - gdat_data.label{end+1}='P_{rad}'; + gdat_data.x(end+1) =gdat_data.x(end)+1; + gdat_data.label{end+1}='P_{rad}'; end end % add tot power @@ -1743,9 +1743,9 @@ elseif strcmp(mapping_for_aug.method,'switchcase') 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); + 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); + 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)'; @@ -1753,24 +1753,24 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.psi_lcfs(it)= psi_it(end); gdat_data.rhopolnorm(:,it) = sqrt(abs((psi_it-gdat_data.psi_axis(it)) ./(gdat_data.psi_lcfs(it)-gdat_data.psi_axis(it)))); if strcmp(DIAG,'EQR'); - % q value has only a few values and from center to edge, assume they are from central rhopol values on - % But they are every other point starting from 3rd - ijk=find(gdat_data.qvalue(:,it)~=0); - if length(ijk)>2 - % now shots have non-zero axis values in eqr - rhoeff=gdat_data.rhopolnorm(ijk,it); - qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above - if ijk(1)>1 - rhoeff = [0.; rhoeff]; - qeff = [qeff(1) ;qeff]; - end - ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); - qfit = zeros(size(gdat_data.rhopolnorm(:,it))); - qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); - else - qfit = zeros(size(gdat_data.rhopolnorm(:,it))); - end - gdat_data.qvalue(:,it) = qfit; + % q value has only a few values and from center to edge, assume they are from central rhopol values on + % But they are every other point starting from 3rd + ijk=find(gdat_data.qvalue(:,it)~=0); + if length(ijk)>2 + % now shots have non-zero axis values in eqr + rhoeff=gdat_data.rhopolnorm(ijk,it); + qeff=gdat_data.qvalue(ijk,it); % radial order was already inverted above + if ijk(1)>1 + rhoeff = [0.; rhoeff]; + qeff = [qeff(1) ;qeff]; + end + ij_nonan=find(~isnan(gdat_data.rhopolnorm(:,it))); + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + qfit(ij_nonan)=interpos(rhoeff,qeff,gdat_data.rhopolnorm(ij_nonan,it),-0.01,[1 0],[0 0],[300; ones(size(qeff(1:end-1)))]); + else + qfit = zeros(size(gdat_data.rhopolnorm(:,it))); + end + gdat_data.qvalue(:,it) = qfit; end % get rhotor values phi_it = phi_tree.value(it,Lpf1:-1:1)'; @@ -1785,13 +1785,13 @@ elseif strcmp(mapping_for_aug.method,'switchcase') q_edge=interpos(gdat_data.rhotornorm(1:end-1,it),qpsi.value(it,Lpf1:-1:2),1,-0.1); gdat_data.qvalue(:,it) = qpsi.value(it,Lpf1:-1:1)'; if abs(gdat_data.qvalue(end,it)) > 1e3 - % assume diverted - gdat_data.qvalue(end,it) = 2. * q_edge; + % assume diverted + gdat_data.qvalue(end,it) = 2. * q_edge; end % $$$ if qpsi.value(ijok(1),1)<0 -% $$$ gdat_data.qvalue(:,it) = max(qpsi.value(it,Lpf1:-1:1)',-max_qValue); +% $$$ 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); +% $$$ gdat_data.qvalue(:,it) = min(qpsi.value(it,Lpf1:-1:1)',max_qValue); % $$$ end end gdat_data.x = gdat_data.rhopolnorm; @@ -1820,11 +1820,11 @@ elseif strcmp(mapping_for_aug.method,'switchcase') [psi_tree,e]=rdaAUG_eff(shot,DIAG,'PFL',exp_name_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); 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 + % 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; + itotransposeback = 0; end psi_tree=adapt_rda(psi_tree,NTIME,ndimrho,itotransposeback); ijnan=find(isnan(psi_tree.value)); @@ -1832,14 +1832,14 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.dim{1} = gdat_data.t; gdat_data.dimunits{1} = 'time [s]'; for it=1:NTIME - Lpf1 =min(Lpf1_t(it),size(psi_tree.value,2)); - 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 + Lpf1 =min(Lpf1_t(it),size(psi_tree.value,2)); + 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]; @@ -1859,35 +1859,35 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.(sources_avail{i}).label=[]; end rap_signals={'ECE_f','Hmode','Ip','P_LH','Phib','beta','cs','exflag','g_f','it','li3','ne_f','res','tcomp','ECEm', ... - 'F','Vp','chie','dk_k','dpsi','dte','g1','g23or','g2','g3','iota','jbs','jec','jnb','jpar','ne','ni', ... - 'p','paux','pbrem','pec','pei','pnbe','poh','prad','psi','q','rhtECE','sdk_k','shear','signeo','sxk_k', ... - 'te','ti','upl','xk_k','ze'}; + 'F','Vp','chie','dk_k','dpsi','dte','g1','g23or','g2','g3','iota','jbs','jec','jnb','jpar','ne','ni', ... + 'p','paux','pbrem','pec','pei','pnbe','poh','prad','psi','q','rhtECE','sdk_k','shear','signeo','sxk_k', ... + 'te','ti','upl','xk_k','ze'}; rap_desc_signals={'rts:Dia/RAPTOR/ECEflag.val','rts:Dia/RAPTOR/Hmode.val','rts:Dia/RAPTOR/Ip.val','rts:Dia/RAPTOR/P_LH.val',... - 'rts:Dia/RAPTOR/Phib.val','rts:Dia/RAPTOR/beta.val','rts:Dia/RAPTOR/conf_state.val', ... - 'rts:Dia/RAPTOR/exitflag.val','rts:Dia/RAPTOR/gflag.val','rts:Dia/RAPTOR/it.val', ... - 'rts:Dia/RAPTOR/li3.val','rts:Dia/RAPTOR/neflag.val','rts:Dia/RAPTOR/res.val', ... - 'rts:Dia/RAPTOR/dtcomp.val','rts:Dia/RAPTOR/ECEmask.val','rts:Dia/RAPTOR/F.val', ... - 'rts:Dia/RAPTOR/Vp.val','rts:Dia/RAPTOR/chie.val','rts:Dia/RAPTOR/dk_k.val', ... - 'rts:Dia/RAPTOR/dpsi.val','rts:Dia/RAPTOR/dte.val','rts:Dia/RAPTOR/g1.val', ... - 'rts:Dia/RAPTOR/g23or.val','rts:Dia/RAPTOR/g2.val','rts:Dia/RAPTOR/g3.val', ... - 'rts:Dia/RAPTOR/iota.val','rts:Dia/RAPTOR/jbs.val','rts:Dia/RAPTOR/jec.val', ... - 'rts:Dia/RAPTOR/jnb.val','rts:Dia/RAPTOR/jpar.val','rts:Dia/RAPTOR/ne.val', ... - 'rts:Dia/RAPTOR/ni.val','rts:Dia/RAPTOR/p.val','rts:Dia/RAPTOR/paux.val', ... - 'rts:Dia/RAPTOR/pbrem.val','rts:Dia/RAPTOR/pec.val','rts:Dia/RAPTOR/pei.val', ... - 'rts:Dia/RAPTOR/pnbe.val','rts:Dia/RAPTOR/poh.val','rts:Dia/RAPTOR/prad.val', ... - 'rts:Dia/RAPTOR/psi.val','rts:Dia/RAPTOR/q.val','rts:Dia/RAPTOR/rhotorECE.val', ... - 'rts:Dia/RAPTOR/sdk_k.val','rts:Dia/RAPTOR/shear.val','rts:Dia/RAPTOR/signeo.val', ... - 'rts:Dia/RAPTOR/sxk_k.val','rts:Dia/RAPTOR/te.val','rts:Dia/RAPTOR/ti.val', ... - 'rts:Dia/RAPTOR/upl.val','rts:Dia/RAPTOR/xk_k.val','rts:Dia/RAPTOR/ze.val'}; + 'rts:Dia/RAPTOR/Phib.val','rts:Dia/RAPTOR/beta.val','rts:Dia/RAPTOR/conf_state.val', ... + 'rts:Dia/RAPTOR/exitflag.val','rts:Dia/RAPTOR/gflag.val','rts:Dia/RAPTOR/it.val', ... + 'rts:Dia/RAPTOR/li3.val','rts:Dia/RAPTOR/neflag.val','rts:Dia/RAPTOR/res.val', ... + 'rts:Dia/RAPTOR/dtcomp.val','rts:Dia/RAPTOR/ECEmask.val','rts:Dia/RAPTOR/F.val', ... + 'rts:Dia/RAPTOR/Vp.val','rts:Dia/RAPTOR/chie.val','rts:Dia/RAPTOR/dk_k.val', ... + 'rts:Dia/RAPTOR/dpsi.val','rts:Dia/RAPTOR/dte.val','rts:Dia/RAPTOR/g1.val', ... + 'rts:Dia/RAPTOR/g23or.val','rts:Dia/RAPTOR/g2.val','rts:Dia/RAPTOR/g3.val', ... + 'rts:Dia/RAPTOR/iota.val','rts:Dia/RAPTOR/jbs.val','rts:Dia/RAPTOR/jec.val', ... + 'rts:Dia/RAPTOR/jnb.val','rts:Dia/RAPTOR/jpar.val','rts:Dia/RAPTOR/ne.val', ... + 'rts:Dia/RAPTOR/ni.val','rts:Dia/RAPTOR/p.val','rts:Dia/RAPTOR/paux.val', ... + 'rts:Dia/RAPTOR/pbrem.val','rts:Dia/RAPTOR/pec.val','rts:Dia/RAPTOR/pei.val', ... + 'rts:Dia/RAPTOR/pnbe.val','rts:Dia/RAPTOR/poh.val','rts:Dia/RAPTOR/prad.val', ... + 'rts:Dia/RAPTOR/psi.val','rts:Dia/RAPTOR/q.val','rts:Dia/RAPTOR/rhotorECE.val', ... + 'rts:Dia/RAPTOR/sdk_k.val','rts:Dia/RAPTOR/shear.val','rts:Dia/RAPTOR/signeo.val', ... + 'rts:Dia/RAPTOR/sxk_k.val','rts:Dia/RAPTOR/te.val','rts:Dia/RAPTOR/ti.val', ... + 'rts:Dia/RAPTOR/upl.val','rts:Dia/RAPTOR/xk_k.val','rts:Dia/RAPTOR/ze.val'}; if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) gdat_data.gdat_params.source = sources_avail; elseif ~iscell(gdat_data.gdat_params.source) if ischar(gdat_data.gdat_params.source) - gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source); - if strcmp(gdat_data.gdat_params.source(1),'o'); gdat_data.gdat_params.source = 'observer'; end - if strcmp(gdat_data.gdat_params.source(1),'p'); gdat_data.gdat_params.source = 'predictive'; end + gdat_data.gdat_params.source = lower(gdat_data.gdat_params.source); + if strcmp(gdat_data.gdat_params.source(1),'o'); gdat_data.gdat_params.source = 'observer'; end + if strcmp(gdat_data.gdat_params.source(1),'p'); gdat_data.gdat_params.source = 'predictive'; end if ~any(strmatch(gdat_data.gdat_params.source,lower(sources_avail))) if (gdat_params.nverbose>=1) warning(['source= ' gdat_data.gdat_params.source ' is not part of the available sources: ' sprintf('''%s'' ',sources_avail{:})]); @@ -1903,8 +1903,8 @@ elseif strcmp(mapping_for_aug.method,'switchcase') else for i=1:length(gdat_data.gdat_params.source) gdat_data.gdat_params.source{i} = lower(gdat_data.gdat_params.source{i}); - if strcmp(gdat_data.gdat_params.source{i}(1),'o'); gdat_data.gdat_params.source{i} = 'observer'; end - if strcmp(gdat_data.gdat_params.source{i}(1),'p'); gdat_data.gdat_params.source{i} = 'predictive'; end + if strcmp(gdat_data.gdat_params.source{i}(1),'o'); gdat_data.gdat_params.source{i} = 'observer'; end + if strcmp(gdat_data.gdat_params.source{i}(1),'p'); gdat_data.gdat_params.source{i} = 'predictive'; end if ~any(strmatch(gdat_data.gdat_params.source{i},lower(sources_avail))) if gdat_data.gdat_params.nverbose>=1 warning(['source = ' gdat_data.gdat_params.source{i} ' not expected with data_request= ' data_request_eff]) @@ -1916,23 +1916,23 @@ elseif strcmp(mapping_for_aug.method,'switchcase') for i=1:length(gdat_data.gdat_params.source) it_def=0; for j=1:length(rap_signals) - gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev; - gdat_params_eff.data_request = {'RAP',[rap_signals{j} '_' gdat_data.gdat_params.source{i}(1)],gdat_params_eff.exp_name}; - abcd = gdat_aug(gdat_data.shot,gdat_params_eff); - gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}) = abcd; - gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}).description = rap_desc_signals{j}; - if strcmp(rap_signals{j},'q') && ~isempty(abcd.data) - % promote q at top level of gdat_data.(gdat_data.gdat_params.source{i}) - gdat_data.(gdat_data.gdat_params.source{i}).data = abcd.data; - gdat_data.(gdat_data.gdat_params.source{i}).t = abcd.t; - gdat_data.(gdat_data.gdat_params.source{i}).x = abcd.x; - gdat_data.(gdat_data.gdat_params.source{i}).units = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho']; - gdat_data.(gdat_data.gdat_params.source{i}).dim = abcd.dim; - gdat_data.(gdat_data.gdat_params.source{i}).dimunits = abcd.dimunits; - gdat_data.(gdat_data.gdat_params.source{i}).dimunits{1} = 'various rhotor'; - gdat_data.(gdat_data.gdat_params.source{i}).data_fullpath = abcd.data_fullpath; - gdat_data.(gdat_data.gdat_params.source{i}).label = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho']; - end + gdat_params_prev = gdat_data.gdat_params; gdat_params_eff = gdat_params_prev; + gdat_params_eff.data_request = {'RAP',[rap_signals{j} '_' gdat_data.gdat_params.source{i}(1)],gdat_params_eff.exp_name}; + abcd = gdat_aug(gdat_data.shot,gdat_params_eff); + gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}) = abcd; + gdat_data.(gdat_data.gdat_params.source{i}).(rap_signals{j}).description = rap_desc_signals{j}; + if strcmp(rap_signals{j},'q') && ~isempty(abcd.data) + % promote q at top level of gdat_data.(gdat_data.gdat_params.source{i}) + gdat_data.(gdat_data.gdat_params.source{i}).data = abcd.data; + gdat_data.(gdat_data.gdat_params.source{i}).t = abcd.t; + gdat_data.(gdat_data.gdat_params.source{i}).x = abcd.x; + gdat_data.(gdat_data.gdat_params.source{i}).units = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho']; + gdat_data.(gdat_data.gdat_params.source{i}).dim = abcd.dim; + gdat_data.(gdat_data.gdat_params.source{i}).dimunits = abcd.dimunits; + gdat_data.(gdat_data.gdat_params.source{i}).dimunits{1} = 'various rhotor'; + gdat_data.(gdat_data.gdat_params.source{i}).data_fullpath = abcd.data_fullpath; + gdat_data.(gdat_data.gdat_params.source{i}).label = ['q RAPTOR_' gdat_data.gdat_params.source{i}(1:3) ' various rho']; + end end end % @@ -1985,14 +1985,14 @@ elseif strcmp(mapping_for_aug.method,'switchcase') switch data_request_eff case {'volume_rho', 'rhovol'} for it=1:NTIME - Lpf1 = Lpf1_t(it); - psi_it=psi_tree.value(it,Lpf1:-1:1)'; - psi_axis(it)= psi_it(1); - psi_lcfs(it)= psi_it(end); - gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); - gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)'; - % gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi - gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it))); + Lpf1 = Lpf1_t(it); + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + psi_axis(it)= psi_it(1); + psi_lcfs(it)= psi_it(end); + gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); + gdat_data.vol(:,it)=Vol.value(it,2*Lpf1-1:-2:1)'; + % gdat_data.dvoldpsi(:,it)=Vol.value(it,2*Lpf1:-2:2)'; % 2nd index are dV/dpsi + gdat_data.rhovolnorm(:,it) = sqrt(abs(gdat_data.vol(:,it) ./ gdat_data.vol(end,it))); end gdat_data.dim{1} = gdat_data.x; gdat_data.dim{2} = gdat_data.t; @@ -2000,46 +2000,46 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.dimunits{2} = 'time [s]'; gdat_data.volume_edge = gdat_data.vol(end,:); if strcmp(data_request_eff,'volume_rho') - gdat_data.data = gdat_data.vol; - gdat_data.units = Vol.units; + gdat_data.data = gdat_data.vol; + gdat_data.units = Vol.units; else strcmp(data_request_eff,'rhovol') - gdat_data.data = gdat_data.rhovolnorm; - gdat_data.units = ' '; + gdat_data.data = gdat_data.rhovolnorm; + gdat_data.units = ' '; end case {'rhotor', 'rhotor_edge', 'rhotor_norm'} b0=gdat(shot,'b0'); gdat_data.b0 = interpos(b0.t,b0.data,gdat_data.t,-1); for it=1:NTIME - Lpf1 = Lpf1_t(it); - gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)'; - gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it))); - gdat_data.rhotor(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.b0(it) ./ pi)); - psi_it=psi_tree.value(it,Lpf1:-1:1)'; - psi_axis(it)= psi_it(1); - psi_lcfs(it)= psi_it(end); - gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); + Lpf1 = Lpf1_t(it); + gdat_data.phi(:,it) = phi_tree.value(it,Lpf1:-1:1)'; + gdat_data.rhotornorm(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.phi(end,it))); + gdat_data.rhotor(:,it) = sqrt(abs(gdat_data.phi(:,it) ./ gdat_data.b0(it) ./ pi)); + psi_it=psi_tree.value(it,Lpf1:-1:1)'; + psi_axis(it)= psi_it(1); + psi_lcfs(it)= psi_it(end); + gdat_data.x(:,it) = sqrt(abs((psi_it-psi_axis(it)) ./(psi_lcfs(it)-psi_axis(it)))); end % always provide rhotor_edge so field always exists gdat_data.rhotor_edge = gdat_data.rhotor(end,:); if strcmp(data_request_eff,'rhotor') - gdat_data.data = gdat_data.rhotor; - gdat_data.units = 'm'; - gdat_data.dim{1} = gdat_data.x; - gdat_data.dim{2} = gdat_data.t; - gdat_data.dimunits{1} = 'rhopolnorm'; - gdat_data.dimunits{2} = 'time [s]'; + gdat_data.data = gdat_data.rhotor; + gdat_data.units = 'm'; + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; elseif strcmp(data_request_eff,'rhotor_edge') - gdat_data.data = gdat_data.rhotor(end,:); - gdat_data.units = 'm'; - gdat_data.dim{1} = gdat_data.t; - gdat_data.dimunits{1} = 'time [s]'; + gdat_data.data = gdat_data.rhotor(end,:); + gdat_data.units = 'm'; + gdat_data.dim{1} = gdat_data.t; + gdat_data.dimunits{1} = 'time [s]'; elseif strcmp(data_request_eff,'rhotor_norm') - gdat_data.data = gdat_data.rhotornorm; - gdat_data.units = ' '; - gdat_data.dim{1} = gdat_data.x; - gdat_data.dim{2} = gdat_data.t; - gdat_data.dimunits{1} = 'rhopolnorm'; - gdat_data.dimunits{2} = 'time [s]'; + gdat_data.data = gdat_data.rhotornorm; + gdat_data.units = ' '; + gdat_data.dim{1} = gdat_data.x; + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits{1} = 'rhopolnorm'; + gdat_data.dimunits{2} = 'time [s]'; end end gdat_data.data_fullpath = [DIAG '/PFL extract in .data: ' data_request_eff]; @@ -2053,7 +2053,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase') if ~isfield(gdat_data.gdat_params,'source') || isempty(gdat_data.gdat_params.source) gdat_data.gdat_params.source = 'G'; elseif length(gdat_data.gdat_params.source)>1 || ... - upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M' + upper(gdat_data.gdat_params.source)<'F' || upper(gdat_data.gdat_params.source)>'M' if gdat_data.gdat_params.nverbose>=1 warning(['source = ' gdat_data.gdat_params.source ' not expected with data_request= ' data_request_eff]) end @@ -2070,7 +2070,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase') dim1_len='nb_of_chords'; else if gdat_data.gdat_params.nverbose>=1 - warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff]) + warning(['camera = ' gdat_data.gdat_params.camera ' not expected with data_request= ' data_request_eff]) end return end @@ -2081,27 +2081,27 @@ elseif strcmp(mapping_for_aug.method,'switchcase') gdat_data.gdat_params.time_interval = []; end [aa,bb]=unix(['ssh ' '-o "StrictHostKeyChecking no" sxaug20.aug.ipp.mpg.de WhichSX ' num2str(shot) ' ' ... - upper(gdat_data.gdat_params.source)]); + upper(gdat_data.gdat_params.source)]); chords_ok_diags=regexpi(bb,'(?<chord>[A-Z]_[0-9]+)\saddress:\s[0-9]+\sDiag:\s(?<diag>[A-Z]+)','names'); if isempty(chords_ok_diags) chords_ok_nb = []; else for i=1:length(chords_ok_diags) - chords_ok_nb(i) = str2num(chords_ok_diags(i).chord(3:end)); + chords_ok_nb(i) = str2num(chords_ok_diags(i).chord(3:end)); end end if isempty(gdat_data.gdat_params.camera); gdat_data.gdat_params.camera = chords_ok_nb; else for i=1:length(gdat_data.gdat_params.camera) - ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i)); - if ~isempty(ij) - chords_ok_diags_tmp(i) = chords_ok_diags(ij); - else - % chord not in use - chords_ok_diags_tmp(i).chord = []; - chords_ok_diags_tmp(i).diag = []; - end + ij=find(chords_ok_nb==gdat_data.gdat_params.camera(i)); + if ~isempty(ij) + chords_ok_diags_tmp(i) = chords_ok_diags(ij); + else + % chord not in use + chords_ok_diags_tmp(i).chord = []; + chords_ok_diags_tmp(i).diag = []; + end end chords_ok_diags = chords_ok_diags_tmp; chords_ok_nb = gdat_data.gdat_params.camera; @@ -2116,49 +2116,49 @@ elseif strcmp(mapping_for_aug.method,'switchcase') for i=1:length(gdat_data.gdat_params.camera) if ischar(gdat_data.gdat_params.freq) && strcmp(gdat_data.gdat_params.freq,'slow'); chords_ok_diags(i).diag = 'SSX'; end if ~isempty(chords_ok_diags(i).diag) && ~isempty(chords_ok_diags(i).chord) - [a,e]=rdaAUG_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff, ... - gdat_data.gdat_params.time_interval,gdat_data.gdat_params.extra_arg_sf2sig); + [a,e]=rdaAUG_eff(shot,chords_ok_diags(i).diag,chords_ok_diags(i).chord,exp_name_eff, ... + gdat_data.gdat_params.time_interval,gdat_data.gdat_params.extra_arg_sf2sig); else - a.data = []; - a.t = []; + a.data = []; + a.t = []; end if ~isempty(a.data) - icount = icount + 1; - if icount == 1 - % first time has data - % assume all chords have same time base even if from different shotfile - % time missing one point - if length(a.value) == length(a.t)+1 - a.t=linspace(a.range(1),a.range(2),length(a.value)); - a_time.index(2) = length(a.value); - end - gdat_data.t = a.t(1:nnth:end); - gdat_data.units = a.units; - if strcmp(dim1_len,'from_chord_nb') - gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t)); - gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central - else - gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t)); - gdat_data.dim{1} = chords_ok_nb; - end - gdat_data.dim{2} = gdat_data.t; - gdat_data.dimunits = [{'chord nb'}; {'s'}]; - gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ... - setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')']; - gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)]; - end - try - if strcmp(dim1_len,'from_chord_nb') - gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end); - else - gdat_data.data(i,:) = a.data(1:nnth:end); - end - catch - if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end + icount = icount + 1; + if icount == 1 + % first time has data + % assume all chords have same time base even if from different shotfile + % time missing one point + if length(a.value) == length(a.t)+1 + a.t=linspace(a.range(1),a.range(2),length(a.value)); + a_time.index(2) = length(a.value); + end + gdat_data.t = a.t(1:nnth:end); + gdat_data.units = a.units; + if strcmp(dim1_len,'from_chord_nb') + gdat_data.data = NaN*ones(max(chords_ok_nb),length(gdat_data.t)); + gdat_data.dim{1} = [1:max(chords_ok_nb)]; % simpler to have index corresponding to chord number, except for central + else + gdat_data.data = NaN*ones(length(chords_ok_nb),length(gdat_data.t)); + gdat_data.dim{1} = chords_ok_nb; + end + gdat_data.dim{2} = gdat_data.t; + gdat_data.dimunits = [{'chord nb'}; {'s'}]; + gdat_data.data_fullpath = ['sxr from source=''' gdat_data.gdat_params.source ' uses shotfiles SX' ... + setdiff(unique(strvcat(chords_ok_diags(:).diag)),'SX')']; + gdat_data.label = ['SXR/' upper(gdat_data.gdat_params.source)]; + end + try + if strcmp(dim1_len,'from_chord_nb') + gdat_data.data(chords_ok_nb(i),:) = a.data(1:nnth:end); + else + gdat_data.data(i,:) = a.data(1:nnth:end); + end + catch + if (gdat_params.nverbose>=1); disp(['problem with ' chords_ok_diags(i).diag ' ; ' chords_ok_diags(i).chord]); end - end + end else - % add fields not yet defined in case all cases have empty data + % add fields not yet defined in case all cases have empty data end end gdat_data.chords = chords_ok_nb; @@ -2178,29 +2178,29 @@ elseif strcmp(mapping_for_aug.method,'switchcase') TRANSP_signals; for i=1:size(transp_sig,1) if strcmp(lower(transp_sig{i,2}),'signal') || strcmp(lower(transp_sig{i,2}),'signal-group') - try - abcd = rdaAUG_eff(shot,diagname,transp_sig{i,1},shotfile_exp_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); - eval(['gdat_data.' transp_sig{i,1} '= abcd;']); - % eval(['[gdat_data.' transp_sig{i,1} ',e]=rdaAUG_eff(shot,diagname,''' transp_sig{i,1} ''',shotfile_exp_eff);']); - catch - eval(['gdat_data.' transp_sig{i,1} '=[];']); - end + try + abcd = rdaAUG_eff(shot,diagname,transp_sig{i,1},shotfile_exp_eff,[],gdat_data.gdat_params.extra_arg_sf2sig); + eval(['gdat_data.' transp_sig{i,1} '= abcd;']); + % eval(['[gdat_data.' transp_sig{i,1} ',e]=rdaAUG_eff(shot,diagname,''' transp_sig{i,1} ''',shotfile_exp_eff);']); + catch + eval(['gdat_data.' transp_sig{i,1} '=[];']); + end elseif strcmp(lower(transp_sig{i,2}),'area-base') - clear adata_area - try - eval(['[adata_area]=sf2ab(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']); - catch - adata_area.value = cell(0); - end - eval(['gdat_data.' transp_sig{i,1} '=adata_area;']); + clear adata_area + try + eval(['[adata_area]=sf2ab(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']); + catch + adata_area.value = cell(0); + end + eval(['gdat_data.' transp_sig{i,1} '=adata_area;']); elseif strcmp(lower(transp_sig{i,2}),'time-base') - clear adata_time - try - eval(['[adata_time]=sf2tb(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff,shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']); - catch - adata_time.value = cell(0); - end - eval(['gdat_data.' transp_sig{i,1} '=adata_time;']); + clear adata_time + try + eval(['[adata_time]=sf2tb(diagname,shot,transp_sig{i,1},''-exp'',shotfile_exp_eff,shotfile_exp_eff' extra_arg_sf2sig_eff_string ');']); + catch + adata_time.value = cell(0); + end + eval(['gdat_data.' transp_sig{i,1} '=adata_time;']); end end % copy TIME to .t