Finish using nodes for liuqe quantities

matlab/TCV_IMAS/tcv_get_ids_core_sources.m

@@ -49,31 +49,40 @@ ohm_gdat = gdat(params_core_sources.shot,params_eff);
 params_eff = params_eff_ref; params_eff.data_request='powers';
 powers_gdat = gdat(params_core_sources.shot,params_eff);
 
-%% load liuqe data from nodes for conversions 
-params_eff = params_eff_ref;
-% normalized sqrt poloidal flux including axis
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
-rho_pol_norm_liu = gdat(shot,params_eff);
-% <1/R^2>
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
-Rm2_fs = gdat(shot,params_eff);
-% Volume
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
-vol= gdat(shot,params_eff);
-% T(\psi) = R*B_tor
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
-T = gdat(shot,params_eff);
-% Q1Q: -dpsi/dV
-%mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:DPSI_VOL';
-Vprime = gdat(shot,params_eff); 
-% I_tor = I_p, total plasma current
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:i_pl';
-Ip = gdat(shot,params_eff);
-
+%% load liuqe data from nodes for conversions
+if params_eff_ref.liuqe == 1
+  nodes_liuqe = 'equil';
+else
+  nodes_liuqe = ['equil',num2str(params_eff_ref.liuqe)];
+end
+if check_nodes_filled(params_core_sources.shot,nodes_liuqe)
+  params_eff = params_eff_ref;
+  % normalized sqrt poloidal flux including axis
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
+  rho_pol_norm_liu = gdat(params_core_sources.shot,params_eff);
+  % <1/R^2>
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
+  Rm2_fs_liu = gdat(params_core_sources.shot,params_eff);
+  % Volume
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
+  vol_liu= gdat(params_core_sources.shot,params_eff);
+  % T(\psi) = R*B_tor
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
+  T_liu = gdat(params_core_sources.shot,params_eff);
+  % safety factor q = -dPhi/dPsi
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:q';
+  q_liu = gdat(params_core_sources.shot,params_eff);
+  % I_tor = I_p, total plasma current
+  params_eff.data_request='\tcv_shot::top.results.equil_1.results:i_pl';
+  Ip_liu = gdat(params_core_sources.shot,params_eff);
+  % define liuqe_time from Ip_liu
+  liuqe_time = Ip_liu.t;
+else
+  error('Liuqe nodes %s not filled. Contact O. Sauter.')
+end
 %% initialize source from template
 % ohm
-ohm_t_grid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_t_grid);
+ohm_tgrid = powers_gdat.ohm.t; ohm_n_t = numel(ohm_tgrid);
 
 main_desc = 'Source from ohmic heating'; production = 'IBS nodes';
 id_ohm.description = sprintf('%s from %s',main_desc,production);
@@ -86,35 +95,24 @@ ohm_data = ohm_gdat.ohm.ohm_data;
 
 % load LIUQE data to convert
 % \tilde{j}_// = <jdotB>/(R0<Bphi/R>) to j_//0 = <jdotB>/B0; Bphi=T(psi)/R
-[L,~,LY] = liuqe(shot,ohm_t_grid,'iterq',50,'ilim',3,'icsint',true);
-params_eff = params_eff_ref;
-% normalized sqrt poloidal flux including axis
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rho';
-rho_pol_norm_liu = gdat(shot,params_eff);
-% <1/R^2>
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
-Rm2_fs = gdat(shot,params_eff);
-% Volume
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:vol';
-vol= gdat(shot,params_eff);
-% T(\psi) = R*B_tor
-params_eff.data_request='\tcv_shot::top.results.equil_1.results:rbtor_rho';
-T = gdat(shot,params_eff);
-%rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
-%T      = LY.TQ;
-%Rm2_fs = LY.Q2Q;
-%vol    = LY.VQ;
+% interpolate liuqe outputs on ohm_tgrid
+T      = interp1(liuqe_time,T_liu.data.',     ohm_tgrid)';
+Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ohm_tgrid)';
+vol    = interp1(liuqe_time,vol_liu.data.',   ohm_tgrid)';
+% get vacuum field data from ids
 R0 = ids_core_sources.vacuum_toroidal_field.r0;
-B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ohm_t_grid);
-jpar_tilde_to_jpar_0 = R0*T.data.*Rm2_fs.data./B0';
+B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ohm_tgrid);
+
+jpar_tilde_to_jpar_0 = R0*T.*Rm2_fs./B0'; 
 
 % map volume and conversion factor on rho_pol grid of ohm_data cd_dens
-jpar_tilde_to_jpar_0_mapped = interp1(rho_pol_norm_liu,jpar_tilde_to_jpar_0,ohm_data.cd_dens.rhopol_norm');
-vol_mapped = interp1(rho_pol_norm_liu,vol,ohm_data.cd_dens.rhopol_norm');
+jpar_tilde_to_jpar_0_mapped = ...
+  interp1(rho_pol_norm_liu.data,jpar_tilde_to_jpar_0,ohm_data.cd_dens.rhopol_norm);
+vol_mapped = interp1(rho_pol_norm_liu.data,vol,ohm_data.cd_dens.rhopol_norm);
 
 for ii = 1:ohm_n_t
   % profiles_1d
-  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = ohm_t_grid(ii);
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = ohm_tgrid(ii);
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
     ohm_data.cd_dens.rhopol_norm';
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
@@ -125,7 +123,7 @@ for ii = 1:ohm_n_t
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
 
   % globals
-  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_t_grid(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = ohm_tgrid(ii);
   ids_core_sources.source{last_index+1}.global_quantities{ii}.power = powers_gdat.ohm.data(ii);
   ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
 end
@@ -135,7 +133,7 @@ last_index = last_index+1;  % add if statement to only increment if ohmic source
 %% bs
 params_eff = params_eff_ref; params_eff.data_request='bs_data';
 bs_gdat = gdat(params_core_sources.shot,params_eff); bs_data = bs_gdat.bs.bs_data;
-bs_t_grid = bs_gdat.bs.t; bs_n_t = numel(bs_t_grid);
+bs_tgrid = bs_gdat.bs.t; bs_n_t = numel(bs_tgrid);
 
 main_desc = 'Bootstrap current'; production = 'IBS nodes';
 id_bs.description = sprintf('%s from %s',main_desc,production);
@@ -147,7 +145,7 @@ ids_core_sources.source{last_index+1}.global_quantities(1:bs_n_t) = {globals_tem
 
 for ii = 1:bs_n_t
   % profiles_1d
-  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = bs_t_grid(ii);
+  ids_core_sources.source{last_index+1}.profiles_1d{ii}.time = bs_tgrid(ii);
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_pol_norm = ...
     bs_data.cd_dens.rhopol_norm';
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.grid.rho_tor_norm = ...
@@ -159,7 +157,7 @@ for ii = 1:bs_n_t
   ids_core_sources.source{last_index+1}.profiles_1d{ii}.current_parallel_inside = integrated_jpar_0_tmp;
 
   % globals
-  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_t_grid(ii);
+  ids_core_sources.source{last_index+1}.global_quantities{ii}.time = bs_tgrid(ii);
   ids_core_sources.source{last_index+1}.global_quantities{ii}.current_parallel = integrated_jpar_0_tmp(end);
 end
 
@@ -187,11 +185,11 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   % find times where EC is on to define time grid with extra time slice just
   % before/after  EC power and at start/end of shot
   itime_ec = find(powers_gdat.ec.data(:,end)>0);
-  if ec_powers_tgrid(itime_ec(end))>=ohm_t_grid(end)
-    i_time_end = iround(ec_powers_tgrid,ohm_t_grid(end));
-    ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
+  if ec_powers_tgrid(itime_ec(end))>=ohm_tgrid(end)
+    i_time_end = iround(ec_powers_tgrid,ohm_tgrid(end));
+    ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:i_time_end)'];
   else
-    ec_tgrid_out = [ohm_t_grid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_t_grid(end)];
+    ec_tgrid_out = [ohm_tgrid(1),ec_powers_tgrid(itime_ec(1)-1:itime_ec(end)+1)',ohm_tgrid(end)];
   end
   nt_ec_out = numel(ec_tgrid_out);
   p_ec_injected = interpos(ec_powers_tgrid,powers_gdat.ec.data(:,end),ec_tgrid_out);
@@ -210,23 +208,21 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
 
   % load LIUQE data to convert
   % \tilde{V,TORAY}_// = dI_\phi/dV = (1/2\pi)<j_\phi/R> to j_//0 = <jdotB>/B0;
-  % conversion using j_//,cd = <j_cd/B> B, thus
-  % j_// = 2\pi <B^2> / (R_0 <B_\phi/R>^2)
-  [L,~,LY] = liuqe(shot,ec_tgrid_toray,'iterq',50,'ilim',3,'icsint',true);
-  rho_pol_norm_liu = L.pQ; % normalized sqrt poloidal flux including axis
+  % conversion using j_//,cd = <j_cd/B> B, thus j_// = 2\pi <B^2> / (B_0 <B_\phi/R>)
   mu0     = 4*pi*10^-7;
-  Rm2_fs  = LY.Q2Q;
-  params_eff = params_eff_ref;
-  params_eff.data_request='\tcv_shot::top.results.equil_1.results:gpsi0_r2_fsa';
-  Rm2_fs = gdat(shot,params_eff);
+  % interpolate liuqe outputs on ohm_tgrid
+  T      = interp1(liuqe_time,T_liu.data.',     ec_tgrid_toray)';
+  Rm2_fs = interp1(liuqe_time,Rm2_fs_liu.data.',ec_tgrid_toray)';
+  vol    = interp1(liuqe_time,vol_liu.data.',   ec_tgrid_toray)';
+  Ip     = interp1(liuqe_time,Ip_liu.data,      ec_tgrid_toray); 
+  q      = interp1(liuqe_time,q_liu.data.',     ec_tgrid_toray)';   
+  mVprime = -2*pi*q./T./Rm2_fs;
+  B2_fs   = mu0*Ip./mVprime + T.^2.*Rm2_fs;
   
-  mVprime = 1./LY.Q1Q; % Q1Q: -dpsi/dV
-  Ip_psi  = LY.ItQ; % Ip(psi) = I_tor(psi)
-  T       = LY.TQ;
-  vol     = LY.VQ;
-  B2_fs   = mu0*Ip_psi./mVprime + T.^2.*Rm2_fs;
+  % get vacuum field data from ids
   R0 = ids_core_sources.vacuum_toroidal_field.r0;
   B0 = interp1(ids_core_sources.time,ids_core_sources.vacuum_toroidal_field.b0,ec_tgrid_toray);
+  
   jtoray_to_jpar0 = 2*pi./B0./(T.*Rm2_fs).*B2_fs;
 
   % interpolate on ec_data rho_pol grid
@@ -234,9 +230,9 @@ if ~isempty(ec_gdat.ec.data) % if EC data available, fill sources
   vol_mapped = nan(size(ec_data.cd_dens.x));
   for ii = 1:nt_ec_toray
     jtoray_to_jpar0_mapped(:,ii) = ...
-      interp1(rho_pol_norm_liu,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+      interp1(rho_pol_norm_liu.data,jtoray_to_jpar0(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
     vol_mapped(:,ii) = ...
-      interp1(rho_pol_norm_liu,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
+      interp1(rho_pol_norm_liu.data,vol(:,ii),ec_data.cd_dens.grids.rho_pol_norm(:,ii));
   end
 
   % load 1d_profiles from ec_data
@@ -373,11 +369,11 @@ if numel(active_nbi)>0
   % find times where NBI is on to define time grid with extra time slice just
   % before/after  NBI power and at start/end of shot
   itime_nbi = find((powers_gdat.(nbi_names{active_nbi(1)}).data>0));
-  if nbi_powers_tgrid(itime_nbi(end))>=ohm_t_grid(end)
-    i_time_end = iround(nbi_powers_tgrid,ohm_t_grid(end));
-    nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
+  if nbi_powers_tgrid(itime_nbi(end))>=ohm_tgrid(end)
+    i_time_end = iround(nbi_powers_tgrid,ohm_tgrid(end));
+    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:i_time_end)'];
   else
-    nbi_tgrid_out = [ohm_t_grid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_t_grid(end)];
+    nbi_tgrid_out = [ohm_tgrid(1),nbi_powers_tgrid(itime_nbi(1)-1:itime_nbi(end)+1)',ohm_tgrid(end)];
   end
   nt_nbi_out = numel(nbi_tgrid_out);
 
@@ -425,11 +421,11 @@ if ~isempty(powers_gdat.dnbi)
   % find times where DNBI is on to define time grid with extra time slice just
   % before/after  DNBI power and at start/end of shot
   itime_dnbi = find((powers_gdat.dnbi.data>0));
-  if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_t_grid(end)
-    i_time_end = iround(dnbi_powers_tgrid,ohm_t_grid(end));
-    dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
+  if dnbi_powers_tgrid(itime_dnbi(end))>=ohm_tgrid(end)
+    i_time_end = iround(dnbi_powers_tgrid,ohm_tgrid(end));
+    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:i_time_end)'];
   else
-    dnbi_tgrid_out = [ohm_t_grid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_t_grid(end)];
+    dnbi_tgrid_out = [ohm_tgrid(1),dnbi_powers_tgrid(itime_dnbi(1)-1:itime_dnbi(end)+1)',ohm_tgrid(end)];
   end
   nt_dnbi_out = numel(dnbi_tgrid_out);