Skip to content
Snippets Groups Projects
Commit fe73c8fe authored by Olivier Sauter's avatar Olivier Sauter
Browse files

fix equil time intervals for calculating rhotor etc for diagnostics like cxrs, ece, thomson 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@5598 d63d8f72-b253-0410-a779-e742ad2e26cf
parent 89b1216d
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
crpptbx/AUG/CEZ_CMZ_fit.m
+ 47
− 32
View file @ fe73c8fe
......@@ -10,18 +10,29 @@ if exist(filename_data,'file')
else
cxrs=gdat(shot,'cxrs_rho','doplot',1);
cmz=gdat(shot,'cxrs_rho','doplot',1,'source','cmz');
eval(['save ' filename_data ' cxrs cmz'])
end
cez_cmz_data = [];
cez_cmz_fit = [];
main_time_base = cxrs.t;
ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
if length(ij)==1
disp('not yet ready')
if isempty(cxrs.data) || isempty(cxrs.t) || ischar(cxrs.data)
disp('no main cxrs')
return
end
if ~exist(filename_data,'file')
eval(['save ' filename_data ' cxrs cmz'])
end
main_time_base = cxrs.t;
if ~exist('time_interval_in') || isempty(time_interval_in)
time_interval_in = [cxrs.t(1) cxrs.t(end)];
ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
else
ij=find(main_time_base>=time_interval_in(1) & main_time_base<=time_interval_in(2));
if length(ij)==1
disp('not yet ready')
return
end
end
cez_cmz_data.time = main_time_base(ij);
for i=1:length(ij)
......@@ -103,32 +114,36 @@ for i=1:length(ij)
end
[rhoeff,irhoeff] = sort(rhotor_data_tofit);
tidata_tofit_eff = tidata_tofit(irhoeff);
vrotdata_tofit_eff = vrotdata_tofit(irhoeff);
tierr_tofit_eff = tierr_tofit(irhoeff);
vroterr_tofit_eff = vroterr_tofit(irhoeff);
cez_cmz_data.ti{i}.rhotornorm = rhoeff;
cez_cmz_data.ti{i}.data = tidata_tofit_eff;
cez_cmz_data.ti{i}.err = tierr_tofit_eff;
cez_cmz_data.vrot{i}.rhotornorm = rhoeff;
cez_cmz_data.vrot{i}.data = vrotdata_tofit_eff;
cez_cmz_data.vrot{i}.err = vroterr_tofit_eff;
xeff = [0. rhoeff];
yeffti = [tidata_tofit_eff(1) tidata_tofit_eff];
yeffvrot = [vrotdata_tofit_eff(1) vrotdata_tofit_eff];
erreffti = [100.*max(tierr_tofit_eff) tierr_tofit_eff];
erreffvrot = [100.*max(vroterr_tofit_eff) vroterr_tofit_eff];
[tifit(:,i),dtifitdrn(:,i)]=interpos(xeff,yeffti,rhotorfit,tension_ti_eff,[1 0],[0 0],erreffti);
[vrotfit(:,i),dvrotfitdrn(:,i)]=interpos(xeff,yeffvrot,rhotorfit,tension_vrot_eff,[1 0],[0 0],erreffvrot);
if doplot
figure(11);clf
errorbar(rhoeff,tidata_tofit_eff,tierr_tofit_eff,'*');
hold all
plot(rhotorfit,tifit(:,i))
figure(12);clf
errorbar(rhoeff,vrotdata_tofit_eff,vroterr_tofit_eff,'*');
hold all
plot(rhotorfit,vrotfit(:,i))
pause(0.01)
if isempty(tidata_tofit_eff)
keyboard
else
vrotdata_tofit_eff = vrotdata_tofit(irhoeff);
tierr_tofit_eff = tierr_tofit(irhoeff);
vroterr_tofit_eff = vroterr_tofit(irhoeff);
cez_cmz_data.ti{i}.rhotornorm = rhoeff;
cez_cmz_data.ti{i}.data = tidata_tofit_eff;
cez_cmz_data.ti{i}.err = tierr_tofit_eff;
cez_cmz_data.vrot{i}.rhotornorm = rhoeff;
cez_cmz_data.vrot{i}.data = vrotdata_tofit_eff;
cez_cmz_data.vrot{i}.err = vroterr_tofit_eff;
xeff = [0. rhoeff];
yeffti = [tidata_tofit_eff(1) tidata_tofit_eff];
yeffvrot = [vrotdata_tofit_eff(1) vrotdata_tofit_eff];
erreffti = [100.*max(tierr_tofit_eff) tierr_tofit_eff];
erreffvrot = [100.*max(vroterr_tofit_eff) vroterr_tofit_eff];
[tifit(:,i),dtifitdrn(:,i)]=interpos(xeff,yeffti,rhotorfit,tension_ti_eff,[1 0],[0 0],erreffti);
[vrotfit(:,i),dvrotfitdrn(:,i)]=interpos(xeff,yeffvrot,rhotorfit,tension_vrot_eff,[1 0],[0 0],erreffvrot);
if doplot
figure(11);clf
errorbar(rhoeff,tidata_tofit_eff,tierr_tofit_eff,'*');
hold all
plot(rhotorfit,tifit(:,i))
figure(12);clf
errorbar(rhoeff,vrotdata_tofit_eff,vroterr_tofit_eff,'*');
hold all
plot(rhotorfit,vrotfit(:,i))
pause(0.01)
end
end
end
......
crpptbx/AUG/gdat_aug.m
+ 9
− 9
View file @ fe73c8fe
......@@ -479,20 +479,20 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
end
gdat_data.gdat_params.equil = params_equil.equil;
gdat_data.equil = equil;
inb_chord_cxrs=size(gdat_data.r,1);
inb_time_cxrs=size(gdat_data.r,2);
inb_chord_cxrs=size(gdat_data.data,1);
inb_time_cxrs=size(gdat_data.data,2);
psi_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhopolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhotornorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
rhovolnorm_out = NaN*ones(inb_chord_cxrs,inb_time_cxrs);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
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));
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);
......@@ -701,12 +701,12 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
rhotornorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
rhovolnorm_out = NaN*ones(inb_chord_ece,inb_time_ece);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-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));
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);
......@@ -1089,12 +1089,12 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
rhotornorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
rhovolnorm_out_edge = NaN*ones(inb_chord_edge,inb_time_edge);
% constructs intervals within which a given equil is used: [time_equil(i),time_equil(i+1)]
time_equil=[1.5*equil.t(1)-0.5*equil.t(2) 0.5.*(equil.t(1:end-1)+equil.t(2:end)) 1.5*equil.t(end)-0.5*equil.t(end-1)];
time_equil=[min(gdat_data.(lower(node_child_nameeff)).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.(lower(node_child_nameeff)).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_core_inequil = find(gdat_data.(lower(node_child_nameeff)).t>=time_equil(itequil) & gdat_data.(lower(node_child_nameeff)).t<=time_equil(itequil+1));
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);
......@@ -1109,7 +1109,7 @@ elseif strcmp(mapping_for_aug.method,'switchcase')
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));
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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment