function plot_extender(varargin)
%PLOT_EXTENDER(mesh,field) Plots the extender fields on a grid
% PLOT_EXTENDER(mesh,field) Produces a plot of the extender fields on the
% phi=0 cutplane. The mesh and field variables are the structure
% returned by READ_EXTENDER after reading the extended_mesh and
% extended_field files respectively.
%
% PLOT_EXTENDER(mesh,field,'phi',phival) Controls the cutplane which is
% plotted by phival (should be 0 to pi/nfp).
%
% PLOT_EXTENDER(mesh,field,'MODB') Produces a plot of |B| on a given
% phi cut plane (default = 0.0).
%
% PLOT_EXTENDER(mesh,field,'BR') Produces a plot of |B| on a given
% phi cut plane (default = 0.0).
%
% PLOT_EXTENDER(mesh,field,'BPHI') Produces a plot of |B| on a given
% phi cut plane (default = 0.0).
%
% PLOT_EXTENDER(mesh,field,'BZ') Produces a plot of |B| on a given
% phi cut plane (default = 0.0).
%
% PLOT_EXTENDER(mesh,field,'2DVECTOR') Produces a plot of B_phi with a
% quiver plot (arrows) of the R and Z vectorfield, on a given phi cut
% plane (default = 0.0).
%
% PLOT_EXTENDER(mesh,field,'GRID') Produces a 3D cutplane view of the
% geometry.
%
% Example:
% extended_mesh=read_extender('extended_mesh.test');
% extended_field=read_extender('extended_field.test');
% plot_extender(extended_mesh,extended_field);
%
% Written by: S.Lazerson (lazerson@pppl.gov)
% Version: 1.0
% Date: 5/11/11
% This code utilzes the PIXPLOT routine available at
% http://www.mathworks.com/matlabcentral/fileexchange/28714-pixplot
% Defaults
phi=0.0;
plottype='2D';
mesh=[];
field=[];
% Handle varargin
i=1;
while i <= nargin
if isstruct(varargin{i})
if isfield(varargin{i},'datatype')
switch varargin{i}.datatype
case 'EXTENDER_MESH'
mesh=varargin{i};
case 'EXTENDER_FIELD'
field=varargin{i};
end
end
elseif ischar(varargin{i})
switch upper(varargin{i})
case '2D'
plottype=varargin{i};
case 'BR'
plottype=varargin{i};
case 'BPHI'
plottype=varargin{i};
case 'BZ'
plottype=varargin{i};
case 'MODB'
plottype=varargin{i};
case '2DVECTOR'
plottype=varargin{i};
case 'GRID'
plottype=varargin{i};
case 'PHI'
i=i+1;
phi=varargin{i};
end
end
i=i+1;
end
% Check for mesh and field
if isempty(mesh) && isempty(field)
disp('ERROR: Must supply mesh and field structures');
return
end
% First get the field data on the full grid
br=zeros(1,mesh.npoints);
bphi=zeros(1,mesh.npoints);
bz=zeros(1,mesh.npoints);
br(mesh.point_address > -1)=field.br;
bphi(mesh.point_address > -1)=field.bphi;
bz(mesh.point_address > -1)=field.bz;
br=reshape(br,[field.nr field.nphi field.nz]);
bphi=reshape(bphi,[field.nr field.nphi field.nz]);
bz=reshape(bz,[field.nr field.nphi field.nz]);
% Find index of phi
nphi=1;
for i=1:mesh.nphi
if (mesh.phi(nphi) <= phi)
nphi=i;
end
end
% Now make plot
switch upper(plottype)
case 'BR'
pixplot(mesh.raxis,mesh.zaxis,squeeze(br(:,nphi,:)));
axis equal; axis tight; colorbar;
title('Radial Magnetic Field (B_R)');
xlabel('R [m]');
ylabel('Z [m]');
case 'BPHI'
pixplot(mesh.raxis,mesh.zaxis,squeeze(bphi(:,nphi,:)));
axis equal; axis tight; colorbar;
title('Toridal Magnetic Field (B_\phi)');
xlabel('R [m]');
ylabel('Z [m]');
case 'BZ'
pixplot(mesh.raxis,mesh.zaxis,squeeze(bz(:,nphi,:)));
axis equal; axis tight; colorbar;
title('Vertical Magnetic Field (B_Z)');
xlabel('R [m]');
ylabel('Z [m]');
case 'MODB'
modb=sqrt(br.*br+bphi.*bphi+bz.*bz);
pixplot(mesh.raxis,mesh.zaxis,squeeze(modb(:,nphi,:)));
axis equal; axis tight; colorbar;
title('Magnetic Field Strength (|B|)');
xlabel('R [m]');
ylabel('Z [m]');
case '2D'
subplot(2,2,1);
pixplot(mesh.raxis,mesh.zaxis,squeeze(br(:,nphi,:)));
axis equal; axis tight;
title('Radial Magnetic Field (B_R)');
subplot(2,2,2);
pixplot(mesh.raxis,mesh.zaxis,squeeze(bphi(:,nphi,:)));
axis equal; axis tight;
title('Toridal Magnetic Field (B_\phi)');
subplot(2,2,3);
pixplot(mesh.raxis,mesh.zaxis,squeeze(bz(:,nphi,:)));
axis equal; axis tight;
title('Vertical Magnetic Field (B_Z)');
subplot(2,2,4);
pixplot(mesh.raxis,mesh.zaxis,squeeze(bphi(:,nphi,:)));
axis equal; axis tight;
hold on
quiver(mesh.raxis,mesh.zaxis,squeeze(br(:,nphi,:))',squeeze(bz(:,nphi,:))','k');
hold off
axis equal; axis tight;
title('Total Field');
case '2DVECTOR'
pixplot(mesh.raxis,mesh.zaxis,squeeze(bphi(:,nphi,:)));
axis equal; axis tight;
hcbar=colorbar;
ylabel(hcbar,'Toroidal Magnetic Field [T]');
hold on
quiver(mesh.raxis,mesh.zaxis,squeeze(br(:,nphi,:))',squeeze(bz(:,nphi,:))','k');
hold off
axis equal; axis tight;
title('Total Field');
xlabel('R [m]');
ylabel('Z [m]');
case 'GRID'% Plot 3D grid cutplanes
% Because the actual grids are very fine we abstract to a
% 20x20 grid for visualization
rtemp=mesh.rmin:(mesh.rmax-mesh.rmin)/19:mesh.rmax;
ztemp=mesh.zmin:(mesh.zmax+mesh.zmax)/19:mesh.zmax;
% First we need to 2d arrays for r and z
raxis2d=repmat(rtemp',[1 20]);
zaxis2d=repmat(ztemp,[20 1]);
% Next we need to map to X(nr,nphi,nz) and Y(nr,nphi,nz)
% Z is the same on each cut plane
x2d=zeros(20,mesh.nphi,20);
y2d=zeros(20,mesh.nphi,20);
for i=1:20
for j=1:mesh.nphi
x2d(i,j,:)=raxis2d(i,:).*cos(mesh.phi(j));
y2d(i,j,:)=raxis2d(i,:).*sin(mesh.phi(j));
end
end
% Now for each surface we create a create the patch object
% The verex list will run over nr
for i=1:mesh.nphi
vertex=[0 0 0];
for j=1:20
for k=1:20
vertex=[vertex; x2d(k,i,j) y2d(k,i,j) zaxis2d(k,j)];
end
end
vertex=vertex(2:size(vertex,1),:);
faces=[1 2 2+20 1+20];
for j=1:19
for k=1:19
index=j+20*(k-1);
faces=[faces; index index+1 index+1+20 index+20];
end
end
faces=faces(2:size(faces,1),:);
hold on
patch('Vertices',vertex,'Faces',faces,'FaceColor','none',...
'EdgeAlpha',0.1)
hold off
end
xlim([0 mesh.rmax]);
ylim([0 mesh.rmax]);
zlim([-mesh.zmax mesh.zmax]);
axis equal; axis tight;
view(3)
otherwise
disp(['ERROR: Unrecognized plot type ' strtrim(plottype)]);
end
return
end
