From e426f98243bd7974ae5c553b3cc8ae25b8dd1630 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Thu, 9 Feb 2023 17:38:30 +0100
Subject: [PATCH] start mv mds/efit files from d3d to private
---
matlab/D3D/efit/private/Contents.m | 50 +
matlab/D3D/efit/private/cc_efit_to_tok.m | 242 ++++
matlab/D3D/efit/private/efit_movie.m | 59 +
matlab/D3D/efit/private/eq_mds.m | 81 ++
matlab/D3D/efit/private/eq_mod.m | 113 ++
matlab/D3D/efit/private/equil_to_gdata.m | 197 +++
matlab/D3D/efit/private/get_efit_data1.m | 274 ++++
matlab/D3D/efit/private/gfile_def.m | 52 +
matlab/D3D/efit/private/inside_plasma.m | 292 ++++
.../D3D/efit/private/make_mhdindat_from_tds.m | 211 +++
matlab/D3D/efit/private/plasma_current.m | 200 +++
matlab/D3D/efit/private/read_mhdin.m | 392 ++++++
matlab/D3D/efit/private/run_efit.m | 93 ++
matlab/D3D/efit/private/shot_from_gfile.m | 125 ++
matlab/D3D/efit/private/std_efit_units.m | 1238 +++++++++++++++++
matlab/D3D/efit/private/trace_boundary.m | 229 +++
matlab/D3D/efit/private/write_afile.m | 634 +++++++++
matlab/D3D/efit/private/write_cc_file.m | 55 +
matlab/D3D/efit/private/write_eqdsk_efit.m | 65 +
matlab/D3D/efit/private/write_gfile.m | 320 +++++
matlab/D3D/efit/private/write_kfile.m | 190 +++
.../D3D/mdsplus_func/private/get_mds_tree.m | 291 ++++
matlab/D3D/mdsplus_func/private/getmds.m | 11 +
.../mdsplus_func/private/make_mdsipmex_east.m | 28 +
matlab/D3D/mdsplus_func/private/mds_eq.m | 298 ++++
.../D3D/mdsplus_func/private/mds_sub_tree.m | 70 +
26 files changed, 5810 insertions(+)
create mode 100644 matlab/D3D/efit/private/Contents.m
create mode 100644 matlab/D3D/efit/private/cc_efit_to_tok.m
create mode 100644 matlab/D3D/efit/private/efit_movie.m
create mode 100644 matlab/D3D/efit/private/eq_mds.m
create mode 100644 matlab/D3D/efit/private/eq_mod.m
create mode 100644 matlab/D3D/efit/private/equil_to_gdata.m
create mode 100644 matlab/D3D/efit/private/get_efit_data1.m
create mode 100644 matlab/D3D/efit/private/gfile_def.m
create mode 100644 matlab/D3D/efit/private/inside_plasma.m
create mode 100644 matlab/D3D/efit/private/make_mhdindat_from_tds.m
create mode 100644 matlab/D3D/efit/private/plasma_current.m
create mode 100644 matlab/D3D/efit/private/read_mhdin.m
create mode 100644 matlab/D3D/efit/private/run_efit.m
create mode 100644 matlab/D3D/efit/private/shot_from_gfile.m
create mode 100644 matlab/D3D/efit/private/std_efit_units.m
create mode 100644 matlab/D3D/efit/private/trace_boundary.m
create mode 100644 matlab/D3D/efit/private/write_afile.m
create mode 100644 matlab/D3D/efit/private/write_cc_file.m
create mode 100644 matlab/D3D/efit/private/write_eqdsk_efit.m
create mode 100644 matlab/D3D/efit/private/write_gfile.m
create mode 100644 matlab/D3D/efit/private/write_kfile.m
create mode 100644 matlab/D3D/mdsplus_func/private/get_mds_tree.m
create mode 100644 matlab/D3D/mdsplus_func/private/getmds.m
create mode 100644 matlab/D3D/mdsplus_func/private/make_mdsipmex_east.m
create mode 100644 matlab/D3D/mdsplus_func/private/mds_eq.m
create mode 100644 matlab/D3D/mdsplus_func/private/mds_sub_tree.m
diff --git a/matlab/D3D/efit/private/Contents.m b/matlab/D3D/efit/private/Contents.m
new file mode 100644
index 00000000..705ea442
--- /dev/null
+++ b/matlab/D3D/efit/private/Contents.m
@@ -0,0 +1,50 @@
+ %
+% Matlab EFIT Area (Plus MDS-EFIT Read)
+%
+% read_afile Reads a0-eqdsk file and puts variables in Matlab env
+% read_gfile Reads g0-eqdsk file and puts variables in Matlab env
+% read_gfile_func Same as read_gfile but function call & structure output
+% read_gfile_tok read gfile structure for particular machines d3d,kstar...
+% read_mfile_func read mfile in netcdf format into struct (ONLY for Matlab 7 V14)
+% plot_efit Plots efit flux contours
+% get_efit_data1 Gets data for lots of EFIT slices and plots it
+% Note: Conflict with daves get_efit_data in /plasma_models 7/03
+% efit_movie Make movie from EFIT flux data
+% run_efit Runs EFIT from inside Matlab much as EFIT runs from unix
+% get_efit_profile Get EFIT current profile parametrization from esave.dat
+% read_esave Read the esave.dat file (fortran mex function)
+% efit_lp Calculate plasma inductance from gfile information
+% plasma_field Calculate plasma flux and field at arbitrary rpt,zpt points
+% shot_from_gfile returns shot,time,dir (or tree,server) based on gfile name
+%
+% New Routines to extract EFIT "eq." structure from MDS: (work in progress)
+%
+% cc_efit_to_tok converts data fetched from mdsplus or gfile to toksys units
+% eq_mds Fetch an EFIT tree out of mdsplus
+% eq_mk_time_last Modifies eq. to put time index as last index in all vectors
+% gdata_to_eq Reads gdata from read_gfile_tok and puts into "eq." structure
+% read_mds_eqdsk Return data like output of read_gfile_tok.m from mdsplus
+% read_mds_eq_func Return data like output of read_gfile_func.m from mdsplus
+% equil_to_gdata Identical output as read_gfile_func but reads from mds+
+%
+% Below eq. routines work in progress:
+% eq_time_lim Reduces time data in sructure eq. within limits; also makes var
+% eq_plasma_cur Generates the plasma current from flux inside eq.
+% plasma_current Generates plasma current and current density from efit psi
+% inside_plasma Function Determines all points inside plasma
+% eq_ga_env Makes standard G (& later A file) variables from eq. into env.
+%
+% Misc EFIT scripts
+%
+% calc_fpcurr Emulates fpcurr in efit (ffprime calculation from polynom.)
+% calc_ppcurr Emulates ppcurr in efit (pprime calculation from polynom.)
+% calc_bsffel Emulates bsffel in efit (iparm'th term of poly representation)
+% calc_bsppel Emulates bsppel in efit (iparm'th term of poly rep. of pprime)
+% gfile_def Generates structure of EFIT gfile variable definitions
+% test_get_efit_profile - Example for get_efit_profile.m,current reconstruction
+%
+% Emmanuel Joffrin's EFIT tools:
+%
+% fluxfun - function to extract profile data from EFIT and make flux mappings
+% fluxmoy - called by fluxfun...
+% isoflux - called by fluxfun...
diff --git a/matlab/D3D/efit/private/cc_efit_to_tok.m b/matlab/D3D/efit/private/cc_efit_to_tok.m
new file mode 100644
index 00000000..66b35a9d
--- /dev/null
+++ b/matlab/D3D/efit/private/cc_efit_to_tok.m
@@ -0,0 +1,242 @@
+function equil_I = cc_efit_to_tok(vac_objs,equil_data)
+ %
+% SYNTAX: equil_I = cc_efit_to_tok(vac_objs,equil_data)
+%
+% PURPOSE: Convert equilibrium conductor currents in internal EFIT
+% representation (contained in equil_data) into a representation
+% compatible with toksys environment.
+%
+% INPUT:
+% vac_objs = vacuum model data object structure
+% equil_data = equilibrium data structure
+% idx_efit_to_tok = optional map of efit pf coil indices to toksys indices, s.t. if
+% I=currents in efit order, then
+%` I(idx_efit_to_tok)=currents in toksys order
+% idxvv_efit_to_tok = " vv indices " " "
+%
+% OUTPUT: (units defined by imks, iterminal in vac_objs)
+% equil_I = structure containing:
+% cc0t = coil current equilibrium vector, toksys conductors
+% cc0e = coil current equilibrium vector, EFIT conductors
+% vc0t = vessel current equilibrium vector, toksys conductors
+% vc0e = vessel current equilibrium vector, EFIT conductors
+% cprojet = mapping: cc0t = cprojet*cc0e
+% cprojte = mapping: cc0e = cprojte*cc0t
+% vprojte = mapping: vc0e = vprojte*vc0t
+% (Equilibrium vessel currents are usually all zero. In those cases, vc0e
+% and vprojte are not produced.)
+
+% RESTRICTIONS:
+%
+% METHOD:
+%
+% WRITTEN BY: Mike Walker ON 8/31/09 (from in-line code in rzrig)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%Derived values:
+ if vac_objs.imks, iscale=1e6; else iscale=1; end
+ iterminal = vac_objs.iterminal;
+ ncc = size(vac_objs.mcc,1);
+
+% Extract Needed EFIT data from equil_data:
+
+ if(~isfield(vac_objs,'ecnturn') | isempty(vac_objs.ecnturn))
+ ecnturn = 0;
+ else
+ ecnturn = vac_objs.ecnturn;
+ end
+ if(ecnturn~=0)
+ ccnturn = [ecnturn; vac_objs.fcnturn];
+ if(~iterminal)
+ ccnturn = ones(size(ccnturn));
+ end
+ else
+ ccnturn = vac_objs.fcnturn;
+ if(~iterminal)
+ ccnturn = ones(size(ccnturn));
+ end
+ end
+% make all "turn" objects into columns, consistent with coil current vectors
+ [s1,s2]=size(ccnturn); if(s1==1), ccnturn = ccnturn'; end
+
+% Define desired currents for calculating equilibrium field.
+% Assumes cc came out of read_gfile always in MA-turns.
+
+if(~isfield(equil_data,'cc') | isempty(equil_data.cc))
+
+ wait('ERROR cc_efit_to_tok: no valid cc data - returning empty matrices')
+ equil_I = struct( ...
+ 'cc0t', [], ...
+ 'cc0e', [], ...
+ 'cprojte', []);
+
+else
+ ccx = iscale*equil_data.cc(:); % ccx now in proper units (MA or Amps)
+
+% Complicated logic required here because either turnfc or fcturn can hold
+% the fraction of current values - depends on how used internally in EFIT.
+
+ if(any(equil_data.fcturn<1))
+ Ifrac = equil_data.fcturn;
+ nturn = equil_data.turnfc;
+ elseif(any(equil_data.turnfc<1))
+ Ifrac = equil_data.turnfc;
+ nturn = equil_data.fcturn;
+ elseif(any(equil_data.turnfc>1))
+ Ifrac = equil_data.fcturn;
+ nturn = equil_data.turnfc;
+ elseif(any(equil_data.fcturn>1))
+ Ifrac = equil_data.turnfc;
+ nturn = equil_data.fcturn;
+ else % all must be = 1
+ Ifrac = equil_data.turnfc;
+ nturn = equil_data.fcturn;
+ end
+% if any ecoils, insert before fcoils:
+
+ idxf = abs(equil_data.fcid); % use abs in case anti-series connection
+ isign = sign(equil_data.fcid);
+ if(~isempty(equil_data.ecid))
+ nee = length(unique(equil_data.ecid));
+ Ifrac = [ones(1,nee) Ifrac];
+ for i=nee:-1:1
+ idx = find(equil_data.ecid==i);
+ neturn = sum(equil_data.ecturn(idx));
+ nturn = [neturn nturn];
+ end
+ idx = [[1:nee] idxf+nee];
+ isign = [ones(1,nee) isign];
+ else
+ idx = idxf;
+ end
+
+ [mm,i1] = min(size(Ifrac)); [mm,i2] = min(size(ccx));
+ if(i1~=i2) Ifrac = Ifrac'; end;
+ cc0 = isign'.*ccx(idx).*Ifrac;
+
+ nn = length(cc0);
+ if(isfield(equil_data,'idx_efit_to_tok'))
+ idx_efit_to_tok = equil_data.idx_efit_to_tok;
+ else
+ idx_efit_to_tok = 1:length(nturn);
+ end
+
+ if(iterminal)
+ [mm,i1] = min(size(nturn)); [mm,i2] = min(size(ccnturn));
+ if(i1~=i2) nturn = nturn'; end;
+% (length(ccnturn) > length(nturn) is allowed. Extra currents are 0.)
+ if(norm(ccnturn(1:length(nturn)) - nturn(idx_efit_to_tok)))
+ fprintf('WARNING cc_efit_to_tok:\n');
+ fprintf(' #turns equil. data=');
+ fprintf('%d ',nturn); fprintf('\n');
+ fprintf(' #turns vacuum objects=');
+ fprintf('%d ',ccnturn(1:length(nturn))); fprintf('\n');
+% wait
+ end
+ cc0 = cc0./nturn;
+ end
+
+% The following accounts for coils assumed to carry no steady state current
+% NOTE that these are always assumed to be listed last.
+ if(length(cc0)<ncc)
+ cc0 = [cc0; zeros(ncc-length(cc0),1)];
+ end
+
+% At this point, cc0 is in correct units, lumped/terminal mode, and represents
+% equilibrium currents assigned to all coils represented by vacuum objects.
+
+cprojet = zeros(length(cc0),length(equil_data.cc));
+cprojte = zeros(length(equil_data.cc),length(cc0));
+if(isfield(vac_objs,'ecnturn'))
+ nec = length(vac_objs.ecnturn);
+ for k=1:nec
+ cprojet(k,k)=1; % ohmic coil
+ cprojte(k,k)=1; % ohmic coil
+ end
+else
+ nec = 0;
+end
+for k=nec+1:size(cprojte,1)
+ idx = find(equil_data.fcid==k-nec);
+ if(~isempty(idx))
+ cprojet(idx+nec,k)=1;
+ cprojte(k,idx(1)+nec)=1;
+ end
+end
+
+if(isfield(equil_data,'idx_efit_to_tok'))
+ nn = length(cc0);
+ ni = max(idx_efit_to_tok);
+ temp = zeros(nn);
+ for k=1:ni
+ temp(idx_efit_to_tok(k),k)=1;
+ end
+ cc0t = zeros(size(cc0));
+ cc0t(1:ni) = cc0(idx_efit_to_tok);
+ cprojte = cprojte*temp;
+ cprojet = temp'*cprojet;
+else
+ cc0t = cc0;
+end
+
+cc0e = cprojte*cc0t;
+
+equil_I = struct( ...
+'cc0t', cc0t, ...
+'cc0e', cc0e, ...
+'cprojet', cprojet, ...
+'cprojte', cprojte);
+
+description = struct( ...
+'cc0t', 'coil currents corresponding to toksys vacuum objects', ...
+'cc0e', 'coil currents corresponding to EFIT fitted coils', ...
+'cprojet', 'mapping: cc0t = cprojet*cc0e', ...
+'cprojte', 'mapping: cc0e = cprojte*cc0t');
+
+end
+
+% If vessel currents in equilibrium, these must be used in rzrig calculation.
+if(isfield(equil_data,'vc') & ~isempty(equil_data.vc))
+ vcx = iscale*equil_data.vc; % vcx now in proper units (MA or Amps)
+ vprojte= proj_turn(equil_data.vvid,equil_data.vvfrac);
+ vc0= vprojte'*vcx;
+
+
+% construct a projection without inverting
+ vprojte= proj_turn(equil_data.vvid,1./equil_data.vvfrac);
+ ii = find(vprojte~=0);
+ vprojtej= zeros(size(vprojte));
+ vprojtej(ii)= 1;
+ num= sum(vprojtej')';
+ vprojtej= (1./num)*ones(1,size(vprojte,2));
+ vprojte= vprojte.*vprojtej;
+
+ % Translate conductors if toksys & efit have different indices or #
+ if(isfield(equil_data,'idxvv_efit_to_tok'))
+ nn = length(vc0);
+ ni = max(equil_data.idxvv_efit_to_tok);
+ temp = zeros(nn,ni);
+ for k=1:ni
+ temp(equil_data.idxvv_efit_to_tok(k),k)=1;
+ end
+ vc0t = zeros(vac_objs.nv,1);
+ vc0t(1:ni) = vc0(equil_data.idxvv_efit_to_tok);
+ vprojte = vprojte*temp;
+ else
+ vc0t = vc0;
+ end
+
+ vc0e = vprojte*vc0t;
+ equil_I.vc0t = vc0t;
+ equil_I.vc0e = vc0e;
+ equil_I.vprojte = vprojte;
+
+ description.vc0t='vessel currents corresponding to toksys vacuum objects from equil_data';
+ description.cc0e='vessel currents corresponding to toksys.def_connect fitted conductors';
+ description.vprojte= 'mapping: vc0e = vprojte*vc0t';
+else
+ equil_I.vc0t = zeros(vac_objs.nv,1);
+ description.vc0t='vessel currents corresponding to toksys vacuum objects';
+end
+
+equil_I.description = description;
diff --git a/matlab/D3D/efit/private/efit_movie.m b/matlab/D3D/efit/private/efit_movie.m
new file mode 100644
index 00000000..fa710f3b
--- /dev/null
+++ b/matlab/D3D/efit/private/efit_movie.m
@@ -0,0 +1,59 @@
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% USAGE: >>efit_movie
+%
+% PURPOSE: Script to display movie of EFIT flux patterns from set of EFITs.
+%
+% INPUTS:
+% eqdir = string giving directory containing EFIT g-files
+% (final character is '/', e.g. '/u/humphrys/efit/')
+% shot = integer shot number
+% t1ms = initial EFIT slice time (msec)
+% t2ms = final EFIT slice time (msec)
+% dtms = EFIT slices time interval (msec)
+%
+% OUTPUTS:
+% EFIT movie displayed.
+%
+% RESTRICTIONS:
+% EFIT g0-files identified by shot, t1ms, t2ms, dtms must exist in eqdir.
+% Input variables must be defined; no defaults provided EXCEPT for
+% eqdir, defaulted to './' (so uses default directory).
+%
+% METHOD:
+
+% WRITTEN BY: Dave Humphreys ON 1/1/97
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Default inputs:
+ if ~exist('eqdir')
+ eqdir = './';
+ end
+
+% Prelims:
+ mu0 = 0.4*pi;
+ clf,hold off %clear current figure
+
+% Derived Values:
+ npts = (t2ms-t1ms)/dtms + 1;
+
+% Acquire data and make movie:
+ for ii=1:npts
+ disp(['ii=',int2str(ii),' of ',int2str(npts)])
+ itime = t1ms + (ii-1)*dtms;
+ filename = [eqdir,'g0',int2str(shot),'.0',int2str(itime)];
+ hold off
+ plot_efit
+ if ii==1, M=moviein(npts); end %initialize M first time
+ M(:,ii) = getframe;
+ end
+
+% Play movie:
+ disp('Playing movie once...')
+ movie(M)
+ disp('Movie complete. To play again N times, do >>movie(M,N)')
+
+
+
diff --git a/matlab/D3D/efit/private/eq_mds.m b/matlab/D3D/efit/private/eq_mds.m
new file mode 100644
index 00000000..b587d8bb
--- /dev/null
+++ b/matlab/D3D/efit/private/eq_mds.m
@@ -0,0 +1,81 @@
+ function eq= eq_mds(shot, tree, server, toupper, verbose)
+%
+% SYNTAX:
+% eq= eq_mds(shot, tree, server, toupper, verbose); % full call
+% eq= eq_mds(shot); % defaults to DIII-D server
+% eq= eq_mds(shot, 'EFIT01', 'NSTX'); % for NSTX
+% eq= eq_mds(shot, 'NB', 'DIII-D'); % Other MDS+ Tree Data -beam
+% Should also Run:
+% eq= eq_mk_time_last(eq); % puts time at end (i.e. (130,1) => (1,130)
+%
+% PURPOSE: Get EFIT equilibrium (GEQDSK AEQDSK MEQDSK) from mdsplus database.
+%
+% INPUT: <default>
+% shot = shot number
+% tree = which efit tree to use <'EFIT01'>
+% server = which MDS+ database to use: defaults to <'DIII-D'>
+% 'DIII-D' opens 'atlas.gat.com'
+% 'NSTX' opens 'europa.pppl.gov:8501'
+% 'OPEN' assumes mdsconnect already called
+% server arbitrary server is called using mdsconnect(server)
+% toupper= 1= all variables made upper case, =-1 all var. made lower case
+% [0]= no change, variables made depending on mds case (typical UC)
+% verbose = set to 1 to get diagnostic prints during execution
+%
+% OUTPUT:
+% eq = structure containing EFIT eqdsk variables for all times
+% structure follows MDS tree structure with 'TOP' replaced with 'eq'
+% See: eq.allnames
+%
+% EXAMPLE Variables (Examples below are for DIII-D -- NSTX is different)
+%
+% Ex: eq.RESULTS.AEQDSK.IPMEAS, eq.RESULTS.GEQDSK.PSIRZ
+% Use: plot(eq.RESULTS.GEQDSK.GTIME,eq.RESULTS.GEQDSK.CPASMA)
+% contour(eq.RESULTS.GEQDSK.R(:,1), eq.RESULTS.GEQDSK.Z(:,1), ...
+% eq.RESULTS.GEQDSK.PSIRZ(:,:,round(length(eq.RESULTS.GEQDSK.GTIME)/2))')
+%
+% Some Extra items added to structure (all lower case)
+% eq.id = sting array of important data identifyer enf
+% eq.shot = shot number
+% eq.server= MDS+ server
+% eq.allnames= list of all variables in structure with full struc. address
+% eq.mdsnames= list of all variables in structure with mds struc. address
+%
+% NOTE: 1)Function returns all TIME data for shot in large arrays.
+% Example:
+% eq.RESULTS.AEQDSK.BETAN(255,1)
+% eq.RESULTS.AEQDSK.ATIME(255,1); % Time data vector for AEQDSK data
+% plot(eq.RESULTS.AEQDSK.ATIME,eq.RESULTS.AEQDSK.BETAN)
+% eq.RESULTS.GEQDSK.GTIME(130,1)
+% eq.RESULTS.GEQDSK.PPRIME(65,130)
+% eq.RESULTS.GEQDSK.PSIRZ(65,65,130)
+% 2) Time is in vector eq...GTIME and is in ms. All other units are as
+% they come from EFIT. psi(Vs/r) I(A) R(m) ...
+%
+% SEE ALSO: eq= eq_mk_time_last(eq); (should be run on eq)
+% eq_time_lim eq_ga_env str_to_ws
+%
+% NOTE: Written very generally and should return any subtree of MDS+
+% Example: ions= eq_mds(shot,'IONS');
+
+% WRITTEN BY: Jim Leuer ON 3/1/05
+% taken from mds_eq.m uses sub-structure to store
+% USES: eq_mod
+% To see MDS structure on HYDRA run traverser
+% tested on DIII-D and NSTX data and should work for JET data but not tested
+% ==========================================================================
+% @(#)eq_mds.m 1.7 06/23/10
+
+wait('WARNING eq_mds: eq_mds is OBSOLETE - use get_mds_tree instead')
+
+if(nargin<2)
+ wait('ERROR eq_mds: requires at least 2 inputs')
+elseif(nargin==2)
+ eq = get_mds_tree(shot, tree);
+elseif(nargin==3)
+ eq = get_mds_tree(shot, tree, server);
+elseif(nargin==4)
+ eq = get_mds_tree(shot, tree, server, toupper);
+elseif(nargin==5)
+ eq = get_mds_tree(shot, tree, server, toupper, verbose);
+end
diff --git a/matlab/D3D/efit/private/eq_mod.m b/matlab/D3D/efit/private/eq_mod.m
new file mode 100644
index 00000000..2644e28c
--- /dev/null
+++ b/matlab/D3D/efit/private/eq_mod.m
@@ -0,0 +1,113 @@
+ function eqmod= eq_mod(eq)
+%
+% SYNTAX:
+% eqmod= eq_mod(eq);
+%
+% PURPOSE: Puts time index as last index of each array.
+%
+% INPUT: <default>
+% eq = structure from mds_eq
+%
+% OUTPUT:
+% eqmod = structure containing EFIT eqdsk variables with time last entry
+% Example: eq.BCENTR(255,1) => eqmod.BCENTER(1,255)
+%
+% see also mds_eq
+%
+% WRITTEN BY: Jim Leuer ON 6/11/03
+% ==========================================================================
+%
+% Defaults
+ if nargin==0
+ disp('% eq_mod needs at least a "eq" argument')
+ help eq_mod
+ return
+ end
+
+% -----------------------------------------------
+ eqmod= eq;
+
+ if isfield(eq,'GTIME') & isfield(eq,'gnames')
+
+ nt= length(eq.GTIME);
+
+% loop over all variables
+ for ii= 1:length(eq.gnames);
+% ii=0; ii=ii+1;
+ nam= deblank(eq.gnames(ii,:));
+
+ str= ['ndim= ndims(eq.' nam ');'];
+ eval(str);
+ if ndim<=2
+ str= ['siz= size(eq.' nam ');'];
+ eval(str);
+ if siz(1)==nt
+ str= ['eqmod.' nam '= eq.' nam ''' ;'];
+ eval(str);
+ end
+ end
+ end
+
+ end
+
+% AFILE
+
+if isfield(eq,'ATIME') & isfield(eq,'anames')
+
+ nt= length(eq.ATIME);
+
+% loop over all variables
+ for ii= 1:length(eq.anames);
+% ii=0; ii=ii+1;
+ nam= deblank(eq.anames(ii,:));
+
+ str= ['ndim= ndims(eq.' nam ');'];
+ eval(str);
+ if ndim<=2
+ str= ['siz= size(eq.' nam ');'];
+ eval(str);
+ if siz(1)==nt
+ str= ['eqmod.' nam '= eq.' nam ''' ;'];
+ eval(str);
+ end
+ end
+ end
+
+ end
+
+% MFILE:
+
+
+if (isfield(eq,'TIME') | isfield(eq,'MTIME')) & isfield(eq,'mnames')
+
+ if isfield(eq,'MTIME')
+ nt= length(eq.MTIME);
+ else
+ nt= length(eq.TIME);
+ end
+
+% loop over all variables
+ for ii= 1:length(eq.mnames);
+% ii=0; ii=ii+1;
+ nam= deblank(eq.mnames(ii,:));
+
+ str= ['ndim= ndims(eq.' nam ');'];
+ eval(str);
+ if ndim<=2
+ str= ['siz= size(eq.' nam ');'];
+ eval(str);
+ if siz(1)==nt
+ str= ['eqmod.' nam '= eq.' nam ''' ;'];
+ eval(str);
+ end
+ end
+ end
+
+ end
+
+ return
+
+% Testing
+% eq= mds_eq(114504, 'EFIT01', 0);
+% eqmod= eq_mod(eq)
+
diff --git a/matlab/D3D/efit/private/equil_to_gdata.m b/matlab/D3D/efit/private/equil_to_gdata.m
new file mode 100644
index 00000000..50978b71
--- /dev/null
+++ b/matlab/D3D/efit/private/equil_to_gdata.m
@@ -0,0 +1,197 @@
+ function [gdata,ireadok]= equil_to_gdata(equil,tms)
+ %
+% USAGE: [gdata,ireadok] = equil_to_gdata(equil, tms);
+% USAGE: [gdata,ireadok] = equil_to_gdata(mds_string); % calls read_mds_func
+%
+%
+% PURPOSE: Convert equil data generated by read_mds_eq_func to read_gfile_func
+%
+% INPUTS: <default>
+% equil= equilibrium data: equil=read_mds_eq_func(shot,tree,server);
+%
+% mds_string= optionally if equil is a string containing MDS+ tree information:
+% format: <Server><.><tree>.shot.time
+% Examples:
+% 'D3D.EFIT09.g131498.02600'
+% 'D3D.EFITRT1.131498.02600'
+% '.g131498.02600' => 'D3D.EFIT01.g131498.02600'
+% tms= time [ms] to use from equil data
+%
+% OUTPUTS:
+% gdata = data structure containing EFIT G-file variables
+% ireadok = flag to report good read of gfile (0=bad, 1=good)
+%
+% RESTRICTIONS:
+% MDS EFIT Tree must exist for call with mds_string string convention
+%
+% NOTE TESTED WITH ANY MDS+ except D3D
+%
+
+% METHOD: Based on read_gfile_func & read_mds_eq_func
+
+% WRITTEN BY: Jim Leuer ON 8Mar2010
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% ================================================================
+% See if input is a 'DOT' format convention and then get data from MDS+
+% ================================================================
+
+ if isstr(equil) | nargin <=1
+ mds_string= equil;
+ id= strfind(mds_string,'.'); % Format: SERVER.TREE.SHOTNUM.TIME
+ if length(id) >=2 % must have 2 dots: .SHOTNUM.TIME
+ [shot,tms,tree,server] = shot_from_gfile(mds_string); % Dot Format
+ disp([' % NOTE equil_to_gdata is reading MDS+ for: ',mds_string]);
+
+ jphi= []; % make sure this is variable and not function
+ [equil,neq,eq] = read_mds_eq_func1(shot,tree,server);
+
+ else % if length(id)<=2
+
+ disp([' %ERROR equil_to_gdata: needs "server.tree.shot.time" ']);
+ gdata=[];
+ ireadok= 0;
+ return
+ end % if ~isempty(id)
+ end % if isstr(equil)
+
+% =========================================================================
+% Get only ONE time
+ if ~isempty(equil)
+ if ~isfield(equil,'tms')
+ if strcmp(upper(server),'NSTX')
+ equil.time= eq.results.geqdsk.gtime;
+ equil.tms= eq.results.geqdsk.gtime*1e+3;
+ elseif strcmp(upper(server),'KSTAR') | findstr(upper(tree),'KSTR') |...
+ strcmp(upper(server),'EAST') | findstr(upper(tree),'EAST')
+ equil.time= eq.results.geqdsk.gtime;
+ equil.tms= eq.results.geqdsk.gtime*1e+3;
+ else % DIIID
+ equil.time= eq.results.geqdsk.gtime*1e-3;
+ equil.tms= eq.results.geqdsk.gtime;
+ end
+ end
+
+ if tms >= equil.tms(1) & tms <= equil.tms(end)
+ idx= interp1(equil.tms,1:length(equil.tms),tms,'nearest');
+ else
+ disp([' %ERROR equil_to_gdata, tms: ' num2str(tms) ...
+ ' not between equil.tms ', ...
+ num2str(equil.tms(1)) ' ' num2str(equil.tms(end))]);
+ ireadok= 0;
+ gdata= [];
+ return
+ end
+% =========================================================================
+
+ gdef = gfile_def; % definitions of all variables read in from G-file
+
+ gdata = struct( ...
+ 'gdef', gdef, ...
+ 'brsp', equil.gdata(idx).brsp, ...
+ 'bzero', equil.gdata(idx).bzero, ...
+ 'cpasma', equil.gdata(idx).cpasma, ...
+ 'ecase', equil.gdata(idx).ecase, ...
+ 'ecurrt', equil.gdata(idx).ecurrt, ...
+ 'ffprim', equil.gdata(idx).ffprim, ...
+ 'fpol', equil.gdata(idx).fpol, ...
+ 'limitr', equil.gdata(idx).limitr, ...
+ 'nbbbs', equil.gdata(idx).nbbbs, ...
+ 'nh', equil.gdata(idx).nh, ...
+ 'nw', equil.gdata(idx).nw, ...
+ 'pcurrt', equil.gdata(idx).pcurrt, ...
+ 'pprime', equil.gdata(idx).pprime, ...
+ 'pres', equil.gdata(idx).pres, ...
+ 'psirz', equil.gdata(idx).psirz, ...
+ 'qpsi', equil.gdata(idx).qpsi, ...
+ 'rbbbs', equil.gdata(idx).rbbbs, ...
+ 'rgrid1', equil.gdata(idx).rgrid1, ...
+ 'rmaxis', equil.gdata(idx).rmaxis, ...
+ 'rzero', equil.gdata(idx).rzero, ...
+ 'ssibry', equil.gdata(idx).ssibry, ...
+ 'ssimag', equil.gdata(idx).ssimag, ...
+ 'xdim', equil.gdata(idx).xdim, ...
+ 'xlim', equil.gdata(idx).xlim, ...
+ 'ylim', equil.gdata(idx).ylim, ...
+ 'zbbbs', equil.gdata(idx).zbbbs, ...
+ 'zdim', equil.gdata(idx).zdim, ...
+ 'zmaxis', equil.gdata(idx).zmaxis, ...
+ 'zmid', equil.gdata(idx).zmid);
+
+ gdata.time= equil.time(idx);
+ gdata.tms= equil.tms(idx);
+ gdata.gdef.time= ['time in seconds; ', ...
+ 'Note: gtime,atime = ms in D3D,EAST,KSTAR but sec in NSTX (TRUTH)'];
+ gdata.gdef.tms= ['time in ms; '];
+ if exist('shot')
+ gdata.shot= shot;
+ gdata.gdef.shot= 'Shot number; used for mds call';
+ end
+ if exist('tms')
+ gdata.tms_target= tms;
+ gdata.gdef.tms_target= ...
+ 'Target time in ms; (.tms and .time are actual mds times found)';
+ end
+ if exist('tree')
+ gdata.tree= tree;
+ gdata.gdef.tree= 'tree; used for mds call';
+ end
+ if exist('server')
+ gdata.server= server;
+ gdata.gdef.server= 'server; used for mds call';
+ end
+ if exist('mds_string')
+ gdata.mds_string= mds_string;
+ gdata.gdef.mds_string= '1st input argument to equil_to_gdata';
+ end
+
+ ireadok= 1;
+ else
+ disp([' %ERROR equil_to_gdata, no gdata generated ']);
+ gdata= [];
+ ireadok= 0;
+ end % ~isempty(equil)
+
+ return
+% ==========================================================
+% testing
+ gfile= '/u/leuer/efit/d3d/shot131498/g131498.02600'
+ g0= read_gfile_func(gfile);
+
+ mds_name='d3d.efit01.g131498.02600';
+ [g1,ireadok] = equil_to_gdata(mds_name)
+ [g2,ireadok] = read_gfile_func(mds_name)
+
+ nms= char(fieldnames(g0));
+ s= char(' '*ones(size(nms,1),1));
+ [int2str((1:size(nms,1))') s nms]
+
+ id= [2:12,14:15 17,19:26,28:30];
+ for ii=id
+ n= deblank(nms(ii,:))
+ s= ['[g0.' n ' g1.' n ']'];
+ eval(s)
+ pause
+ end
+
+ contour(g0.pcurrt - g1.pcurrt)
+
+ for ii=1:g0.nw
+ [g0.pcurrt(:,ii) g1.pcurrt(:,ii)]
+ pause
+ end
+
+ for ii=1:g0.nw
+ [g0.psirz(:,ii) g1.psirz(:,ii)]
+ pause
+ end
+
+ [g0.rbbbs(1:g0.nbbbs) g1.rbbbs(1:g1.nbbbs)]
+
+% test read_gfile_func
+ [g2,iok]= read_gfile_func(filename);
+ [g3,iok]= read_gfile_func(gfilename);
+
+
+
diff --git a/matlab/D3D/efit/private/get_efit_data1.m b/matlab/D3D/efit/private/get_efit_data1.m
new file mode 100644
index 00000000..f9201aa0
--- /dev/null
+++ b/matlab/D3D/efit/private/get_efit_data1.m
@@ -0,0 +1,274 @@
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SYNTAX: get_efit_data
+%
+% PURPOSE: Script to extract many slices of EFIT data from set of a0-eqdsk
+% files assumed present in default directory and plot selected data.
+%
+% INPUT:
+% eqdir = directory for EFIT data files (default='./');
+% eqdir string must end in "/" (eg '/users/humphrys/matlab/')
+% shot = shot number
+% iplot = 1 if want to plot many frames on one plot (default)
+% = 2 if want to plot Z,R,Ip(t) on different plots
+% = 0 to skip plotting
+% EFIT time slices are specified using EITHER of the following methods:
+% tmiefit = 1st efit time slice (msec)
+% dt = delta time (msec)
+% nslices = # of efit slices
+% OR:
+% tefit = vector of efit times (msec) to read from a0 files.
+% If tefit is specified, input variables tmiefit, dt, nslices are ignored.
+% ipltref = present figure if want to start plotting after present fig
+% (default=0) e.g. if ipltref=2, EFIT plots start at fig 3
+% Currently used only for iplot=1
+%
+% OUTPUT:
+% Extracts data from a0 files and plots things.
+
+% RESTRICTIONS:
+% A0-files specified by inputs must exist in default directory.
+% In present version, efit times must all have 4 digits (ie t must be
+% > 999, and < 9999)...
+%
+% METHOD:
+% Uses read_afile to read A0-eqdsk file for series of EFIT slices.
+%
+% WRITTEN BY: Dave Humphreys ON 6/25/95 Adapted from compare_efitz.m
+%
+% MODIFICATION HISTORY:
+% Finally completed to work properly. DAH 2/96
+% Modified to accept inputs from Matlab environment instead of definitions
+% at beginning of script. MLW 5/15/97
+% Plotting logic modified. DAH 5/16/97
+%
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Initial definitions:
+fprintf('shot = %d\n',shot)
+if ~exist('iplot'), iplot = 1; end
+fprintf('iplot = %d\n',iplot)
+if(exist('tefit'))
+ fprintf('time slices specified by tefit \n')
+else
+ fprintf('tmiefit = %d\n',tmiefit)
+ fprintf('dt = %d\n',dt)
+ fprintf('nslices = %d\n',nslices)
+ tefit = linspace(tmiefit,tmiefit+dt*(nslices-1),nslices)'; %time vec in msec
+end
+
+% Default for EFIT data directory:
+ if ~exist('eqdir')
+ eqdir = './'; %default to default directory
+ end
+
+%Derived Quantities:
+ nt = length(tefit);
+ shotno=int2str(shot);
+
+% Data vectors:
+ t=tefit*1e-3; %time vector in secs
+ t1=t(1); t2=max(t);
+ zmagt=t; rmagt=t; zcurt=t; rcurt=t; routt=t; ipt=t; kappat=t;
+ aoutt=t; zsep1t=t; rsep1t=t; zsep2t=t; rsep2t=t; zoutt=t;
+ gaptopt=t; gapint=t; gapoutt=t;
+ betapt=t; lit=t; qet=t;
+ cmpr2t = [];
+
+for ii=1:nt
+ if shot<10000
+ stmp='00'
+ elseif shot<100000
+ stmp='0'
+ else
+ stmp='';
+ end
+ if tefit(ii) < 1e4
+ padstr = '.0';
+ if tefit(ii) < 1e3
+ padstr = '.00';
+ if tefit(ii) < 1e2
+ padstr = '.000';
+ if tefit(ii) < 1e1
+ padstr = '.0000';
+ end
+ end
+ end
+ end
+ filename = [eqdir,'a',stmp,shotno,padstr,int2str(tefit(ii))];
+ read_afile;
+ zmagt(ii) = zmagx*0.01; %cm->m
+ rmagt(ii) = rmagx*0.01;
+ zcurt(ii) = zcurrt*0.01;
+ rcurt(ii) = rcurrt*0.01;
+ zoutt(ii) = zout*0.01;
+ routt(ii) = rout*0.01;
+ aoutt(ii) = aout*0.01;
+ kappat(ii) = eout;
+ zsep1t(ii) = zseps1*0.01;
+ rsep1t(ii) = rseps1*0.01;
+ zsep2t(ii) = zseps2*0.01;
+ rsep2t(ii) = rseps2*0.01;
+ gaptopt(ii) = otop*0.01;
+ gapoutt(ii) = oright*0.01;
+ gapint(ii) = oleft*0.01;
+ ipt(ii) = pasmat*1e-6; %A->MA
+ betapt(ii) = betap;
+ lit(ii) = ali;
+ qet(ii) = qpsib;
+ cmpr2t = [cmpr2t; cmpr2'];
+ cpasmat(ii) = cpasma;
+end
+
+%Derived and others:
+ zlim = -1.3675; %floor limiter vertical position [m]
+ zrelt = (zcurt-zlim); %vert posit rel to floor limiter
+ dzdipt = (zrelt(2:nt)-zrelt(1:nt-1))./(ipt(2:nt)-ipt(1:nt-1));
+
+%Plot things:
+if iplot==1
+if ~exist('ipltref'), ipltref=0; end
+figure(1+ipltref),clf,hold off
+nplots=4;
+subplot(nplots,1,1)
+ plot(t,rcurt)
+ hold on
+ plot(t,rmagt,'m--')
+ plot(t,routt,'c-.')
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Rcur,mag,out [m]')
+title([shotno,': EFIT Ovrvu 1 '])
+subplot(nplots,1,2)
+ plot(t,zcurt)
+ hold on
+ plot(t,zmagt,'m--')
+ plot(t,zoutt,'c-.')
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Zcur,mag,out [m]')
+subplot(nplots,1,3)
+ plot(t,aoutt)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('a [m]')
+subplot(nplots,1,4)
+ plot(t,kappat)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Kappa')
+xlabel('t [sec]')
+
+figure(2+ipltref),clf,hold off
+nplots=4;
+subplot(nplots,1,1)
+ plot(t,rsep1t)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Rx1 [m]')
+title([shotno,': EFIT Ovrvu 2'])
+subplot(nplots,1,2)
+ plot(t,zsep1t)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Zx1 [m]')
+subplot(nplots,1,3)
+ plot(t,rsep2t)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Rx2 [m]')
+subplot(nplots,1,4)
+ plot(t,zsep2t)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Zx2 [m]')
+xlabel('t [sec]')
+
+figure(3+ipltref),clf,hold off
+nplots=4;
+subplot(nplots,1,1)
+ plot(t,zsep1t)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Zx1 [m]')
+title([shotno,': EFIT Ovrvu 3'])
+subplot(nplots,1,2)
+ plot(t,gapoutt)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Gout [m]')
+subplot(nplots,1,3)
+ plot(t,gaptopt)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Gtop [m]')
+subplot(nplots,1,4)
+ plot(t,gapint)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Gin [m]')
+xlabel('t [sec]')
+
+figure(4+ipltref),clf,hold off
+nplots=4;
+subplot(nplots,1,1)
+ plot(t,ipt)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Ip [MA]')
+title([shotno,': EFIT Ovrvu 4'])
+subplot(nplots,1,2)
+ plot(t,lit)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('li')
+subplot(nplots,1,3)
+ plot(t,betapt)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('Betap')
+subplot(nplots,1,4)
+ plot(t,qet)
+ axis0=axis; axis([t1 t2 axis0(3:4)])
+ grid
+ ylabel('q_{edge}')
+xlabel('t [sec]')
+
+end %end iplot==1
+
+
+if iplot==2
+ figure(1),clf,hold off
+ plot(t,zrelt)
+ hold on
+ plot(t,zmagt-zlim,'m--') %vert posit of mag axis rel to floor
+ xlabel('t [sec]')
+ ylabel('Zrel [m]')
+ title([shotno,': Z Relative to Floor'])
+ grid
+
+figure(2),clf,hold off
+ plot(t,rcurt)
+ hold on
+ plot(t,rmagt,'m--')
+ xlabel('t [sec]')
+ ylabel('Rcur [m]')
+ title([shotno,': Rcur '])
+ grid
+
+figure(3),clf,hold off
+ plot(t,ipt)
+ xlabel('t [sec]')
+ ylabel('Ip_meas [MA]')
+ title([shotno,': Plasma Current'])
+ grid
+
+%figure(3),clf,hold off
+% plot(t(1:nt-1),dzdipt)
+% hold on
+% plot(t(1:nt-1),smoothdt(dzdipt,t,4e-3),'m--')
+% xlabel('t [sec]')
+% ylabel('dz/dIp [m/MA]')
+% title([shotno,': dZ/dIp During VDE'])
+% grid
+
+end %end iplot==2
+
diff --git a/matlab/D3D/efit/private/gfile_def.m b/matlab/D3D/efit/private/gfile_def.m
new file mode 100644
index 00000000..9809de13
--- /dev/null
+++ b/matlab/D3D/efit/private/gfile_def.m
@@ -0,0 +1,52 @@
+ function gdef = gfile_def()
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SYNTAX: gdef = gfile_def
+%
+% PURPOSE: Produces G-file variable definition structure for inclusion in
+% gfile structure
+%
+% INPUT: (none)
+%
+% OUTPUT:
+% gdef = data structure containing EFIT G-file variable definitions
+
+% RESTRICTIONS:
+%
+% METHOD:
+%
+% WRITTEN BY: Jim Leuer ON 01/05/2005
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% From efit_definitions.xls (sheet gfile)
+
+ gdef = struct( ...
+ 'brsp', 'OPTIONAL F-coil current vector; BRSP(1:NFCOIL); UNITS: depend on fcturn in mhdin.dat: fcturn()=1 brsp=[A-turns]; fcturn()=#turns brsp= [A]; brsp.*fcturn() is always [A-t], brsp.*fcturn().*turnfc() is always [A] ', ...
+ 'bzero', 'Vacuum toroidal magnetic field corresponding to RZERO, [T]', ...
+ 'cpasma', 'Plamsma Current [A]', ...
+ 'ecase', 'Identification character string; (1:48 characters of 1st line)', ...
+ 'ecurrt', 'OPTIONAL E-coil current vector; ECURRT(1:NESUM); [A]', ...
+ 'ffprim', 'Poloidal current function derivitive; FFPRIM= f*df/dPsi; ffprim= -sign(cpasma)*ffprime(in file); See FPOL; Jt (Amp/m2) = R*PPRIME + FFPRIM/(R*muo); FFPRIM(1:NW); [(Tm)^2/(Vs/Rad)]', ...
+ 'fpol', 'Poloidal current function; F= R*Bt on flux grid; Jt (Amp/m2) = R dP/dPsi + FdF/dPsi/(mu0*R); fpol(1:NW); [Tm]', ...
+ 'limitr', 'Number of limiter points; See: xlim,ylim; []', ...
+ 'nbbbs', 'Number of plasma boundary points; []', ...
+ 'nh', 'Number of vertical Z grid points ', ...
+ 'nw', 'Number of horizontal R grid points ', ...
+ 'pcurrt', 'OPTIONAL Plasma current on each grid point; PCURRT(1:NH;1:NW); [A]', ...
+ 'pprime', 'Plasma presure function derivitive; PPRIME= dp/dPsi; ppprime= -sign(cpasma)*pprime(in file); See PRES; Jt (Amp/m2) = R*PPRIME + FFPRIM/(R*muo); PPRIME(1:NW); [A/m^3 => N/(m^2)/(Vs/rad)]', ...
+ 'pres', 'Plasma pressure function; PRES= pressure on flux grid, Jt (Amp/m2) = R dP/dPsi + FdF/dpsi/(R*mu0); PRES(1:NW); [N/m^2]', ...
+ 'psirz', 'Poloidal flux on plasma grid in radial direction 1st order; PSIRZ(1:NW,1:NH); [Vs/rad]', ...
+ 'qpsi', 'Safety factor, q function; q(1)= axis; qpsi(1:nw); []', ...
+ 'rbbbs', 'Radial locaion of boundary; RBBBS(1:NBBBS); [m]', ...
+ 'rgrid1', 'Starting Radius of Grid; =rgrid(1); [m]', ...
+ 'rmaxis', 'R of magnetic axis; [m]', ...
+ 'rzero', 'Radius Used in Current Profile Parameterization for Ro; Radial location where BZERO is measured; [m] ', ...
+ 'ssibry', 'Poloidal flux at plasma boundary; [Vs/rad]', ...
+ 'ssimag', 'Poloidal flux at magnetic axis; [Vs/rad]', ...
+ 'xdim', 'Width of Grid; =rgrid(nw)-rgrid(1); [m]', ...
+ 'xlim', 'Radial locaion of limiter points; XLIM(1:LIMITR); [m]', ...
+ 'ylim', 'Z locaion of limiter points; YLIM(1:LIMITR); [m]', ...
+ 'zbbbs', 'Z location of boundary; ZBBBS(1:NBBBS); [m]', ...
+ 'zdim', 'Height of Grid; =zgrid(nh)-zgrid(1); [m]', ...
+ 'zmaxis', 'Z of magnetic axis; [m]', ...
+ 'zmid', 'Z of center of computational box; =(zgrid(1)+zgrid(nh))/2.0; [m] ' ...
+ );
+ return
diff --git a/matlab/D3D/efit/private/inside_plasma.m b/matlab/D3D/efit/private/inside_plasma.m
new file mode 100644
index 00000000..36c20e89
--- /dev/null
+++ b/matlab/D3D/efit/private/inside_plasma.m
@@ -0,0 +1,292 @@
+ function inside= inside_plasma(psizr, rgefit,zgefit, rbbbs,zbbbs,nbbbs,...
+ psimag,psibnd, iplot,dofast)
+%
+% inside_plasma sets matrix "inside" to 0's & 1's for outside & inside plasma
+%
+% SYNTAX:
+% inside= inside_plasma(psizr, rgefit,zgefit, rbbbs,zbbbs,...
+% psimag,psibnd, iplot,dofast);
+% inside= inside_plasma(psirz, rgefit,zgefit, rbbbs,zbbbs); % minimum input
+%
+% PURPOSE: set matrix "inside" with 0's or 1's depending if inside or outside
+% plasma. Uses same input as typically generated by read_gfile
+%
+% INPUTS: <default> [optional]
+%
+% psizr= Flux on plasma grid in (nz,nr) array format (Wb)
+% rgefit= Radius vector of plasma grid (m)
+% zgefit= Z vector of plasma grid (m)
+% rbbbs= Radius vector of closed polygon plasma boundary (m)
+% zbbbs= Z vector of closed polygon plasma boundary (m)
+% nbbbs= number of elements to use in rbbbs,zbbbs
+% [psimag]= Flux on axis (Wb) <max(max(psizr))>
+% [psibnd]= Flux on boundary (Wb) <average of psi on boundary>
+% [iplot]= 1; plot results <0>, 2= outside plasma only
+% [dofast]= 1; % increases algorithm speed at risk of including some
+% private flux in solution. Only extreme shapes like
+% crescents and beans are expected to need dofast=0.
+% Default is dofast=0 since speed is still very fast. <0>
+%
+% OUTPUTS:
+%
+% inside= Matrix of 1's & 0's for inside or outside plasma (nz,nr)
+%
+% RESTRICTIONS:
+% rb,zb adjacent elements cant be identical.
+% Not tested for coarse plasma boundary polygon mesh.
+%
+% METHOD:
+% uses limits to do most screening and finally FE triangular formula
+
+% WRITTEN BY: Jim Leuer ON 6-3-03
+%
+% MODIFICATION HISTORY:
+%
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Input Checks:
+
+ if nargin==9
+ dofast= 0;
+ elseif nargin==8
+ iplot=1;
+ dofast= 0;
+ elseif nargin==7
+ psibnd= mean(interp2(rgefit,zgefit,psizr,rbbbs,zbbbs));
+ iplot=0;
+ dofast= 0
+ elseif nargin==6
+ psimag= max(max(psizr));
+ psibnd= mean(interp2(rgefit,zgefit,psizr,rbbbs,zbbbs));
+ iplot=0;
+ dofast= 0;
+ elseif nargin==5
+ nbbb= length(rbbb);
+ elseif nargin<=4
+ disp('%ERROR inside_plasma: Needs atleast 4 arguments to work')
+ return
+ end
+
+% ========================================================================
+% INITIAL SETUP PARAMTERS (boundary, grid, plasma)
+
+ inside= zeros(size(psizr)); % 0's mean plasma 1's means no plasma
+
+% ---------------------------------------------------------------------
+% boundary polygon manipuliation to get in CCW and assending angle order
+
+ rbb= rbbbs(1:nbbbs); zbb= zbbbs(1:nbbbs);
+ if(size(rbb,1)==1), rbb=rbb'; end; % make sure its a column vector
+ if(size(zbb,1)==1), zbb=zbb'; end; % make sure its a column vector
+ rb0= mean(rbb); zb0= mean(zbb);
+ if rbb(1)==rbb(end) & zbb(1)==zbb(end) % remove duplicate last point
+ rbb= rbb(1:end-1); zbb= zbb(1:end-1);
+ end
+
+ drb= rbb-rb0;
+ dzb= zbb-zb0;
+% assending angular order from -pi makes this Counter Clockwise CCW
+ angb= atan2(dzb,drb);
+ [angb,id]= sort(angb);
+% make sure -pi to +pi is covered in assending order
+ rbb= rbb([id(end); id; id(1)]); % add points at front and back for 2pi
+ zbb= zbb([id(end); id; id(1)]); % add points at front and back for 2pi
+ angb= angb([id(end); id; id(1)]); % add points at front and back for 2pi
+ angb(1)= angb(1)-2*pi; % make 1 and end correct
+ angb(end)= angb(end)+2*pi; %
+ drb= drb([id(end); id; id(1)]);
+ dzb= dzb([id(end); id; id(1)]);
+ radb= sqrt(drb.^2+dzb.^2);
+ radbmn= min(radb);
+
+ drbmx= max(drb);
+ drbmn= min(drb);
+ dzbmx= max(dzb);
+ dzbmn= min(dzb);
+
+% plot(drb,dzb,'.b'), hold on, plot(drb([1,end]),dzb([1,end]),'rx'), grid on
+
+% Geometrical setup for particular grid
+
+ r2= ones(length(zgefit),1)*rgefit'; % same size as psizr
+ z2= zgefit*ones(1,length(rgefit));
+ dr2= r2-rb0;
+ dz2= z2-zb0;
+
+% plasma boundary psi limits clipped at grid limits
+
+ psibmx= max(max(psizr));
+ psibmn= min(min(psizr));
+
+ psibmx= min([psibmx, max(psibnd,psimag)]); % insure clip at grid values
+ psibmn= max([psibmn, min(psibnd,psimag)]); % insure clip at grid values
+
+% ======================================
+% Main Selection Part
+% ======================================
+
+% ======================================
+% 1st sort on center to boundary flux
+% figure(3)
+% clf
+% plot(rbbbs,zbbbs,'k')
+% hold on
+% axis equal
+
+ id1= find(psizr >= psibmn & psizr <= psibmx); % find points within min:max psi
+ inside(id1)= 1;
+
+% idplt= find(inside);
+% plot(r2(idplt),z2(idplt),'r.')
+
+
+% ======================================
+% 2nd sort on min/max Delta dimensions
+
+ ii= ones(size(id1));
+ id2= find( dr2(id1)<drbmn | dr2(id1)>drbmx | dz2(id1)<dzbmn | dz2(id1)>dzbmx);
+ ii(id2)= 0;
+
+ id2no= id1(find(~ii)); % find 0's which are outside plasma
+ inside(id2no)= 0; % set points outside plasma to zero
+ id22= id1(find(ii)); % this is remaining set that have ones
+
+% idplt= find(inside);
+% plot(r2(idplt),z2(idplt),'b*')
+
+% ===========================================================================
+% For most plasma we really dont have to do the below detailed calculations
+ if ~dofast
+
+% ======================================
+% 3rd eliminate points close to center based on smallest radius
+
+ rad= sqrt(dr2(id22).^2 + dz2(id22).^2);
+
+ ii= ones(size(id22));
+ id3= find(rad <= radbmn);
+ ii(id3)= 0; % this are GOOD point but dont have to do angle; inside already 1
+ id23= find(ii);
+ id33= id22(id23);
+% plot(r2(id33),z2(id33),'go') % now we must work on these points with angles
+
+% reduce set to work on:
+
+ rad= rad(id23); % reduced set to work on
+ ang= atan2(dz2(id33),dr2(id33));
+
+% ======================================
+% 4th elimination based on angle
+
+% pause
+
+% loop over boundary lines
+ for ibm=1:length(drb)-1 %
+% ibm=5; ibm= ibm+1
+
+ ibp= ibm+1; % boundary line end point id
+% plot(rbb([ibm,ibp]),zbb([ibm,ibp]),'co')
+% plot(rbb([ibm,ibp]),zbb([ibm,ibp]),'c')
+% plot([rbb(ibm),rb0],[zbb(ibm),zb0],'c')
+% plot([rbb(ibp),rb0],[zbb(ibp),zb0],'c')
+% plot(rb0,zb0,'ro')
+
+ id= find(ang>= angb(ibm) & ang<= angb(ibp)); % all points to process
+% plot(r2(id33(id)),z2(id33(id)),'b.')
+
+ if ~isempty(id)
+ radbm= min(radb(ibm:ibp));
+ id3= find(rad(id) > radbm); % Eliminates most points inside radius
+
+% Final detailed calculation based on local proximity to boundary
+% new method uses FE rotation dircetion of triangle see segerlind, pg29
+% assumes boundary is traversed in CCW direction
+% i=>j is boundary line, k is plasma point (use delta coordinates)
+ if ~isempty(id3) % work on all points which outside minimum bound rad
+ xi= drb(ibm);
+ yi= dzb(ibm);
+ xj= drb(ibp);
+ yj= dzb(ibp);
+
+ for ii=1:length(id3) % loop over remaining points to see if inside
+% ii=0, ii=ii+1;
+ id4= id33(id(id3(ii))); % pointer to main arrays
+% plot(r2(id4),z2(id4),'ro')
+ xk= dr2(id4);
+ yk= dz2(id4);
+ ccw= (xj-xi)*(yk-yj)-(xk-xj)*(yj-yi); % inside= + ccw around triang
+
+% plot([xi,xj,xk,xi]+rb0,[yi,yj,yk,yi]+zb0,'r')
+% arrow_head([xi;xj]+rb0,[yi;yj]+zb0,.1,'r');
+% arrow_head([xj;xk]+rb0,[yj;yk]+zb0,.1,'r');
+
+ if ccw < 0
+ inside(id4)=0; % this is where final points eliminated
+% plot([rb0,xx+rb0],[zb0,yy+zb0],'r')
+ end
+
+% pause(3)
+% disp('paused')
+% pause
+ end % for ii
+
+ end % ~isempty(id3)
+
+% need to reduce rad and ang by points processed
+ jj= ones(size(rad)); % all indexes in rad & ang;
+ jj(id)= 0; % eliminate all point processed
+ id23= find(jj);
+ id33= id33(id23); % reduce main pointer count
+ rad= rad(id23); % reduced set to work on
+ ang= ang(id23);
+% pause
+
+ end % ~isempty(id)
+
+ end % for ibm
+ end % dofast
+
+% ======================================
+% plot results
+
+ if iplot
+% plot(rgefit([1;end;end;1;1]),zgefit([1;1;end;end;1]),'r')
+ plot(rbbbs(1:nbbbs),zbbbs(1:nbbbs),'k')
+ hold on
+ axis equal
+% jal ERROR WORK AROUND FOR MAC/X11/ssh/THOR/MATLAB13 4/6/2005
+ if size(r2,1)*size(z2,2) >= 2145
+ inside1= inside(1:2:size(inside,1),1:2:size(inside,2));
+ r2= r2(1:2:size(r2,1),1:2:size(r2,2));
+ z2= z2(1:2:size(z2,1),1:2:size(z2,2));
+ else
+ inside1= inside;
+ end
+ idd= find(~inside1);
+ plot(r2(idd),z2(idd),'b.')
+
+ if iplot==1
+ idd= find(inside1);
+ plot(r2(idd),z2(idd),'r.')
+ end
+ end
+
+% ======================================
+
+ return
+
+
+% test: see also inside_plasma_old.m for different different algorithm
+
+% filename= '/users/walker/NSTX/g107167.00223'; % looks like BAD current data
+% read_gfile_nstx
+% nargin=8;
+% inside= inside_plasma(psizr,rgefit,zgefit,rbbbs,zbbbs,nbbbs,psimag,psibnd,1,0);
+% clear
+
+% filename= '/u/leuer/efit/diiid/g110214.01740';
+% read_gfile
+% inside= inside_plasma(psizr,rgefit,zgefit,rbbbs,zbbbs,nbbbs,psimag,psibnd,1,0);
+
+
+
diff --git a/matlab/D3D/efit/private/make_mhdindat_from_tds.m b/matlab/D3D/efit/private/make_mhdindat_from_tds.m
new file mode 100644
index 00000000..737f069e
--- /dev/null
+++ b/matlab/D3D/efit/private/make_mhdindat_from_tds.m
@@ -0,0 +1,211 @@
+function make_mhdindat_from_tds(tok_data_struct,options)
+% SYNTAX: make_mhdindat_from_tds(tok_data_struct,options)
+%
+% PURPOSE: Creates EFUND input file mhdin.dat from Toksys tok_data_struct
+%
+% INPUT:
+% tok_data_struct
+% options = structure
+% fcid = Connections between fcoils (see cccirc in build_tokamak_system) [1:nc]
+% use_vv = 1:Add VV elements as coils 0:Don't [0]
+% vvid = VV element conenctions (see vvid in build_tokamamk_system) [1:nv]
+%
+% OUTPUT: mhdin.dat file
+%
+% RESTRICTIONS:
+
+% WRITTEN BY: Nick Eidietis ON 2016-05-18
+% based upon Ander Welanders make_file_mhdin_dat4east
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ if exist('options','var')
+ struct_to_ws(options);
+ end
+ struct_to_ws(tok_data_struct);
+ if ~exist('use_vv','var')
+ use_vv = 0;
+ end
+
+
+ bpnam=' ';
+ for j=1:size(bpnames,1)
+ bpnam = [bpnam '''' deblank(bpnames(j,:)) ''','];
+ if j/5==round(j/5), bpnam=[bpnam 10 32]; end
+ end
+ xmp2=' ';
+ for j=1:size(bpnames,1)
+ s = num2str(bpdata(2,j)); while length(s)<5, s=[s '0']; end
+ xmp2 = [xmp2 s ', '];
+ if j/6==round(j/6), xmp2=[xmp2 10 32 32]; end
+ end
+ ymp2=' ';
+ for j=1:size(bpnames,1)
+ s = num2str(bpdata(1,j));
+ if bpdata(1,j)==0, s='0.';end
+ if bpdata(1,j)>=0, s=[' ' s];end
+ while length(s)<6, s=[s '0']; end
+ ymp2 = [ymp2 s ','];
+ if j/6==round(j/6), ymp2=[ymp2 10 32]; end
+ end
+ amp2='';
+ for j=1:size(bpnames,1)
+ s = [num2str(bpdata(3,j)) ','];
+ if bpdata(3,j)>=0, s=[' ' s];end
+ if abs(bpdata(3,j))<10, s=[' ' s];end
+ if abs(bpdata(3,j))<100, s=[' ' s];end
+ while length(s)<7, s=[s ' ']; end
+ amp2 = [amp2 s];
+ if j/6==round(j/6), amp2=[amp2 10 32]; end
+ end
+ smp2='';
+ for j=1:size(bpnames,1)
+ s = [num2str(bpdata(4,j)) ','];
+ if bpdata(4,j)>=0, s=[' ' s];end
+ if abs(bpdata(4,j))<10, s=[' ' s];end
+ while length(s)<7, s=[s ' ']; end
+ smp2 = [smp2 s];
+ if j/6==round(j/6), smp2=[smp2 10]; end
+ end
+ patmp2=[num2str(length(tok_data_struct.bpnames)),'*0.0'];
+
+ flnam=' ';
+ for j=1:size(flnames,1)
+ flnam = [flnam '''' deblank(flnames(j,:)) ''','];
+ if j/5==round(j/5), flnam=[flnam 10 32]; end
+ end
+ rsi=' ';
+ for j=1:size(flnames,1)
+ s = num2str(fldata(2,j)); while length(s)<5, s=[s '0']; end
+ rsi = [rsi s ', '];
+ if j/6==round(j/6), rsi=[rsi 10 32 32]; end
+ end
+ zsi=' ';
+ for j=1:size(flnames,1)
+ s = num2str(fldata(1,j));
+ if fldata(1,j)==0, s='0.';end
+ if fldata(1,j)>=0, s=[' ' s];end
+ while length(s)<6, s=[s '0']; end
+ zsi = [zsi s ','];
+ if j/6==round(j/6), zsi=[zsi 10 32]; end
+ end
+
+
+ fcdat = '';
+ for ii=[1:size(fcdata,2)]
+ fcdat = strvcat(fcdat,sprintf('%-12.4f%-12.4f%-12.4f%-12.4f%-12.4f%-12.4f\n',fcdata([2 1 4 3 5 6],ii)));
+ end
+ if use_vv
+ for ii=[1:nv]
+ fcdat = strvcat(fcdat,sprintf('%-12.4f%-12.4f%-12.4f%-12.4f%-12.4f%-12.4f\n',vvdata([2 1 4 3 5 6],ii)));
+ end
+ end
+
+
+ % a holds the characters that go in the file
+ a = [ ...
+ ' $IN3' 10 ...
+ ' RLEFT=',num2str(min(rg)),' RRIGHT=',num2str(max(rg)) 10 ...
+ ' ZBOTTO=',num2str(min(zg)),' ZTOP=',num2str(max(zg)), ...
+ ' NSMP2=25' 10 ...
+ ' IGRID=1 IFCOIL=1 IVESEL=0 ISLPFC=1 IECOIL=0 IACOIL=0' 10];
+
+ if ~exist('fcid','var')
+ fcid = [1:tok_data_struct.nc];
+ end
+ if ~use_vv
+ vvid = [];
+ end
+ if ~exist('vvid','var')
+ vvid = [1:tok_data_struct.nv];
+ end
+
+ [FCID, FCTURN, TURNFC] = make_fc_strings(tok_data_struct,fcid,vvid);
+
+
+ a = [a, ...
+ FCID ...
+ FCTURN ...
+ TURNFC ...
+ ' LPNAME=' 10 ...
+ flnam ...
+ '' 10 ...
+ ' MPNAM2=' 10 ...
+ bpnam ...
+ '' 10 ...
+ ' XMP2=' 10 ...
+ xmp2 ...
+ '' 10 ...
+ ' YMP2=' 10 ...
+ ymp2 ...
+ '' 10 ...
+ ' RSI=' 10 ...
+ rsi ...
+ '' 10 ...
+ ' ZSI=' 10 ...
+ zsi ...
+ '' 10 ...
+ '' 10 ...
+ ' AMP2=' 10 ...
+ amp2 ...
+ '' 10 ...
+ ' SMP2= ' 10 ...
+ smp2 ...
+ '' 10 ...
+ [' PATMP2 = ' patmp2 10] ...
+ '' 10 ...
+ '' 10 ...
+ ' $' 10 ...
+ ];
+
+
+ f=fopen('mhdin.dat','w');
+ fwrite(f,a);
+ fwrite(f,fcdat');
+ fclose(f);
+ disp('Wrote file mhdin.dat')
+end
+
+
+function [fcidStr, fcturnStr, turnfcStr] = make_fc_strings(tds,fcid,vvid)
+ allid = [fcid max(fcid)+vvid];
+ allTurn = [tds.ccnturn' ones(size(vvid))];
+ %uniqueGroups = unique(abs(fcid));
+ uniqueGroups = unique(abs(allid));
+ nGroup = length(uniqueGroups);
+ nTurnGroup = zeros(1,nGroup);
+ for ii=1:nGroup
+ idxGroup = find(abs(allid) == uniqueGroups(ii));
+ nTurnGroup(ii) = sum(allTurn(idxGroup));
+ end
+ for ii=1:length(allid)
+ fcturn(ii) = sign(allid(ii))*allTurn(ii)/nTurnGroup(abs(allid(ii)));
+ end
+
+ fcidStr = ' FCID=';
+ for ii=allid(:)'
+ fcidStr = [fcidStr num2str(ii) '.,'];
+ end
+ fcidStr = [fcidStr 10];
+
+ fcturnStr = ' FCTURN=';
+ for ii=1:length(fcturn)
+ if fcturn(ii) < 1
+ fcturnStr = [fcturnStr num2str(abs(fcturn(ii))) ','];
+ else
+ fcturnStr = [fcturnStr num2str(abs(fcturn(ii))) '.,'];
+ end
+ end
+ fcturnStr = [fcturnStr 10];
+
+ turnfcStr = ' TURNFC=';
+ for ii=1:nGroup
+ if ismember(uniqueGroups(ii),fcid)
+ turnfcStr = [turnfcStr num2str(nTurnGroup(ii)) '.,'];
+ else
+ turnfcStr = [turnfcStr '1.,'];
+ end
+ end
+ turnfcStr = [turnfcStr 10];
+end
+
+
diff --git a/matlab/D3D/efit/private/plasma_current.m b/matlab/D3D/efit/private/plasma_current.m
new file mode 100644
index 00000000..4971c41a
--- /dev/null
+++ b/matlab/D3D/efit/private/plasma_current.m
@@ -0,0 +1,200 @@
+ function [pcur,jphi,cpasma]= plasma_current(psizr, pprime, ffprim,...
+ psimag,psibnd,...
+ rgefit,zgefit,...
+ rbbbs,zbbbs,nbbbs,...
+ iplot,dofast,cpasma)
+% Function generates plasma current (pcurrt) from flux (psizr)
+%
+% SYNTAX:
+% [pcur,jphi,cpasma]= plasma_current(psizr, pprime, ffprim,...
+% psimag,psibnd,...
+% rgefit,zgefit,...
+% rbbbs,zbbbs,nbbbs,...
+% iplot,dofast)
+%
+% PURPOSE: Generate plasma current from flux and current profile information
+% Input is same as read_gfile input or mds_eq => eq_ga_env
+% pcur is almost identical to that obtained using iecurr=2 in EFIT
+% (error <~ 0.2% and max near boundary)
+%
+% INPUTS: [default] <optional>
+%
+% psizr= true flux on grid (Wb)
+% pprime= p' vector (pprime from read_gfile) (nt/m^2)/(Vs/rad)
+% ffprim= ff' vector (ffprim in read_gfile) (mT)^2/(Vs/rad)
+% psimag= Flux on axis (Wb) [max(max(psizr))]
+% psibnd= Flux on boundary (Wb) [average of flux on boundary]
+% rgefit= Radius vector of grid (m)
+% zgefit= Z vector of grid (m)
+% rbbbs= Radius vector of boundary (m)
+% zbbbs= Z vector of boundary (m)
+% nbbbs= number of elements in rbbbs zbbbs [length(rbbbs]
+% iplot= 2; % 1= plot boundary, 2= contour jphi; 3= clabel [0]
+% dofast= 1; % increases algorithm speed at risk of including some
+% private flux in solution. Only extreme shapes like
+% crescents and beans are expected to need dofast=0.
+% Default is dofast=0 since speed is still very fast. [0]
+%
+% <cpasma>= Plasma Current (A)
+%
+% OUTPUTS:
+% pcur = plasma current on nr,nz grid [A]
+% jphi= plasma current density on nr,nz grid [MA/m^2]
+% cpasma= plasma current = sum(sum(pcur))
+%
+% Warning messages displayed if sum(sum(pcur)) and existing cpasma
+% are widely different. cpasma then overwritten.
+%
+% NOTE: ALL UNITS ARE MKS (except jphi is MA/m^2)
+% CAUTION: Overwrites pcur, jphi & cpasma if read_gfile read from gfile
+%
+% RESTRICTIONS:
+%
+% SEE ALSO: inside_plasma
+
+% WRITTEN BY: Jim Leuer ON 6-3-03
+%
+% Routines Used: inside_plasma
+%
+% MODIFICATION HISTORY:
+% jal30sep10: +cpasma normalization, spline fit; garea per efit
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% @(#)plasma_current.m 1.4 07/15/09
+
+% Default:
+% Note: cpasma is 12 argument
+
+ if nargin==11
+ dofast= 0;
+ elseif nargin==10
+ iplot=0;
+ dofast= 0;
+ elseif nargin==9
+ nbbb= length(rbbb);
+ elseif nargin<=8
+ disp('%ERROR plasma_current: Needs at least 9 arguments to work')
+ return
+ end
+
+ if isempty(psimag),
+ disp('% CAUTION: plasma_current: psimag is empty using max on grid')
+ psimag= max(max(psizr)),
+ end
+
+ if isempty(psibnd) % average on boundary
+ disp('% CAUTION: plasma_current: psibnd is empty using ave. on boundary')
+ temp = interp2(rgefit,zgefit,psizr,rbbbs,zbbbs);
+ idx = find(isnan(temp));
+ idx = setdiff(1:length(temp),idx);
+ temp = temp(idx);
+ psibnd= mean(temp);
+ end
+ if exist('iplot')~=1, iplot=0; end % no plotting
+ if exist('dofast')~=1, dofast=0; end % 100% guaranteed mapping
+
+% ====================================================
+% Find Plasma Grid that is "INSIDE" plasma Boundary:
+
+ inside= inside_plasma(psizr,rgefit,zgefit,rbbbs,zbbbs,nbbbs, psimag,psibnd,...
+ iplot,dofast);
+
+% ==================================================================
+% Compute (see GA-A14490, pg8): Jphi= r*Pprime + F*Fprime/(amu0*r)
+
+ garea= (rgefit(2)-rgefit(1))*(zgefit(2)-zgefit(1)); % EFIT area def
+ areaomu= garea/(4*pi*1e-7);
+ [nz,nr]= size(psizr);
+ jphi= zeros(nz,nr); % initialize to zero current
+ pcur= jphi; % Plasma current on grid (initial 0)
+ rg= ones(nz,1)*rgefit'; % R on grid (same size as psizr)
+ id= find(inside); % work only on inside plasma points
+ if(isempty(id)) % handle badly converged EFIT
+ cpasma = 0;
+ return;
+ end
+ rr= rg(id); % reduced set of R's Inside plasma
+ psi= linspace(psimag,psibnd,length(pprime))'; % from center to edge
+ pp= interp1(psi,pprime,psizr(id),'spline'); % Pprime
+ ffp= interp1(psi,ffprim,psizr(id),'spline'); % F*Fprime
+ pcur(id)= garea*rr.*pp + areaomu*ffp./rr; % Toroidal current density (inside)
+ cpas= sum(pcur(id));
+ if exist('cpasma')==1 & ~isempty(cpasma)
+ if abs((cpas-cpasma)/cpasma) > 0.01
+ disp(['%CAUTION: plasma_current sum(pcur)-cpasma>1.0% want,got: ',...
+ num2str(cpasma),' ',num2str(cpas)])
+ disp(['% Re-Normalizing current distribution to cpasma: ', ...
+ num2str(cpasma)])
+ end
+ pcur(id)= pcur(id)*cpasma/cpas;
+ cpas= cpasma;
+ end
+
+ jphi(id)= pcur(id)/(garea*1e+6); % Convert to MA/m^2 (daves Convention)
+ cpasma= cpas; % set plasma current to sum of grid current
+
+ if iplot>=2
+ hold on
+ [c,h]= contour(rgefit,zgefit,jphi);
+ if iplot>=3
+ clabel(c,h)
+ end
+ hold on
+ plot(rbbbs(1:nbbbs),zbbbs(1:nbbbs),'k');
+ hold on
+ axis image
+ xlabel('R [m]')
+ ylabel('Z [m]')
+ title([' Plasma Current Density [MA/m^2], Ip= ',...
+ num2str(1e-6*cpasma), ' MA']);
+ end
+
+ return
+
+% =============================================================================
+% Testing (FOR MAIN TESTING SEE Leuer's efit area: test_plasma_current)
+
+% filename= '/u/leuer/efit/diiid/g110214.01740'; % DIIID file with: iecurr=2
+ filename= '/u/leuer/efit/diiid/shot113697/g113697.03000'; %
+ read_gfile
+ cpasma0= cpasma;
+ if exist('pcurrt')
+ pcurrt0= pcurrt; % save gfile currents for later error analysis
+ end
+
+ iplot=3;
+ dofast=0;
+
+% do plasma_current
+
+ figure(3)
+ clf
+ iplot=3;
+ [pcur,jphi,cpasma]= plasma_current(psizr, pprime, ffprim,...
+ psimag,psibnd,...
+ rgefit,zgefit,...
+ rbbbs,zbbbs,nbbbs,...
+ iplot,dofast);
+% do some checks to see error
+ if exist('pcurrt0')==1
+ figure(7)
+ clf
+ colormap('default');
+ cmap= colormap;
+ cmap= [cmap(1:end-1,:);[1 1 1]];
+ colormap(cmap);
+ err1= zeros(size(pcurrt0));
+ id= find(pcurrt0~=0);
+ err1(id)= (pcur(id)-pcurrt0(id))./pcurrt0(id);
+ contourf(rgefit,zgefit,100*err1)
+ hold on
+ plot(rbbbs(1:nbbbs),zbbbs(1:nbbbs),'k');
+ hold on
+ axis image
+ title(['%Error= 100*(pcur-pcurrt0)/pcurrt0, Max value= ',...
+ num2str(100*max(max(abs(err1)))),' %'])
+ colorbar('vertical')
+
+% Results: Max error for 98549.04000: 0.11% toward edge Most error <.04%
+ end
+
diff --git a/matlab/D3D/efit/private/read_mhdin.m b/matlab/D3D/efit/private/read_mhdin.m
new file mode 100644
index 00000000..2096cbd0
--- /dev/null
+++ b/matlab/D3D/efit/private/read_mhdin.m
@@ -0,0 +1,392 @@
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% USAGE: >>read_mhdin
+%
+% PURPOSE: reads efit file mhdin.dat. All variables are added in UPPER CASE
+%
+% INPUTS: [default]
+% filename = ['mhdin.dat']; % EFIT mhdin.dat file name Can contain directory
+%
+% Note: below dimensions are needed only if geometry is read below namelist
+% defaults for D3D
+%
+% nfcoil = [18]; % # of F-coils; must be specified for non-D3D (Opt.)
+% defaults to D3D values of nfcoil=18
+% necoil = [122];% # number of E-coils elements specify for non-D3D (Opt.)
+% defaults to D3D value of nesum=6
+% nsilop = [41]; % # flux loops, needed if not read in namelist (opt.)
+% nvesel = [24]; % # vessels, needed if not read in namelist (opt.)
+%
+% OUTPUTS:
+% all name list items: IGRID, RLEFT, RRIGHT, ... FCID, FCTURN, ...
+% PF coil geometry if nfcoil>0,
+% E-coil geometry if nesum>0
+% ireadok = flag to report good read of gfile (0=bad, 1=good)
+%
+% CAUTION: Namelist is read with to_upper=+1 which makes all variables UPPER
+% CASE.
+% New read_namelist in /thor/leuer/matlab/util/read_namelist.m
+% allows change of case. We want to read and make all variables UC
+%
+% To change UPPER CASE to lower case use:
+% mhdin_list= mk_uc_lc_var; % makes all lower case variables (nothing to UPCASE)
+% mhdin_list= mk_uc_lc_var('-remove_upper','*'); % remove upper case
+
+% Generated by Jim Leuer 03/25/2004
+% taken from frotran of efund.for
+% updated to all UPPER CASE 6/2/2004
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Prelims:
+
+% if exist('nfcoil')~=1, nfcoil=18, end %default to D3D values...
+% if exist('necoil')~=1, necoil=122, end
+% if exist('nsilop')~=1, nsilop=41, end
+% if exist('nvesel')~=1, nvesel=24, end
+
+ to_upper= +1; % this is now default read into all UPPER CASE
+
+% -------------------------- Open File: filename
+
+ if exist('filename')~=1
+ filename= 'mhdin.dat'
+ end
+
+ if ~isstr(filename)
+ disp([' %ERROR Filename must be a string ',filename]);
+ ireadok= 0;
+ return
+ end
+
+ if exist(filename)~=2
+ disp([' %ERROR Filename doesnt exist ',filename]);
+ ireadok= 0;
+ return
+ end
+
+ fid= fopen(filename, 'r');
+
+ if fid == (-1)
+ disp([' %ERROR read_mhdin: Couldn''t open file ',filename]);
+ ireadok= 0;
+ return
+ end
+
+ ireadok= 1;
+ fclose(fid);
+
+% =============================================================
+% Start
+ read_namelist; % new version works for most namelist types
+
+% NOT USED SINCE TOUPPER in read_namelist CONTROLS CASE
+% Convert all upper case to lower case
+% disp('% NOTE: We are making lower case variables out of ALL UPPER CASE')
+% disp('% NOTE: We are REMOVING ALL UPPER CASE VARIABLES')
+
+% mhdin_list= mk_uc_lc_var; % makes all lower case variables (nothing to UPCASE)
+% mhdin_list= mk_uc_lc_var('-remove_upper','*'); % remove upper case
+
+%
+% =============================================================
+% Now start reading data at end of namelist
+
+% if (rf(1).lt.0.0)
+% .read (nin,10000) (rf(i),zf(i),wf(i),hf(i),af(i),af2(i),
+% . i=1,nfcoil)
+% 10000 format (6e12.6)
+
+
+% open filename and find end of namelist
+ fid= fopen(filename, 'r');
+
+
+% find start of namelist
+ nlcnt = 0;
+ line1 = fgets(fid);
+ totstr = '';
+ if(~isempty(findstr(line1,'$')) | ~isempty(findstr(line1,'&')))
+% fprintf('Begin processing %s\n',line1)
+ nlcnt = nlcnt+1;
+ temp = deblank(line1);
+ if(size(nlnames,2)~=0)
+ if(length(temp)>size(nlnames,2))
+ bb = blanks(length(temp)-size(nlnames,2));
+ nlnames = [nlnames bb*ones(size(nlnames,1))];
+ end
+ end
+ nlnames(nlcnt,1:length(temp)) = temp;
+ while line1~=-1
+ line1 = fgets(fid);
+
+% find end of namelist
+ idx = union(findstr(line1,'$'),findstr(line1,'&'));
+ if(~isempty(idx))
+ break
+ end;
+ if(length(line1)>1 & strcmp(line1(2),'/'))
+ break
+ end;
+ end
+ end
+
+ if isempty(findstr(upper(nlnames(1,:)),'IN3'))
+ disp('%WARNING: couldnt find namelist IN3 in the input file')
+ end
+
+% ===========================
+% Read F-coil
+% ===========================
+% next lines should be F-coil
+% if (rf(1).lt.0.0)
+% .read (nin,10000) (rf(i),zf(i),wf(i),hf(i),af(i),af2(i),
+% . i=1,nfcoil)
+% 10000 format (6e12.6)
+
+ disp('% Looking after namelist for other data: F-coil, Flux, B-probe, E, VV')
+
+% see if these were written in name list and make all dummy variables not read
+ if exist('RF')==1 & exist('HF')==1
+ if exist('WF')~=1
+ WF= zeros(size(RF));
+ end
+ if exist('HF')~=1
+ HF= zeros(size(RF));
+ end
+ if exist('AF')~=1
+ AF= zeros(size(RF));
+ end
+ if exist('AF2')~=1
+ AF2= zeros(size(RF));
+ end
+ end
+
+ if exist('RF')~=1 RF= -1; end
+ if RF(1) < 0.0
+ disp(['% reading ',int2str(nfcoil),' PF coils: rf,zf,wf,hf,af,af2'])
+ clear RF ZF
+ RF= zeros(nfcoil,1);
+ ZF= RF;
+ HF= RF;
+ AF= RF;
+ AF2= RF;
+ for ii=1:nfcoil
+ str= fgetl(fid);
+ dum= sscanf(str,'%f');
+ len= length(dum);
+ RF(ii,1)= dum(1);
+ ZF(ii,1)= dum(2);
+ if len>=3 WF(ii,1)= dum(3); end
+ if len>=4 HF(ii,1)= dum(4); end
+ if len>=5 AF(ii,1)= dum(5); end
+ if len>=6 AF2(ii,1)= dum(6); end
+ end % do
+ disp('F-COIL RF,ZF,WF,HF,AF,AF2')
+ disp([RF,ZF,WF,HF,AF,AF2])
+ end
+
+% ===========================
+% Read Flux Loops
+% ===========================
+% if (rsi(1).lt.0.0)
+% .read (nin,10000) (rsi(i),zsi(i),wsi(i),hsi(i),as(i),as2(i),
+% . i=1,nsilop)
+
+% see if UC was read in namelist
+ if exist('RSI')==1 & exist('ZSI')==1
+ if exist('WSI')~=1
+ WSI= zeros(size(RSI));
+ end
+ if exist('HSI')~=1
+ HSI= zeros(size(RSI));
+ end
+ if exist('AS')~=1
+ AS= zeros(size(RSI));
+ end
+ if exist('AS2')~=1
+ AS2= zeros(size(RSI));
+ end
+ if exist('LPNAME')~=1
+ LPNAME= int2str((1:length(RSI))');
+ end
+ end
+
+ if exist('RSI')~=1 RSI(1)=-1; end
+ if RSI(1) < 0.0
+ disp(['% reading ',int2str(nsilop),' flux loops: rsi,zsi,wsi,hsi,as,as2'])
+ clear RSI ZSI
+ RSI= zeros(nsilop,1);
+ ZSI= RSI;
+ WSI= RSI;
+ HSI= RSI;
+ AS= RSI;
+ AS2= RSI;
+ for ii=1:nsilop
+ str= fgetl(fid);
+ dum= sscanf(str,'%f');
+ len= length(dum);
+ RSI(ii,1)= dum(1);
+ ZSI(ii,1)= dum(2);
+ if len>=3 WSI(ii,1)= dum(3); end
+ if len>=4 HSI(ii,1)= dum(4); end
+ if len>=5 AS(ii,1)= dum(5); end
+ if len>=6 AS2(ii,1)= dum(6); end
+ end % do
+ disp('flux loops: RSI,ZSI,WSI,HSI,AS,AS2')
+ disp([RSI,ZSI,WSI,HSI,AS,AS2])
+ end
+
+
+% ===========================
+% Read Ecoil & Vessel
+% ===========================
+% if ((iecoil.gt.0).or.(ivesel.gt.0)) then
+% if (re(1).lt.0.0)
+% .read (nin,10020) (re(i),ze(i),we(i),he(i),ecid(i),i=1,necoil)
+
+% if (ivesel.gt.0.and.rvs(1).lt.0.0) then
+% if (wvs(1).lt.0.0) then
+% read (nin,10000) (rvs(i),zvs(i),wvs(i),hvs(i),avs(i),avs2(i),
+% . i=1,nvesel)
+% else
+% do i=1,nvesel
+% read (nin,*) rvs(i),zvs(i)
+% enddo
+% endif
+% endif
+% endif
+
+ if exist('IECOIL')~=1 IECOIL= 0; end
+
+ if IECOIL>=1
+ if exist('RE')~=1 RE= -1; end % must read in ecoil
+ if RE(1) < 0
+ disp(['% Reading ',int2str(necoil),' Ecoil: re ze we he ecid'])
+ clear RE ZE
+ WE= zeros(necoil,1);
+ HE= WE;
+ ECID= (1:length(necoil))';
+ for ii=1:necoil
+ str= fgetl(fid);
+ dum= sscanf(str,'%f');
+ len= length(dum);
+ RE(ii,1)= dum(1);
+ ZE(ii,1)= dum(2);
+ if len>=3 WE(ii,1)= dum(3); end
+ if len>=4 HE(ii,1)= dum(4); end
+ if len>=5 ECID(ii,1)= dum(5); end
+ end % end do
+ disp('Ecoil: RE, ZE, WE, HE, ECID')
+ disp([RE, ZE, WE, HE, ECID])
+ end
+ end
+
+ if exist('IVESEL')~=1 IVESEL= 0; end
+
+ if IVESEL>=1
+ if exist('RVS')~=1
+ RVS= -1; % must read in vessel
+ WVS= -1; % parallogram
+ end
+ if RVS(1) < 0
+ disp(['% Reading ',int2str(nvesel),' Vessels: rvs zvs wvs hvs avs avs2'])
+ clear RVS ZVS
+ WVS= zeros(nvesel,1);
+ HVS= WVS;
+ AVS= WVS;
+ AVS2= WVS;
+ for ii=1:nvesel
+ str= fgetl(fid);
+ dum= sscanf(str,'%f');
+ len= length(dum);
+ RVS(ii,1)= dum(1);
+ ZVS(ii,1)= dum(2);
+ if len>=3 WVS(ii,1)= dum(3); end
+ if len>=4 HVS(ii,1)= dum(4); end
+ if len>=5 AVS(ii,1)= dum(5); end
+ if len>=6 AVS2(ii,1)= dum(6); end
+ end % end do
+ disp('Vessels: RVS,ZVS,WVS,HVS,AVS,AVS2')
+ disp([RVS,ZVS,WVS,HVS,AVS,AVS2])
+ end
+ end
+
+% ========================================================
+% Do Some ERROR Checking
+ if exist('RSI')~=1
+ disp('% ERROR: No Flux loops in mhdin.dat file')
+ end
+ if exist('XMP2')~=1
+ disp('% ERROR: No Bprobes in mhdin.dat file')
+ end
+ if exist('RF')~=1
+ disp('% ERROR: No PF Coils in mhdin.dat file')
+ end
+
+% ========================================================
+% Do some Plotting:
+ if exist('RF') & exist('AF2')
+ plot_box(RF,ZF,WF,HF,'r',AF,AF2)
+ hold on
+ end
+ if exist('RE') & exist('HE')
+ plot_box(RE,ZE,WE,HE,'k')
+ hold on
+ end
+ if exist('RVS') & exist('AVS')
+ plot_box(RVS,ZVS,WVS,HVS,'r',AVS,AVS2)
+ hold on
+ end
+ if exist('RSI') & exist('ZSI')
+ plot(RSI,ZSI,'xb')
+ hold on
+ end
+ if exist('XMP2') & exist('YMP2')
+ plot(XMP2,YMP2,'og')
+ hold on
+ end
+ if exist('RLEFT') & exist('ZBOTTO')
+ plot([RLEFT RRIGHT RRIGHT RLEFT RLEFT],[ZBOTTO ZBOTTO ZTOP ZTOP ZBOTTO],'k')
+ hold on
+ end
+ axis image
+ grid
+ title([fix_undscr(filename)])
+ ylabel('R')
+ xlabel('Z')
+
+ return
+% ========================================================
+
+% ========================================================
+% NOW generate some zero arrays if they dont exist
+ if exist('RSI')~=1
+ if exist('HSI')~=1 HSI= zeros(size(RSI)); end
+ if exist('WSI')~=1 WSI= zeros(size(RSI)); end
+ if exist('AS')~=1 AS= zeros(size(RSI)); end
+ if exist('AS2')~=1 AS2= zeros(size(RSI)); end
+ if exist('HSI')~=1 HSI= zeros(size(RSI)); end
+ end
+
+ if exist('XMP2')~=1 XMP2= zeros(nbprobe,1); end
+ if exist('AMP2')~=1 AMP2= zeros(size(XMP2)); end
+ if exist('SMP2')~=1 SMP2= zeros(size(XMP2)); end
+ if exist('MPNAM2')~=1 MPNAM2= int2str((1:length(XMP2))'); end
+ if exist('PATMP2')~=1 PATMP2= zeros(size(XMP2)); end
+ if exist('MPNAM2')~=1 MPNAM2= int2str((1:length(XMP2))'); end
+
+% Probably will need to add more default arrays here
+
+ return
+
+% =====================================================================
+% testing:
+% tested for east ok
+% Not tested for any other mhdin.dat but should be reasonable
+ filename='/u/leuer/efit/east/green/run/ef20040325.est'
+ nfcoil=14
+ read_mhdin
+
+
+
+
+
diff --git a/matlab/D3D/efit/private/run_efit.m b/matlab/D3D/efit/private/run_efit.m
new file mode 100644
index 00000000..47ad68b8
--- /dev/null
+++ b/matlab/D3D/efit/private/run_efit.m
@@ -0,0 +1,93 @@
+
+ function run_efit(shot,tims,dtms,nslices,snap,snapidx,efitver,efitdir)
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% SYNTAX: run_efit(shot,tims,dtms,nslices,snap,snapidx,efitver,efitdir)
+%
+% PURPOSE:
+% Run EFIT from inside Matlab: uses basic input format just as in
+% standard invokation from command line (so will perform many slices
+% as usual with starting time, delta time, nslices).
+%
+% INPUTS:
+% shot = shot number
+% tims = starting time (in ms) at which to calculate EFIT
+% dtms = delta time (in ms) between each slice
+% nslices = # of slices to calculate
+% snap = (optional) string specifying:
+% if snapidx=2: data file to read input data from
+% if snapidx=3: use snap file in default directory (or in /link)
+% if snapidx=7: snap file extension (eg 'jta_t')
+% snapidx = (optional) EFIT input mode (2=input file, 3=snap file)
+% Default=3 (used if nargin<6). Note if snapidx=2, the shot,
+% tims,dtms,nslices inputs are ignored, and snap specifies the
+% input file to use.
+% efitver = (optional) string specifying EFIT version (eg 'efitd65yd' or 'efitd6565d')
+% efitdir = (optional) string specifying directory in which EFIT version
+% is to be found (default = /link/efit/ or /d/linux/efit/)
+%
+% OUTPUTS:
+% None
+%
+% RESTRICTIONS:
+%
+% METHOD:
+% Uses Matlab unix() command to execute command string built in
+% this function for desired inputs.
+
+% WRITTEN BY: Dave Humphreys ON 11/1/00
+%
+% VERSION @(#)run_efit.m 1.1 08/27/10
+%
+% MODIFICATION HISTORY: ASW 4/30/04 runs on either thor or uscws9 transparently
+% ASW 8/27/10 changed from int2str to num2str
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Define defaults:
+ if nargin<5
+ snap=''; %no default: use snap file in default dir...
+ snapidx = 3; %select snap option "3" in efit call (use file in def dir)
+ end
+ if nargin==5
+ snapidx = 7; %select snap option "7" in efit call (use snap file ext)
+ end
+ if nargin<7, efitver='efitd65yd'; end %default version (recommend efitd6565d for 65x65 grid)
+ if nargin<8
+ efitdir='/link/efit/'; % default efit directory for hp
+ % Check to see if we are on Thor (a linux machine)
+ thordisk = getenv('THORDISK');
+ if length(thordisk)==5
+ if thordisk=='/home'
+ efitdir='/d/linux/efit/';
+ end
+ end
+ end
+
+% Prelims:
+ mu0 = 0.4*pi;
+
+% Derived Values:
+
+% Build command string:
+ if snapidx==3 %use snap file input, use snap file "efit_snap.dat"
+ efitcom = ['echo ',' "','3\n', ...
+ num2str(shot),',',num2str(tims),',',num2str(dtms),',', ...
+ num2str(nslices),'" ','| ',efitdir,efitver];
+ end
+
+ if snapidx==7 %use snap file input, use snap file with extension snap
+ efitcom = ['echo ',' "','7\n', ...
+ snap,'\n', ...
+ num2str(shot),',',num2str(tims),',',num2str(dtms),',', ...
+ num2str(nslices),'" ','| ',efitdir,efitver];
+ end
+
+ if snapidx==2 %use input file input, input file specified by snap variable
+ efitcom = ['echo ',' "','2\n', ...
+ '1\n',snap,'\n','" ','| ',efitdir,efitver];
+ end
+
+% Execute EFIT:
+ unix(efitcom);
+
+
diff --git a/matlab/D3D/efit/private/shot_from_gfile.m b/matlab/D3D/efit/private/shot_from_gfile.m
new file mode 100644
index 00000000..24bf7df0
--- /dev/null
+++ b/matlab/D3D/efit/private/shot_from_gfile.m
@@ -0,0 +1,125 @@
+ function [shot,tms,dr,server] = shot_from_gfile(gfile)
+ %
+% SYNTAX:
+% [shot,tms,dr] = shot_from_gfile(gfile) % pointer to file
+% [shot,tms,tree,server] = shot_from_gfile(gfile) % pointer to MDS+ tree
+%
+% PURPOSE: look at string in gfile name and extract shot, tms, etc
+%
+% INPUT: <default>
+% gfile = string contining gfile name or "dot" pointer to MDS tree ex:
+% '/u/leuer/efit/d3d/shot131498/g131498.02000' % file pointer
+% 'D3D.EFITRT1.g131498.02000' % MDS+ pointer
+%
+% OUTPUT: <default>
+% shot = shot number
+% tms = shot time [ms]
+% dr = directory if pointer to gfile
+% tree = EFIT Tree if pointer to MDS tree <'EFIT01'>
+% server= MDS server for efit data <'D3D'>
+
+% NOTE:
+% FULL Pointer to efit trees is: '\D3D::TOP.MHD.EFIT'
+
+% WRITTEN BY: Jim Leuer ON 17Mar2010
+% ==========================================================================
+
+ if nargin==0
+ disp('% shot_from_gfile needs at least a "gfile" argument')
+ help shot_from_gfile
+ return
+ end
+
+ id= strfind(gfile,'.'); % "dot" Format: SERVER.TREE.SHOTNUM.TIME
+
+ if isempty(id)
+ disp(['%ERROR: shot_from_gfile: minimum form "shot.tms" ',gfile])
+ shot=[]; tms=[]; dr=[]; server=[];
+
+ else % isempty(id)
+
+ [dr, nam, ext]= parse_filename(gfile);
+
+% ---------------------
+% 1) Time in ms
+ id1= strfind(ext,'.');
+ if length(ext)~= 6 | isempty(id1)
+ disp([' %ERROR shot_from_gfile file extention needs format ".02600" ',ext])
+ end
+ tms= str2num(ext(2:end));
+ if tms<0 | tms>99999
+ disp([' %ERROR shot_from_gfile:TIME extention error ',ext]);
+ tms=[];
+ end
+
+% ---------------------
+% 2) Shot Number
+ sshot= nam(end-5:end); % string of shotnum last 6 digits of nam
+ shot= sscanf(sshot,'%d');
+ if shot<0 | shot>999999
+ disp([' %ERROR shot_from_gfile:SHOT number error ',sshot]);
+ shot=[];
+ end
+
+ if length(id)<=1 % gfile points to file directory
+
+% ---------------------
+% 3a) Directory
+ server=[];
+ if exist(gfile) ~= 2
+ disp(['%WARNING: shot_from_gfile: gfile doesnt exist ',gfile])
+ end
+
+ else % if length(id)<=1
+
+% ---------------------
+% 3b) tree = optional <EFIT01>
+ if length(id)<=2 id0= 0; else id0= id(end-2); end
+ idd= id0+1:id(end-1)-1;
+ if ~isempty(idd)
+ tree= gfile(idd); % this is TREE
+ else
+ tree= 'EFIT01';
+ disp([' % NOTE shot_from_gfile:is setting TREE to ',tree]);
+ end
+ dr= tree;
+
+% ---------------------
+% 4) server = optional <d3d>
+ if length(id)<=2
+ server= 'D3D';
+ else
+ if length(id)<=3 id0= 0; else id0= id(end-3); end
+ server= gfile(id0+1:id(end-2)-1);
+ end
+% ----------------------
+% Could check HERE to see if tree exists in MDS+
+% Can do for d3d using get_mdsefit_nms
+% fln=[server '.' tree '.' int2str(shot) '.' stms];
+ if strmatch(upper(server),'D3D') & ~isempty(shot) & ...
+ exist('get_mdsefit_nms')==2
+ [efit_nms] = get_mdsefit_nms(shot);
+ if isempty(strmatch(upper(tree),efit_nms))
+ disp([' % CAUTION shot_from_gfile cant find D3D tree in mds+ ' tree])
+ efit_nms
+ end
+ end
+ end % if length(id)<=1
+
+ end % isempty(id)
+
+ return
+% =========================================================
+% testing
+ gfile = '/u/leuer/efit/d3d/shot131498/g131498.02600'
+ gfile = '.g123216.03000'
+ gfile = '.D3D.EFITRT1.131498.02600'
+ gfile = '.D3D.dum.131498.02600'
+ gfile = '131498.02600'
+ gfile = '.xxx.EFITRT1.131498.02600'
+
+ [dr, nam, ext]= parse_filename(gfile)
+
+ [shot,tms,dr] = shot_from_gfile(gfile) % pointer to file
+ [shot,tms,tree,server] = shot_from_gfile(gfile) % pointer to file
+
diff --git a/matlab/D3D/efit/private/std_efit_units.m b/matlab/D3D/efit/private/std_efit_units.m
new file mode 100644
index 00000000..784f1ce8
--- /dev/null
+++ b/matlab/D3D/efit/private/std_efit_units.m
@@ -0,0 +1,1238 @@
+function gfile_data = std_efit_units(gdata,tokamak,options)
+ %
+% SYNTAX: gfile_data = std_efit_units(gdata,tokamak)
+%
+% PURPOSE: Do device-dependent conversions of eqdsk data to a common
+% set of units and conventions.
+%
+% INPUT:
+% gdata = structure containing g eqdsk data
+% tokamak = name of device
+% options = parameters needed to customize logic for a particular device
+%
+% OUTPUT:
+% gfile_data = structure containing g eqdsk data with standard units (i.e., MA-turns)
+
+% RESTRICTIONS:
+%
+% METHOD:
+%
+% WRITTEN BY: Mike Walker ON 8/25/11
+% (extracted from read_gfile_tok.m)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% @(#)std_efit_units.m 1.16 07/12/12
+
+ gfile_data = gdata;
+ if(nargin>2)
+ struct_to_ws(options);
+ end
+ if(nargin<3 | ~isfield(options,'old_efit'))
+ old_efit = 0;
+ end
+
+ switch upper(tokamak)
+
+% ===========================================================================
+ case {'DIII-D','D3D','DIIID'}
+% ===========================================================================
+ if(~exist('nesum','var'))
+ nesum = 6;
+ end
+
+% include d3d build parameters (see build area => must be same)
+ gfile_data.fcturn = ones(1,18);
+ gfile_data.turnfc = [58 58 58 58 58 55 55 58 55 58 58 58 58 58 55 55 58 55];
+ gfile_data.fcid = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18];
+% NOTE: EFIT IS CONSTRUCTED with ecid having 1,2,3,4,5,6 numbers. Even are 1's
+% odds hook to 2
+ gfile_data.ecid = [ ... % must be same as ecdata(5,:)
+ 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 ...
+ 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 ...
+ 1 2 2 1 1 2 2 1 1 1 1 1 2 2 2 2 2 2 2 1 ...
+ 1 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 ...
+ 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 2 2 1 1 ...
+ 2 2 1 1 2 2 1 1 2 2 2 2 2 1 1 1 1 1 1 1 ...
+ 2 2];
+ gfile_data.ecturn = ones(122,1); % required if ecid is not empty
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+
+% Construct E-coil Number of turns to make output MA-turns:
+ dum= min(2,nesum); % if nesum=1 then use ONLY 1 ECOIL current
+
+ if ~isempty(gdata.brsp) & ~isempty(gdata.ecurrt)
+% Construct E-coil Number of turns to make output MA-turns:
+ dum= min(2,nesum); % if nesum=1 then use ONLY 1 ECOIL current
+% find reducted ecoil set data
+ nee = length(unique(gfile_data.ecid));
+ clear neturn
+ for ii=1:nee
+ idx = find(gfile_data.ecid==ii);
+ neturn(ii,1) = sum(gfile_data.ecturn(idx));
+ end
+
+ %cc-vector with 2 e-coil segments (MA-t):
+ cc2 = [gdata.ecurrt(1:dum).*neturn(1:dum);gdata.brsp]*1.e-6;
+
+% CAUTION: we are not using the 6segment e-coil anymore so this is crippled
+%cc-vector with all e-coil segs(2 or 6) (MA-t):
+% 6 ecoils: ECOILA ECOILB E567UP E567DN E89DN E89UP
+% (When only 2 ecoils, ECOILA = E567UP = E89UP, ECOILB = E567DN = E89DN.)
+% cc = [gdata.ecurrt;gdata.brsp]*1.e-6;
+ cc= cc2;
+
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+
+ gfile_data.gdef.cc = 'E/F coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'E/F coil currents in MA-turns (for 2-segment E-coil). CAUTION: cc2 has min(2,nesum) E-coil elements (could be 1 vs std: 2)';
+ end
+
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(~isfield(gdata,'gtime') & isfield(gdata,'time'))
+ gdata.gtime = gdata.time;
+ end
+ gfile_data.time=gdata.gtime*1e-3;
+ gfile_data.tms= gdata.gtime;
+ end
+
+% ===========================================================================
+ case {'EAST'}
+% ===========================================================================
+
+% Construct coil current vector:
+ ncadd=2; % fix cc so it has extra control coils
+ ncadd=0; % fix cc so it has extra control coils
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp*1e-6;zeros(ncadd,1)];
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+% gfile_data.fcturn = [1 1 1 0.17741935 0.82258065 1 1 ...
+% 1 1 1 0.17741935 0.82258065 1 1]';
+ gfile_data.fcturn = [1 1 1 44/(44+204) 204/(44+204) 1 1 ...
+ 1 1 1 44/(44+204) 204/(44+204) 1 1 1/2 -1/2]';
+ gfile_data.turnfc = [140 140 140 44 204 60 32 140 140 140 44 204 60 32 2 2]';
+ gfile_data.fcid = [1 2 3 4 4 5 6 7 8 9 10 10 11 12 13 13];
+ gfile_data.ecid = [];
+
+% EFIT order of coils is not the same as toksys order:
+ gfile_data.idx_efit_to_tok = [1 8 2 9 3 10 4 11 5 12 6 13 7 14 15 16];
+
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(~isfield(gdata,'gtime') & isfield(gdata,'time'))
+ gdata.gtime = gdata.time;
+ end
+ gfile_data.time = gdata.gtime;
+ gfile_data.tms = gdata.gtime*1e+3;
+ end
+
+% ===========================================================================
+ case {'KSTAR'}
+% ===========================================================================
+
+% select case based on gfile header ecase date which is from EFITD fortran
+ ecase= gdata.ecase;
+ icase = 0;
+ if ~isempty(findstr('03/16/2004',ecase))
+ icase = 0; % Leuer's Original KSTAR efit's
+ elseif ~isempty(findstr('01/23/2002',ecase)) % this is SABBAGH EFIT
+ dum= sscanf(ecase(27:38),'%f');
+ if ~isempty(dum)
+ shot= dum(1);
+ if length(dum)>=2
+ time= dum(2);
+ end
+ if shot<=3000 % ? guess where is shot transition from 2009 to 2010?
+ icase = 1; % Sabbagh efit 2008/2009
+ nfcoil= 14;
+ nves= 56;
+ nesum= 1;
+ else
+ icase = 3; % Sabbagh efit 2010
+ nfcoil= 14;
+ nves= 60;
+ nesum= 1;
+ end % if shot
+ end % if ~isempty(dum)
+ elseif (gdata.nh==33 & gdata.nw==33 & length(gdata.brsp)==32) % rtEFIT v7b (2011)
+ icase = 5;
+ nfcoil = 18;
+ nves = 70; % Like sabbagh, vv is listed as Fcoil
+ nesum = 1;
+ elseif (gdata.nh==33 & gdata.nw==33 & length(gdata.brsp)==42) % rtEFIT v7d2 (2012 active version w/ more passive els)
+ icase= 6;
+ nfcoil = 18;
+ nves = 180;
+ elseif (gdata.nh==33 & gdata.nw==33 & length(gdata.brsp)==37) % rtefit_v7e (2012 active version, series coils and slightly shifted divertor indices)
+ icase = 7;
+ nfcoil = 13;
+ nves = 180;
+ end
+
+% =============
+ switch icase
+% =============
+ case 0 % Original Leuer EFIT
+ ncadd=4; % fix cc so it has extra control coils
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp*1e-6;zeros(ncadd,1)];
+ cc2= cc;
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns ';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+ gfile_data.fcturn = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]';
+ gfile_data.turnfc = [180 144 72 108 208 128 72 180 144 72 108 208 128 72 6 4 6 4];
+ gfile_data.fcid = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18];
+ gfile_data.ecid = [];
+
+
+% =============
+ case 1 % Sabbagh efit 2009 => should use tok_sys.config_name= sab09
+% % PF Order is: % f1u:7u,7l,6l,5l...1l, must add 4 zero current IC's IC1U..
+ ncadd=4; % add in extra ic coils with zero current
+ gfile_data.fcid= 1:(nfcoil+4); % add in 4 IC's
+ gfile_data.turnfc= [180 144 72 108 208 128 72 72 128 208 108 72 144 180 6 4 6 4];
+ gfile_data.fcturn= ones(1,nfcoil+4);
+
+ % VV has 56elements broken into 6 cur groups. We expand 6 actual currents to 56 elements
+ % load /home/leuer/tokamaks/kstar/make/sab/vvid.mat
+ % disp(int2str(vvid'))
+ % load /home/leuer/tokamaks/kstar/make/sab/vvfrac.mat
+ vvid= [...
+ 1 1 2 2 2 2 2 2 3 3 3 3 3 4 4 5 5 5 5 5 6 6 6 6 6 6 1 1 ...
+ 1 1 2 2 2 2 2 2 3 3 3 3 3 4 4 5 5 5 5 5 6 6 6 6 6 6 1 1]';
+ gfile_data.vvfrac= ones(nves,1);
+
+ if ~isempty(gdata.brsp)
+ % EFIT brsp is terminal current (A) => Changing cc and cc2 are Mamp-turns
+ cc_terminal = [gdata.brsp(1:nfcoil)*1.e-6; zeros(ncadd,1)] ; % 1st 14 brsp => F-coils [A]
+ cc= cc_terminal.*gfile_data.turnfc'; % convert to MA-turns
+
+ % VV has 56elements broken into 6 cur groups. We expand 6 actual currents to 56 elements
+ proj= proj_turn(vvid,gfile_data.vvfrac);
+ vc6_terminal= gdata.brsp(nfcoil + (1:6)); % 6 currents after nfcoil are VV groups [A]
+ vc= proj'*vc6_terminal*1e-6; % expand 6 => 56 and change to MA-turn
+
+ if any( abs(gdata.brsp((nfcoil+6+1):end)) > 1e-5*max(abs(gdata.brsp)) )
+ wait(' % CAUTION: some Sabbagh EFIT brsp(21:60) have substantial current???')
+ end
+
+ gfile_data.cc= cc;
+ gfile_data.cc2= cc;
+ gfile_data.gdef.cc= 'PF coil currents in MA-turns; sab09; 14F & 4IC=0 ';
+ gfile_data.gdef.cc2= 'PF coil currents in MA-turns ';
+ gfile_data.vc= vc;
+ gfile_data.gdef.vc = 'vessel currents in MA; sab09; 6 groups => 56elm; 1 turn per elm.';
+ end
+ gfile_data.vvid= 1:nves; % NOTE: we are expanding 6=>56 so keep all VV elements
+ gfile_data.gdef.vvid = 'VV identifier; Note: sab09 has 6 VV groups but we maintain 56elm';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation; sab09 WILL CHANGE LATER ';
+
+% =============
+ case 2 % KI You efit 2009
+ ncadd=4; % fix cc so it has extra control coils holding place of IC
+ if ~isempty(gdata.brsp)
+ gfile_data.cc = [gdata.brsp*1e-6;zeros(ncadd,1)];
+ gfile_data.cc2 = gfile_data.cc;
+ end
+
+% =============
+ case 3 % Sabbagh efit 2010
+ % PF Order is: % f1u:7u,7l,6l,5l...1l, must add 4 zero current IC's IC1U..
+ ncadd=4; % add in extra ic coils with zero current
+ % Need to map Sab which is clockwise from inner midplane to toksysKSTAR which is CW top and CCW bottom
+
+ id= [1 2 3 4 5 6 7 14 13 12 11 10 9 8]; % EFIT/Sab PF coil convert to sab10.mat
+
+ gfile_data.fcid= 1:(nfcoil+4); % add in 4 IC's
+
+ % NOTE: Some how the conventions for turnfc and fcturn seems opposite of EFIT!!!!!!
+ % gfile order gfile_data.turnfc= [180 144 72 108 208 128 72 72 128 208 108 72 144 180 6 4 6 4]';
+ gfile_data.turnfc= [180 144 72 108 208 128 72 180 144 72 108 208 128 72 6 4 6 4]';
+ gfile_data.fcturn= ones(1,length(gfile_data.fcid));
+
+ % VV see /home/leuer/tokamaks/kstar/efit/sab10/
+ vvid0= [...
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ... % Outer VV
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ]'; % Inner VV
+ vvid= [vvid0; (13:36)']; % add in 10 passive plates and 14 inner limiter = 24 elements
+ gfile_data.vvid= vvid;
+
+ vvfrac0= [...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.2500 0.2500 0.1250 0.1250 0.1250 0.1250 0.2500 ...
+ 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 0.1250 ...
+ 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 ...
+ 0.2500 0.2500]'; % this is Outer VV
+
+ vvfrac0= [vvfrac0; vvfrac0]; % add in Inner VV which is same
+
+ vvfrac= [vvfrac0; ones(24,1)]; % add in 10 PP + 14 IL
+
+ if ~isempty(gdata.brsp)
+ % EFIT brsp in SABBAGH EFIT is terminal current (A) => Changing cc and cc2 are Mamp-turns
+
+ id= [1 2 3 4 5 6 7 14 13 12 11 10 9 8]; % EFIT/Sab PF coil convert to sab10.mat
+ cc_terminal = [gdata.brsp(id)*1.e-6; zeros(ncadd,1)] ; % 1st 14 brsp => F-coils [MA]
+ cc= cc_terminal.*gfile_data.turnfc; % convert to MA => MA-turns
+
+ % VV has 60elements broken into 12 cur groups. We expand 12 actual currents to 60 elements
+ proj= proj_turn(vvid0,vvfrac0);
+ vc6_terminal= gdata.brsp(nfcoil + (1:12)); % 12 currents after nfcoil are VV groups [A]
+ vc0= proj'*vc6_terminal*1e-6; % expand 12 => 60 and change to MA-turn
+ vc= [vc0; zeros(24,1)]; % add in zeros for 10PP and 12IL
+ vvid= (1:length(vvid))'; % since we have expaned up to all elements keep all
+ gfile_data.gdef.vvid = 'VV id; We expand 12 groups in 60 elemnts +24 extra so 1:84';
+ vvfrac= ones(size(vvfrac)); % since we have expaned up fractons are one
+ gfile_data.gdef.vvid = 'VV id; We expand 12 groups in 60 elemnts +24 extra so 1:84';
+
+ if any( abs(gdata.brsp((nfcoil+12+1):end)) > 1e-5*max(abs(gdata.brsp)) )
+ wait(' % CAUTION: some Sabbagh EFIT brsp(15:74) have substantial current???')
+ end
+
+ gfile_data.cc= cc;
+ gfile_data.cc2= cc;
+ gfile_data.gdef.cc= 'PF coil currents in MA-turns; sab09; 14F & 4IC=0 ';
+ gfile_data.gdef.cc2= 'PF coil currents in MA-turns ';
+ gfile_data.vc= vc; % example of vv is vc(1)= brsp(15)*1e-6*[vvturn(1) == 0.25]
+ gfile_data.gdef.vc = 'vessel currents in MA; sab10; 12 groups => 60elm;';
+
+ vvid= (1:length(vvid))'; % since we have expaned up to all elements keep all
+ gfile_data.gdef.vvid = 'VV id; We expand 12 groups in 60 elemnts +24 extra so 1:84';
+
+ end
+ gfile_data.vvfrac= vvfrac;
+ gfile_data.vvid= vvid; % NOTE: we are expanding 12=>56+24 so keep all VV elements
+ gfile_data.gdef.vvid = 'VV identifier; Note: sab10 has 12 VV groups but we maintain 60elm';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation; sab10 WILL CHANGE LATER ';
+
+ case 4 % KI You efit 2010
+
+ case 5 % 2011 rtEFIT
+% This EFIT seems to be different than SAB2010 in that:
+% 1) PF ordering is like DIII-D so we do not have to reorder to get to TokSys (DIII-D) ordering
+% 2) Inside Coils (IC) are included with the coils (14PF + 4IC = 18 currents)
+% 3) Only the passive plates (10 elements) are added VV; BUT NOT Center Limiter (14 elements)
+% Below"if 0" allows VC ouput of EFIT (reduced) size elements or EFUNT (full size) elements
+
+ ncadd = 0;
+ id = [1:18]; % Fcoils
+ % PF coils
+ gfile_data.fcturn = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]';
+ gfile_data.turnfc = [180 144 72 108 208 128 72 180 144 72 108 208 128 72 6 4 6 4];
+ gfile_data.fcid = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18];
+ gfile_data.ecid = [];
+ gfile_data.fcturn= ones(1,length(gfile_data.fcid));
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp(id)*1e-6;zeros(ncadd,1)];
+ cc2= cc;
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns ';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+
+ % Vessel = 60 VV -> 12 sections + 10 PP -> 2 sections
+ vvid= [...
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ... % Outer VV
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ... % Inner VV
+ 13 13 13 14 14 14 13 13 14 14]'; % Passive Plates
+ vvfrac0= [...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.2500 0.2500 0.1250 0.1250 0.1250 0.1250 0.2500 ...
+ 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 0.1250 ...
+ 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 ...
+ 0.2500 0.2500]'; % this is Outer VV
+ vvfrac0= [vvfrac0; vvfrac0]; % add in Inner VV which is same
+ vvfrac= [vvfrac0; 0.2*ones(10,1)]; % add in 10 PP elements
+ vvfrac = [vvfrac;
+ [0.250, 0.250, 0.250, 0.250, 0.200, 0.200, 0.200, 0.200, 0.200, 0.200,...
+ 0.200, 0.200, 0.200, 0.200, 0.050, 0.050, 0.050, 0.050, 0.050, 0.050,...
+ 0.050, 0.050, 0.050, 0.050, 0.050, 0.050, 0.050, 0.050, 0.050, 0.050,...
+ 0.050, 0.050, 0.050, 0.050, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034,...
+ 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034,...
+ 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034, 0.034,...
+ 0.034, 0.034, 0.034, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062,...
+ 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.062, 0.040,...
+ 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040,...
+ 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040, 0.040,...
+ 0.040, 0.040, 0.040, 0.040]'];
+
+ vvid(71:174) = 0; % These elements do not exist in 2011 rtefit_v7b
+
+ vc14_terminal= gdata.brsp(nfcoil + (1:14)); % 14 currents after nfcoil are VV groups [A]
+ vc = vc14_terminal*1e-6; % Keep EFIT size (14) vc=> 1:14
+
+ gfile_data.vc= vc; % example of vv is vc(1)= brsp(15)*1e-6*[vvturn(1) == 0.25]
+ gfile_data.vvid = vvid; % NOTE: we are expanding 14=>70 so keep all VV elements
+ gfile_data.vvfrac = vvfrac; % Expanded out, so all ones
+
+ gfile_data.gdef.vc = 'vessel currents in MA; 14 groups => 70elm;';
+ gfile_data.gdef.vvid = 'VV identifier';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation';
+
+ case 6 % rtEFIT v7d2 (2012 intermediate version)
+ % Similar to previous rtEFIT, except includes inner limiter, cryo, and additional
+ % anti-parallel circuit to model eddy currents in passive plate. This means
+ % 180 vessel elements instead of 70.
+ ncadd = 0;
+ id = [1:18]; % Fcoils
+ % PF coils
+ gfile_data.fcturn = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]';
+ gfile_data.turnfc = [180 144 72 108 208 128 72 180 144 72 108 208 128 72 6 4 6 4];
+ gfile_data.fcid = [1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18];
+ gfile_data.ecid = [];
+ gfile_data.fcturn= ones(1,length(gfile_data.fcid));
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp(id)*1e-6;zeros(ncadd,1)];
+ cc2= cc;
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns ';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+
+ % Vessel = 180 VV -> 12 vessel sections + 2 inner lim + 2 diverter + 2 symmetrical
+ % passive plates + 4 cryo sections + 2 antiparallel pp eddy circuits = 24 sections
+ vvid= [...
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 ...
+ 8 9 9 10 10 11 11 11 12 12 1 1 2 2 2 3 3 4 4 5 ...
+ 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ...
+ 13 13 13 14 14 14 15 15 16 16 15 15 16 16 15 15 16 16 17 17 ...
+ 17 18 18 18 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 ...
+ 19 19 19 19 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 ...
+ 20 20 20 20 20 20 20 20 20 20 20 20 20 21 21 21 21 21 21 21 ...
+ 21 21 21 21 21 21 21 21 21 22 22 22 22 22 22 22 22 22 22 22 ...
+ 22 22 22 22 22 22 22 22 22 22 22 22 22 22 23 23 23 24 24 24]'; %rtefit_v7e
+ vvfrac= [...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 0.2500 0.2500 0.1250 ...
+ 0.1250 0.1250 0.1250 0.2500 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 ...
+ 0.1250 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 0.2500 0.2500 0.1250 ...
+ 0.1250 0.1250 0.1250 0.2500 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 ...
+ 0.1250 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.1667 0.1667 0.1667 0.1667 ...
+ 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.3333 0.3333 ...
+ 0.3333 0.3333 0.3333 0.3333 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 ...
+ 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 ...
+ 0.0500 0.0500 0.0500 0.0500 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 ...
+ 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.500 0.500 -1.000 -1.000 0.500 0.500]';
+ idxvv_efit_to_tok = [1:length(vvfrac)-6]; % toksys model does not have anti-series PP
+
+ vc = gdata.brsp(nfcoil + (1:24))*1e-6; % 24 currents after nfcoil are VV groups [MA-turn]
+
+ gfile_data.vc= vc; % example of vv is vc(1)= brsp(15)*1e-6*[vvturn(1) == 0.25]
+ gfile_data.vvid = vvid;
+ gfile_data.vvfrac = vvfrac;
+ gfile_data.idxvv_efit_to_tok = idxvv_efit_to_tok;
+
+ gfile_data.gdef.vc = 'vessel currents in MA; 24 groups => 180 elm;';
+ gfile_data.gdef.vvid = 'VV identifier';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation';
+
+ case 7 % rtEFIT v7e ()
+ % Similar to rtEFIT_v7d, but with series PF coils and slightly shifted divertor indices
+ ncadd = 0;
+ id = [1:13]; % Fcoils
+ % PF coils
+ gfile_data.fcid = [ 1 2 3 4 5 6 7 1 2 8 9 10 11 7 12 13 12 13];
+ gfile_data.fcturn = [.5 .5 1 1 1 1 .5 .5 .5 1 1 1 1 .5 .5 .5 -.5 .5]';
+ gfile_data.turnfc = [360 288 72 108 208 128 144 72 108 208 128 12 8];
+
+ gfile_data.ecid = [];
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp(id)*1e-6;zeros(ncadd,1)];
+ cc2= cc;
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns ';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+
+ % Vessel = 180 VV -> 12 vessel sections + 2 inner lim + 2 diverter + 2 symmetrical
+ % passive plates + 4 cryo sections + 2 antiparallel pp eddy circuits = 24 sections
+ vvid= [...
+ 1 1 2 2 2 3 3 4 4 5 5 5 5 6 6 7 7 8 8 8 ...
+ 8 9 9 10 10 11 11 11 12 12 1 1 2 2 2 3 3 4 4 5 ...
+ 5 5 5 6 6 7 7 8 8 8 8 9 9 10 10 11 11 11 12 12 ...
+ 13 13 13 14 14 14 15 15 16 16 15 15 16 16 15 15 16 16 17 17 ...
+ 17 18 18 18 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 19 ...
+ 19 19 19 19 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 20 ...
+ 20 20 20 20 20 20 20 20 20 20 20 20 20 21 21 21 21 21 21 21 ...
+ 21 21 21 21 21 21 21 21 21 22 22 22 22 22 22 22 22 22 22 22 ...
+ 22 22 22 22 22 22 22 22 22 22 22 22 22 22 23 23 23 24 24 24]'; %rtefit_v7e
+ vvfrac= [...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 0.2500 0.2500 0.1250 ...
+ 0.1250 0.1250 0.1250 0.2500 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 ...
+ 0.1250 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 0.2500 0.2500 0.1250 ...
+ 0.1250 0.1250 0.1250 0.2500 0.2500 0.2500 0.2500 0.1250 0.1250 0.1250 ...
+ 0.1250 0.2500 0.2500 0.2500 0.2500 0.1667 0.1667 0.1667 0.2500 0.2500 ...
+ 0.3333 0.3333 0.3333 0.3333 0.3333 0.3333 0.1667 0.1667 0.1667 0.1667 ...
+ 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.3333 0.3333 ...
+ 0.3333 0.3333 0.3333 0.3333 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 ...
+ 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 0.0500 ...
+ 0.0500 0.0500 0.0500 0.0500 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 0.0345 ...
+ 0.0345 0.0345 0.0345 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 ...
+ 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 0.0400 ...
+ 0.0400 0.0400 0.0400 0.0400 0.500 0.500 -1.000 -1.000 0.500 0.500]';
+
+ idxvv_efit_to_tok = [1:length(vvfrac)-6]; % toksys model does not have anti-series PP
+
+ vc = gdata.brsp(nfcoil + (1:24))*1e-6; % 24 currents after nfcoil are VV groups [MA-turn]
+
+ gfile_data.vc= vc; % example of vv is vc(1)= brsp(15)*1e-6*[vvturn(1) == 0.25]
+ gfile_data.vvid = vvid;
+ gfile_data.vvfrac = vvfrac;
+ gfile_data.idxvv_efit_to_tok = idxvv_efit_to_tok;
+
+ gfile_data.gdef.vc = 'vessel currents in MA; 24 groups => 180 elm;';
+ gfile_data.gdef.vvid = 'VV identifier';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation';
+
+ otherwise
+ wait(['ERROR read_gfile_tok: unsupported KSTAR efit format ' gfile_data.ecase])
+ return;
+ end
+
+ gfile_data.ecid = [];
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(~isfield(gdata,'gtime') & isfield(gdata,'time'))
+ gdata.gtime = gdata.time;
+ end
+ gfile_data.time = gdata.gtime;
+ gfile_data.tms = gdata.gtime*1e+3;
+ end
+
+% ===========================================================================
+ case {'NSTXU'}
+% ===========================================================================
+
+ if(~exist('nfcoil','var'))
+ nfcoil = [];
+ end
+ if(~exist('nesum','var'))
+ nesum = [];
+ end
+ if(~exist('nves','var'))
+ nves = [];
+ end
+
+ icase = 'rtefit_v7a';
+
+ switch icase
+
+ case 'lrdfit' %used for initial model development
+
+ if isempty(nfcoil), nfcoil=14; end % default NSTXU values
+ if isempty(nves), nves=12; end % default NSTXU values
+ if isempty(nesum), nesum=1; end % default NSTXU values
+ nec=6; %number of ohmic coil elements
+
+ % Believe turnfc should be tok_data_struct.fcnturn
+ if(~isfield(gfile_data,'turnfc'))
+ %This turnfc consistent with Menard Nucl. Fusion 52 (2012) 083015
+ gfile_data.turnfc=[64 32 20 14 14 7 8 7 8 4 5 8 4 5 ...
+ 8 12 12 12 12 7 8 7 8 14 14 20 32 64];
+ %This turnfc tweaked so build_tokamak_system doesn't complain
+ %about PF3U,L coil currents being different (line 562)
+ %Change is to increase turns to 32 instead of 30
+ %Eg turnfc=8 8 8 8 instaed of 7 8 7 8
+ % (there must be a better way to do this)
+ %gfile_data.turnfc=[64 32 20 14 14 8 8 8 8 4 5 8 4 5 ...
+ % 8 12 12 12 12 8 8 8 8 14 14 20 32 64];
+ end
+
+ % Group 28 VF coils into 14 circuits
+ if(~isfield(gfile_data,'fcid'))
+ gfile_data.fcid=[1 2 3 4 4 5 5 5 5 6 6 6 7 7 7 8 8 9 9 10 10 10 10 11 11 12 13 14];
+ end
+ % Specify fraction of current in each coil for the specified fcid
+ % Hence, if you do k = find(fcid==2) % where 2 is an example
+ % then sum(fcturn(k)) % should be 1
+ if(~isfield(gfile_data,'fcturn'))
+ for k=1:nfcoil
+ idfc = gfile_data.fcid(k);
+ tmp = sum(gfile_data.turnfc(gfile_data.fcid==idfc)); % get total turns according to fcid
+ gfile_data.fcturn(k)=gfile_data.turnfc(k)/tmp; % apply fcturn definition
+ end
+ end
+
+ % vv
+ gfile_data.vvid = [1:nves];
+ gfile_data.vvfrac = ones(nves,1);
+
+ % ec
+ % To make the object have a consistent interpretation, this needs to be the
+ % fraction of current carried by each defined group????
+ % gfile_data.ecturn = gfile_data.ecturn/sum(gfile_data.ecturn);
+ if(~isfield(gfile_data,'ecid'))
+ gfile_data.ecid=ones(59,1);
+ end
+ if(~isfield(gfile_data,'ecturn'))
+ gfile_data.ecturn=repmat(15,[59 1]);
+ end
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+
+ % Changing cc and cc2 so they are amp-turns
+ % brsp is terminal current (A) NOT Amp-turns because NSTX uses groups and turns
+ % for each group.
+
+ % for LRDFIT equilibria, brsp and ecurrt are NaN.
+
+ if length(gdata.brsp)==1 & isnan(gdata.brsp)
+ fprintf('WARNING std_efit_units: brsp is NaN, coil currents cannot be processed\n')
+ gfile_data.vc = [];
+ gfile_data.cc = [];
+ gfile_data.cc2 = [];
+ elseif length(gdata.ecurrt)==1 & isnan(gdata.ecurrt)
+ fprintf('WARNING std_efit_units: ecurrt is NaN, coil currents cannot be processed\n')
+ gfile_data.vc = [];
+ gfile_data.cc = [];
+ gfile_data.cc2 = [];
+ elseif (~isempty(gdata.brsp) & ~isempty(gdata.ecurrt) & length(gdata.brsp)>=(nfcoil+nves))
+ gfile_data.vc = gdata.brsp(nfcoil+1:nfcoil+nves)*1.e-6;
+ cc_terminal = [gdata.ecurrt; gdata.brsp(1:nfcoil)]*1.e-6;
+
+ turns = zeros(2,1); %force to be column
+ for k=1:nesum
+ turns(k) = sum(gfile_data.ecturn);
+ end
+ for k=1:nfcoil
+ idx = find(gfile_data.fcid==k);
+ turns(k+nesum) = sum(gfile_data.turnfc(idx));
+ end
+
+ gfile_data.cc = cc_terminal.*turns;
+ gfile_data.cc2 = gfile_data.cc;
+
+ gfile_data.gdef.cc = 'Fcoil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'Fcoil currents in MA-turns ';
+ gfile_data.gdef.vc = 'vessel currents in MA';
+ end
+
+
+ case 'rtefit_v7a' %first rtefit for NSTXU based on Sabbagh mhdin.dat dated 1/15/15
+
+ if isempty(nfcoil), nfcoil=14; end % default NSTXU values
+ if isempty(nves), nves=40; end % default NSTXU values
+ if isempty(nesum), nesum=1; end % default NSTXU values
+ nec=8; %number of ohmic coil elements
+
+ % 22 sub coils now so can't use with LRDFIT models
+ % Order is from rtefit_NSTX_diagnostics.h
+ % "PF1AU","PF1BU","PF1CU", "PF2U","PF3U","PF4U","PF5U","PF5L","PF4L","PF3L","PF2L","PF1CL","PF1BL","PF1AL",
+ gfile_data.fcid = [1.00000000E+00, 2.00000000E+00, 3.00000000E+00, 4.00000000E+00,...
+ 4.00000000E+00, 5.00000000E+00, 5.00000000E+00, 6.00000000E+00,...
+ 6.00000000E+00, 7.00000000E+00, 7.00000000E+00, 8.00000000E+00,...
+ 8.00000000E+00, 9.00000000E+00, 9.00000000E+00, 1.00000000E+01,...
+ 1.00000000E+01, 1.10000000E+01, 1.10000000E+01, 1.20000000E+01,...
+ 1.30000000E+01, 1.40000000E+01];
+ gfile_data.turnfc = [6.40000000E+01, 3.20000000E+01, 2.00000000E+01, 1.40000000E+01,...
+ 1.40000000E+01, 1.50000000E+01, 1.50000000E+01, 9.00000000E+00,...
+ 8.00000000E+00, 1.20000000E+01, 1.20000000E+01, 1.20000000E+01,...
+ 1.20000000E+01, 9.00000000E+00, 8.00000000E+00, 1.50000000E+01,...
+ 1.50000000E+01, 1.40000000E+01, 1.40000000E+01, 2.00000000E+01,...
+ 3.20000000E+01, 6.40000000E+01];
+ gfile_data.fcturn = [1, 1, 1,...
+ 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,...
+ 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,...
+ 1, 1, 1];
+ gfile_data.vvid = [1, 1, 2, 2, 3, 4, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7,...
+ 8, 8, 8, 8, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 9, 9, 9, 9, 9, 9,...
+ 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10,...
+ 11, 11, 11, 11, 11, 12, 12, 13, 13, 13, 14, 15, 16, 17, 18, 18, 18, 19, 19, 20,...
+ 20, 20, 20, 20, 21, 21, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 22,...
+ 22, 22, 22, 22, 22, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23,...
+ 24, 24, 24, 24, 25, 25, 25, 25, 25, 23, 23, 23, 23, 23, 23, 26, 26, 26, 27, 28,...
+ 29, 29, 30, 30, 31, 31, 32, 32, 33, 34, 35, 36, 37, 38, 39, 40];
+ gfile_data.vvfrac = [0.5000, 0.5000, 0.5000, 0.5000, 1.0000, 1.0000, 0.3333, 0.3333, 0.3333, 0.0625,...
+ 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.2000, 0.2000, 0.2000, 0.2000, 0.2000,...
+ 0.2500, 0.2500, 0.2500, 0.2500, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625,...
+ 0.0625, 0.0625, 0.0625, 0.0625, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417,...
+ 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417,...
+ 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.5000, 0.5000,...
+ 0.2000, 0.2000, 0.2000, 0.2000, 0.2000, 0.5000, 0.5000, 0.3333, 0.3333, 0.3333,...
+ 1.0000, 1.0000, 1.0000, 1.0000, 0.3333, 0.3333, 0.3333, 0.5000, 0.5000, 0.2000,...
+ 0.2000, 0.2000, 0.2000, 0.2000, 0.5000, 0.5000, 0.0417, 0.0417, 0.0417, 0.0417,...
+ 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417,...
+ 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417, 0.0417,...
+ 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625,...
+ 0.2500, 0.2500, 0.2500, 0.2500, 0.2000, 0.2000, 0.2000, 0.2000, 0.2000, 0.0625,...
+ 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.3333, 0.3333, 0.3333, 1.0000, 1.0000,...
+ 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 0.5000, 1.0000, 1.0000,...
+ 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000];
+ gfile_data.ecid = [1, 1, 1, 1, 1, 1, 1, 1];
+ gfile_data.ecturn = [112.0, 110.0, 109.5, 108.5, 108.5, 109.5, 110.0, 112.0];
+
+ turns = zeros(2,1); %force to be column
+ for k=1:nesum
+ turns(k) = sum(gfile_data.ecturn);
+ end
+ for k=1:nfcoil
+ idx = find(gfile_data.fcid==k);
+ turns(k+nesum) = sum(gfile_data.turnfc(idx));
+ end
+
+ %brsp has also the vessel currents
+ cc_terminal = [gdata.ecurrt; gdata.brsp(1:nfcoil)]*1.e-6;
+ gfile_data.cc = cc_terminal.*turns;
+ gfile_data.cc2 = gfile_data.cc;
+
+
+ gfile_data.gdef.fcid = 'Fcoil circuit id';
+ gfile_data.gdef.fcturn = 'Multiplier used in generation of EFIT greens tables; if =1 brsp is in A-turn; fraction of total turns';
+ gfile_data.gdef.turnfc = 'Multiplier of EFIT input currents, if =1 brsp is in terminal [A]';
+ gfile_data.gdef.vvid = 'VV identifier; Note: sab09 has 6 VV groups but we maintain 56elm';
+ gfile_data.gdef.vvfrac = 'VV fraction used in Green generation; sab09 WILL CHANGE LATER ';
+ gfile_data.gdef.ecid = 'Ohmic coil circuit id';
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+ gfile_data.gdef.cc = 'Fcoil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'Fcoil currents in MA-turns ';
+ gfile_data.gdef.vc = 'vessel currents in MA';
+
+
+ end %icase switch
+
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(~isfield(gdata,'gtime') & isfield(gdata,'time'))
+ gdata.gtime = gdata.time;
+ end
+ gfile_data.time = gdata.gtime;
+ gfile_data.tms = gdata.gtime*1e+3;
+ end
+
+% ===========================================================================
+ case {'NSTX'}
+% ===========================================================================
+
+ if(~exist('nfcoil','var'))
+ nfcoil = [];
+ end
+ if(~exist('nesum','var'))
+ nesum = [];
+ end
+ if(~exist('nves','var'))
+ nves = [];
+ end
+
+ if old_efit==0 % NEW Apr 05 EFIT Config: 04202005Av1.0
+
+ if isempty(nfcoil), nfcoil=17; end % default NSTX values
+ if isempty(nves), nves=35; end % default NSTX values
+ if isempty(nesum), nesum=1; end % default NSTX values
+
+% Coil Note except for changes to pf1aul,
+% below I add a 2nd vector which adds from the Old to the New(Apr05)
+ gfile_data.fcid= [[1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 11] ...
+ [12 13 14 15 16 17]];
+ gfile_data.turnfc=[[20 14 14 15 15 9 8 12 12 12 12 9 8 15 15 14 14 20 32]...
+ [1 1 1 1 48 48]];
+
+ gfile_data.fcturn = ...
+ [[1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1 1]...
+ [1 1 1 1 1 1]];
+
+% vv has some differences internal so new from "make"
+ gfile_data.vvid = [ ...
+ 1 2 3 4 4 5 6 6 7 7 7 8 8 8 8 8 8 8 8 9 9 9 9 9 10 10 11 11 11 12 ...
+ 13 14 15 16 16 16 17 17 18 18 18 18 18 19 19 19 19 19 20 20 21 ...
+ 22 22 23 24 25 26 26 27 27 28 29 30 31 32 33 34 35];
+
+ gfile_data.vvfrac = [ ...
+ 1 1 1 0.5 0.5 1 ...
+ 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 ...
+ 0.25 0.25 0.25 0.25 ...
+ 0.33333 0.33333 0.33333 0.33333 0.33333 0.33333 ...
+ 0.5 0.5 0.5 0.5 ...
+ 0.33333 0.33333 0.33333 0.33333 0.33333 0.33333 ...
+ 0.25 0.25 0.25 0.25 ...
+ 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 ...
+ 1 0.5 0.5 1 1 1 0.5 0.5 0.5 0.5 1 1 1 1 1 1 1 1];
+
+ gfile_data.vvfrac = [ ...
+ 1.0000 1.0000 1.0000 0.5000 0.5000 1.0000 0.5000 ...
+ 0.5000 0.3333333 0.3333333 0.3333333 0.1250 0.1250 0.1250 ...
+ 0.1250 0.1250 0.1250 0.1250 0.1250 0.2000 0.2000 ...
+ 0.2000 0.2000 0.2000 0.5000 0.5000 0.3333333 0.3333333 ...
+ 0.3333333 1.0000 1.0000 1.0000 1.0000 0.3333333 0.3333333 ...
+ 0.3333333 0.5000 0.5000 0.2000 0.2000 0.2000 0.2000 ...
+ 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.5000 ...
+ 0.5000 1.0000 0.5000 0.5000 1.0000 1.0000 1.0000 ...
+ 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 ...
+ 1.0000 1.0000 1.0000 1.0000 1.0000];
+
+% ecoil - no change old to new
+ gfile_data.ecid = ones(240,1);
+ gfile_data.ecturn = [ ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ];
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+
+% To make the object have a consistent interpretation, this needs to be the
+% fraction of current carried by each defined group????
+% gfile_data.ecturn = gfile_data.ecturn/sum(gfile_data.ecturn);
+% gfile_data.turnec = ?; % WHAT IS THIS??????
+
+% Changing cc and cc2 so they are amp-turns
+% brsp is terminal current (A) NOT Amp-turns because NSTX uses groups and turns
+% for each group.
+
+% for LRDFIT equilibria, brsp and ecurrt are NaN.
+
+ if length(gdata.brsp)==1 & isnan(gdata.brsp)
+ fprintf('WARNING std_efit_units: brsp is NaN, coil currents cannot be processed\n')
+ gfile_data.vc = [];
+ gfile_data.cc = [];
+ gfile_data.cc2 = [];
+ elseif length(gdata.ecurrt)==1 & isnan(gdata.ecurrt)
+ fprintf('WARNING std_efit_units: ecurrt is NaN, coil currents cannot be processed\n')
+ gfile_data.vc = [];
+ gfile_data.cc = [];
+ gfile_data.cc2 = [];
+ elseif ~isempty(gdata.brsp) & ~isempty(gdata.ecurrt)
+ gfile_data.vc = gdata.brsp(nfcoil+1:nfcoil+nves)*1.e-6;
+ cc_terminal = [gdata.ecurrt; gdata.brsp(1:nfcoil)]*1.e-6;
+
+ turns = zeros(2,1); %force to be column
+ for k=1:nesum
+ turns(k) = sum(gfile_data.ecturn);
+ end
+ for k=1:nfcoil
+ idx = find(gfile_data.fcid==k);
+ turns(k+nesum) = sum(gfile_data.turnfc(idx));
+ end
+
+ gfile_data.cc = cc_terminal.*turns;
+ gfile_data.cc2 = gfile_data.cc;
+
+ gfile_data.gdef.cc = 'Fcoil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'Fcoil currents in MA-turns ';
+ gfile_data.gdef.vc = 'vessel currents in MA';
+ end
+
+% ----------------------------------------------------------------------------
+ else % if old_efit==1 OLD Feb 2002 EFIT BELOW Config_name= '02072002Av1.0';
+% ----------------------------------------------------------------------------
+
+ if isempty(nfcoil), nfcoil=11; end % default NSTX values
+ if isempty(nves), nves=30; end % default NSTX values
+ if isempty(nesum), nesum=1; end % default NSTX values
+
+ gfile_data.fcid = [1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 11];
+ gfile_data.turnfc = [48 14 14 15 15 9 8 12 12 12 12 9 8 15 15 14 14 48 32];
+ gfile_data.fcturn = ...
+ [1 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 1 1];
+
+ gfile_data.vvid = [ ...
+ 1 2 3 4 4 5 6 6 6 6 6 6 6 7 7 7 7 8 8 8 9 9 9 10 10 11 11 12 12 12 ...
+ 13 13 13 14 14 14 14 15 15 15 15 15 15 15 16 17 17 18 19 20 21 21 ...
+ 22 22 23 24 25 26 27 28 29 30];
+ gfile_data.vvfrac = [ ...
+ 1 1 1 0.5 0.5 1 ...
+ 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 ...
+ 0.25 0.25 0.25 0.25 ...
+ 0.33333 0.33333 0.33333 0.33333 0.33333 0.33333 ...
+ 0.5 0.5 0.5 0.5 ...
+ 0.33333 0.33333 0.33333 0.33333 0.33333 0.33333 ...
+ 0.25 0.25 0.25 0.25 ...
+ 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 0.14286 ...
+ 1 0.5 0.5 1 1 1 0.5 0.5 0.5 0.5 1 1 1 1 1 1 1 1];
+
+ gfile_data.ecid = ones(240,1);
+ gfile_data.ecturn = [ ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 4.0875 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 4.033333 4.033333 4.033333 4.033333 4.033333 4.033333 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.995833 3.995833 3.995833 3.995833 3.995833 3.995833 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ...
+ 3.933333 3.933333 3.933333 3.933333 3.933333 3.933333 ];
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+
+% To make the object have a consistent interpretation, this needs to be the
+% fraction of current carried by each defined group????
+% gfile_data.ecturn = gfile_data.ecturn/sum(gfile_data.ecturn);
+% gfile_data.turnec = ?; % WHAT IS THIS??????
+
+% Changing cc and cc2 so they are amp-turns
+% brsp is terminal current (A) NOT Amp-turns because NSTX uses groups and turns
+% for each group.
+
+ if ~isempty(gdata.brsp) & ~isempty(gdata.ecurrt)
+ gfile_data.vc = gdata.brsp(nfcoil+1:nfcoil+nves)*1.e-6;
+ cc_terminal = [gdata.ecurrt; gdata.brsp(1:nfcoil)]*1.e-6;
+
+ turns = zeros(2,1); %force to be column
+ for k=1:nesum
+ turns(k) = sum(gfile_data.ecturn);
+ end
+ for k=1:nfcoil
+ idx = find(gfile_data.fcid==k);
+ turns(k+nesum) = sum(gfile_data.turnfc(idx));
+ end
+
+ gfile_data.cc = cc_terminal.*turns;
+ gfile_data.cc2 = gfile_data.cc;
+
+ gfile_data.gdef.cc = 'Fcoil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'Fcoil currents in MA-turns ';
+ gfile_data.gdef.vc = 'vessel currents in MA';
+ end
+% --------------------------------------------------
+ end % if old_efit
+% --------------------------------------------------
+
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(isfield(gfile_data,'gtime'))
+ gfile_data.time = gdata.gtime;
+ gfile_data.tms = gdata.gtime*1e+3;
+ else
+ gfile_data.time = gdata.time;
+ gfile_data.tms = gdata.time*1e+3;
+ end
+ end
+
+% ===========================================================================
+ case {'ITER'}
+% ===========================================================================
+
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp*1e-6];
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+ gfile_data.fcid = 1:12;
+
+% use date of file to switch between turns in system: iter07 Vs 2010, See /u/leuer/efit/iter/
+ date2010= datenum('01-Jan-2010'); % Runs prior to 2010 use old turns
+ d= dir(filename);
+ if datenum(d.date) < datenum('01-Jan-2010') % Runs prior to 2010 use old turns
+ disp('% NOTE: read_gfile using OLD 2007 ITER TURNS: 249 106 185 169 217 425 548 548..')
+ gfile_data.turnfc = [249 106 185 169 217 425 548 548 548 548 548 548];
+ else
+ disp('% NOTE: read_gfile using 2010 ITER TURNS: 249 115 186 170 217 459 553 553..')
+ gfile_data.turnfc = [249 115 186 170 217 459 553 553 553 553 553 553];
+ end
+ gfile_data.fcturn = ones(12,1);
+
+ gfile_data.ecid = [];
+
+% disp('Using Build for ITER07')
+
+% ===========================================================================
+ case {'CTF'}
+% ===========================================================================
+
+%cc-vector with NO-ECOIL (MA-t):
+ if ~isempty(gdata.brsp)
+ cc = gdata.brsp*1.e-6;
+
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc;
+
+ gfile_data.gdef.cc = '14 F coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = '14 F coil currents in MA-turns ';
+ end
+ gfile_data.fcturn = ones(1,nfcoil); % used to build efit greens tables
+ gfile_data.turnfc = 50*ones(1,nfcoil); % Ipf*turnfc => efit greens table
+ gfile_data.fcid = 1:14;
+ gfile_data.ecid = [];
+% disp(['read_gfile_tok done reading CTF gfile: ',filename])
+
+% ===========================================================================
+ case {'FDF'}
+% ===========================================================================
+
+%cc-vector with NO-ECOIL (MA-t):
+ if ~isempty(gdata.brsp)
+
+% special reordering of brsp for symmetric EFIT runs:
+ if size(gdata.brsp,1)==22
+ sm= 0;
+ mx= 0;
+ for ii= 1:11
+ sm= (gdata.brsp(ii)-gdata.brsp(ii+11)).^2 + sm;
+ mx= max([mx,abs(gdata.brsp(ii)),abs(gdata.brsp(ii+1))]);
+ end
+ err= sqrt(sm/11)/mx; % diff is typically of order 6e-7 for symmetric
+ if err <= 1e-5
+ gfile_data.brsp= gfile_data.brsp([1:11 22:-1:12]); % reorder for sym efit
+ end
+ end
+
+ cc = gdata.brsp*1.e-6;
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc;
+
+ gfile_data.gdef.cc = 'F coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'F coil currents in MA-turns ';
+ end
+ gfile_data.fcturn = ones(1,nfcoil); % used to build efit greens tables
+ %Apparently turnfc must be ones, even though MHDIN.DAT claims = 50;
+ % ==> Logic in rzrig or here may be broken... DAH 5/31/08
+% Below depends if MATLAB is build with fcnturn=1 or 50 JAL 6may2011
+ gfile_data.turnfc = turnfc0*ones(1,nfcoil); % Ipf*turnfc => efit greens table
+ gfile_data.fcid = 1:nfcoil;
+ gfile_data.ecid = [];
+% disp(['read_gfile_tok done reading FDF gfile: ',filename])
+
+
+% ===========================================================================
+ case {'CFETR'}
+% ===========================================================================
+
+ if ~isempty(gdata.brsp)
+ cc = [gdata.brsp*1e-6];
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc;
+ gfile_data.gdef.cc = 'PF coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+ nfc=14;
+ gfile_data.fcid = 1:nfc;
+ gfile_data.turnfc = ones(nfc,1);
+ gfile_data.fcturn = ones(nfc,1);
+ gfile_data.ecid = [];
+
+ disp('cfetr')
+
+% ===========================================================================
+ case {'PEGASUS'}
+% ===========================================================================
+
+ if ~isempty(gdata.brsp)
+% ncadd=4; % fix cc so it has extra control coils
+% cc = [gdata.brsp*1e-6;zeros(ncadd,1)];
+% cc2= cc;
+
+% gfile_data.cc = cc;
+% gfile_data.cc2 = cc2;
+% gfile_data.gdef.cc = 'PF coil currents in MA-turns ';
+% gfile_data.gdef.cc2 = 'PF coil currents in MA-turns ';
+ end
+
+% ??? Do we care that length(cc)~=length(turnfc) ???
+
+% gfile_data.fcturn = [1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]';
+ gfile_data.turnfc = [2 3 5 5 5 5 3 2 5 5 280 4 4 1 1 1 1];
+ gfile_data.fcid = [1 2 3 4 4 3 2 5 6 7 8 9 10 11 12 13 14];
+ gfile_data.ecid = [];
+
+ gfile_data.fcturn = [1 0.5 0.5 0.5 0.5 0.5 0.5 1 1 1 1 1 1 1 1 1 1]';
+
+% ===========================================================================
+ case {'SST'}
+% ===========================================================================
+
+ gfile_data.fcturn = ones(1,14);
+ gfile_data.turnfc = [40 40 192 40 40 8 1 40 40 192 40 40 8 1];
+ gfile_data.fcid = [1 2 3 4 5 6 7 8 9 10 11 12 13 14];
+ gfile_data.ecid = [ ... % must be same as ecdata(5,:)
+ 1 1 1 1 1 ];
+ gfile_data.ecturn = [692 40 8 40 8]; % required if ecid is not empty
+ gfile_data.gdef.ecturn = 'individual turns in ecoil(s)';
+
+% Construct E-coil Number of turns to make output MA-turns:
+ dum= min(2,nesum); % if nesum=1 then use ONLY 1 ECOIL current
+
+ if ~isempty(gdata.brsp) & ~isempty(gdata.ecurrt)
+% Construct E-coil Number of turns to make output MA-turns:
+ dum= min(2,nesum); % if nesum=1 then use ONLY 1 ECOIL current
+% find reduced ecoil set data
+ nee = length(unique(gfile_data.ecid));
+ clear neturn
+ for ii=1:nee
+ idx = find(gfile_data.ecid==ii);
+ neturn(ii,1) = sum(gfile_data.ecturn(idx));
+ end
+
+ %cc-vector with 2 e-coil segments (MA-t):
+ cc2 = [gdata.ecurrt(1:dum).*neturn(1:dum);gdata.brsp]*1.e-6;
+
+% CAUTION: we are not using the 6segment e-coil anymore so this is crippled
+%cc-vector with all e-coil segs(2 or 6) (MA-t):
+ cc= cc2;
+
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+
+ gfile_data.gdef.cc = 'E/F coil currents in MA-turns; convert to toksys cc0 using cc_efit_to_tok';
+ gfile_data.gdef.cc2 = 'E/F coil currents in MA-turns ';
+ end
+% ===========================================================================
+ case {'HL2M'} % Jim Leuer 2jul2015
+% ===========================================================================
+ if(~exist('nesum','var'))
+ nesum = 1;
+ end
+
+% include d3d build parameters (see build area => must be same)
+ gfile_data.fcturn = [28 28 28 28 28 27 28 28 28 28 28 28 28 27 28 28];
+ gfile_data.turnfc = ones(1,16);
+ gfile_data.fcid = 1:16;
+ gfile_data.ecid = [ 1 1];
+ gfile_data.ecturn = [48 48]; % ? NOT SURE IF THIS IS 1 or 48?
+ gfile_data.gdef.ecturn = 'individual turns in upper and lower ecoils';
+
+
+ if ~isempty(gdata.brsp) & ~isempty(gdata.ecurrt)
+% Construct E-coil Number of turns to make output MA-turns:
+ dum= 1; % ONLY 1 ECOIL current
+% find reducted ecoil set data
+ nee = length(unique(gfile_data.ecid));
+ clear neturn
+ for ii=1:nee
+ idx = find(gfile_data.ecid==ii);
+ neturn(ii,1) = sum(gfile_data.ecturn(idx));
+ end
+
+ % cc-vector with 1 e-coil segment is in Amps, brsp is in Amps so conversion (MA-t):
+ cc2 = [gdata.ecurrt(1:dum).*neturn(1:dum);...
+ gdata.brsp.*gfile_data.fcturn']*1.e-6; % this is now in MA-t
+
+ cc= cc2;
+
+ gfile_data.cc = cc;
+ gfile_data.cc2 = cc2;
+
+ gfile_data.gdef.cc = 'E/F currents in MA-t; convert to toksys cc0 using std_efit_units';
+ gfile_data.gdef.cc2 = 'E/F currents in MA-turns';
+ end
+
+ if(isfield(gdata,'gtime') | isfield(gdata,'time'))
+ if(~isfield(gdata,'gtime') & isfield(gdata,'time'))
+ gdata.gtime = gdata.time;
+ end
+ gfile_data.time=gdata.gtime*1e-3;
+ gfile_data.tms= gdata.gtime;
+ end
+% Validated using test_build jal 7jul2014
+
+% ===========================================================================
+ otherwise
+% ===========================================================================
+
+ disp(['%ERROR: std_efit_units does not recognize tokamak:' tokamak])
+
+end % switch
+
+ gfile_data.gdef.fcturn = 'Multiplier used in generation of EFIT greens tables; if =1 brsp is in A-turn';
+ gfile_data.gdef.turnfc = 'Multiplier of EFIT input currents, if =1 brsp is in terminal [A]';
+ gfile_data.gdef.fcid = 'Coil identifier used in construction of greens tables';
+ gfile_data.gdef.ecid = 'Coil identifier used in construction of greens tables';
+
+% ===========================================================================
+% DERIVED DATA
+% ===========================================================================
+
+% Derive efit grid:
+
+ rg= linspace(gdata.rgrid1, gdata.rgrid1+gdata.xdim, gdata.nw)';
+ zg= linspace(gdata.zmid-gdata.zdim/2,gdata.zmid+gdata.zdim/2,gdata.nh)';
+ dr= (rg(end)-rg(1))/(gdata.nw-1);
+ dz= (zg(end)-zg(1))/(gdata.nh-1);
+
+% Create useful fluxes (even for iecurr not = 2)
+% Convert flux objects to REAL units:
+
+ psizr = -gdata.psirz'*2*pi*sign(gdata.cpasma);
+ psimag = -gdata.ssimag*2*pi*sign(gdata.cpasma);
+ psibry = -gdata.ssibry*2*pi*sign(gdata.cpasma);
+ psibnd = -gdata.ssibry*2*pi*sign(gdata.cpasma);
+
+% Convert pcurrt to nh x nw array and scale to MA/m^2:
+
+ if ~isempty(gdata.pcurrt) & ~isfield(gfile_data,'jphi')
+ jphi = gdata.pcurrt*1.e-6/(dz*dr);
+ gfile_data.jphi = jphi;
+ gfile_data.gdef.jphi = 'current density on grid MA/m^2';
+ end
+
+% Construct output object:
+
+ gfile_data.rg= rg;
+ gfile_data.zg= zg;
+ gfile_data.dr= dr;
+ gfile_data.dz= dz;
+ gfile_data.psizr = psizr;
+ gfile_data.psimag = psimag;
+ gfile_data.psibry = psibry;
+ gfile_data.psibnd = psibnd;
+
+ gfile_data.gdef.rg= 'radial coordinates of plasma grid';
+ gfile_data.gdef.zg= 'vertical coordinates of plasma grid';
+ gfile_data.gdef.dr= 'radial distance between plasma grid points';
+ gfile_data.gdef.dz= 'vertical distance between plasma grid points';
+ gfile_data.gdef.psizr = 'true total flux on grid in Wb';
+ gfile_data.gdef.psimag = 'axis flux in true Wb';
+ gfile_data.gdef.psibry = 'flux on plasma boundary';
+ gfile_data.gdef.psibnd = 'boundary flux in true Wb (psibry also defined same)';
+
diff --git a/matlab/D3D/efit/private/trace_boundary.m b/matlab/D3D/efit/private/trace_boundary.m
new file mode 100644
index 00000000..a2bb4343
--- /dev/null
+++ b/matlab/D3D/efit/private/trace_boundary.m
@@ -0,0 +1,229 @@
+function [rb,zb,rx,zx,ra,za,r0,z0,ilimited,psimag,psibry]=trace_boundary(psizr,rg,zg,limdata)
+%
+% USAGE: [rb,zb,rx,zx,ra,za,r0,z0,ilimited,psimag,psibry]=trace_boundary(psizr,rg,zg,limdata)
+%
+% PURPOSE: Find boundary
+%
+% INPUTS: psizr
+%
+% OUTPUTS: Coordinates of boundary: rb, zb
+% Coordinates of x or touch point: r0,z0
+% Coordinates of axis: ra, za
+%
+% RESTRICTIONS:
+%
+% METHOD:
+%
+
+% VERSION @(#)trace_boundary.m 1.1 01/12/11
+%
+% WRITTEN BY: Anders Welander ON 3/31/10
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % Prelims
+ mu0 = .4e-6*pi;
+ twopi = 2*pi;
+ nr = length(rg); nz = length(zg);
+ dr = rg(2)-rg(1); dz = zg(2)-zg(1);
+ rb=0; zb=0; rx=0; zx=0; ra=0; za=0; r0=0; z0=0; ilimited=0; % Always return values to avoid error
+
+ if exist('limdata','var')
+ if size(limdata,2)<size(limdata,1), limdata=limdata'; end
+ nlim = size(limdata,2);
+ rlm = interp1(1:nlim,limdata(2,:),1:.01:nlim);
+ zlm = interp1(1:nlim,limdata(1,:),1:.01:nlim);
+ psilim = interp2(rg,zg,psizr,rlm,zlm,'spline');
+ end
+
+ % Determine whether the axis is at max or min flux
+ [psimax, imax] = max(psizr); zmax = zg(imax); [psimax, imax] = max(psimax); rmax = rg(imax); zmax=zmax(imax);
+ [psimin, imin] = min(psizr); zmin = zg(imin); [psimin, imin] = min(psimin); rmin = rg(imin); zmin=zmin(imin);
+ % Pick the extremum closest to the middle of the grid.
+ if (rmax-mean(rg))^2+(zmax-mean(zg))^2 < (rmin-mean(rg))^2+(zmin-mean(zg))^2
+ ra = rmax; za = zmax; % First stab at ra, za
+ else
+ ra = rmin; za = zmin;
+ psizr = -psizr; % This way the flux always has max at axis
+ end
+
+ % Find weights in grid points to calculate value in a point using cubic Hermite spline (ref. wikipedia)
+ mx = [0 2 0 0;-1 0 1 0;2 -5 4 -1;-1 3 -3 1]/2;
+ neighbors = reshape([-1-nz -1 -1+nz -1+2*nz;-nz 0 nz 2*nz;1-nz 1 1+nz 1+2*nz;2-nz 2 2+nz 2+2*nz],1,16);
+
+ % Magnetic axis
+ for j=1:50
+ iz = max(find(zg<za)); ir = max(find(rg<ra)); ia = iz+nz*(ir-1);
+ if length(ia)==0
+ disp('trace_boundary: Could not find axis.')
+ return
+ end
+ iia = ia+neighbors; % iia indexes 16 points around magnetic axis
+ if min(iia)<1 | max(iia)>nr*nz
+ disp('trace_boundary: Could not find axis.')
+ return
+ end
+ pp = psizr(iia);
+ tr = (ra-rg(ir))/dr; tz = (za-zg(iz))/dz;
+ wa = reshape(([1 tz tz^2 tz^3]*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ war = reshape(([1 tz tz^2 tz^3]*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ waz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ warr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 0 2 6*tr]/dr^2*mx,1,16);
+ wazz = reshape(([0 0 2 6*tz]/dz^2*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ warz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ ashift_Ba = -inv([pp*warr' pp*warz';pp*warz' pp*wazz']);
+ c = ashift_Ba*[pp*war'; pp*waz']; ra = ra+c(1)/5; za = za+c(2)/5;
+ end
+ psimag = pp*wa';
+
+ % Find x-point
+ na = 2*(nr-1); psix = -1e99; rho = linspace(0,sqrt((zg(end)-zg(1))^2+(rg(end)-rg(1))^2)/2,100);
+ for ia = 1:na
+ r = ra+rho*cos(ia*2*pi/na); z = za+rho*sin(ia*2*pi/na);
+ p = interp2(rg,zg,psizr,r,z,'spline');
+ j = 5;
+ while p(j+1)<p(j) & j<99, j=j+1; end
+ if p(j) > psix
+ psix = p(j);
+ rx = r(j); zx=z(j);
+ end
+ end
+
+ % Now we are close
+ for j=1:30
+ if rx<rg(3) | rx>rg(end-3) | zx<zg(3) | zx>zg(end-3)
+ if exist('limdata','var')
+ ilimited = 1; j=99;
+ else
+ disp('trace_boundary: Could not find x-point.')
+ return
+ end
+ end
+ if ~ilimited
+ iz = max(find(zg<zx)); ir = max(find(rg<rx)); ix = iz+nz*ir;
+ iix = ix+neighbors;
+ pp = psizr(iix);
+ tr = (rx-rg(ir))/dr; tz = (zx-zg(iz))/dz;
+ wx = reshape(([1 tz tz^2 tz^3]*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ wxz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 0 2 6*tr]/dr^2*mx,1,16);
+ wxzz = reshape(([0 0 2 6*tz]/dz^2*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ xshift_Bx = -inv([pp*wxrr' pp*wxrz';pp*wxrz' pp*wxzz']);
+ c = xshift_Bx*[pp*wxr'; pp*wxz']; rx = rx+c(1)/5; zx = zx+c(2)/5;
+ end
+ end
+ if ilimited
+ [psibry, ilim] = max(psilim);
+ r0 = rlm(ilim); z0 = zlm(ilim);
+ else
+ psibry = pp*wx';
+ w0 = wx; r0 = rx; z0 = zx;
+ end
+
+ th0 = angle(r0-ra+sqrt(-1)*(z0-za)); dth = 2*pi/na;
+ % Trace boundary
+ rb([1 na+1]) = r0; zb([1 na+1]) = z0;
+ for ia = 2:na
+ th = th0+dth*(ia-1);
+ r = ra+rho*cos(th); z = za+rho*sin(th);
+ p = interp2(rg,zg,psizr,r,z,'spline');
+ j = 5;
+ while p(j+1)<p(j) & j<99 & p(j)>psibry, j=j+1; end
+ if p(j) > psibry
+ if j<99 % In this case we are close to the other x-point
+ rx2 = r(j); zx2 = z(j);
+ rx2a = r(min(100,j+3)); zx2a = z(min(100,j+3));
+ rx2b = r(max(1,j-3)); zx2b = z(max(1,j-3));
+ for k=1:9
+ iz = max(find(zg<zx2)); ir = max(find(rg<rx2)); ix = iz+nz*ir;
+ if length(ix)==0 | iz<3 | iz>nz-3 | ir<3 | ir>nr-3 % Found no x-point so set point like normal
+ rx2 = spline(p(3:j),r(3:j),psibry);
+ zx2 = spline(p(3:j),z(3:j),psibry);
+ k = 9;
+ else
+ iix = ix+neighbors;
+ pp = psizr(iix);
+ tr = (rx2-rg(ir))/dr; tz = (zx2-zg(iz))/dz;
+ wx = reshape(([1 tz tz^2 tz^3]*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ wxz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 0 2 6*tr]/dr^2*mx,1,16);
+ wxzz = reshape(([0 0 2 6*tz]/dz^2*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ xshift_Bx = -inv([pp*wxrr' pp*wxrz';pp*wxrz' pp*wxzz']);
+ c = xshift_Bx*[pp*wxr'; pp*wxz']; rx2 = rx2+c(1); zx2 = zx2+c(2);
+ end
+ end
+ c = sort([rx2a rx2 rx2b]); rb(ia) = c(2);
+ c = sort([zx2a zx2 zx2b]); zb(ia) = c(2);
+ else
+ disp('trace_boundary: error in tracing boundary. Could not find flux value below psibry.')
+ end
+ else
+ rb(ia) = spline(p(3:j),r(3:j),psibry); zb(ia) = spline(p(3:j),z(3:j),psibry);
+ end
+ end
+ if ~ilimited & exist('limdata','var')
+ k = find(zlm<max(zb) & zlm>min(zb));
+ if max(psilim(k))>psibry
+ ilimited = 1;
+ [psibry, ilim] = max(psilim(k));
+ r0 = rlm(k(ilim)); z0 = zlm(k(ilim));
+ th0 = angle(r0-ra+sqrt(-1)*(z0-za)); dth = 2*pi/na;
+ % Trace boundary again because it was limited
+ rb([1 na+1]) = r0; zb([1 na+1]) = z0;
+ for ia = 2:na
+ th = th0+dth*(ia-1);
+ r = ra+rho*cos(th); z = za+rho*sin(th);
+ p = interp2(rg,zg,psizr,r,z,'spline');
+ j = 5;
+ while p(j+1)<p(j) & j<99 & p(j)>psibry, j=j+1; end
+ if p(j) > psibry
+ if j<99 % In this case we are close to the other x-point
+ rx2 = r(j); zx2 = z(j);
+ rx2a = r(min(100,j+3)); zx2a = z(min(100,j+3));
+ rx2b = r(max(1,j-3)); zx2b = z(max(1,j-3));
+ for k=1:9
+ iz = max(find(zg<zx2)); ir = max(find(rg<rx2)); ix = iz+nz*ir;
+ iix = ix+neighbors;
+ pp = psizr(iix);
+ tr = (rx2-rg(ir))/dr; tz = (zx2-zg(iz))/dz;
+ wx = reshape(([1 tz tz^2 tz^3]*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ wxz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrr = reshape(([1 tz tz^2 tz^3]*mx)'*[0 0 2 6*tr]/dr^2*mx,1,16);
+ wxzz = reshape(([0 0 2 6*tz]/dz^2*mx)'*[1 tr tr^2 tr^3]*mx,1,16);
+ wxrz = reshape(([0 1 2*tz 3*tz^2]/dz*mx)'*[0 1 2*tr 3*tr^2]/dr*mx,1,16);
+ xshift_Bx = -inv([pp*wxrr' pp*wxrz';pp*wxrz' pp*wxzz']);
+ c = xshift_Bx*[pp*wxr'; pp*wxz']; rx2 = rx2+c(1); zx2 = zx2+c(2);
+ end
+ c = sort([rx2a rx2 rx2b]); rb(ia) = c(2);
+ c = sort([zx2a zx2 zx2b]); zb(ia) = c(2);
+ else
+ disp('trace_boundary: error in tracing boundary. Could not find flux value below psibry.')
+ end
+ else
+ rb(ia) = spline(p(3:j),r(3:j),psibry); zb(ia) = spline(p(3:j),z(3:j),psibry);
+ end
+ end
+ end
+ end
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/matlab/D3D/efit/private/write_afile.m b/matlab/D3D/efit/private/write_afile.m
new file mode 100644
index 00000000..2a4d8031
--- /dev/null
+++ b/matlab/D3D/efit/private/write_afile.m
@@ -0,0 +1,634 @@
+function write_afile(equil,eqversion,aday)
+% USAGE: write_afile(equil,eqversion,aday)
+%
+% PURPOSE: write afiles
+%
+% INPUTS: equil: equilibrium on the toksys format (that is returned by read_mds_eqdsk)
+% eqversion: (optional) source of equilibrium (default is TOKSYS)
+% aday: (optional) the date the equilibrium was created
+%
+% OUTPUTS: afiles to disk
+%
+% RESTRICTIONS:
+%
+
+% VERSION @(#)write_afile.m 1.3 10/09/14
+%
+% WRITTEN BY: Anders Welander ON 6/1/11
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ xlim = 0; ylim = 0; % fixes matlab bug
+ mu0 = 4e-7*pi;
+
+ struct_to_ws(equil);
+ if exist('adata','var') & isstruct(adata)
+ struct_to_ws(adata);
+ end
+ if ~exist('shotnum','var')
+ shotnum = 0;
+ end
+ if ~exist('dr','var') & exist('dw','var')
+ dr = dw;
+ end
+ if ~exist('dz','var') & exist('dh','var')
+ dz = dh;
+ end
+
+ if ~exist('uday','var')
+ uday = datestr(date,'dd-mmm-yyyy');
+ uday = uday([1:7 10:11]);
+ end
+ while length(uday)<10, uday = [' ' uday]; end
+ if ~exist('mfvers','var')
+ mfvers = datestr(date,23);
+ end
+
+ if ~exist('tims','var')
+ if exist('time','var')
+ tims = round(1000000*time)/1000;
+ else
+ tims = 0;
+ end
+ end
+ ss = num2str(shotnum); while length(ss)<6, ss = ['0' ss]; end
+ tt = num2str(fix(tims)); while length(tt)<5, tt = ['0' tt]; end
+ if mod(tims,1), tt = [tt '_' num2str(mod(tims,1))]; end
+ filename = ['a' ss '.' tt];
+
+ fid = fopen(filename,'w');
+
+% 1 -------------------------- Start writing file
+ fprintf(fid,'%s %s\n',uday,mfvers);
+
+% 2 --------------------------------------------
+ if ~exist('ktime','var'), ktime = 1; end % ktime is total number of time slices
+ if length(num2str(shotnum))>5
+ fprintf(fid,'%7d',shotnum);
+ else
+ fprintf(fid,'%6d',shotnum);
+ end
+ fprintf(fid,'%16d\n',ktime);
+
+% 3 --------------------------------------------
+ fprintf(fid,' % 11.9E\n',tims); % time in ms
+
+% 4 --------------------------------------------
+ if ~exist('jflag','var'), jflag = 1; end
+ if ~exist('lflag','var'), lflag = 0; end % lflag is an efit error flag
+ if ~exist('limloc','var'), limloc = 'SNB'; end % Type of shape: IN , OUT, TOP, BOT, DN , SNT, SNB, MAR
+ if ~exist('mco2v','var'), mco2v = 3; end % Number of co2v items
+ if ~exist('mco2r','var'), mco2r = 2; end % Number of co2r items
+ if ~exist('fwtqa','var'), fwtqa = 0; end % Fit weight for assumed q on axis
+ if ~exist('qvfit','var'), qvfit = 0; end % A parameter for q-profile model used by fitting logic
+ if fwtqa > 0 & qvfit > 0
+ qmflag='FIX';
+ else
+ qmflag='CLC';
+ end
+ if ~exist('n1old','var'), n1old = 39; end
+ if ~exist('n1new','var'), n1new = 40; end
+ fprintf(fid,'*%8.3f',tims);
+ fprintf(fid,'%14d',jflag);
+ fprintf(fid,'%16d',lflag);
+ fprintf(fid,'%4s',limloc);
+ fprintf(fid,'%4d%4d',mco2v,mco2r);
+ fprintf(fid,'%4s',qmflag);
+ fprintf(fid,'%6d',n1old);
+ fprintf(fid,'%5d',n1new);
+ fprintf(fid,'\n ');
+
+% 5 --------------------------------------------
+ if ~exist('rcentr','var')
+ rcentr = rzero;
+ end
+ if ~exist('tsaisq','var'), tsaisq = 0; end % saisq for 1 or several time slices
+ if ~exist('rcencm','var'), rcencm = rzero*100; end % rcentr*100 for 1 or several time slices
+ if ~exist('bcentr','var'), bcentr = bzero; end % Bt at rcentr for 1 or several time slices
+ if ~exist('pasmat','var'), pasmat = cpasma; end % Measured? plasma current for 1 or several time slices
+ skriv(fid,tsaisq);
+ skriv(fid,rcencm);
+ skriv(fid,bcentr);
+ skriv(fid,pasmat);
+ fprintf(fid,'\n ');
+
+% 6 --------------------------------------------
+ if ~exist('rout','var'), rout = 100*(max(rbbbs(1:nbbbs))+min(rbbbs(1:nbbbs)))/2; end
+ if ~exist('zout','var'), zout = 100*(max(zbbbs(1:nbbbs))+min(zbbbs(1:nbbbs)))/2; end
+ if ~exist('aout','var'), aout = 100*(max(rbbbs(1:nbbbs))-min(rbbbs(1:nbbbs)))/2; end
+ skriv(fid,cpasma);
+ skriv(fid,rout);
+ skriv(fid,zout);
+ skriv(fid,aout);
+ fprintf(fid,'\n ');
+
+% 7 --------------------------------------------
+ if ~exist('eout','var')
+ eout = 100*(max(zbbbs(1:nbbbs))-min(zbbbs(1:nbbbs)))/2/aout;
+ end
+ if ~exist('doutu','var')
+ [ztop, k] = max(zbbbs(1:nbbbs));
+ rtop = 100*rbbbs(k);
+ doutu = (rout-rtop)/aout;
+ end
+ if ~exist('doutl','var')
+ [zbot, k] = min(zbbbs(1:nbbbs));
+ rbot = 100*rbbbs(k);
+ doutl = (rout-rbot)/aout;
+ end
+ if ~exist('vout','var')
+ vout = 0;
+ for j = 1:nh
+ for k = 1:nw
+ if jphi(j,k)~=0
+ vout = vout+2*pi*rg(k)*dr*dz*1e6;
+ end
+ end
+ end
+ end
+ skriv(fid,eout);
+ skriv(fid,doutu);
+ skriv(fid,doutl);
+ skriv(fid,vout);
+ fprintf(fid,'\n ');
+
+% 8 --------------------------------------------
+ if ~exist('rcurrt','var')
+ rcurrt = 0;
+ for j = 1:nh
+ for k = 1:nw
+ rcurrt = rcurrt+jphi(j,k)*rg(k)*100;
+ end
+ end
+ rcurrt = rcurrt/sum(sum(jphi));
+ end
+ if ~exist('zcurrt','var')
+ zcurrt = 0;
+ for j = 1:nh
+ for k = 1:nw
+ zcurrt = zcurrt+jphi(j,k)*zg(j)*100;
+ end
+ end
+ zcurrt = zcurrt/sum(sum(jphi));
+ end
+ aspect = rout/aout;
+ if ~exist('qsta','var')
+ qsta = abs(rcentr*bcentr/mu0/cpasma*(eout^2+1)/2/aspect^2);
+ end
+ if ~exist('betat','var'), betat = 0; end
+ skriv(fid,rcurrt); % R of current centroid? in cm
+ skriv(fid,zcurrt); % Z of current centroid? in cm
+ skriv(fid,qsta);
+ skriv(fid,betat);
+ fprintf(fid,'\n ');
+
+% 9 --------------------------------------------
+ if ~exist('betap','var'), betap = 0; end
+ if ~exist('ali','var')
+ if exist('li','var')
+ ali = li;
+ else
+ ali = 0;
+ end
+ end
+ if ~exist('oleft','var') % Closest distance to limiter
+ oleft = 1e10;
+ end
+ if ~exist('oright','var'), oright = 1e10; end
+ skriv(fid,betap);
+ skriv(fid,ali);
+ skriv(fid,oleft);
+ skriv(fid,oright);
+ fprintf(fid,'\n ');
+
+% 10 --------------------------------------------
+ if ~exist('otop','var'), otop = 1e10; end
+ if ~exist('obott','var'), obott = 1e10; end
+ if ~exist('qpsib','var')
+ psibar = linspace(0,1,nw);
+ qpsib = spline(psibar(1:nw-1),qpsi(1:nw-1),0.95);
+ end
+ if ~exist('vertn','var'), vertn = 0; end
+ skriv(fid,otop);
+ skriv(fid,obott);
+ skriv(fid,qpsib);
+ skriv(fid,vertn);
+ fprintf(fid,'\n ');
+
+% 11 --------------------------------------------
+ if ~exist('rco2v','var'), rco2v = zeros(1,mco2v); end
+ if ~exist('dco2v','var'), dco2v = zeros(1,mco2v); end
+ for j = 1:mco2v
+ skriv(fid,rco2v(j));
+ end
+ fprintf(fid,'\n ');
+
+% 12 --------------------------------------------
+ for j = 1:mco2v
+ skriv(fid,dco2v(j));
+ end
+ fprintf(fid,'\n ');
+
+% 13 --------------------------------------------
+ if ~exist('rco2r','var'), rco2r = zeros(1,mco2r); end
+ if ~exist('dco2r','var'), dco2r = zeros(1,mco2r); end
+ for j = 1:mco2r
+ skriv(fid,rco2r(j));
+ end
+ fprintf(fid,'\n ');
+
+% 14 --------------------------------------------
+ for j = 1:mco2r
+ skriv(fid,dco2r(j));
+ end
+ fprintf(fid,'\n ');
+
+% 15 --------------------------------------------
+ if ~exist('shearb','var'), shearb = 0; end
+ if ~exist('bpolav','var'), bpolav = 0; end
+ if ~exist('s1','var'), s1 = 0; end
+ if ~exist('s2','var'), s2 = 0; end
+ skriv(fid,shearb);
+ skriv(fid,bpolav);
+ skriv(fid,s1);
+ skriv(fid,s2);
+ fprintf(fid,'\n ');
+
+% 16 --------------------------------------------
+ if ~exist('s3','var'), s3 = 0; end
+ if ~exist('qout','var'), qout = 0; end
+ if ~exist('olefs','var'), olefs = 0; end
+ if ~exist('orighs','var'), orighs = 0; end
+ skriv(fid,s3);
+ skriv(fid,qout);
+ skriv(fid,olefs);
+ skriv(fid,orighs);
+ fprintf(fid,'\n ');
+
+% 17 --------------------------------------------
+ if ~exist('otops','var'), otops = 0; end
+ if ~exist('sibdry','var'), sibdry = 0; end
+ if ~exist('areao','var'), % Area of plasma in cm^2
+ areao = 0;
+ for j = 1:nh
+ for k = 1:nw
+ if jphi(j,k) ~= 0
+ areao = areao+dr*dz*1e4;
+ end
+ end
+ end
+ end
+ if ~exist('wplasm','var'), wplasm = 0; end
+ skriv(fid,otops);
+ skriv(fid,sibdry);
+ skriv(fid,areao);
+ skriv(fid,wplasm);
+ fprintf(fid,'\n ');
+
+% 18 --------------------------------------------
+ if ~exist('terror','var'), terror = 0; end
+ if ~exist('elongm','var'), elongm = 0; end
+ if ~exist('qqmagx','var'), qqmagx = qpsi(1); end
+ if ~exist('cdflux','var'), cdflux = 0; end
+ skriv(fid,terror);
+ skriv(fid,elongm);
+ skriv(fid,qqmagx);
+ skriv(fid,cdflux);
+ fprintf(fid,'\n ');
+
+% 19 --------------------------------------------
+ if ~exist('alpha','var'), alpha = 0; end
+ if ~exist('rttt','var'), rttt = 0; end
+ if ~exist('psiref','var'), psiref = 0; end
+ if ~exist('xndnt','var'), xndnt = 0; end
+ skriv(fid,alpha);
+ skriv(fid,rttt);
+ skriv(fid,psiref);
+ skriv(fid,xndnt);
+ fprintf(fid,'\n ');
+
+% 20 --------------------------------------------
+ if ~exist('rseps1','var'), rseps1 = 0; end
+ if ~exist('zseps1','var'), zseps1 = 0; end
+ if ~exist('rseps2','var'), rseps2 = 0; end
+ if ~exist('zseps2','var'), zseps2 = 0; end
+% THIS NEEDS TO BE FIXED to correctly handle empty objects.
+ if isempty(rseps1), rseps1=0; end
+ if isempty(zseps1), zseps1=0; end
+ if isempty(rseps2), rseps2=0; end
+ if isempty(zseps2), zseps2=0; end
+ skriv(fid,rseps1);
+ skriv(fid,zseps1);
+ skriv(fid,rseps2);
+ skriv(fid,zseps2);
+ fprintf(fid,'\n ');
+
+% 21 --------------------------------------------
+ if ~exist('sepexp','var'), sepexp = 0; end
+ if ~exist('obots','var'), obots = 0; end
+ if ~exist('btaxp','var'), btaxp = 0; end
+ if ~exist('btaxv','var'), btaxv = 0; end
+ skriv(fid,sepexp);
+ skriv(fid,obots);
+ skriv(fid,btaxp);
+ skriv(fid,btaxv);
+ fprintf(fid,'\n ');
+
+% 22 --------------------------------------------
+ if ~exist('aaq1','var'), aaq1 = 0; end
+ if ~exist('aaq2','var'), aaq2 = 0; end
+ if ~exist('aaq3','var'), aaq3 = 0; end
+ if ~exist('seplim','var'), seplim = 0; end
+ skriv(fid,aaq1);
+ skriv(fid,aaq2);
+ skriv(fid,aaq3);
+ skriv(fid,seplim);
+ fprintf(fid,'\n ');
+
+% 23 --------------------------------------------
+ if ~exist('rmagx','var'), rmagx = 100*rmaxis; end
+ if ~exist('zmagx','var'), zmagx = 100*zmaxis; end
+ if ~exist('simagx','var'), simagx = psimag/2/pi; end
+ if ~exist('taumhd','var'), taumhd = 0; end
+ skriv(fid,rmagx);
+ skriv(fid,zmagx);
+ skriv(fid,simagx);
+ skriv(fid,taumhd);
+ fprintf(fid,'\n ');
+
+% 24 --------------------------------------------
+ if ~exist('betapd','var'), betapd = 0; end
+ if ~exist('betatd','var'), betatd = 0; end
+ if ~exist('wplasmd','var'), wplasmd = 0; end
+ if ~exist('fluxx','var'), fluxx = 0; end
+ skriv(fid,betapd);
+ skriv(fid,betatd);
+ skriv(fid,wplasmd);
+ skriv(fid,fluxx);
+ fprintf(fid,'\n ');
+
+% 25 --------------------------------------------
+ if ~exist('vloopt','var'), vloopt = 0; end
+ if ~exist('taudia','var'), taudia = 0; end
+ if ~exist('qmerci','var'), qmerci = 0; end
+ if ~exist('tavem','var'), tavem = 0; end
+ skriv(fid,vloopt);
+ skriv(fid,taudia);
+ skriv(fid,qmerci);
+ skriv(fid,tavem);
+ fprintf(fid,'\n ');
+
+% 26 --------------------------------------------
+ if ~exist('nsilop','var'), nsilop = 1; end
+ if ~exist('magpri','var'), magpri = 1; end
+ if ~exist('nfcoil','var'), nfcoil = 1; end
+ if ~exist('nesum','var')
+ if exist('ecurrt','var')
+ nesum = length(ecurrt);
+ if nesum == 0
+ nesum = 1;
+ ecurrt = 0;
+ end
+ else
+ nesum = 1;
+ end
+ end
+ fprintf(fid,'%5i',nsilop);
+ fprintf(fid,'%5i',magpri);
+ fprintf(fid,'%5i',nfcoil);
+ fprintf(fid,'%5i\n ',nesum);
+
+% 27 --------------------------------------------
+ if ~exist('csilop','var'), csilop = zeros(nsilop,1); end
+ if ~exist('cmpr2','var'), cmpr2 = zeros(magpri,1); end
+ if ~exist('ccbrsp','var'), ccbrsp = zeros(nfcoil,1); end
+ if ~exist('ecurrt','var'), ecurrt = zeros(nesum,1); end
+ for j = 1:nsilop+magpri
+ if j<=nsilop
+ skriv(fid,csilop(j));
+ else
+ skriv(fid,cmpr2(j-nsilop));
+ end
+ if mod(j,4) == 0
+ fprintf(fid,'\n ');
+ end
+ end
+ if mod(j,4) ~= 0
+ fprintf(fid,'\n ');
+ end
+
+% 27+ --------------------------------------------
+ for j = 1:nfcoil
+ skriv(fid,ccbrsp(j));
+ if mod(j,4) == 0
+ fprintf(fid,'\n ');
+ end
+ end
+ if mod(j,4) ~= 0
+ fprintf(fid,'\n ');
+ end
+
+% 28+ --------------------------------------------
+ for j = 1:nesum
+ skriv(fid,ecurrt(j));
+ if mod(j,4) == 0
+ fprintf(fid,'\n ');
+ end
+ end
+ if mod(j,4) ~= 0
+ fprintf(fid,'\n ');
+ end
+
+% 29+ --------------------------------------------
+ if ~exist('pbinj','var'), pbinj = 0; end
+ if ~exist('rvsin','var'), rvsin = 0; end
+ if ~exist('zvsin','var'), zvsin = 0; end
+ if ~exist('rvsout','var'), rvsout = 0; end
+ skriv(fid,pbinj);
+ skriv(fid,rvsin);
+ skriv(fid,zvsin);
+ skriv(fid,rvsout);
+ fprintf(fid,'\n ');
+
+ if ~exist('zvsout','var'), zvsout = 0; end
+ if ~exist('vsurfa','var'), vsurfa = 0; end
+ if ~exist('wpdot','var'), wpdot = 0; end
+ if ~exist('wbdot','var'), wbdot = 0; end
+ skriv(fid,zvsout);
+ skriv(fid,vsurfa);
+ skriv(fid,wpdot);
+ skriv(fid,wbdot);
+ fprintf(fid,'\n ');
+
+ if ~exist('slantu','var'), slantu = 0; end
+ if ~exist('slantl','var'), slantl = 0; end
+ if ~exist('zuperts','var'), zuperts = 0; end
+ if ~exist('chipre','var'), chipre = 0; end
+ skriv(fid,slantu);
+ skriv(fid,slantl);
+ skriv(fid,zuperts);
+ skriv(fid,chipre);
+ fprintf(fid,'\n ');
+
+ if ~exist('cjor95','var'), cjor95 = 0; end
+ if ~exist('pp95','var'), pp95 = 0; end
+ if ~exist('ssep','var'), ssep = 0; end
+ if ~exist('yyy2','var'), yyy2 = 0; end
+ skriv(fid,cjor95);
+ skriv(fid,pp95);
+ skriv(fid,ssep);
+ skriv(fid,yyy2);
+ fprintf(fid,'\n ');
+
+ if ~exist('xnnc','var'), xnnc = 0; end
+ if ~exist('cprof','var'), cprof = 0; end
+ if ~exist('oring','var'), oring = 0; end
+ if ~exist('cjor0','var'), cjor0 = 0; end
+ skriv(fid,xnnc);
+ skriv(fid,cprof);
+ skriv(fid,oring);
+ skriv(fid,cjor0);
+ fprintf(fid,'\n ');
+
+ if ~exist('fexpan','var'), fexpan = 0; end
+ if ~exist('qqmin','var'), qqmin = min(qpsi); end
+ if ~exist('chigamt','var'), chigamt = 0; end
+ if ~exist('ssi01','var'), ssi01 = 0; end
+ skriv(fid,fexpan);
+ skriv(fid,qqmin);
+ skriv(fid,chigamt);
+ skriv(fid,ssi01);
+ fprintf(fid,'\n ');
+
+ if ~exist('fexpvs','var'), fexpvs = 0; end
+ if ~exist('sepnose','var'), sepnose = 0; end
+ if ~exist('ssi95','var'), ssi95 = 0; end
+ if ~exist('rqqmin','var'), rqqmin = 0; end
+ skriv(fid,fexpvs);
+ skriv(fid,sepnose);
+ skriv(fid,ssi95);
+ skriv(fid,rqqmin);
+ fprintf(fid,'\n ');
+
+ if ~exist('cjor99','var'), cjor99 = 0; end
+ if ~exist('cjlave','var'), cjlave = 0; end
+ if ~exist('rmidin','var') % min(rbbbs) at zbbbs=0
+ k = find(zbbbs(1:nbbbs-2).*zbbbs(2:nbbbs-1)<=0);
+ [dum, kk] = sort(rbbbs(k)); k = k(kk(1));
+ rmidin = interp1(zbbbs(k+[0 1]),rbbbs(k+[0 1]),0);
+ end
+ if ~exist('rmidout','var')
+ k = find(zbbbs(1:nbbbs-2).*zbbbs(2:nbbbs-1)<=0);
+ [dum, kk] = sort(rbbbs(k)); k = k(kk(end));
+ rmidout = interp1(zbbbs(k+[0 1]),rbbbs(k+[0 1]),0);
+ end
+ skriv(fid,cjor99);
+ skriv(fid,cjlave);
+ skriv(fid,rmidin);
+ skriv(fid,rmidout);
+ fprintf(fid,'\n ');
+
+ if ~exist('psurfa','var'), psurfa = 0; end
+ if ~exist('peak','var'), peak = 0; end
+ if ~exist('dminux','var'), dminux = 0; end
+ if ~exist('dminlx','var'), dminlx = 0; end
+ skriv(fid,psurfa);
+ skriv(fid,peak);
+ skriv(fid,dminux);
+ skriv(fid,dminlx);
+ fprintf(fid,'\n ');
+
+ if ~exist('dolubaf','var'), dolubaf = 0; end
+ if ~exist('dolubafm','var'), dolubafm = 0; end
+ if ~exist('diludom','var'), diludom = 0; end
+ if ~exist('diludomm','var'), diludomm = 0; end
+ skriv(fid,dolubaf);
+ skriv(fid,dolubafm);
+ skriv(fid,diludom);
+ skriv(fid,diludomm);
+ fprintf(fid,'\n ');
+
+ if ~exist('ratsol','var'), ratsol = 0; end
+ if ~exist('rvsiu','var'), rvsiu = 0; end
+ if ~exist('zvsiu','var'), zvsiu = 0; end
+ if ~exist('rvsid','var'), rvsid = 0; end
+ skriv(fid,ratsol);
+ skriv(fid,rvsiu);
+ skriv(fid,zvsiu);
+ skriv(fid,rvsid);
+ fprintf(fid,'\n ');
+
+ if ~exist('zvsid','var'), zvsid = 0; end
+ if ~exist('rvsou','var'), rvsou = 0; end
+ if ~exist('zvsou','var'), zvsou = 0; end
+ if ~exist('rvsod','var'), rvsod = 0; end
+ skriv(fid,zvsid);
+ skriv(fid,rvsou);
+ skriv(fid,zvsou);
+ skriv(fid,rvsod);
+ fprintf(fid,'\n ');
+
+ if ~exist('zvsod','var'), zvsod = 0; end
+ if ~exist('condno','var'), condno = 0; end
+ if ~exist('psin32','var'), psin32 = 0; end
+ if ~exist('psin21','var'), psin21 = 0; end
+ skriv(fid,zvsod);
+ skriv(fid,condno);
+ skriv(fid,psin32);
+ skriv(fid,psin21);
+ fprintf(fid,'\n ');
+
+ if ~exist('rq32in','var'), rq32in = 0; end
+ if ~exist('rq21top','var'), rq21top = 0; end
+ if ~exist('chilibt','var'), chilibt = 0; end
+ if ~exist('ali3','var'), ali3 = 0; end
+ skriv(fid,rq32in);
+ skriv(fid,rq21top);
+ skriv(fid,chilibt);
+ skriv(fid,ali3);
+ fprintf(fid,'\n ');
+
+ fprintf(fid,'\n ');
+
+ if ~exist('header','var'), header = blanks(42); end
+ if ~exist('fit_type','var'), fit_type = blanks(3); end
+ fprintf(fid,' %42s',header);
+ fprintf(fid,' %3s',fit_type);
+
+ fclose(fid);
+ return
+
+ % function to print to file EXACTLY like Fortran's write(fid,1040) does
+ function skriv(fid,val)
+ if val >= 0
+ fprintf(fid,' %11.9E',val);
+ else
+ fprintf(fid,'%11.9E',val);
+ end
+ return
+
+ % function copied from efit code
+ function dismin = dslant(x,y,np,xmin,xmax,ymin,ymax,x1,y1,x2,y2)
+ dismin=inf;
+ delx=x2-x1;
+ dely=y2-y1;
+ dels=sqrt(delx^2+dely^2);
+ nn=dels/0.002;
+ nn=max(5,nn);
+ delx=delx/(nn-1);
+ dely=dely/(nn-1);
+ for j=1:nn
+ xw=x1+delx *(j-1);
+ yw=y1+dely *(j-1);
+ for m=1:np
+ if x(m) > xmin & x(m)<xmax & y(m) > ymin & y(m) < ymax
+ disw=sqrt((xw-x(m))^2+(yw-y(m))^2);
+ dismin=min(dismin,disw);
+ end
+ end
+ end
+ dismin=dismin*100;
+ return
+
diff --git a/matlab/D3D/efit/private/write_cc_file.m b/matlab/D3D/efit/private/write_cc_file.m
new file mode 100644
index 00000000..ebdb778b
--- /dev/null
+++ b/matlab/D3D/efit/private/write_cc_file.m
@@ -0,0 +1,55 @@
+function [cc_filename,ier] = write_cc_file(eq);
+% USAGE: [cc_filename,ier] = write_cc_file(equilibria);
+%
+% PURPOSE: return absolute path to newly generated cc_file
+%
+% INPUTS: filename:name of geqdsk file
+%
+% OUTPUTS: cc_filename: string containing name of cc_file
+%
+% RESTRICTIONS:
+%
+% METHOD:
+%
+% cc_file= coil current file name [cc=load(cc_file)] or coil current vector
+% use when coil currents not available from efit (or to override)
+%
+% ohmic coil current is first, then poloidal shaping currents
+% (units = MA-turns)
+%
+% WRITTEN BY: Matthew J. Lanctot
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ %inits
+ cc_filename='';
+ ier = 0;
+
+ cc_filename = ['cc' num2str(eq.shotnum) '.' sprintf('%05g',eq.time*1e3)];
+
+ if isfield(eq,'brsp')
+ brsp = eq.brsp; %have to fit PF coil, even with small weight
+ else
+ ier=1;
+ end
+
+ if isfield(eq,'ecurrt')
+ ecurrt = eq.ecurrt; %have to fit OH coil, even with small weight
+ else
+ ier=2;
+ end
+
+ switch ier
+ case 1
+ warning('gfile does not contain brsp data. were PF coils fit?')
+ return
+ case 2
+ warning('gfile does not contain ecurrt data. was OH coil fit?')
+ return
+ otherwise %write cc file
+ nf = length(unique( eq.fcid));
+ data = [ecurrt brsp(1:nf)']';
+ dlmwrite(cc_filename,data,'\t',0,0,14);
+ end
+
+
diff --git a/matlab/D3D/efit/private/write_eqdsk_efit.m b/matlab/D3D/efit/private/write_eqdsk_efit.m
new file mode 100644
index 00000000..6ddb7334
--- /dev/null
+++ b/matlab/D3D/efit/private/write_eqdsk_efit.m
@@ -0,0 +1,65 @@
+function [equilibria, neq, eqraw, cc_filename] = write_eqdsk(shotnum,tefit,efit_source,tokamak,write_cc)
+% USAGE: [equilibria, neq, eqraw] = write_eqdsk(shotnum,tefit,efit_source,tokamak,write_cc)
+%
+% PURPOSE: Retrieves eqdsk data from MDS and saves g,a files to dir
+%
+% INPUTS:
+%
+% shotnum: Shotnumber for equilibrium eg. 146970
+% tefit: time of efit, can be two element array specifying time range as in read_eq
+% efit_source: EFIT source eg. 'EFIT01'
+% tokamak: device name eg. 'd3d'
+%
+% OUTPUTS: eq = equilibrium data on toksys form, that can be used with
+% cc_efit_to_tok and build_tokamak_system
+% neq = number of equilibria returned
+% eqraw = 'raw' eq data on original form before conversion to toksys convention
+%
+% RESTRICTIONS:
+%
+% METHOD:
+%
+%
+% WRITTEN BY: Matthew J. Lanctot
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+[equilibria, neq, eqraw] = read_eq(shotnum,tefit,efit_source,tokamak);
+
+cc_filename='';
+if write_cc
+ if neq==1
+ [cc_filename,ier] = write_cc_file(equilibria);
+ else
+ for k=1:neq
+ [cc_filename{k},ier] = write_cc_file(equilibria.gdata(k));
+ end
+ end
+end
+
+if isfield(equilibria,'gdata')
+ % to simply gfile names, do not support sub-ms time resolution if times are unique
+ uttmp = unique(equilibria.time-mod(equilibria.time,1e-3));
+ if length(uttmp)==length(equilibria.time)
+ equilibria.time = equilibria.time-mod(equilibria.time,1e-3);
+ fix_gtime = 1;
+ else
+ fix_gtime = 0;
+ end
+
+ %write files
+ for k=1:length(equilibria.time)
+ if fix_gtime==1
+ equilibria.gdata(k).time=equilibria.gdata(k).time-mod(equilibria.gdata(k).time,1e-3);
+ end
+ write_gfile(equilibria.gdata(k));
+ % afile from read_eq doesn't have rzero
+ %if isfield(equilibria,'adata'), write_afile(equilibria.adata(k));end
+ end
+else
+ write_gfile(equilibria);
+end
diff --git a/matlab/D3D/efit/private/write_gfile.m b/matlab/D3D/efit/private/write_gfile.m
new file mode 100644
index 00000000..1c7c9045
--- /dev/null
+++ b/matlab/D3D/efit/private/write_gfile.m
@@ -0,0 +1,320 @@
+function write_gfile(eq)
+% USAGE: write_gfile(eq)
+%
+% PURPOSE: write gfiles
+%
+% INPUTS: eq: equilibrium on the toksys format (that is returned by read_mds_eqdsk)
+%
+% OUTPUTS: gfiles to disk
+%
+% RESTRICTIONS:
+%
+
+% VERSION @(#)write_gfile.m 1.3 09/17/12
+%
+% WRITTEN BY: Anders Welander ON 6/1/11
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ xlim = 0; ylim = 0; % fixes matlab bug
+
+ struct_to_ws(eq);
+
+ if ~exist('pasmat','var')
+ pasmat = cpasma; % pasmat is perhaps measured plasma current
+ end
+
+ if ~exist('ecurrt','var')
+ ecurrt = [];
+ end
+
+ tims = round(1000000*time)/1000;
+ ss = num2str(shotnum); while length(ss)<6, ss = ['0' ss]; end
+ tt = num2str(fix(tims)); while length(tt)<5, tt = ['0' tt]; end
+ if mod(tims,1), tt = [tt '_' num2str(mod(tims,1))]; end
+ filename = ['g' ss '.' tt];
+
+ fid = fopen(filename,'w');
+
+ fprintf(fid,' EFITD ');
+ fprintf(fid,datestr(now,23));
+
+ str = num2str(shotnum);
+ while length(str)<6, str = [' ' str]; end
+ str = [' #' str];
+ fprintf(fid,str);
+
+ str = [num2str(tims) 'ms'];
+ while length(str)<8, str = [' ' str]; end
+ str = [str blanks(8)];
+ fprintf(fid,str);
+
+ fprintf(fid,'%4d%4d%4d\n',3,nw,nh);
+
+ if ~exist('xdim','var')
+ xdim = rg(end)-rg(1);
+ end
+ if ~exist('zdim','var')
+ zdim = zg(end)-zg(1);
+ end
+ if ~exist('zmid','var')
+ zmid = (zg(end)+zg(1))/2;
+ end
+ fprintf(fid,'% 11.9E',xdim);
+ fprintf(fid,'% 11.9E',zdim);
+ fprintf(fid,'% 11.9E',rzero);
+ fprintf(fid,'% 11.9E',rg(1));
+ fprintf(fid,'% 11.9E\n',zmid);
+
+ fprintf(fid,'% 11.9E',rmaxis);
+ fprintf(fid,'% 11.9E',zmaxis);
+ fprintf(fid,'% 11.9E',ssimag);
+ fprintf(fid,'% 11.9E',ssibry);
+ fprintf(fid,'% 11.9E\n',bzero); % bcentr = bzero?
+
+ xdum = 0; % introducing xdum...
+
+ fprintf(fid,'% 11.9E',cpasma);
+ fprintf(fid,'% 11.9E',ssimag);
+ fprintf(fid,'% 11.9E',xdum);
+ fprintf(fid,'% 11.9E',rmaxis);
+ fprintf(fid,'% 11.9E\n',xdum);
+
+ fprintf(fid,'% 11.9E',zmaxis);
+ fprintf(fid,'% 11.9E',xdum);
+ fprintf(fid,'% 11.9E',ssibry);
+ fprintf(fid,'% 11.9E',xdum);
+ fprintf(fid,'% 11.9E\n',xdum);
+
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',fpol(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',pres(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',-sign(pasmat)*ffprim(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',-sign(pasmat)*pprime(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ for j = 1:nh
+ for k = 1:nw
+ fprintf(fid,'% 11.9E',psirz(k,j));
+ if mod((j-1)*nw+k,5) == 0, fprintf(fid,'\n'); end
+ end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',qpsi(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+
+ if ~exist('limitr','var')
+ limitr = length(xlim);
+ end
+ fprintf(fid,'% 5d',nbbbs);
+ fprintf(fid,'% 5d\n',limitr);
+ for j = 1:nbbbs
+ fprintf(fid,'% 11.9E',rbbbs(j));
+ if mod(2*j-1,5) == 0, fprintf(fid,'\n'); end
+ fprintf(fid,'% 11.9E',zbbbs(j));
+ if mod(2*j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(2*nbbbs,5)~=0
+ fprintf(fid,'\n');
+ end
+ for j = 1:limitr
+ fprintf(fid,'% 11.9E',xlim(j));
+ if mod(2*j-1,5) == 0, fprintf(fid,'\n'); end
+ fprintf(fid,'% 11.9E',ylim(j));
+ if mod(2*j,5) == 0, fprintf(fid,'\n'); end
+ end
+
+ % Rotation information
+ if ~exist('kvtor','var')
+ kvtor = 0;
+ end
+ if ~exist('rvtor','var')
+ rvtor = 1.7;
+ end
+ if ~exist('nmass','var')
+ nmass = 0;
+ end
+ if mod(2*limitr,5)~=0
+ fprintf(fid,'\n');
+ end
+
+ fprintf(fid,'% 5d',kvtor);
+ fprintf(fid,'% 11.9E',rvtor);
+ fprintf(fid,'% 5d\n',nmass);
+
+ % write out rotation information
+ if kvtor>0
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',pressw(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ fprintf(fid,'\n');
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',-sign(pasmat)*pwprim(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ fprintf(fid,'\n');
+ end
+
+ % write out ion mass density profile if available
+ if nmass>0
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',dmion(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ fprintf(fid,'\n');
+ end
+
+ if ~exist('rhovn','var')
+ rhovn(1:nw) = 0;
+ end
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',rhovn(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+
+ if ~exist('keecur','var')
+ keecur = 0;
+ end
+ fprintf(fid,'% 5d\n',keecur);
+ if keecur>0
+ for j = 1:nw
+ fprintf(fid,'% 11.9E',-sign(pasmat)*epoten(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ end
+
+ if exist('jphi','var') % if true then iplcout was 1
+ fprintf(fid,'% 5d',nw);
+ fprintf(fid,'% 5d',nh);
+ fprintf(fid,'% 5d',shotnum);
+ fprintf(fid,'% 5d\n',tims);
+
+ fprintf(fid,'% 11.9E',rg(1));
+ fprintf(fid,'% 11.9E',rg(end));
+ fprintf(fid,'% 11.9E',zg(1));
+ fprintf(fid,'% 11.9E\n',zg(end));
+
+ if ~exist('nfcoil','var')
+ nfcoil = length(brsp);
+ end
+ for j = 1:nfcoil
+ fprintf(fid,'% 11.9E',brsp(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+
+ if ~exist('nesum','var')
+ nesum = length(ecurrt);
+ end
+ for j = 1:nesum
+ fprintf(fid,'% 11.9E',ecurrt(j));
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+
+ for j = 1:nw*nh
+ fprintf(fid,'% 11.9E',pcurrt(j)); % Need to check r,z order
+ if mod(j,5) == 0, fprintf(fid,'\n'); end
+ end
+ if mod(j,5), fprintf(fid,'\n'); end
+ end
+
+ out1{ 1,1} = 'ishot';
+ out1{ 1,2} = 'shot';
+ out1{ 2,1} = 'itime';
+ out1{ 2,2} = 'tims';
+ out1{ 3,1} = 'betap0';
+ out1{ 3,2} = 'betap';
+ out1{ 4,1} = 'btor';
+ out1{ 4,2} = 'bzero';
+ out1{ 5,1} = 'rcentr';
+ out1{ 5,2} = 'rzero';
+ out1{ 6,1} = 'rbdry';
+ out1{ 6,2} = 'rbbbs';
+ out1{ 7,1} = 'zbdry';
+ out1{ 7,2} = 'zbbbs';
+ out1{ 8,1} = 'nbdry';
+ out1{ 8,2} = 'nbbbs';
+ out1{ 9,1} = 'mw';
+ out1{ 9,2} = 'nw';
+ out1{ 9,3} = 'nr';
+ out1{10,1} = 'mh';
+ out1{10,2} = 'nh';
+ out1{10,3} = 'nz';
+ out1{11,1} = 'psirz';
+ out1{12,1} = 'xlim';
+ out1{13,1} = 'ylim';
+ out1{14,1} = 'limitr';
+ out1{15,1} = 'brsp';
+ out1{15,2} = 'cc';
+ out1{15,3} = 'ic';
+
+ fprintf(fid,[' &OUT1' 10]);
+ for i = 1:size(out1,1)
+ S = upper(out1{i,1});
+ for j = 1:size(out1,2)
+ s = out1{i,j};
+ if ~isempty(s) & exist(s,'var')
+ dum = eval(s);
+ n = min(size(dum));
+ for k = 1:numel(dum)
+ if k == 1
+ V = [' ' S ' = '];
+ else
+ V = 32+zeros(1,length(S)+3);
+ end
+ if dum(k) > 0
+ x = ' ';
+ else
+ x = '';
+ end
+ m = num2str(dum(k));
+ if isa(dum,'logical')
+ if dum(k)
+ m = 'T';
+ else
+ m = 'F';
+ end
+ end
+ if k/n == round(k/n) % Add new line at the end
+ if (k-1)/n == round((k-1)/n) % Add V at beginning of line
+ fprintf(fid,[V x m ',' 10]);
+ else
+ fprintf(fid,[x m ', ' 10]);
+ end
+ else
+ if (k-1)/n == round((k-1)/n) % Add V at beginning of line
+ fprintf(fid,[V x m ',']);
+ else
+ fprintf(fid,[x m ', ']);
+ end
+ end
+ end
+ break % Don't look for alternatives
+ end
+ end
+ end
+ fprintf(fid,[' /' 10]);
+
+ fclose(fid);
diff --git a/matlab/D3D/efit/private/write_kfile.m b/matlab/D3D/efit/private/write_kfile.m
new file mode 100644
index 00000000..2cecd1c8
--- /dev/null
+++ b/matlab/D3D/efit/private/write_kfile.m
@@ -0,0 +1,190 @@
+function write_kfile(IN1,INWANT,INS,efitin)
+% USAGE: write_kfile(eq) % Creates k file from TokSys equilibrium
+% Alternatively, create k file with exact content of variables:
+% write_kfile(IN1,INWANT,INS,efitin)
+%
+% PURPOSE: write k-files
+%
+% INPUTS: eq, TokSys equilibrium with fields
+% List shows: field name, alternative names, description
+% Note that some alternatives aren't exactly the same things
+%shot, a shot number
+%time, a time in seconds
+%plasma, rog, cpasma, MEASURED total plasma current [A]
+%expmp2, bp MEASURED signals from magnetic probes
+%coils, fl, MEASURED signals from flux loops
+%btor, bzero, toroidal field at rcentr, rzero
+%rcentr, rzero, radius where btor, bzero is measured
+%psibit, physics units per digitizer bit for flux loops
+%bitmpi, physics units per digitizer bit for probes
+%bitip, physics units per digitizer bit for rogowski
+%bitfc, physics units per digitizer bit for F coil currents
+%bitec, physics units per digitizer bit for E coil currents
+%tgamma, tangent for measured pitch angles
+%sgamma, uncertainties in tgamma
+%fitdelz, boolean flag
+
+
+%ISHOT, a shot number
+%ITIME, a time in ms
+%PLASMA, measured total plasma current [A]
+%EXPMP2, measured signals from magnetic probes
+%COILS, measured signals from flux loops
+%BTOR, toroidal field at RCENTR
+%RCENTR, radius where BTOR measured
+%PSIBIT, physics units per digitizer bit for flux loops
+%BITMPI, physics units per digitizer bit for probes
+%BITIP, physics units per digitizer bit for rogowski
+%BITFC, physics units per digitizer bit for F coil currents
+%BITEC, physics units per digitizer bit for E coil currents
+%tgamma, tangent for measured pitch angles
+%sgamma, uncertainties in tgamma
+%fitdelz, boolean flag
+
+%MORE VARIABLES
+%timeu, some other time, default 0
+%qvfit, ?, default 0
+%fwtsi, fitting weights for flux loops
+%fwtcur, fitting weight for rogowski?, default 1
+%limitr, ?, default -nw?
+%fwtmp2, fitting weights for magnetic probes
+%kffcur, number of knots in spline for ffprim
+%kppcur, number of knots in spline for pprime
+%fwtqa, ?, default 0
+%ierchk, ?, default 1
+%fwtbp, ?, default 0
+%serror, ?, default 3e-2
+%nextra, ?, default 2
+%scrape, ?, default 4e-2
+%itrace, ?, default 1
+%xltype, ?, default 0
+%rcentr, rzero, radius where btor, bzero is measured
+%
+%
+
+% OUTPUTS: kfiles to disk
+%
+% RESTRICTIONS:
+%
+
+% WRITTEN BY: Anders Welander ON 2016-05-31
+%
+% MODIFICATION HISTORY:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if nargin == 1
+ struct_to_ws(IN1);
+
+ if ~exist('ishot','var') & exist('shot','var')
+ ishot = shot;
+ end
+
+ tims = round(1000000*time)/1000;
+ ss = num2str(ishot); while length(ss)<6, ss = ['0' ss]; end
+ tt = num2str(fix(tims)); while length(tt)<5, tt = ['0' tt]; end
+ if mod(tims,1), tt = [tt '_' num2str(mod(tims,1))]; end
+ filename = ['k' ss '.' tt];
+
+ in1{ 1,1} = 'ishot';
+ in1{ 1,2} = 'shot';
+ in1{ 2,1} = 'itime';
+ in1{ 2,2} = 'tims';
+ in1{ 3,1} = 'plasma';
+ in1{ 3,2} = 'rog';
+ in1{ 3,3} = 'cpasma';
+ in1{ 4,1} = 'expmp2';
+ in1{ 4,2} = 'bp';
+ in1{ 5,1} = 'coils';
+ in1{ 5,2} = 'fl';
+ in1{ 6,1} = 'btor';
+ in1{ 6,2} = 'bzero';
+ in1{ 7,1} = 'rcentr';
+ in1{ 7,2} = 'rzero';
+ in1{ 8,1} = 'psibit';
+ in1{ 9,1} = 'bitmpi';
+ in1{10,1} = 'bitip';
+ in1{11,1} = 'bitfc';
+ in1{12,1} = 'bitec';
+ in1{13,1} = 'tgamma';
+ in1{14,1} = 'sgamma';
+
+ inwant{1,1} = 'fitdelz';
+
+
+ fid = fopen(filename,'w');
+
+ fprintf(fid,['&IN1 ' 10]);
+ for i = 1:size(in1,1)
+ S = upper(in1{i,1});
+ for j = 1:size(in1,2)
+ s = in1{i,j};
+ if ~isempty(s) & exist(s,'var')
+ dum = eval(s);
+ for k = 1:numel(dum)
+ if k == 1
+ V = [S ' = '];
+ else
+ V = 32+zeros(1,length(S)+3);
+ end
+ if dum(k) > 0
+ x = ' ';
+ else
+ x = '';
+ end
+ m = num2str(dum(k));
+ if isa(dum,'logical')
+ if dum(k)
+ m = 'T';
+ else
+ m = 'F';
+ end
+ end
+ fprintf(fid,[V x m ',' 10]);
+ end
+ break % Don't look for alternatives
+ end
+ end
+ end
+ fprintf(fid,['/' 10]);
+
+ fprintf(fid,['&INWANT ' 10]);
+ for i = 1:size(inwant,1)
+ S = upper(inwant{i,1});
+ for j = 1:size(inwant,2)
+ s = inwant{i,j};
+ if ~isempty(s) & exist(s,'var')
+ dum = eval(s);
+ for k = 1:numel(dum)
+ if k == 1
+ V = [S ' = '];
+ else
+ V = 32+zeros(1,length(S)+3);
+ end
+ if dum(k) > 0
+ x = ' ';
+ else
+ x = '';
+ end
+ m = num2str(dum(k));
+ if isa(dum,'logical')
+ if dum(k)
+ m = 'T';
+ else
+ m = 'F';
+ end
+ end
+ fprintf(fid,[V x m ',' 10]);
+ end
+ break % Don't look for alternatives
+ end
+ end
+ end
+ fprintf(fid,['/' 10]);
+
+ fclose(fid);
+
+else
+ % Create namelist files for variables
+
+end
diff --git a/matlab/D3D/mdsplus_func/private/get_mds_tree.m b/matlab/D3D/mdsplus_func/private/get_mds_tree.m
new file mode 100644
index 00000000..4bf12703
--- /dev/null
+++ b/matlab/D3D/mdsplus_func/private/get_mds_tree.m
@@ -0,0 +1,291 @@
+function [data,ier] = get_mds_tree(shot, tree, server, toupper, verbose)
+%
+% SYNTAX:
+% data= get_mds_tree(shot, tree, server, toupper,verbose); % full call
+% data= get_mds_tree(shot); % defaults to DIII-D
+% data= get_mds_tree(shot, 'EFIT01', 'NSTX'); % for NSTX
+% data= get_mds_tree(shot, 'NB', 'DIII-D'); % example other tree
+%
+% PURPOSE: Get entire selected mds tree from mdsplus database.
+%
+% INPUT: <default>
+% shot = shot number
+% tree = tree to use <'EFIT01'>
+% server = MDS+ database to use: 'DIII-D'(default),'NSTX','EAST',
+% 'THOR', 'OPEN'(assumes mdsconnect already called). Other
+% inputs invoke mdsconnect(server) to connect to server.
+% toupper = 1= all variables made upper case, =-1 all var. made lower case
+% [0]= no change, variables made depending on mds case (typical UC)
+% verbose = set to 1 to get diagnostic prints during execution
+% <1> server= 'NSTX' or 'EAST'; (since can take some time to get)
+% <0> server= 'DIII-D' & all others;
+%
+% OUTPUT:
+% data = structure containing all data in MDS+ tree, with same tree
+% structure except that 'TOP' replaced by 'tree'
+% See: data.allnames for list of all variables in full structure
+% Some Extra items added to structure (all lower case)
+% data.id = sting array of important data identifyer enf
+% data.shot = shot number
+% data.server= MDS+ server
+% data.allnames= list of variables in structure with path relative to "TOP"
+% data.mdsnames= list of variables in structure with full mds path
+% ier = error code
+%
+% WRITTEN BY: Jim Leuer ON 3/1/05 (original name get_mds_tree)
+% taken from get_mds_tree.m uses sub-structure to store
+%
+% USES: eq_mod
+% To see MDS structure on HYDRA run traverser
+% tested on DIII-D and NSTX data and should work for JET data but not tested
+%
+% CHANGE LOG: SMF 20140923 - Changed getnci call to use nid_numbers due to
+% fullpath not accepting wildcards.
+%
+% ==========================================================================
+
+ if nargin==0
+ disp('% get_mds_tree needs at least a "shot" argument')
+ help get_mds_tree
+ return
+ end
+ if nargin < 4
+ toupper=0;
+ end
+ if nargin < 3
+ server='DIII-D';
+ end
+ if nargin < 2
+ tree= 'EFIT01';
+ end
+ if nargin < 5
+ if strcmp(server,'NSTX') | strcmp(server,'EAST')
+ verbose=1;
+ else
+ verbose=0;
+ end
+ end
+ server = upper(server);
+
+ if isempty(toupper) toupper= 0; end
+ if isempty(server) server= 'DIII-D'; end
+ if isempty(tree) tree= 'EFIT01'; end
+ if isempty(verbose)
+ if strcmp(server,'NSTX') | strcmp(server,'EAST')
+ verbose=1;
+ else
+ verbose=0;
+ end
+ end
+
+% -----------------------------------------------
+% Open and check conneciton to MDSPLUS data base:
+ tic
+ ier= 0; status=1;
+ mdsisopen=0;
+% [shoto,status]=mdsopen('atlas.gat.com::EFIT01',shot)
+ if strcmp(server,'DIII-D') | strcmp(server,'DIIID') | strcmp(server,'D3D')
+ status = mdsconnect('atlas.gat.com');
+ elseif strcmp(server,'EAST')
+ status = mdsconnect('202.127.204.12'); % NOT SURE THIS WORKS
+ elseif strcmp(server,'NSTX')
+ status = mdsconnect('skylark.pppl.gov:8501');
+ elseif strcmp(server,'KSTAR')
+ if(strcmp(getenv('HOST'),'datagrid') | strcmp(getenv('HOST'),'ksim2')) % On site at NFRI
+ mds_server = '172.17.100.200:8300';
+ else % Offsite
+% mds_server = '203.230.125.212:8005'; % not able to connect to NFRI server directly
+ mds_server = 'kd';
+ end
+ status = mdsconnect(mds_server);
+ elseif strcmp(server,'THOR')
+ status = mdsconnect('thor');
+ elseif strcmp(server,'VIDAR')
+ status = mdsconnect('vidar');
+ elseif strcmp(server,'OPEN')
+ disp(['get_mds_tree: Assuming MDSCONNECT already called and MDS is ',server])
+ status = 1;
+ else
+ disp(['get_mds_tree: Attempting to connect to server = ',server])
+ status = mdsconnect(server);
+ end
+
+ if ~mod(status,2)
+ ier=1;
+ disp(['ERROR get_mds_tree: unable to connect to ' server])
+ data=[];
+ return;
+ end
+
+% Calling this twice seems to work better for NSTX (who knows why?)
+ [shoto,status]=mdsopen(tree,shot);
+ if(~mod(status,2))
+ [shoto,status]=mdsopen(tree,shot);
+ end
+ if ~mod(status,2)
+ ier=1;
+ disp(['ERROR get_mds_tree: unable to open ' tree ' for shot ' ...
+ int2str(shot)])
+ data=[];
+ status=mdsdisconnect;
+ return;
+ else
+ if(verbose)
+ disp(['% get_mds_tree opened tree, shot: ' tree ' ' int2str(shot)])
+ end
+ end
+
+% add identifier string to structure:
+ data.id= str2mat('mds_efit ', int2str(shot), tree, date);
+ data.shot= shot;
+ data.tree= tree;
+ data.server= server;
+ data.creator= 'get_mds_tree';
+ data.allnames= char([]);
+ data.mdsnames= char([]);
+
+% ===============================================================
+% Get List of tree node names using recursive algorithm & store in mds_all
+% ===============================================================
+ mds_all= char([]); % storage of all mds path names below TOP
+ mds_nam0= ['\' tree '::TOP']; % Top tree name
+ mds_nam= mds_nam0; % starting tree name
+ top_num= length(mds_nam);
+ mds_ot= mds_nam;
+ while ~isempty(mds_ot)
+ mds_all= strvcat(mds_all,mds_ot);
+ mds_ot= mds_sub_tree(mds_ot);
+ end
+
+% ===============================================================
+% Process each name in mds_all for all variables present
+% ===============================================================
+ for kk=1:size(mds_all,1)
+% kk=0; kk=kk+1
+% kk
+ mds_nam = deblank(mds_all(kk,:));
+ mdscmd = ['getnci("\' mds_nam ':*","NID_NUMBER")'];
+ [mds_nids,mstatus] = mdsvalue(mdscmd);
+ if ~mod(mstatus,2) % handle nstx different format
+ mdscmd = ['getnci("\\\' mds_nam ':*","NID_NUMBER")'];
+ [mds_nids,mstatus] = mdsvalue(mdscmd);
+ end
+ varnames = [];
+ if mod(mstatus,2)
+ num_nids = length(mds_nids);
+ varnames = cell(1,num_nids);
+ for j=1:num_nids
+ mdscmd = ['getnci(' num2str(mds_nids(j)) ',"FULLPATH")'];
+ varnames{j} = char(mdsvalue(mdscmd));
+ end
+ else
+% fprintf('WARNING get_mds_tree: Bad mds status for command = %s\n',mdscmd);
+ end
+
+ varnamesc = char(varnames);
+ idgood = strmatch(upper(mds_nam0),upper(varnamesc));
+ varnamesc = varnamesc(idgood,:);
+ data_namesc = varnamesc(:,top_num+2:end);
+ numvar = size(varnamesc,1);
+
+% loop over all variables
+ if numvar>=30
+ tic
+ for ii= 1:numvar
+% ii=0; ii=ii+1;
+ fullnam= deblank(varnamesc(ii,:));
+ data.mdsnames= strvcat(data.mdsnames,fullnam);
+ dat= mdsvalue(fullnam); % Actual value of data
+
+ subnam= deblank(data_namesc(ii,:)); %
+ id= findstr(subnam,':');
+ subnam(id)= '.';
+ if toupper==1
+ subnam= upper(subnam);
+ elseif toupper==-1
+ subnam= lower(subnam);
+ end
+ totnam= deblank(subnam); % now relative address rather than absolute address allow for data name chg
+ if toc > 5 | ii==1
+ if(verbose)
+ disp(['% get_mds_tree Reading ' int2str(ii) '/', int2str(numvar),' variable: ', totnam]);
+ end
+ tic
+ end % if toc
+ data.allnames= strvcat(data.allnames,totnam);
+ str= ['data.' totnam '= dat;'];
+ eval(str);
+ end % for ii
+
+ elseif numvar>=1 % if numvar >=1
+ for ii= 1:numvar
+ fullnam= deblank(varnamesc(ii,:));
+ data.mdsnames= strvcat(data.mdsnames,fullnam);
+ dat= mdsvalue(fullnam); % Actual value of data
+
+ subnam= deblank(data_namesc(ii,:)); %
+ id= findstr(subnam,':');
+ subnam(id)= '.';
+ if toupper==1
+ subnam= upper(subnam);
+ elseif toupper==-1
+ subnam= lower(subnam);
+ end
+ totnam= deblank(subnam);
+ if ii==1
+ if(verbose)
+ disp(['% get_mds_tree Reading ' int2str(ii) '/',int2str(numvar),' variable: ', totnam]);
+ end
+ end % if ii
+ data.allnames= strvcat(data.allnames,totnam);
+ str= ['data.' totnam '= dat;'];
+ eval(str);
+ end % for ii
+ end % if numvar
+ end % for kk
+
+ if strcmp(server,'NSTX') | strcmp(server,'EAST')
+ status=mdsdisconnect; % Exit MDS+ if conneced to remote server
+ end
+
+ return
+
+% ========================================================
+% Testing
+% ========================================================
+
+% Testing SEE test_get_mds_tree
+% (WATCH OUT 114504 has problems use 98549 Ferron High Performance)
+
+% data= get_mds_tree(114504);
+% data= get_mds_tree(98549, 'EFIT01', 'DIII-D');
+% data= get_mds_tree(98549, 'd3d', 'DIII-D');
+
+ data= get_mds_tree(98549, [], [], -1);
+ data= eq_mk_time_last(data); % puts time at end (i.e. (130,1) => (1,130)
+
+ shot=98549; tree='EFIT01'; server='DIII-D'; toupper=-1; % make all lc variables
+% tree='d3d'; tree='IONS'
+ data= get_mds_tree(shot, tree, server, toupper);
+ data= get_mds_tree(98549, [], [], -1);
+ data= get_mds_tree(shot, tree, 'OPEN', toupper);
+
+ filename='/u/leuer/efit/diiid/s98549/g098549.04000' %
+ read_gfile
+
+ shot=98549; tree='IONS'; server='DIII-D'; toupper=-1;
+ ions= get_mds_tree(shot, tree, server, toupper);
+ shot=98549; tree='NB'; server='DIII-D'; toupper=-1;
+ nb= get_mds_tree(shot, tree, server, toupper);
+
+% ========================================================
+% Check NSTX:
+% shot=110843; tree='EFIT01'; server='NSTX';
+
+ data= get_mds_tree(113363, 'EFIT01','NSTX');
+ data= eq_mk_time_last(data); % puts time at end (i.e. (130,1) => (1,130)
+
+ shot=113363; tree='EFIT01'; server='NSTX';
+ clear data
+ data= get_mds_tree(shot, tree, server);
+
diff --git a/matlab/D3D/mdsplus_func/private/getmds.m b/matlab/D3D/mdsplus_func/private/getmds.m
new file mode 100644
index 00000000..78e3cb2e
--- /dev/null
+++ b/matlab/D3D/mdsplus_func/private/getmds.m
@@ -0,0 +1,11 @@
+% MDSplus data manipulation functions.
+%
+% getmds - get data from the MDSplus database
+% get_mds_tree - get entire tree out of mdsplus
+% mdsclose - close an MDSplus tree
+% mdsconnect - connect to MDSplus server
+% mdsdisconnect - disconnect from MDSplus server
+% mdsopen - open an MDSplus tree
+% mdsput - put data in MDSplus database
+% mds_sub_tree -
+% mdsvalue - get data from MDSplus database
diff --git a/matlab/D3D/mdsplus_func/private/make_mdsipmex_east.m b/matlab/D3D/mdsplus_func/private/make_mdsipmex_east.m
new file mode 100644
index 00000000..e0c62231
--- /dev/null
+++ b/matlab/D3D/mdsplus_func/private/make_mdsipmex_east.m
@@ -0,0 +1,28 @@
+% make_mdsipmex.m
+% make sure ../mdsplus points to current version of mdsplus distribution
+% to make the debug version, debug=1;make_mdsipmex;
+% to make the shared version, shared=1;make_mdsipmex;
+% Note; if shared, install MdsIpShr library or point $LD_LIBRARY_PATH to shareables
+% Basil P. DUVAL, May 2000
+
+%MDSPLUS=getenv('MDSPLUS');
+%if(length(MDSPLUS)==0)
+% disp('shell variable MDSPLUS must point to /f/mdsplus/linux/mdsplus before compilation');end
+%MDSPLUS='/f/mdsplus/linux/mdsplus'
+
+MDSPLUS = '/usr/local/mdsplus'
+
+if(exist('debug','var'));DEBUG = '-DDEBUG';else;DEBUG='';end
+
+if(~strcmp(computer,'VMS'));PASSWD = '-DPASSWD';else;PASSWD='mdsipmex.opt';end
+
+if(~exist('shared','var'))
+comm = sprintf('mex -v %s %s mdsipmex.c %s/mdstcpip/mdsipshr.c %s/mdstcpip/mdsiputil.c %s/mdstcpip/mdsip_socket_io.o -DNOCOMPRESSION -I%s/include -I%s/mdstcpip',...
+ DEBUG,PASSWD,MDSPLUS,MDSPLUS,MDSPLUS,MDSPLUS);
+else
+comm = sprintf('mex -v %s %s mdsipmex.c -DNOCOMPRESSION -I%s/include -I%s/mdstcpip -L%s/lib -lMdsIpShr',...
+ DEBUG,PASSWD,MDSPLUS,MDSPLUS,MDSPLUS);
+end
+disp(comm);
+comm
+eval(comm);
diff --git a/matlab/D3D/mdsplus_func/private/mds_eq.m b/matlab/D3D/mdsplus_func/private/mds_eq.m
new file mode 100644
index 00000000..6dccbd10
--- /dev/null
+++ b/matlab/D3D/mdsplus_func/private/mds_eq.m
@@ -0,0 +1,298 @@
+ function eq= mds_eq(shot, tree, gam, mk_var, server)
+%
+% SYNTAX:
+% eq= mds_eq(shot, tree, gam, mk_var); % defaults to DIII-D server
+% mds_eq(shot, 'EFIT01, 'gam', 1, 'NSTX'); % makes vars in workspace
+%
+% PURPOSE: Get EFIT equilibrium (GEQDSK AEQDSK MEQDSK) from mdsplus database.
+%
+% INPUT: <default>
+% shot = shot number
+% tree = which efit tree to use <'EFIT01'>
+% gam = ['g']; get only g file, a= only a file, m= mfile, <gam>
+% mk_var = [1]; Will make efit variables in calling workspace, <0> makes none
+% 2 will make lower case, 3 will make upper case variable
+% note default variables in MDS are upper case
+% NOTE: see eq_to_env for way to make tstart<=t<tend arrays
+% server = which data base to use: defaults to [DIII-D] also NSTX
+%
+% OUTPUT:
+% eq = structure containing EFIT eqdsk variables
+% See: eq.gnames .anames .mnames for g,a,m file names
+% All EFIT variables are as in MDSPLUS => Upper Case
+% Ex: eq.BCENTR, eq.BDRY, eq.PSIRZ,
+% Use: plot(eq.GTIME, eq.CPASMA),
+% contour(eq.R, eq.Z, eq.PSIRZ(:,:,round(length(eq.GTIME)/2))')
+%
+% Extra items added to structure (all lower case)
+% eq.id = sting array of important data identifyer enf
+% eq.gnames= sting array of all gfile names collected (eq to see all)
+% eq.anames= sting array of all afile names collected (eq to see all)
+% eq.mnames= sting array of all mfile names collected (eq to see all)
+% eq.shot = shot number
+% eq.tree = mds tree used to get data (ex. EFIT01 EFIT02 (=MSD));
+%
+% if mk_var>0 it makes variables in base workspace
+% BCENTR, BDRY, CPASMA, PSIRZ, ...
+%
+% NOTE: 1)Function returns all TIME data for shot in large arrays.
+% Example:
+% eq.BCENTER(255,1)
+% eq.GTIME(255,1); % Time data vector for all .GNAME arrays 255 times
+% eq.ATIME(251,1); % Time data vector for all .ANAME arrays 251 times
+% eq.MTIME or eq.TIME % Note old "m" uses TIME new "m" uses MTIME
+% eq.PPRIME(65,255),PSIRZ(65,65,225)
+%
+% Similar arrays (without eq. strucutre) are generated in workspace
+% if mk_var>=1 : BCENTER(255,1), GTIME(255,1) or LC if mk_var==2
+%
+% 2) Time is in vector eq.GTIME and is in ms. All other units are as
+% they come from EFIT. psi(Vs/r) I(A) R(m) ...
+%
+% CAUTION: All time vectors are last index of vector of array.
+% CAUTION: GTIME ATIME & (MTIME or TIME) are not necesaarly the same.
+% CAUTION: Current Bad Variable in d3dMDSPLUS ARE: CASE, HEADER, ZGRID (use Z)
+% CAUTION: THIS ROUTINE GETS ALL MDS EQDSK DATA FOR SHOT AND EQ CAN BE HUGE
+% CAUTION: Afile read 1st then Gfile: variables in common: BCENTR, CPASMA, ...
+%
+% SEE ALSO: eq_time_lim eq_ga_env
+%
+
+% WRITTEN BY: Jim Leuer ON 6/11/03
+% USES: eq_mod
+% To see MDS structure on HYDRA run traverser
+% should work for JET data but not tested
+% ==========================================================================
+%
+% Defaults
+ wait('%CAUTION: mds_eq is OBSOLETE use eq_mds instead')
+ return;
+
+ if nargin==4
+ server='DIII-D';
+ elseif nargin==3
+ server='DIII-D';
+ mk_var= 0;
+ elseif nargin==2
+ server='DIII-D';
+ mk_var= 0;
+ gam= 'gam';
+ elseif nargin==1
+ server='DIII-D';
+ mk_var= 0;
+ gam= 'gam';
+ tree= 'EFIT01';
+ elseif nargin==0
+ disp('% mds_eq needs at least a "shot" argument')
+ help mds_eq
+ return
+ end
+
+ if isempty(server) server= 'DIII-D'; end
+ if isempty(tree) tree= 'EFIT01'; end
+ if isempty(gam) gam= 'gam'; end
+ if isempty(mk_var) mk_var= 0; end
+
+% -----------------------------------------------
+% Open and check conneciton to MDSPLUS data base:
+% -----------------------------------------------
+ tic
+ ier= 0;
+% [shoto,status]=mdsopen('atlas.gat.com::EFIT01',shot)
+ if strcmp(upper(server),'DIII-D')
+ mdsconnect('atlas.gat.com');
+ elseif strcmp(upper(server),'NSTX')
+ mdsconnect('europa.pppl.gov:8501');
+ else
+ disp(['%ERROR: Couldnt recognize parameter server =',server])
+ end
+ [shoto,status]=mdsopen(tree,shot);
+ if ~mod(status,2)
+ ier=1;
+ fprintf(['%ERROR mds_eq: unable to open ' tree ' for shot ' ...
+ int2str(shot) '\n'])
+ eq=[];
+ status=mdsdisconnect;
+ return;
+ end
+
+% add identifier string to structure:
+ id= str2mat('mds_eq ', int2str(shot), tree, date);
+ eq.id= id;
+ eq.shot= shot;
+ eq.tree= tree;
+ eq.creator= 'mds_eq';
+
+% ===============================================================
+% START READ OF AEQDSK DATA
+% ===============================================================
+
+ if ~isempty(strfind(gam,'a'))
+
+ eq.anamesmds= mdsvalue('getnci(".RESULTS.AEQDSK:*","FULLPATH")');
+
+ if isempty(eq.anamesmds)
+ disp('%ERROR mds_eq cant read AEQDSK Names, returning')
+ status=mdsdisconnect;
+ return
+ end
+
+ fullnames= char(eq.anamesmds);
+ id= strfind(fullnames(1,:),'.RESULTS.AEQDSK:');
+ eq.anames= fullnames(:,id+16:end);
+
+% loop over all variables
+ for ii= 1:length(eq.anames);
+% ii=0; ii=ii+1;
+ fullnam= char(eq.anamesmds(ii));
+ nam= deblank(eq.anames(ii,:));
+ dat= mdsvalue(fullnam); %
+ str= ['eq.' nam '= dat;'];
+ eval(str);
+ if mk_var
+ switch mk_var
+ case 1
+ assignin('base',nam,dat);
+ case 2
+ assignin('base',lower(nam),dat);
+ case 3
+ assignin('base',upper(nam),dat);
+ end
+ end
+ if toc > 10
+ disp(['%mds_eq still reading at ' int2str(ii) ' of ',...
+ int2str(length(eq.anames)),' at variable: ', nam]);
+ tic
+ end
+ end
+ end
+
+% ===============================================================
+% START READ OF MEQDSK DATA which is in measurments area time in eq.TIME
+% ===============================================================
+
+ if ~isempty(strfind(gam,'m'))
+
+ eq.mnamesmds= mdsvalue('getnci(".MEASUREMENTS:*","FULLPATH")');
+
+ if ~isempty(eq.mnamesmds)
+ fullnames= char(eq.mnamesmds);
+ id= strfind(fullnames(1,:),'.MEASUREMENTS:');
+ eq.mnames= fullnames(:,id+14:end);
+
+% loop over all variables
+ for ii= 1:length(eq.mnames);
+% ii=0; ii=ii+1;
+ fullnam= char(eq.mnamesmds(ii));
+ nam= deblank(eq.mnames(ii,:));
+ dat= mdsvalue(fullnam); %
+ str= ['eq.' nam '= dat;'];
+ eval(str);
+ if mk_var
+ switch mk_var
+ case 1
+ assignin('base',nam,dat);
+ case 2
+ assignin('base',lower(nam),dat);
+ case 3
+ assignin('base',upper(nam),dat);
+ end
+ end
+ if toc > 10
+ disp(['%mds_eq still reading at ' int2str(ii) ' of ',...
+ int2str(length(eq.mnames)),' at variable: ', nam]);
+ tic
+ end
+ end
+ else
+ disp('%mds_eq COULD NOT FIND ANY EFIT Measurement data: mnamesmds=[]');
+ end
+
+ end
+
+% ===============================================================
+% START READ OF GEQDSK DATA
+% ===============================================================
+
+ if ~isempty(strfind(gam,'g'))
+
+ eq.gnamesmds= mdsvalue('getnci(".RESULTS.GEQDSK:*","FULLPATH")');
+
+ if ~isempty(eq.gnamesmds)
+ fullnames= char(eq.gnamesmds);
+ id= strfind(fullnames(1,:),'.RESULTS.GEQDSK:');
+ eq.gnames= fullnames(:,id+16:end);
+
+% loop over all variables
+ for ii= 1:length(eq.gnames);
+% ii=0; ii=ii+1;
+ fullnam= char(eq.gnamesmds(ii));
+ nam= deblank(eq.gnames(ii,:));
+ dat= mdsvalue(fullnam); %
+ if isfield(eq,nam);
+ rmfield(eq,nam);
+ disp(['%Note: mds_eq is overwriting structure field: ',nam]);
+ end
+ str= ['eq.' nam '= dat;'];
+ eval(str);
+ if mk_var
+ switch mk_var
+ case 1
+ assignin('base',nam,dat);
+ case 2
+ assignin('base',lower(nam),dat);
+ case 3
+ assignin('base',upper(nam),dat);
+ end
+ end
+ if toc > 10
+ disp(['%mds_eq still reading at ' int2str(ii) ' of ',...
+ int2str(length(eq.gnames)),' at variable: ', nam]);
+ tic
+ end
+ end
+ else
+ disp('%mds_eq COULD NOT FIND ANY EFIT GEQDSK Data: gnamesmds=[]');
+ end
+
+ end
+
+% ========================================================
+
+ status=mdsdisconnect;
+
+ eq= eq_mod(eq); % fix time so end of vectors
+
+ return
+
+% Testing SEE test_mds_eq
+% (WATCH OUT 114504 has problems use 98549 Ferron High Performance)
+
+% eqg= mds_eq(114504, 'EFIT01', 'g', 0);
+% eqa= mds_eq(114504, 'EFIT01', 'a', 0);
+% eqm= mds_eq(98549, 'EFIT01', 'm', 0);
+
+ shot=98549; tree='EFIT01'; gam='gam'; mk_var= 0;
+ filename='/u/leuer/efit/diiid/s98549/g098549.04000' %
+ read_gfile
+ eq= mds_eq(shot, tree, 'gam', 0);
+ eq = mds_eq(shot);
+ eq= mds_eq(shot, tree, 'gm', 0 );
+ eq= eq_time_lim(eq,4,4); % make only 4s time
+ eq_ga_env(eq,[],[],'g');
+
+% mds_eq(114504, 'EFIT01',[], 1);
+% mds_eq(114504, 'EFIT02',[], 2);
+
+% eqg= mds_eq(114504, 'EFIT01', 'g', 0);
+% id= round(length(eqg.GTIME)/2);
+% contour(eqg.R, eqg.Z, eqg.PSIRZ(:,:,id)'); hold on, axis equal
+% plot(eqg.RBBBS(1:eqg.NBBBS(id),id),eqg.ZBBBS(1:eqg.NBBBS(id),id),'k')
+
+% ========================================================
+% Check NSTX:
+ shot=110843; tree='EFIT01'; gam='gam'; mk_var= 0; server='NSTX'
+ clear eq
+ eq= mds_eq(shot, tree, 'gam', 0, server);
+ plot(eq.ATIME,eq.IPMHD)
+
diff --git a/matlab/D3D/mdsplus_func/private/mds_sub_tree.m b/matlab/D3D/mdsplus_func/private/mds_sub_tree.m
new file mode 100644
index 00000000..b82ef479
--- /dev/null
+++ b/matlab/D3D/mdsplus_func/private/mds_sub_tree.m
@@ -0,0 +1,70 @@
+function mds_ot = mds_sub_tree(mds_in)
+%
+% SYNTAX:
+% mds_ot= mds_sub_tree(mds_in)
+%
+% PURPOSE:gets mds subtrees for next level down.
+%
+% INPUT: <default>
+% mds_in = Input mds tree
+%
+% OUTPUT:
+% mds_ot = Output mds subtree names
+%
+% WRITTEN BY: Jim Leuer ON 3/3/05
+%
+% CHANGE LOG: SMF 20140923 - Changed getnci call to use nid_numbers due to
+% fullpath not accepting wildcards.
+%
+% ==========================================================================
+% @(#)mds_sub_tree.m 1.3 07/09/09
+
+ % Initialize return variable
+ mds_ot = char([]);
+
+ % For all mds paths provided as input
+ for i = 1:size(mds_in,1)
+
+ % Strip trailing blanks from the mds path
+ mds_inn = deblank(mds_in(i,:));
+
+ % Call GETNCI() to get the NID_NUMBERS
+ mdscmd = ['getnci("\' mds_inn '.*","NID_NUMBER")'];
+ [mds_nids,status] = mdsvalue(mdscmd);
+ if ~status
+ mdscmd = ['getnci("\\\' mds_inn '.*","NID_NUMBER")'];
+ [mds_nids,status] = mdsvalue(mdscmd);
+ end
+
+ % Call GETNCI() to get the FULLPATH for all NIDs
+ if status
+ mds_fpath = [];
+ num_nids = length(mds_nids);
+ mds_fpath = cell(1,num_nids);
+ for j=1:num_nids
+ mdscmd = ['getnci(' num2str(mds_nids(j)) ',"FULLPATH")'];
+ mds_fpath{j} = char(mdsvalue(mdscmd));
+ end
+
+ % ?
+ mds_char = char(mds_fpath);
+ id = strmatch(upper(mds_inn), upper(mds_char));
+ if ~isempty(id)
+ mds_ot = strvcat(mds_ot,mds_char(id,:));
+ end
+ end
+
+ end
+
+ return
+
+% testing
+
+ mds_in = '\EFIT01::TOP';
+ mds_ot = mds_sub_tree(mds_in)
+ mds_in = mds_ot
+ mds_ot = mds_sub_tree(mds_in)
+ mds_in = mds_ot
+ mds_ot = mds_sub_tree(mds_in)
+ mds_in = mds_ot
+
--
GitLab