Home

This is the wiki for the matlabVMEC group of subroutines. These subroutines help the user read, process and import data from a common set of STELLARATOR codes, the principal one being VMEC. Below you can find the basic instructions for setting up and using these tools. As always, the user may type help <subroutine> in Matlab for a full list.

---

%%  Installation Instructions
%  This file will walk you through using the matlabVMEC routines.  In 
%  order to effectively use these routines they'll need to be added to 
%  your 'path'.  
%
%  1) Make sure you cloned the matlabVMEC repo to a permanent location.
%  2) The easiest way to setup access to these files is to use the 
%     addpath command and place the command in your startup.m file.
%     you can find this file by the 'which startup.m' command.
%  3) Add each sudirectory like so:
addpath('/Users/username/src/matlabVMEC/VMEC');
%
%  Now the matlabVMEC routines can be called from any working directory.
%
%  NETCDF Note:
%  If the user wishes to read netCDF files (*.nc) they will need to
%  download the java class file from the NCAR website.  Please note you'll
%  need the netcdfALL jar package which can now be found here:
%       https://downloads.unidata.ucar.edu/netcdf-java/
%  This can be added to your startup.m file with a line like:
javaaddpath('/Users/lazerson/username/src/java/netcdfAll-5.5.2.jar');


%%  VMECplot
%   The VMECplot routine provides users with the ability to interactively
%   plot the contents of VMEC 'wout' files.  There are a few way to envoke
%   this routine but the simplest is simply to go to a directory containing
%   'wout' files and invoke it
    VMECplot
%   It can take a few seconds to load up.  The code reads the wout files in
%   the current working directory, loads the first one it finds, converts
%   the Fourier space arrays to real space, and then you get the interface.
%   The conversion to real space can take a few minutes on slower machines
%   so have patience.

%%  Using VMECplot
%   The VMECplot interface should be fairly straightforward.  There is a
%   menu which lists all the 'wout' files in the current working directory.
%   A menu for the various parameters that can be plotted.  Buttons which
%   control what type of plots to make.  A textfield for outputting plots
%   to a file (format is determined by file suffix).  Sliders for
%   controlling plots.  And two textboxes which control the number of
%   gridpoints used in theta and phi for plotting.  Status messages are
%   displayed in the lower left corner.  Note that some plots can only be
%   displayed in 1D or 2D or 3D.  Those plots will turn off sliders and
%   plotting options as they are selected.

%%  VMECedit
%   The VMECedit routine is an GUI which constructs a VMEC input namelist
%   file.  It can also be used to modify existing VMEC input namelist
%   files.  Its main usefulness is it's ability to dynamically plot the
%   boundary flux surface as the user enters the RBC and ZBS values.  It
%   was built around a VMEC v6.9 so it's still a beta but has been added to
%   the matlabVMEC package for users to provide feedback.
    VMECedit
    
%%  Using matlabVMEC (nuts and bolts)
%   The matlabVMEC package contains various functions for reading in data
%   files, doing transformations, and plotting results.  All routines have
%   help test associated with them so you can alway type: help read_vmec
%   for example.  I'll now break down some basic operations the user may
%   routinely want to use to make more complex plots.

%%  Reading a wout file
%   A dummy variable will be used here for the file name.  In general this
%   can be passed as a string directly to read_vmec like:
%   read_vmec('wout.test')
    filename='';
    data=read_vmec(filename);
%   Some notes should be made regarding field names.  All fields for
%   quantities in Fourier space end with a 'c' or 's' to denote that they
%   are even or odd (require cosine or sine transforms).  Older version of
%   VMEC supported non-axisymmetric runs but didn't necessarily output all
%   variables, so support for non-axisymmetric runs isn't really there.
%   Some quantities are only in 1D and an attempt to map everything to full
%   VMEC grid (from the half-grid) has been made.
    
%%  Getting R and Z in real space (axisymmetric)
%   Here we define the points in space (theta, phi) at which we wish to
%   know the R and Z coordinates.  A note should be make regarding the
%   choice of phi.  These functions differ from previous versions as we now
%   use the mn by ns arrays which are vectorized over the Fourier modes.
%   Doing this provides the code a substantial speedup (the old functions
%   have been renamed sfunct_old and cfunct_old for those who still wish to
%   use them)
    theta=0:2*pi/(359):2*pi;    % 36 points define theta [0,6.28]
    phi=0:2*pi/(359):2*pi;     % 361 points define phi [0,6.28]
    r=cfunct(theta,phi,data.rmnc,data.xm,data.xn);  % Cos transform
    z=sfunct(theta,phi,data.zmns,data.xm,data.xn);  % Sin transform
%   Note we have not taken advantage of the VMEC representation of the
%   poloidal coordinate (lambda).  An example can be seen later in this
%   tutorial file.
 
%%  2D flux surface plots
%   Here we only note that r and z are indexed (s,theta,phi)
    phidex=1;   % Index at which to plot the flux surfaces
    plot(r(data.ns,:,phidex),z(data.ns,:,phidex),'r');  % Boundary
    hold on
    for i=2:data.ns-1
        plot(r(i,:,phidex),z(i,:,phidex),'k');          % Inner flux
    end
    plot(r(1,1,phidex),z(1,1,phidex),'+k');             % Magnetic Axis
    hold off
    xlabel('R [m]');
    ylabel('Z [m]');
    title(['Flux Surfaces at phi=' num2str(phi(phidex),'%4.2f')]);
    
%%  2D Contour plots (torocont)
%   Here we'll plot |B| (modB).  First we'll map it into real space
    modb=cfunct(theta,phi,data.bmnc,data.xm_nyq,data.xn_nyq);    % Cos Transform
    torocont(r,z,modb,phidex);
    xlabel('R [m]');
    ylabel('Z [m]');
    title(['|B| at phi=' num2str(phi(phidex),'%4.2f')]);

%%  3D Flux surface plots (isotoro)
    isotoro(r,z,phi,data.ns);  % Just plot a flux surface
    colorbar on;
    xlabel('x [m]');
    ylabel('y [m]');
    zlabel('z [m]');
    title(['Flux Surface (ns=' num2str(data.ns,'%3d') ')']);
%   We can plot quantities on the surface using (modb for example)
    isotoro(r,z,phi,data.ns,modb);
%   We can plot segments
    isotoro(r,z,phi(0:180),data.ns);
%   We can also plot multiple surfaces
    isotoro(r,z,phi(0:180),1:data.ns);
    
%%  Vectors
%   VMEC currently (as of 8.47) outputs the fourier coefficients for the
%   covariant (b_s,b_u,b_v) and contravariant (b^u,b^v, b^s=0) components
%   of the magnetic field.  It is a straightforward matter to transform the
%   contravariant components to (bu and bv structure elements) cylindrical
%   components.  In general we recognize:
%              dR           dR
%   B_R= B^u ------ + B^v ------
%              du           dv
%
%   B_phi= B^u R
%
%              dZ           dZ
%   B_Z= B^u ------ + B^v ------
%              du           dv
%
%   The derivative terms are simply the products between the R and Z
%   coefficients and m,n, and nfp.  This is already handled by the
%   read_data function.  So to create the vectors:
    data=read_vmec(filename);
    theta=0:2*pi/(35):2*pi;    % 36 points define theta [0,6.28]
    phi=0:2*pi/(360):2*pi;     % 361 points define phi [0,6.28]
    r=cfunct(theta,phi,data.rmnc,data.xm,data.xn);  % Cos transform
    z=sfunct(theta,phi,data.zmns,data.xm,data.xn);  % Sin transform
    drdu=sfunct(theta,phi,data.rumns,data.xm,data.xn);
    drdv=sfunct(theta,phi,data.rvmns,data.xm,data.xn);
    dzdu=cfunct(theta,phi,data.zumnc,data.xm,data.xn);
    dzdv=cfunct(theta,phi,data.zvmnc,data.xm,data.xn);
    bu=cfunct(theta,phi,data.bsupumnc,data.xm_nyq,data.xn_nyq);
    bv=cfunct(theta,phi,data.bsupvmnc,data.xm_nyq,data.xn_nyq);
    brho=bu.*drdu+bv.*drdv;    % Make B_R
    bphi=bv.*r;                % Make B_phi
    bz=bu.*dzdu+bv.*dzdv;      % Make B_Z
    phidex=1;
    torocont(r,z,bphi,phidex);
    hold on
    quiver(r(:,:,phidex),z(:,:,phidex),brho(:,:,phidex),bz(:,:,phidex));
    hold off
    xlabel('R [m]');
    ylabel('Z [m]');
    title(['Magnetic Field at phi=' num2str(phi(phidex),'%4.2f')]);
%   It is worthwhile to note that for non-axisymmertic configurations the
%   vectors will appear to have a normal component.  This is due to the
%   twisting of the surface in 3D.  In essence the normal to the surface
%   does not lie in the (R-Z) plane.  
    
%%  EXTRA Routines
%   There are other routines and the user is encourage to use help and look
%   at their source to use them (in general they are all in beta):
%
%   read_vmec_input:    Used by VMECedit to read the VMEC input namelist
%   toroslice:          Automates the plotting 2D flux surface
%   read_vessel:        Reads a piecewise description of the vessel
%                       boundary.  VMECplot will plot this vessel if the
%                       returned data structure is passed to VMECplot.
%                       (eg. VMECplot(ves_data) )
%   write_vmec:         Used by VMECedit to write the VMEC input namelist
%   
%   VMECedit makes use of the fortran_namelist matlab package and the
%   read_netcdf package (which are included in the extras subdirectories).