function output = VMECcomp( filename,datatype )
%VMECcomp(filename,datatype) Creates plots comparing VMEC equilibria.
% The VMECcomp function plots allows the user to create comparrision
% plots of specific quantities between VMEC equilibria. This is
% accomplished by passing either a string with wildcards as the first
% parameter (or a cell array of file names) and a string inicating the
% type of value to plot as the second parameter. Available options are:
% 'iota' Rotational Transform profile
% 'pressure' Pressure profile
% 'current' Toroidal current profile
% 'flux' Flux surfaces at phi=0
% 'fluxpi2' Flux surfaces at quarter field period
% 'fluxpi' Flux surfaces at half field period
% 'magaxis' 3D plot of magnetic axis trajectory
%
% Example usage
% haxis=VMEComp('wout*','iota'); % Comparrision plot of iota
%
% Maintained by: Samuel Lazerson (lazerson@pppl.gov)
% Version: 1.0
if ~iscell(filename)
file_struct=dir(filename);
nfiles=length(file_struct);
filename=cell(1,nfiles);
vmec_data=cell(1,nfiles);
for i=1:nfiles
filename{i}=file_struct(i).name;
vmec_data{i}=read_vmec(filename{i});
end
else
nfiles=max(size(filename));
vmec_data=cell(1,nfiles);
for i=1:nfiles
vmec_data{i}=read_vmec(filename{i});
end
end
switch datatype
case 'iota'
hold on
for i=1:nfiles
ns=vmec_data{i}.ns;
iota=vmec_data{i}.iotaf;
if (i == 1)
plot3(0:1/(ns-1):1,i.*ones(1,ns),iota,'b','LineWidth',2.0);
elseif (i == nfiles)
plot3(0:1/(ns-1):1,i.*ones(1,ns),iota,'r','LineWidth',2.0);
else
plot3(0:1/(ns-1):1,i.*ones(1,ns),iota,'k');
end
end
xlabel('Normalized Flux');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Rotational Transform');
view(3);
output=gca;
case 'pressure'
hold on
for i=1:nfiles
ns=vmec_data{i}.ns;
presf=vmec_data{i}.presf;
if (i == 1)
plot3(0:1/(ns-1):1,i.*ones(1,ns),presf,'b','LineWidth',2.0);
elseif (i == nfiles)
plot3(0:1/(ns-1):1,i.*ones(1,ns),presf,'r','LineWidth',2.0);
else
plot3(0:1/(ns-1):1,i.*ones(1,ns),presf,'k');
end
end
xlabel('Normalized Flux');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Pressure');
view(3);
output=gca;
case 'current'
hold on
for i=1:nfiles
ns=vmec_data{i}.ns;
jcurv=vmec_data{i}.jcurv;
if (i == 1)
plot3(0:1/(ns-1):1,i.*ones(1,ns),jcurv,'b','LineWidth',2.0);
elseif (i == nfiles)
plot3(0:1/(ns-1):1,i.*ones(1,ns),jcurv,'r','LineWidth',2.0);
else
plot3(0:1/(ns-1):1,i.*ones(1,ns),jcurv,'k');
end
end
xlabel('Normalized Flux');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Toroidal Current');
view(3);
output=gca;
case {'flux0','flux'}
ntheta=90;
hold on
for i=1:nfiles
ns=vmec_data{i}.ns;
rmnc=vmec_data{i}.rmnc;
zmns=vmec_data{i}.zmns;
xm=vmec_data{i}.xm;
xn=vmec_data{i}.xn;
r=cfunct(0:2*pi/(ntheta-1):2*pi,0,rmnc,xm,xn);
z=sfunct(0:2*pi/(ntheta-1):2*pi,0,zmns,xm,xn);
for j=10:10:ns
plot3(r(j,:),i.*ones(1,ntheta),z(j,:),'k');
end
plot3(r(1,1),i,z(1,1),'+k');
plot3(r(ns,:),i.*ones(1,ntheta),z(ns,:),'k');
end
xlabel('R [m]');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Z [m]');
view(3);
axis equal
output=gca;
case 'fluxpi2'
ntheta=90;
hold on
for i=1:nfiles
nfp=vmec_data{i}.nfp;
ns=vmec_data{i}.ns;
rmnc=vmec_data{i}.rmnc;
zmns=vmec_data{i}.zmns;
xm=vmec_data{i}.xm;
xn=vmec_data{i}.xn;
r=cfunct(0:2*pi/(ntheta-1):2*pi,pi/2/nfp,rmnc,xm,xn);
z=sfunct(0:2*pi/(ntheta-1):2*pi,pi/2/nfp,zmns,xm,xn);
for j=10:10:ns
plot3(r(j,:),i.*ones(1,ntheta),z(j,:),'k');
end
plot3(r(1,1),i,z(1,1),'+k');
plot3(r(ns,:),i.*ones(1,ntheta),z(ns,:),'k');
end
xlabel('R [m]');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Z [m]');
view(3);
axis equal
output=gca;
case 'fluxpi'
ntheta=90;
hold on
for i=1:nfiles
nfp=vmec_data{i}.nfp;
ns=vmec_data{i}.ns;
rmnc=vmec_data{i}.rmnc;
zmns=vmec_data{i}.zmns;
xm=vmec_data{i}.xm;
xn=vmec_data{i}.xn;
r=cfunct(0:2*pi/(ntheta-1):2*pi,pi/nfp,rmnc,xm,xn);
z=sfunct(0:2*pi/(ntheta-1):2*pi,pi/nfp,zmns,xm,xn);
for j=10:10:ns
plot3(r(j,:),i.*ones(1,ntheta),z(j,:),'k');
end
plot3(r(1,1),i,z(1,1),'+k');
plot3(r(ns,:),i.*ones(1,ntheta),z(ns,:),'k');
end
xlabel('R [m]');
set(gca,'YTick',1:nfiles);
set(gca,'YTickLabel',filename);
zlabel('Z [m]');
view(3);
axis equal
output=gca;
case 'magaxis'
nzeta=360;
cosph=cos(0:2*pi/(nzeta-1):2*pi);
sinph=sin(0:2*pi/(nzeta-1):2*pi);
hold on
for i=1:nfiles
nfp=vmec_data{i}.nfp;
ns=vmec_data{i}.ns;
rmnc=vmec_data{i}.rmnc;
zmns=vmec_data{i}.zmns;
xm=vmec_data{i}.xm;
xn=vmec_data{i}.xn;
r=cfunct(0,0:2*pi/(nzeta-1):2*pi,rmnc,xm,xn);
z=sfunct(0,0:2*pi/(nzeta-1):2*pi,zmns,xm,xn);
plot3(squeeze(r(1,1,:)).*cosph',squeeze(r(1,1,:)).*sinph',squeeze(z(1,1,:)),'k');
end
xlabel('X [m]');
ylabel('Y [m]');
zlabel('Z [m]');
view(3);
axis equal
output=gca;
otherwise
disp(['Error: Datatype ' strtrim(datatype) ' is not supported']);
output = -1;
end
return
end
