diff --git a/crpptbx/AUG/CEZ_CMZ_fit.m b/crpptbx/AUG/CEZ_CMZ_fit.m index 43b7814fb6b1c204bfbe3aeeca76e8962427b99b..ace671e8304a315ea9efa1c685d65c3e4dc7dba0 100644 --- a/crpptbx/AUG/CEZ_CMZ_fit.m +++ b/crpptbx/AUG/CEZ_CMZ_fit.m @@ -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 diff --git a/crpptbx/AUG/gdat_aug.m b/crpptbx/AUG/gdat_aug.m index a9f3c7dd6c3131b972d5d2a83bb04e342f02be5f..708413db8e5dde115970db8dca32d7b5fc40476e 100644 --- a/crpptbx/AUG/gdat_aug.m +++ b/crpptbx/AUG/gdat_aug.m @@ -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);