From 5909d74d5a332acd9eb44fd76a410a164c5cfbb5 Mon Sep 17 00:00:00 2001
From: Olivier Sauter <Olivier.Sauter@epfl.ch>
Date: Thu, 9 Feb 2023 17:36:25 +0100
Subject: [PATCH] start mv mds/efit files from d3d to private
---
matlab/D3D/efit/Contents.m | 50 -
matlab/D3D/efit/cc_efit_to_tok.m | 242 ----
matlab/D3D/efit/efit_movie.m | 59 -
matlab/D3D/efit/eq_mds.m | 81 --
matlab/D3D/efit/eq_mod.m | 113 --
matlab/D3D/efit/equil_to_gdata.m | 197 ---
matlab/D3D/efit/get_efit_data1.m | 274 ----
matlab/D3D/efit/gfile_def.m | 52 -
matlab/D3D/efit/inside_plasma.m | 292 -----
matlab/D3D/efit/make_mhdindat_from_tds.m | 211 ---
matlab/D3D/efit/plasma_current.m | 200 ---
matlab/D3D/efit/read_mhdin.m | 392 ------
matlab/D3D/efit/run_efit.m | 93 --
matlab/D3D/efit/shot_from_gfile.m | 125 --
matlab/D3D/efit/std_efit_units.m | 1238 ------------------
matlab/D3D/efit/trace_boundary.m | 229 ----
matlab/D3D/efit/write_afile.m | 634 ---------
matlab/D3D/efit/write_cc_file.m | 55 -
matlab/D3D/efit/write_eqdsk_efit.m | 65 -
matlab/D3D/efit/write_gfile.m | 320 -----
matlab/D3D/efit/write_kfile.m | 190 ---
matlab/D3D/mdsplus_func/Contents.m | 11 -
matlab/D3D/mdsplus_func/get_mds_tree.m | 291 ----
matlab/D3D/mdsplus_func/getmds.m | 465 -------
matlab/D3D/mdsplus_func/make_mdsipmex.m | 25 -
matlab/D3D/mdsplus_func/make_mdsipmex_east.m | 28 -
matlab/D3D/mdsplus_func/mds_eq.m | 298 -----
matlab/D3D/mdsplus_func/mds_sub_tree.m | 70 -
28 files changed, 6300 deletions(-)
delete mode 100644 matlab/D3D/efit/Contents.m
delete mode 100644 matlab/D3D/efit/cc_efit_to_tok.m
delete mode 100644 matlab/D3D/efit/efit_movie.m
delete mode 100644 matlab/D3D/efit/eq_mds.m
delete mode 100644 matlab/D3D/efit/eq_mod.m
delete mode 100644 matlab/D3D/efit/equil_to_gdata.m
delete mode 100644 matlab/D3D/efit/get_efit_data1.m
delete mode 100644 matlab/D3D/efit/gfile_def.m
delete mode 100644 matlab/D3D/efit/inside_plasma.m
delete mode 100644 matlab/D3D/efit/make_mhdindat_from_tds.m
delete mode 100644 matlab/D3D/efit/plasma_current.m
delete mode 100644 matlab/D3D/efit/read_mhdin.m
delete mode 100644 matlab/D3D/efit/run_efit.m
delete mode 100644 matlab/D3D/efit/shot_from_gfile.m
delete mode 100644 matlab/D3D/efit/std_efit_units.m
delete mode 100644 matlab/D3D/efit/trace_boundary.m
delete mode 100644 matlab/D3D/efit/write_afile.m
delete mode 100644 matlab/D3D/efit/write_cc_file.m
delete mode 100644 matlab/D3D/efit/write_eqdsk_efit.m
delete mode 100644 matlab/D3D/efit/write_gfile.m
delete mode 100644 matlab/D3D/efit/write_kfile.m
delete mode 100644 matlab/D3D/mdsplus_func/Contents.m
delete mode 100644 matlab/D3D/mdsplus_func/get_mds_tree.m
delete mode 100644 matlab/D3D/mdsplus_func/getmds.m
delete mode 100644 matlab/D3D/mdsplus_func/make_mdsipmex.m
delete mode 100644 matlab/D3D/mdsplus_func/make_mdsipmex_east.m
delete mode 100644 matlab/D3D/mdsplus_func/mds_eq.m
delete mode 100644 matlab/D3D/mdsplus_func/mds_sub_tree.m
diff --git a/matlab/D3D/efit/Contents.m b/matlab/D3D/efit/Contents.m
deleted file mode 100644
index 705ea442..00000000
--- a/matlab/D3D/efit/Contents.m
+++ /dev/null
@@ -1,50 +0,0 @@
- %
-% 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/cc_efit_to_tok.m b/matlab/D3D/efit/cc_efit_to_tok.m
deleted file mode 100644
index 66b35a9d..00000000
--- a/matlab/D3D/efit/cc_efit_to_tok.m
+++ /dev/null
@@ -1,242 +0,0 @@
-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/efit_movie.m b/matlab/D3D/efit/efit_movie.m
deleted file mode 100644
index fa710f3b..00000000
--- a/matlab/D3D/efit/efit_movie.m
+++ /dev/null
@@ -1,59 +0,0 @@
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% 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/eq_mds.m b/matlab/D3D/efit/eq_mds.m
deleted file mode 100644
index b587d8bb..00000000
--- a/matlab/D3D/efit/eq_mds.m
+++ /dev/null
@@ -1,81 +0,0 @@
- 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/eq_mod.m b/matlab/D3D/efit/eq_mod.m
deleted file mode 100644
index 2644e28c..00000000
--- a/matlab/D3D/efit/eq_mod.m
+++ /dev/null
@@ -1,113 +0,0 @@
- 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/equil_to_gdata.m b/matlab/D3D/efit/equil_to_gdata.m
deleted file mode 100644
index 50978b71..00000000
--- a/matlab/D3D/efit/equil_to_gdata.m
+++ /dev/null
@@ -1,197 +0,0 @@
- 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/get_efit_data1.m b/matlab/D3D/efit/get_efit_data1.m
deleted file mode 100644
index f9201aa0..00000000
--- a/matlab/D3D/efit/get_efit_data1.m
+++ /dev/null
@@ -1,274 +0,0 @@
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% 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/gfile_def.m b/matlab/D3D/efit/gfile_def.m
deleted file mode 100644
index 9809de13..00000000
--- a/matlab/D3D/efit/gfile_def.m
+++ /dev/null
@@ -1,52 +0,0 @@
- 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/inside_plasma.m b/matlab/D3D/efit/inside_plasma.m
deleted file mode 100644
index 36c20e89..00000000
--- a/matlab/D3D/efit/inside_plasma.m
+++ /dev/null
@@ -1,292 +0,0 @@
- 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/make_mhdindat_from_tds.m b/matlab/D3D/efit/make_mhdindat_from_tds.m
deleted file mode 100644
index 737f069e..00000000
--- a/matlab/D3D/efit/make_mhdindat_from_tds.m
+++ /dev/null
@@ -1,211 +0,0 @@
-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/plasma_current.m b/matlab/D3D/efit/plasma_current.m
deleted file mode 100644
index 4971c41a..00000000
--- a/matlab/D3D/efit/plasma_current.m
+++ /dev/null
@@ -1,200 +0,0 @@
- 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/read_mhdin.m b/matlab/D3D/efit/read_mhdin.m
deleted file mode 100644
index 2096cbd0..00000000
--- a/matlab/D3D/efit/read_mhdin.m
+++ /dev/null
@@ -1,392 +0,0 @@
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% 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/run_efit.m b/matlab/D3D/efit/run_efit.m
deleted file mode 100644
index 47ad68b8..00000000
--- a/matlab/D3D/efit/run_efit.m
+++ /dev/null
@@ -1,93 +0,0 @@
-
- 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/shot_from_gfile.m b/matlab/D3D/efit/shot_from_gfile.m
deleted file mode 100644
index 24bf7df0..00000000
--- a/matlab/D3D/efit/shot_from_gfile.m
+++ /dev/null
@@ -1,125 +0,0 @@
- 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/std_efit_units.m b/matlab/D3D/efit/std_efit_units.m
deleted file mode 100644
index 784f1ce8..00000000
--- a/matlab/D3D/efit/std_efit_units.m
+++ /dev/null
@@ -1,1238 +0,0 @@
-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/trace_boundary.m b/matlab/D3D/efit/trace_boundary.m
deleted file mode 100644
index a2bb4343..00000000
--- a/matlab/D3D/efit/trace_boundary.m
+++ /dev/null
@@ -1,229 +0,0 @@
-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/write_afile.m b/matlab/D3D/efit/write_afile.m
deleted file mode 100644
index 2a4d8031..00000000
--- a/matlab/D3D/efit/write_afile.m
+++ /dev/null
@@ -1,634 +0,0 @@
-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/write_cc_file.m b/matlab/D3D/efit/write_cc_file.m
deleted file mode 100644
index ebdb778b..00000000
--- a/matlab/D3D/efit/write_cc_file.m
+++ /dev/null
@@ -1,55 +0,0 @@
-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/write_eqdsk_efit.m b/matlab/D3D/efit/write_eqdsk_efit.m
deleted file mode 100644
index 6ddb7334..00000000
--- a/matlab/D3D/efit/write_eqdsk_efit.m
+++ /dev/null
@@ -1,65 +0,0 @@
-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/write_gfile.m b/matlab/D3D/efit/write_gfile.m
deleted file mode 100644
index 1c7c9045..00000000
--- a/matlab/D3D/efit/write_gfile.m
+++ /dev/null
@@ -1,320 +0,0 @@
-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/write_kfile.m b/matlab/D3D/efit/write_kfile.m
deleted file mode 100644
index 2cecd1c8..00000000
--- a/matlab/D3D/efit/write_kfile.m
+++ /dev/null
@@ -1,190 +0,0 @@
-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/Contents.m b/matlab/D3D/mdsplus_func/Contents.m
deleted file mode 100644
index 78e3cb2e..00000000
--- a/matlab/D3D/mdsplus_func/Contents.m
+++ /dev/null
@@ -1,11 +0,0 @@
-% 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/get_mds_tree.m b/matlab/D3D/mdsplus_func/get_mds_tree.m
deleted file mode 100644
index 4bf12703..00000000
--- a/matlab/D3D/mdsplus_func/get_mds_tree.m
+++ /dev/null
@@ -1,291 +0,0 @@
-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/getmds.m b/matlab/D3D/mdsplus_func/getmds.m
deleted file mode 100644
index 942451cf..00000000
--- a/matlab/D3D/mdsplus_func/getmds.m
+++ /dev/null
@@ -1,465 +0,0 @@
-function [data,varargout] = getmds(shot,name,range_min,range_max,tree,server,ical);
- %
-% SYNTAX OPTIONS:
-% [data,tvec,ier] = getmds(shot,name,range_min,range_max,tree,server,ical)
-% [data,xvec,tvec,ier] = getmds(shot,name,range_min,range_max,tree,server,ical)
-%
-% PURPOSE: Get mdsplus data from the MDS database.
-%
-% INPUT:
-% shot = shot number
-% name = name of data - must be in single quotes
-% range_min = minimum value of independent variable(s) - use NaN to specify no limit:
-% if scalar, this is min time (s)
-% if vector, last entry is min time (s), other entries are minima
-% for other independent variables
-% range_max = maximum value of independent variable(s) - use NaN to specify no limit:
-% if scalar, this is max time (s)
-% if vector, last entry is max time (s), other entries are maxima
-% for other independent variables
-% tree = which tree to use (optional, default = 'EFIT01')
-% server = which server to use (optional). Options: 'DIII-D'(default),
-% 'JET', 'NSTX', 'THOR'
-% ical = return data in units of digitizer counts (ical=0), volts into
-% digitizer (ical=2), or physics units (ical=1). Only used for
-% data from 'THOR', 'data-server3' (EAST), and 'pcs_kstar' (KSTAR)
-%
-% OUTPUT:
-% data = data for name defined by input (dimensions correspond to dimensions
-% specified in range_min)
-% (xvec..tvec) = variable number of arguments containing vectors of independent
-% variables specified in range_min and range_max
-% ier = error code = 0 if OK, else > 0
-%
-%
-% EXAMPLES
-% getmds(106649,'rmaxis',1.5,5,'EFIT01','DIII-D')
-% getmds(57035,'MAGN/IPLA',1.0,3.0,'PPF','JET')
-% getmds(109070,'ip',1.0,3.0,'wf','NSTX')
-% getmds(999991,'PCIP',0.0,8.0,'PCS_KSTAR','THOR');
-%
-% NEED 3D examples
-%
-% RESTRICTIONS:
-% (1) Only works for range_min/max up to length 2 at the moment (make a request for more).
-% (2) Use the lower level routines starting with "mds" to manipulate the MDS
-% data directly.
-
-% METHOD:
-%
-% WRITTEN BY: Mike Walker ON 6/4/01
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% @(#)getmds.m 1.12 08/08/13
-
-% Error checking:
-
-if (nargin < 4)
- disp('ERROR getmds: Must have at least 4 input arguments')
- return;
-end;
-
-ndim = length(range_min);
-if(ndim~=length(range_max))
- disp('ERROR getmds: range_min and range_max must be equal length')
- return;
-end
-
-% test for existence of optional arguments
-
-if (nargin < 5)
- tree = 'EFIT01';
-end;
-if (nargin < 6)
- server = 'DIII-D';
-end;
-
-if(ndim>1)
- xmin = range_min(1); tmin = range_min(2);
- xmax = range_max(1); tmax = range_max(2);
-else
- tmin = range_min; tmax = range_max;
-end
-
-if(isempty(tmin))
- tmin=-1e+6;
-end
-if(isempty(tmax))
- tmax=1e+6;
-end
-if(nargin < 7)
- ical=1;
-end
-data=[];tvec=[];ier=0;
-if(ical~=0 & ical~=1 & ical~=2)
- fprintf('ERROR getmds: invalid ical = %d\n',ical);
- ier=1; return;
-end
-
-upcase_server = upper(server);
-
-% process default DIII-D server case
-if ((nargin < 6) | ((nargin >= 6) & (strcmp(upcase_server,'DIII-D'))))
-
- status = mdsconnect('atlas.gat.com'); % returns status=1 if success, else -1
- if status ~= 1
- disp(['ERROR getmds: unable to open server ' server])
- ier=1;
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- ier=0;
- [shoto,status] = mdsopen(tree,shot); % returns status=1 on success, else 0
- if ~mod(status,2)
- ier=1;
- fprintf(['ERROR getmds: unable to open ' tree ' for shot ' ...
- int2str(shot) '\n'])
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get data vector:
- dname = ['\' name];
- [data,status] = mdsvalue(dname);
- if ~mod(status,2)
- ier=2;
- fprintf('ERROR getmds: unable to get data %s from %s\n',dname,tree)
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get time vector:
- if(ndim==1)
- tname = ['dim_of(\' name ')'];
- else
- tname = ['dim_of(\' name ',1)'];
- end
- [tvec,status] = mdsvalue(tname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get time data from ' tree '\n'])
- tvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
- if(min(tvec)>10 | max(tvec)>100) % time units must be milliseconds
- tvec = tvec*1e-3; % convert to seconds
- end
-
- % get x vector:
- if(ndim>1)
- xname = ['dim_of(\' name ',0)'];
- [xvec,status] = mdsvalue(xname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get x data from ' tree '\n'])
- xvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,tvec,[]);
- return;
- end
- end
-
-end;
-
-% process non DIII-D server case
-if ((nargin >= 6) & (~strcmp(upcase_server,'DIII-D')))
-
- if (strcmp(upcase_server,'JET'))
- status = mdsconnect('mdsplus.jet.efda.org');
- if status ~= 1
- disp(['ERROR getmds: unable to open server ' server])
- ier=1;
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
- cmdString = '_s = JET("';
- cmdString = strcat('_s = JET("',tree,'/',name,'",',int2str(shot),')');
- data = mdsvalue(cmdString);
- tvec = mdsvalue('DIM_OF(_s)');
- elseif (strcmp(upcase_server,'NSTX'))
- status = mdsconnect('birch.pppl.gov:8501');
- if status ~= 1
- disp(['ERROR getmds: unable to open server ' server])
- ier=1;
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- ier=0;
- [shoto,status] = mdsopen(tree,shot);
- if ~mod(status,2)
- [shoto,status] = mdsopen(tree,shot); % seems to work better to do it twice
- end
- if ~mod(status,2)
- ier=1;
- fprintf(['ERROR getmds: unable to open ' tree ' for shot ' ...
- int2str(shot) '\n'])
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get data vector:
- dname0=name;
- if(name(1:1)=='.')
- dname = name; % relative reference
- elseif(strcmp(name(1:1),'\'))
- dname = name;
- else
- dname = ['\' name];
- end
- [data,status] = mdsvalue(dname);
- if ~mod(status,2)
- ier=2;
- fprintf('ERROR getmds: unable to get data %s from %s\n',...
- deblank(dname),tree);
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get time vector:
- if(ndim==1)
- tname = ['dim_of(' dname ')'];
- else
- tname = ['dim_of(\' name ',1)'];
- end
- [tvec,status] = mdsvalue(tname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get time data from ' tree '\n'])
- tvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- if(ndim>1) % get x vector
- xname = ['dim_of(\' name ',0)'];
- [xvec,status] = mdsvalue(xname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get x data from ' tree '\n'])
- xvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,tvec,[]);
- return;
- end
- end
-
- elseif (strcmp(upcase_server,'THOR') | ...
- strcmp(upcase_server,'202.127.205.8') | ...
- strcmp(upcase_server,'203.230.125.212:8005') | ...
- (strcmp(upcase_server,'202.127.205.59') & strcmp(upper(tree),'PCS_EAST')))
-
- status = mdsconnect(upcase_server);
- if status ~= 1
- disp(['ERROR getmds: unable to open server ' server])
- ier=1;
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- ier=0;
- [shoto,status] = mdsopen(tree,shot);
- if ~mod(status,2)
- ier=1;
- fprintf(['ERROR getmds: unable to open ' tree ' for shot ' ...
- int2str(shot) '\n'])
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- if(strcmp(upper(tree),'PCS_EAST') | ... % EAST PCS tree
- strcmp(upper(tree),'PCS_KSTAR') | ... % KSTAR PCS tree
- strcmp(upper(tree),'RDATA')) % KSTAR non-PCS tree
- dname = ['\' name];
- else
- dname = name;
- end
-
- % get data vector:
- dname0 = upper(dname);
- if(ical==1)
- dname = upper(dname);
- elseif(ical==2)
- dname = [upper(dname) ':RAW'];
- elseif(ical==0)
- dname = [upper(dname) ':RAW'];
- end
- [data,status] = mdsvalue(dname);
- if ~mod(status,2)
- ier=2;
- fprintf(['ERROR getmds: unable to get data from ' tree '\n'])
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get time vector:
- if(ndim==1)
- tname = ['dim_of(' dname0 ')'];
- else
- tname = ['dim_of(' dname0 ',1)'];
- end
- [tvec,status] = mdsvalue(tname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get time data from ' tree '\n'])
- tvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- if(ndim>1) % get x vector
- xname = ['dim_of(\' name ',0)'];
- [xvec,status] = mdsvalue(xname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get x data from ' tree '\n'])
- xvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,tvec,[]);
- return;
- end
- end
-
- if(ical==2)
- attr_name = [upper(name) ':OFFSET'];
- [offset,status] = mdsvalue(attr_name);
-
- attr_name = [upper(name) ':FBITS'];
- [counts_to_volts,status] = mdsvalue(attr_name);
-
- data = (data+offset) * counts_to_volts;
-
- attr_name = [upper(name) ':SCALE'];
- [scale,status] = mdsvalue(attr_name);
-
- end
-
-
- else % unknown server
- status = mdsconnect(server);
- if status ~= 1
- disp(['ERROR getmds: unable to open server ' server])
- ier=1;
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- ier=0;
- [shoto,status] = mdsopen(tree,shot);
- if ~mod(status,2)
- ier=1;
- fprintf(['ERROR getmds: unable to open ' tree ' for shot ' ...
- int2str(shot) '\n'])
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get data vector:
- dname = ['\' name];
- [data,status] = mdsvalue(dname);
- if ~mod(status,2)
- ier=2;
- fprintf(['ERROR getmds: unable to get data from ' tree '\n'])
- disp(data)
- data = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- % get time vector:
- if(ndim==1)
- tname = ['dim_of(\' name ')'];
- else
- tname = ['dim_of(\' name ',1)'];
- end
- [tvec,status] = mdsvalue(tname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get time data from ' tree '\n'])
- tvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,[],[]);
- return;
- end
-
- if(ndim>1) % get x vector
- xname = ['dim_of(\' name ',0)'];
- [xvec,status] = mdsvalue(xname);
- if ~mod(status,2)
- ier=3;
- fprintf(['ERROR getmds: unable to get x data from ' tree '\n'])
- xvec = [];
- varargout = getmds_varargout(nargout,ndim,ier,tvec,[]);
- return;
- end
- end
-
-
- end;
-end
-
-% return only the data requested:
-if(isempty(tvec) | tvec(end) < tmin | tvec(1) > tmax)
- ier = 4;
- data=[]; xvec=[]; tvec=[];
- fprintf(['ERROR getmds: no data between tmin and tmax \n'])
-elseif((ndim>1) & (isempty(xvec) | max(xvec) < xmin | min(xvec) > xmax))
- ier = 4;
- data=[]; xvec=[]; tvec=[];
- fprintf(['ERROR getmds: no data between xmin and xmax \n'])
-elseif((ndim>1) & (isnan(xmin) || isnan(xmax))) % Don't truncate x-series
- [mm,m1]=min(abs(tvec-tmin));
- [mm,m2]=min(abs(tvec-tmax));
- tvec = tvec(m1:m2);
- data = data(:,m1:m2);
-elseif (isnan(tmin) || isnan(tmax)) % Don't truncate t-series
- if(ndim>1)
- [nn,n1]=min(abs(xvec-xmin));
- [nn,n2]=min(abs(xvec-xmax));
- xvec = xvec(n1:n2);
- data = data(n1:n2,:);
- end
-else
- [mm,m1]=min(abs(tvec-tmin));
- [mm,m2]=min(abs(tvec-tmax));
- tvec = tvec(m1:m2);
- if(ndim>1)
- [nn,n1]=min(abs(xvec-xmin));
- [nn,n2]=min(abs(xvec-xmax));
- xvec = xvec(n1:n2);
- data = data(n1:n2,m1:m2);
- else
- data = data(m1:m2);
- end
-end
-
-if(ndim>1)
- varargout = getmds_varargout(nargout,ndim,ier,tvec,xvec);
-else
- varargout = getmds_varargout(nargout,ndim,ier,tvec);
-end
-
-mdsdisconnect;
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-function out = getmds_varargout(nargs,ndim,ier,tvec,xvec)
-
-out = cell(1);
-if(ndim>1)
- out{1} = xvec;
- out{2} = tvec;
-else
- out{1} = tvec;
-end
-if(nargs > ndim+1)
- out{nargs-1} = ier;
-end
diff --git a/matlab/D3D/mdsplus_func/make_mdsipmex.m b/matlab/D3D/mdsplus_func/make_mdsipmex.m
deleted file mode 100644
index 7e1fdb20..00000000
--- a/matlab/D3D/mdsplus_func/make_mdsipmex.m
+++ /dev/null
@@ -1,25 +0,0 @@
-% 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
-
-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/make_mdsipmex_east.m b/matlab/D3D/mdsplus_func/make_mdsipmex_east.m
deleted file mode 100644
index e0c62231..00000000
--- a/matlab/D3D/mdsplus_func/make_mdsipmex_east.m
+++ /dev/null
@@ -1,28 +0,0 @@
-% 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/mds_eq.m b/matlab/D3D/mdsplus_func/mds_eq.m
deleted file mode 100644
index 6dccbd10..00000000
--- a/matlab/D3D/mdsplus_func/mds_eq.m
+++ /dev/null
@@ -1,298 +0,0 @@
- 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/mds_sub_tree.m b/matlab/D3D/mdsplus_func/mds_sub_tree.m
deleted file mode 100644
index b82ef479..00000000
--- a/matlab/D3D/mdsplus_func/mds_sub_tree.m
+++ /dev/null
@@ -1,70 +0,0 @@
-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