finished implementing and testing torbeam_prepare_inputs_and_run and related...

finished implementing and testing torbeam_prepare_inputs_and_run and related functions, will test reading TORBEAM angles from shotfile from Giovanni

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@11874 d63d8f72-b253-0410-a779-e742ad2e26cf

crpptbx/AUG/torbeam_prepare_inputs_and_run.m

@@ -3,6 +3,11 @@ function [torbeam_out,torbeam_in,varargout] = torbeam_prepare_inputs_and_run(sho
 % [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)
+%              leave then times_in empty to rerun at all times in torbeam_in.inbeam.t_some_beam_on or select times points from that array
+%
+% Examples:
+%       [torbeam_out,torbeam_in] = torbeam_prepare_inputs_and_run(30594,linpace(0,10,101));
+%       [torbeam_out2] = torbeam_prepare_inputs_and_run(30594,torbeam_in.inbeam.t_some_beam_on(1:30:end),torbeam_in); % to rerun a few cases
 %
 
 if nargin==0 || isempty(shot)
@@ -32,8 +37,8 @@ else
 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;
+if doload && (~exist('times_for_torbeam_eff') || isempty(times_for_torbeam_eff)),
+  times_for_torbeam_eff = torbeam_in.nete.fit.t(torbeam_in.nete.fit.t>0);
 end
 % load beam releated data including theta,phi in torbeam reference
 if doload, 
@@ -47,12 +52,16 @@ 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;
+  if ~isempty(times_for_torbeam_eff)
+    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
 end
 if doload,
   torbeam_in.eqd_all=gdat(shot,'eqdsk','time',torbeam_in.inbeam.t_some_beam_on,'write',0);
+  torbeam_in.dir_torbeam = ['/tmp/' getenv('USER') '/' num2str(shot) '_torbeam'];
+  unix(['mkdir ' torbeam_in.dir_torbeam]);
 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);')
@@ -65,8 +74,6 @@ 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);
@@ -77,7 +84,7 @@ if doload,
   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']);
+    torbeam_in.fname{it}.te = fullfile(torbeam_in.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]);
@@ -87,7 +94,7 @@ if doload,
     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']);
+    torbeam_in.fname{it}.ne = fullfile(torbeam_in.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]);
@@ -99,7 +106,7 @@ if doload,
     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))]);
+    torbeam_in.fname{it}.topfile = fullfile(torbeam_in.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));
@@ -133,7 +140,7 @@ if doload,
     
     % 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']);
+      torbeam_in.fname{it}.inbeam{ibeam} = fullfile(torbeam_in.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));
@@ -154,48 +161,52 @@ end
 
 % run TORBEAM
 current_dir = pwd;
-eval(['cd ' dir_torbeam])
+eval(['cd ' torbeam_in.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);
+it_inb = iround_os(torbeam_in.inbeam.t,times_for_torbeam_eff);
 torbeam_out.t = torbeam_in.inbeam.t_some_beam_on(it_inb_beam_on);
 tic
+icount = 0;
 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)))',
+  disp(['   t=' num2str(torbeam_in.inbeam.t(it_inb(it_eff))) '   beam=' num2str(find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))','%d')]);
+  for ibeam=find(torbeam_in.inbeam.beam_on(:,it_inb(it_eff)))',
     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');
+    icount = icount + 1;
     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);
+      torbeam_out.run_torbeam(ibeam,it_eff) = true;
+      torbeam_out.pabs(ibeam,it_eff) = 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);
+      torbeam_out.icd(ibeam,it_eff) = 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);
+      torbeam_out.pabs_roamax(ibeam,it_eff) = sscanf(b(imaxat(1)+7:end),'%f',1);
+      torbeam_out.icd_roamax(ibeam,it_eff) = 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}]));
+        torbeam_out.(tb_out_files_to_cp_varname{j}){ibeam,it_eff} = load(fullfile(torbeam_in.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);
+      torbeam_out.cntpr.tbegabs(ibeam,it_eff) = nnn(1);
+      torbeam_out.cntpr.dtabs(ibeam,it_eff) = nnn(2);
+      torbeam_out.pdens_jcd_rho(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,1);
+      torbeam_out.pdens(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,2);
+      torbeam_out.jcd(:,ibeam,it_eff) = torbeam_out.t2_new_LIB{ibeam,it_eff}(:,3);
     else
       disp(['TORBEAM did not run properly, no MW result'])
     end
   end
 end
-disp(['time for running all TORBEAMs: ' num2str(toc)])
+aaa=toc;
+disp(['time for running all TORBEAMs: ' num2str(aaa) ' for ' num2str(icount) ' runs so ' num2str(aaa/icount) 's per TORBEAM run'])
 for ibeam=find(any(torbeam_out.run_torbeam,2))'
   figure;
   subplot(2,1,1)
@@ -211,16 +222,16 @@ end
 
 figure
 subplot(3,1,1)
-plot(torbeam_out.t,torbeam_out.pabs/1e6)
+plot(torbeam_out.t,torbeam_out.pabs/1e6,'*-')
 ylabel('Pabs [MW]')
 
 subplot(3,1,2)
-plot(torbeam_out.t,torbeam_out.icd/1e6)
+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)')
+plot(torbeam_out.t,torbeam_out.pabs_roamax,'*-')
+ylabel('r(max Pdens) / a')
 xlabel('t [s]')