From 21d76a808634bd8496746f9c2c02a4e03ed07c42 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <olivier.sauter@epfl.ch>
Date: Mon, 13 May 2019 07:23:34 +0000
Subject: [PATCH] add torbeam_prepare_inputs_and_run.m to prepare all input
files for various times and run torbeam, as well as read_ech_py to get angles
git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11872 d63d8f72-b253-0410-a779-e742ad2e26cf
---
crpptbx/AUG/read_ech_py.m | 160 +++++++++++++
crpptbx/AUG/torbeam_plot_output_singlerun.m | 158 +++++++++++++
crpptbx/AUG/torbeam_prepare_inputs_and_run.m | 227 +++++++++++++++++++
3 files changed, 545 insertions(+)
create mode 100644 crpptbx/AUG/read_ech_py.m
create mode 100644 crpptbx/AUG/torbeam_plot_output_singlerun.m
create mode 100644 crpptbx/AUG/torbeam_prepare_inputs_and_run.m
diff --git a/crpptbx/AUG/read_ech_py.m b/crpptbx/AUG/read_ech_py.m
new file mode 100644
index 00000000..5945b5db
--- /dev/null
+++ b/crpptbx/AUG/read_ech_py.m
@@ -0,0 +1,160 @@
+function [torbeam] = read_ech_py(shot,varargin);
+%
+% [torbeam] = read_ech_py(shot,varargin);
+%
+% loads angles and powers from ECS/ECN for given shot
+% transforms ECS/ECN into theta,phi of Torbeam
+%
+% varargin{1}: time array for Torbeam relevant data output
+%
+torbeam.rfmod = - ones(1,8); % 1 for O-mode, -1 for X-mode
+torbeam.tb_par.('xxb') = [238., 238., 231.1, 231.1, 236.4, 236.4, 236.4 ,236.4];
+torbeam.tb_par.('xzb') = [0., 0., 0., 0., 32., 32., -32., -32.];
+torbeam.tb_par.('xryyb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4]; % 85.0 in GT's py, 85.4 in tb_gui saved files
+torbeam.tb_par.('xrzzb') = [129.4,129.4,129.4, 129.4, 85.4, 85.4, 85.4, 85.4];
+torbeam.tb_par.('xwyyb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
+torbeam.tb_par.('xwzzb') = [2.97, 2.97, 2.97, 2.97, 2.55, 2.55, 2.55, 2.55];
+
+torbeam.tr2as = { ...
+ 'XECECH', 'RR_n'; 'ZECECH', 'ZZ_n'; ...
+ 'ECB_WIDTH_HOR', 'xwyyb'; 'ECB_WIDTH_VERT', 'xwzzb'; ...
+ 'ECB_CURV_HOR' , 'xryyb'; 'ECB_CURV_VERT' , 'xrzzb' };
+
+%aa=gdat(shot,'pgyro'); % loads ECS and ECN stuff
+
+tec_min = NaN;
+tec_max = NaN;
+
+% following python script:
+
+% WARNING seems alpha, beta inverted in read_ech.py script, so here use alpha for poloidal related angle and beta for tor
+% Thus use alpha_pol and beta_tor to avoid confusion (inverted with respect to .py alpha, beta script)
+
+ecsystems=[4,4];
+nbgy = sum(ecsystems);
+igy_seq=[1:ecsystems(1);ecsystems(1)+1:ecsystems(1)+ecsystems(2)]';
+for isys=1:length(ecsystems)
+ isyseff = isys;
+ if isys==1 && shot > 33725; isyseff = 3; end
+ for igy=1:ecsystems(isys)
+ ps_source = ['P_sy' num2str(isyseff) '_g' num2str(igy)];
+ aa=gdat(shot,{'ECS',ps_source,'','param:GPolPos'},'machine','aug');
+ torbeam.alpha_pol(igy_seq(igy,isys),1) = aa.data;
+ aa=gdat(shot,{'ECS',ps_source,'','param:GTorPos'},'machine','aug');
+ torbeam.beta_tor(igy_seq(igy,isys),1) = aa.data;
+ aa=gdat(shot,{'ECS',ps_source,'','param:gyr_freq'},'machine','aug');
+ if aa.data > 1000
+ torbeam.freq(igy_seq(igy,isys),1) = aa.data;
+ else
+ torbeam.freq(igy_seq(igy,isys),1) = 1.398765E+11;
+ end
+ if isyseff == 1
+ aaon=gdat(shot,{'ECS',ps_source,'','param:gy_on'},'machine','aug');
+ else
+ aaon=gdat(shot,{'ECS',ps_source,'','param:gy_HV_on'},'machine','aug');
+ end
+ if isempty(aaon.data)
+ warning(['problem with detecting if gy_on for isyseff = ' num2str(isyseff) ', igy = ',num2str(igy)]);
+ end
+ torbeam.is_gyro_on(igy_seq(igy,isys)) = false;
+ if aaon.data; torbeam.is_gyro_on(igy_seq(igy,isys)) = true; end
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ if isys==1
+ aa=gdat(shot,{'ECS',['PG' num2str(igy)]},'machine','aug');
+ else
+ aa=gdat(shot,{'ECS',['PG' num2str(igy) 'N']},'machine','aug');
+ end
+ torbeam.pow{igy_seq(igy,isys)}.data = aa.data;
+ torbeam.pow{igy_seq(igy,isys)}.t = aa.t;
+ tec_min = min(min(torbeam.pow{igy_seq(igy,isys)}.t),tec_min);
+ tec_max = max(max(torbeam.pow{igy_seq(igy,isys)}.t),tec_max);
+ else
+ torbeam.pow{igy_seq(igy,isys)}.data = [];
+ torbeam.pow{igy_seq(igy,isys)}.t = [];
+ end
+ end
+end
+
+% Compute theta,phi for TORBEAM
+torbeam.theta = NaN(nbgy,1);
+torbeam.phi = NaN(nbgy,1);
+torbeam.err_status = -ones(nbgy,1);
+torbeam.isystunit = NaN(nbgy,1);
+torbeam.datum = NaN(nbgy,1);
+for isys=1:length(ecsystems)
+ for igy=1:ecsystems(isys)
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ [torbeam.theta(igy_seq(igy,isys),1),torbeam.phi(igy_seq(igy,isys),1),torbeam.err_status(igy_seq(igy,isys),1),torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
+ aug_setval2tp_mex(torbeam.alpha_pol(igy_seq(igy,isys),1),torbeam.beta_tor(igy_seq(igy,isys),1),igy_1_8,shot);
+ end
+ end
+end
+
+% get time dependent data
+
+if shot < 23187
+ % Before ECRH2 installed
+ return
+end
+
+if ~isfinite(tec_min) || ~isfinite(tec_max)
+ warning('no good time interval')
+ tec_min
+ tec_max
+ return
+end
+
+tec_min = max(tec_min, 0.);
+if tec_max <= tec_min
+ warning('no good time interval')
+ tec_min
+ tec_max
+ return
+end
+
+isys=length(ecsystems);
+isyseff = isys;
+diagname = 'ECN';
+if nargin>=2 && ~isempty(varargin{1})
+ time_for_torbeam = varargin{1};
+else
+ time_for_torbeam = linspace(tec_min,tec_max,round(tec_max-tec_min)*3e2); % aim at 3ms
+end
+torbeam.t = time_for_torbeam;
+torbeam.theta_t = NaN(nbgy,length(time_for_torbeam));
+torbeam.phi_t = NaN(nbgy,length(time_for_torbeam));
+torbeam.err_status_t = NaN(nbgy,length(time_for_torbeam));
+% 1st system angles fixed so use above values
+isys = 1;
+for igy=1:ecsystems(isys)
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ torbeam.theta_t(igy_1_8,:) = torbeam.theta(igy_seq(igy,isys),1);
+ torbeam.phi_t(igy_1_8,:) = torbeam.phi(igy_seq(igy,isys),1);
+ torbeam.err_status_t(igy_1_8,:) = torbeam.err_status(igy_seq(igy,isys),1);
+ % linear interpolation for power to avoid over/under shoot
+ torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
+ end
+end
+isys = 2;
+for igy=1:ecsystems(isys)
+ torbeam.alpha_polN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'POL']},'machine','aug');
+ torbeam.alpha_polN{igy}.data = 0.01 .* torbeam.alpha_polN{igy}.data;
+ torbeam.alpha_pol_t = interpos(torbeam.alpha_polN{igy}.t,torbeam.alpha_polN{igy}.data,torbeam.t,-1e2);
+% $$$ torbeam.beta_torN{igy} = gdat(shot,{diagname,['G' num2str(igy) 'TOR']},'machine','aug');
+ igy_1_8 = igy_seq(igy,isys);
+ if torbeam.is_gyro_on(igy_seq(igy,isys))
+ [theta_t,phi_t,err_status_t,torbeam.isystunit(igy_seq(igy,isys),1),torbeam.datum(igy_seq(igy,isys),1)] = ...
+ aug_setval2tp_mex(torbeam.alpha_pol_t,torbeam.beta_tor(igy_seq(igy,isys),1).*ones(size(torbeam.alpha_pol_t)),igy_1_8,shot);
+ torbeam.theta_t(igy_1_8,:) = theta_t';
+ torbeam.phi_t(igy_1_8,:) = phi_t';
+ torbeam.err_status_t(igy_1_8,:) = err_status_t';
+ torbeam.power_t(igy_1_8,:) = interpos(-21,torbeam.pow{igy_1_8}.t,torbeam.pow{igy_1_8}.data,torbeam.t);
+ end
+end
+
+min_power_on = 10;
+torbeam.beam_on = false.*ones(size(torbeam.power_t));
+torbeam.beam_on(torbeam.power_t >= min_power_on) = true;
+torbeam.t_some_beam_on = torbeam.t(any(torbeam.beam_on,1));
diff --git a/crpptbx/AUG/torbeam_plot_output_singlerun.m b/crpptbx/AUG/torbeam_plot_output_singlerun.m
new file mode 100644
index 00000000..b7a7661e
--- /dev/null
+++ b/crpptbx/AUG/torbeam_plot_output_singlerun.m
@@ -0,0 +1,158 @@
+function [torbeam_std_files_struct] = torbeam_plot_output_singlerun;
+%
+% [torbeam_std_files_struct] = torbeam_plot_output_singlerun;
+%
+% reads local topfile, Te.dat, ne.dat input files
+% as well as outputs files 't1_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat','cntpr.dat'
+% and plot results
+%
+
+% read topfile
+fid=fopen('topfile','r');
+% fscanf(fid,'%s %d %.6f\n',' Number of radial and vertical grid points in ',shot,time_for_torbeam(it));
+dummy=fgetl(fid);
+dummy=fgetl(fid);
+dummyaaa = sscanf(dummy,'%d',2);
+zzz.nb_rmesh = dummyaaa(1);
+zzz.nb_zmesh = dummyaaa(2);
+% fscanf(fid,'%s\n','Inside and Outside radius and psi_sep');
+dummy=fgetl(fid);
+dummy=fgetl(fid);
+dummyaaa = sscanf(dummy,'%f',3);
+zzz.rmesh_min = dummyaaa(1);
+zzz.rmesh_max = dummyaaa(2);
+zzz.psi_edge = dummyaaa(3);
+% fscanf(fid,'%s\n','Radial grid coordinates');
+dummy=fgetl(fid);
+zzz.rmesh = fscanf(fid,'%f',zzz.nb_rmesh);
+% fscanf(fid,'%s\n','Vertical grid coordinates');
+dummy=fgetl(fid);
+if isempty(deblank(dummy)); dummy=fgetl(fid); end
+zzz.zmesh = fscanf(fid,'%f',zzz.nb_zmesh);
+% fscanf(fid,'%s\n','B_r values');
+dummy=fgetl(fid);
+if isempty(deblank(dummy)); dummy=fgetl(fid); end
+zzz.BR = fscanf(fid,'%f',[zzz.nb_rmesh,zzz.nb_zmesh]);
+% fscanf(fid,'%s\n','B_t values');
+dummy=fgetl(fid);
+if isempty(deblank(dummy)); dummy=fgetl(fid); end
+zzz.Bphi = fscanf(fid,'%f',[zzz.nb_rmesh,zzz.nb_zmesh]);
+% fscanf(fid,'%s\n','B_z values');
+dummy=fgetl(fid);
+if isempty(deblank(dummy)); dummy=fgetl(fid); end
+zzz.BZ= fscanf(fid,'%f',[zzz.nb_rmesh,zzz.nb_zmesh]);
+% fscanf(fid,'%s\n','psi values');
+dummy=fgetl(fid);
+if isempty(deblank(dummy)); dummy=fgetl(fid); end
+zzz.psi= fscanf(fid,'%f',[zzz.nb_rmesh,zzz.nb_zmesh]);
+fclose(fid);
+
+zzz.prof_to_cp = {'Te.dat','ne.dat'};
+zzz.prof_to_cp_var = lower(strrep(zzz.prof_to_cp,'.dat',''));
+for i=1:length(zzz.prof_to_cp_var)
+ zzz.(zzz.prof_to_cp_var{i}) = [];
+ try
+ [zzz.(zzz.prof_to_cp_var{i})(:,1),zzz.(zzz.prof_to_cp_var{i})(:,2)] = textread(zzz.prof_to_cp{i},'%f%f','headerlines',1);
+ catch
+ end
+end
+
+zzz.tb_out_files_to_cp = {'t1_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat'};
+zzz.tb_out_files_to_cp_var = strrep(zzz.tb_out_files_to_cp,'.dat','');
+for i=1:length(zzz.tb_out_files_to_cp_var)
+ zzz.(zzz.tb_out_files_to_cp_var{i}) = [];
+ try
+ zzz.(zzz.tb_out_files_to_cp_var{i}) = load(zzz.tb_out_files_to_cp{i});
+ catch
+ end
+end
+zzz.tbeg_abs=load('cntpr.dat');
+zzz.tbeg_abs(2) = zzz.tbeg_abs(2) + 1; % counted from 0?
+
+% plots
+
+
+figure
+contour(zzz.rmesh,zzz.zmesh,zzz.BR',100);
+hold on
+contour(zzz.rmesh,zzz.zmesh,zzz.psi',[zzz.psi_edge zzz.psi_edge],'linewidth',2,'linecolor','k');
+for i=1:2:size(zzz.t1_LIB,2)
+ hpl=plot(zzz.t1_LIB(:,i)/1e2,zzz.t1_LIB(:,i+1)/1e2);
+ plotos(zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i)/1e2,zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i+1)/1e2,'-', ...
+ [7 0],[],'k');%get(hpl,'color'));
+end
+axis equal
+colorbar
+title('BR')
+
+figure
+contour(zzz.rmesh,zzz.zmesh,zzz.BZ',100);
+hold on
+contour(zzz.rmesh,zzz.zmesh,zzz.psi',[zzz.psi_edge zzz.psi_edge],'linewidth',2,'linecolor','k');
+for i=1:2:size(zzz.t1_LIB,2)
+ hpl=plot(zzz.t1_LIB(:,i)/1e2,zzz.t1_LIB(:,i+1)/1e2);
+ plotos(zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i)/1e2,zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i+1)/1e2,'-', ...
+ [7 0],[],'k');%get(hpl,'color'));
+end
+axis equal
+colorbar
+title('BZ')
+
+figure
+contour(zzz.rmesh,zzz.zmesh,zzz.Bphi',100);
+hold on
+contour(zzz.rmesh,zzz.zmesh,zzz.psi',[zzz.psi_edge zzz.psi_edge],'linewidth',2,'linecolor','k');
+for i=1:2:size(zzz.t1_LIB,2)
+ hpl=plot(zzz.t1_LIB(:,i)/1e2,zzz.t1_LIB(:,i+1)/1e2);
+ ixx= plotos(zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i)/1e2,zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i+1)/1e2,'-', ...
+ [7 0],[],'k');%get(hpl,'color'));
+end
+axis equal
+colorbar
+title('Bphi')
+
+figure;
+subplot(2,1,1)
+plot(zzz.te(:,1),zzz.te(:,2))
+ylabel('Te')
+subplot(2,1,2)
+plot(zzz.ne(:,1),zzz.ne(:,2))
+ylabel('ne')
+
+figure;
+subplot(2,1,1)
+plot(zzz.t2_new_LIB(:,1),zzz.t2_new_LIB(:,2))
+ylabel('Pdens')
+subplot(2,1,2)
+plot(zzz.t2_new_LIB(:,1),zzz.t2_new_LIB(:,3))
+ylabel('jcd')
+
+figure
+theta=linspace(0,2.*pi,300);
+Rmintheta = (0.3+zzz.rmesh_min).*cos(theta);
+Ymintheta = (0.3+zzz.rmesh_min).*sin(theta);
+Rmaxtheta = (0.3+zzz.rmesh_max).*cos(theta);
+Ymaxtheta = (0.3+zzz.rmesh_max).*sin(theta);
+plot(Rmintheta,Ymintheta,'k')
+hold on
+plot(Rmaxtheta,Ymaxtheta,'k')
+for i=1:2:size(zzz.t1_LIB,2)
+ hpl=plot(zzz.t1tor_LIB(:,i)/1e2,zzz.t1tor_LIB(:,i+1)/1e2);
+ plotos(zzz.t1tor_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i)/1e2,zzz.t1tor_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i+1)/1e2,'-', ...
+ [7 0],[],'k');%get(hpl,'color'));
+end
+
+figure
+contour(zzz.rmesh,zzz.zmesh,zzz.psi',100);
+hold on
+contour(zzz.rmesh,zzz.zmesh,zzz.psi',[zzz.psi_edge zzz.psi_edge],'linewidth',2,'linecolor','k');
+for i=1:2:size(zzz.t1_LIB,2)
+ hpl=plot(zzz.t1_LIB(:,i)/1e2,zzz.t1_LIB(:,i+1)/1e2);
+ plotos(zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i)/1e2,zzz.t1_LIB(zzz.tbeg_abs(2):sum(zzz.tbeg_abs),i+1)/1e2,'-', ...
+ [7 0],[],'k');%get(hpl,'color'));
+end
+axis equal
+colorbar
+title('psi')
+
+torbeam_std_files_struct = zzz;
diff --git a/crpptbx/AUG/torbeam_prepare_inputs_and_run.m b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
new file mode 100644
index 00000000..9925badf
--- /dev/null
+++ b/crpptbx/AUG/torbeam_prepare_inputs_and_run.m
@@ -0,0 +1,227 @@
+function [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(shot,times_in,varargin);
+%
+% [torbeam_out,torbeam_in,varargout] = aaa(shot,times_in,varargin);
+%
+% varargin{1}: torbeam_in (from a previous run, so avoids to reload all the data)
+%
+
+if nargin==0 || isempty(shot)
+ error('requires shot number')
+ return
+end
+
+times_for_torbeam_eff = [];
+if nargin>=2 && ~isempty(times_in)
+ times_for_torbeam_eff = times_in;
+end
+
+doload = 1;
+if nargin>=3 && ~isempty(varargin{1})
+ torbeam_in = varargin{1};
+ doload = 0;
+end
+
+if doload
+ tic
+ torbeam_in.nete=gdat(shot,'nete_rho','machine','aug');
+ % qrho_and_grid=gdat(shot,'q_rho','machine','aug');
+else
+ if ~isfield(torbeam_in,'nete') || isempty(torbeam_in.nete)
+ error('expects nete as field of torbeam_in (varargin{1}) from gdat(shot,''nete_rho'')')
+ end
+end
+
+% base all times on Te data (most important? although ne as well), note take Thomson by default to avoid cut-off
+if ~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff),
+ times_for_torbeam_eff = torbeam_in.nete.fit.t;
+end
+% load beam releated data including theta,phi in torbeam reference
+if doload,
+ [torbeam_in.inbeam] = read_ech_py(shot,times_for_torbeam_eff);
+else
+ if ~isfield(torbeam_in,'inbeam') || isempty(torbeam_in.inbeam)
+ error('expects inbeam as field of torbeam_in (varargin{1}) from read_ech_py(shot,times_for_torbeam)')
+ end
+end
+
+% if torbeam_in provided as input (with its own time input sequence) and times_in provided as well: need to adapt times_in to match available inputs
+% (otherwise would take closest eqdsk)
+if doload==0
+ times_for_torbeam_eff = intersect(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
+else
+ times_for_torbeam_eff = torbeam_in.inbeam.t_some_beam_on;
+end
+if doload,
+ torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0);
+else
+ if ~isfield(torbeam_in,'eqd_all') || isempty(torbeam_in.eqd_all)
+ error('expects eqd_all as field of torbeam_in (varargin{1}) from gdat(shot,''eqdsk'',''time'',torbeam_in.inbeam.t_some_beam_on,''write'',0);')
+ end
+end
+
+if doload; disp(['time for reading data: ' num2str(toc)]); end
+
+if doload,
+ nb_points_prof = 46;
+ rhofit = linspace(0,1,nb_points_prof);
+
+ dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
+ unix(['mkdir ' dir_torbeam]);
+ tb_par_fields=fieldnames(torbeam_in.inbeam.tb_par);
+
+ itime_te = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
+ itime_ne = iround_os(torbeam_in.nete.fit.t,times_for_torbeam_eff);
+ it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
+ cocos_out = 17;
+ tic
+ for it=1:length(torbeam_in.inbeam.t_some_beam_on)
+ % Te
+ it_te = itime_te(it);
+ torbeam_in.fname{it}.te = fullfile(dir_torbeam,['Te' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
+ fid=fopen(torbeam_in.fname{it}.te,'w'); % te.dat file etc per time
+ fprintf(fid,'%d\n',nb_points_prof);
+ % tefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_te) , torbeam_in.nete.fit.te.data(:,it_te),rhofit,-1,[1 2],[0 0]);
+ tefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_te) , torbeam_in.nete.fit.te.data(:,it_te),rhofit,-0.1,[1 0],[0 0]);
+ aa=[rhofit(end:-1:1)' , abs(tefit(end:-1:1))/1e3]; % in keV
+ fprintf(fid,'%.7f %.7f\n',aa');
+ fclose(fid);
+ % ne
+ it_ne = itime_ne(it);
+ torbeam_in.fname{it}.ne = fullfile(dir_torbeam,['ne' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) '.dat']);
+ fid=fopen(torbeam_in.fname{it}.ne,'w'); % ne.dat file etc per time
+ fprintf(fid,'%d\n',nb_points_prof);
+ nefit = interpos(torbeam_in.nete.fit.rhopolnorm(:,it_ne) , torbeam_in.nete.fit.ne.data(:,it_ne),rhofit,-0.1,[1 0],[0 0]);
+ aa=[rhofit(end:-1:1)' , abs(nefit(end:-1:1))/1e19]; % in 1e19
+ fprintf(fid,'%.6f %.6f\n',aa');
+ fclose(fid);
+ % topfile
+ % eqd_with_b = read_eqdsk(torbeam_in.eqd_all.eqdsk{it},17,1);
+ eqd_with_b2=eqdsk_cocos_transform(torbeam_in.eqd_all.eqdsk{it},[17 cocos_out]);
+ eqd_with_b1=eqdsk_transform(eqd_with_b2,65,129,-0.1);
+ eqd_with_b = read_eqdsk(eqd_with_b1,cocos_out,1); % to compute BR, etc
+ torbeam_in.fname{it}.topfile = fullfile(dir_torbeam,['topfile' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it))]);
+ fid=fopen(torbeam_in.fname{it}.topfile,'w');
+ fprintf(fid,'%s %d %.6f\n',' Number of radial and vertical grid points in ',shot,torbeam_in.inbeam.t_some_beam_on(it));
+ fprintf(fid,'%d %d\n',length(eqd_with_b.rmesh),length(eqd_with_b.zmesh));
+ fprintf(fid,'%s\n','Inside and Outside radius and psi_sep');
+ psi_shift = 0.; % or eqd_with_b.psiedge for 0 value at edge
+ fprintf(fid,'%.6f %.6f %.6f\n',min(eqd_with_b.rmesh),max(eqd_with_b.rmesh),eqd_with_b.psiedge-psi_shift);
+ fprintf(fid,'%s\n','Radial grid coordinates');
+ nb_rmesh = length(eqd_with_b.rmesh);
+ nb_zmesh = length(eqd_with_b.zmesh);
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.rmesh(1:nb_rmesh));
+ if mod(nb_rmesh,8) ~= 0; fprintf(fid,'\n'); end
+ fprintf(fid,'%s\n','Vertical grid coordinates');
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.zmesh(1:nb_zmesh));
+ if mod(nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
+ fprintf(fid,'%s\n','B_r values');
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.BR(1:nb_rmesh,1:nb_zmesh));
+ if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
+ fprintf(fid,'%s\n','B_t values');
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.Bphi(1:nb_rmesh,1:nb_zmesh));
+ if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
+ fprintf(fid,'%s\n','B_z values');
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.BZ(1:nb_rmesh,1:nb_zmesh));
+ if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
+% $$$ fprintf(fid,'%s\n','psi norm values');
+% $$$ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',(eqd_with_b.psi(1:nb_rmesh,1:nb_zmesh)-eqd_with_b.psiaxis) ...
+% $$$ ./ (eqd_with_b.psiedge-eqd_with_b.psiaxis));
+ fprintf(fid,'%s\n',['psi cocos_out=' num2str(cocos_out)]);
+ fprintf(fid,'%14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f %14.6f\n',eqd_with_b.psi(1:nb_rmesh,1:nb_zmesh));
+ if mod(nb_rmesh*nb_zmesh,8) ~= 0; fprintf(fid,'\n'); end
+ fclose(fid);
+
+ % extract beam info at given time
+ for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it)))',
+ torbeam_in.fname{it}.inbeam{ibeam} = fullfile(dir_torbeam,['inbeam' '_t' num2str(torbeam_in.inbeam.t_some_beam_on(it)) 'g' num2str(ibeam) '.dat']);
+ [aa,bb]=unix(['cp ~/TORBEAM/inbeam_default.dat ' torbeam_in.fname{it}.inbeam{ibeam}]);
+ fid=fopen(torbeam_in.fname{it}.inbeam{ibeam},'a');
+ fprintf(fid,' xf = %f, ! wave frequency om=2*pi*xf\n',torbeam_in.inbeam.freq(ibeam));
+ fprintf(fid,' nmod = %d, ! mode selection: O-mode (1), X-mode (-1)\n',torbeam_in.inbeam.rfmod(ibeam));
+ fprintf(fid,' xtordeg = %f, ! geom. optics injection angle\n',torbeam_in.inbeam.phi_t(ibeam,it_inb(it)));
+ fprintf(fid,' xpoldeg = %f, ! geom. optics injection angle\n',torbeam_in.inbeam.theta_t(ibeam,it_inb(it)));
+ % may be position and curvature, but get defaults from tb_gui. Check with tb_gui how several launchers are given
+ for i=1:length(tb_par_fields)
+ fprintf(fid,[' ' tb_par_fields{i} ' = %f,\n'],torbeam_in.inbeam.tb_par.(tb_par_fields{i})(ibeam));
+ end
+ fprintf(fid,' xpw0 = %f,\n',torbeam_in.inbeam.power_t(ibeam,it_inb(it))/1.e6);
+ fprintf(fid,'/\n');
+ fclose(fid);
+ end
+ end
+ disp(['time for creating all input files: ' num2str(toc)])
+end
+
+% run TORBEAM
+current_dir = pwd;
+eval(['cd ' dir_torbeam])
+% namelist, ol proj 3 rays(R,z)*100, not used now, tor proj 3 rays(R,Y), vol() dP/dV(r/a) j(r/a)
+%tb_out_files_to_cp = {'echo.dat', 't1_LIB.dat', 't2_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat', 'cntpr.dat'};
+tb_out_files_to_cp = {'t1_LIB.dat', 't1tor_LIB.dat', 'volumes.dat', 't2_new_LIB.dat'};
+tb_out_files_to_cp_varname = strrep(tb_out_files_to_cp,'.dat','');
+torbeam_out.run_torbeam = false .* ones(size(torbeam_in.inbeam.power_t,1),length(times_for_torbeam_eff));
+it_inb_beam_on = iround_os(torbeam_in.inbeam.t_some_beam_on,times_for_torbeam_eff);
+torbeam_out.t = torbeam_in.inbeam.t_some_beam_on(it_inb_beam_on);
+tic
+for it_eff=1:length(times_for_torbeam_eff)
+ it = it_inb_beam_on(it_eff);
+ unix(['cp ' torbeam_in.fname{it}.te ' Te.dat']);
+ unix(['cp ' torbeam_in.fname{it}.ne ' ne.dat']);
+ unix(['cp ' torbeam_in.fname{it}.topfile ' topfile']);
+ disp([' t=' num2str(torbeam_in.inbeam.t(it_inb(it))) ' beam=' num2str(find(torbeam_in.inbeam.beam_on(:,it_inb(it)))','%d')]);
+ for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it)))',
+ unix(['cp ' torbeam_in.fname{it}.inbeam{ibeam} ' inbeam.dat']);
+ [a,b]=unix('/afs/ipp-garching.mpg.de/home/o/osauter/TORBEAM/lib-OUT/a.out');
+ imw = findstr('(MW)',b);
+ if ~isempty(imw)
+ torbeam_out.run_torbeam(ibeam,it) = true;
+ torbeam_out.pabs(ibeam,it) = 1e6 * sscanf(b(imw+4:end),'%f',1);
+ ima = findstr('(MA)',b);
+ torbeam_out.icd(ibeam,it) = 1e6 * sscanf(b(ima+4:end),'%f',1);
+ imaxat = findstr('max. at',b);
+ torbeam_out.pabs_roamax(ibeam,it) = sscanf(b(imaxat(1)+7:end),'%f',1);
+ torbeam_out.icd_roamax(ibeam,it) = sscanf(b(imaxat(2)+7:end),'%f',1);
+ for j=1:length(tb_out_files_to_cp)
+ torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it} = load(fullfile(dir_torbeam,[tb_out_files_to_cp{j}]));
+ end
+ nnn = load('cntpr.dat');
+ torbeam_out.cntpr.tbegabs(ibeam,it) = nnn(1);
+ torbeam_out.cntpr.dtabs(ibeam,it) = nnn(2);
+ torbeam_out.pdens_jcd_rho(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,1);
+ torbeam_out.pdens(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,2);
+ torbeam_out.jcd(:,ibeam,it) = torbeam_out.t2_new_LIB{ibeam,it}(:,3);
+ else
+ disp(['TORBEAM did not run properly, no MW result'])
+ end
+ end
+end
+disp(['time for running all TORBEAMs: ' num2str(toc)])
+for ibeam=find(any(torbeam_out.run_torbeam,2))'
+ figure;
+ subplot(2,1,1)
+ plot(squeeze(torbeam_out.pdens_jcd_rho(:,ibeam,:)),squeeze(torbeam_out.pdens(:,ibeam,:)))
+ title(['for beam = ' num2str(ibeam)])
+ ylabel('Pdens')
+ xlabel('r/a')
+ subplot(2,1,2)
+ plot(squeeze(torbeam_out.pdens_jcd_rho(:,ibeam,:)),squeeze(torbeam_out.jcd(:,ibeam,:)))
+ ylabel('jcd')
+ xlabel('r/a')
+end
+
+figure
+subplot(3,1,1)
+plot(torbeam_out.t,torbeam_out.pabs/1e6)
+ylabel('Pabs [MW]')
+
+subplot(3,1,2)
+plot(torbeam_out.t,torbeam_out.icd/1e6)
+ylabel('Icd [MA]')
+
+subplot(3,1,3)
+plot(torbeam_out.t,torbeam_out.pabs_roamax/1e6)
+ylabel('r/a (max Pdens)')
+xlabel('t [s]')
+
+
+eval(['cd ' current_dir])
--
GitLab