clear line

% specify the r,z, and phi starting points for the field lines
srvals=[0.29:0.007:0.38];
sphivals=(2*pi/32)*[16:0.5:20];
szvals=-0.29;

clear sr sz sphi
%generate starting point arrays
for jnr=1:length(srvals)
    for jnphi=1:length(sphivals)
        for jnz=1:length(szvals)
            jj=jnr+(jnphi-1)*(length(srvals))+(jnz-1)*(length(srvals)*length(sphivals));
sr(jj)=srvals(jnr);
sz(jj)=szvals(jnz);
sphi(jj)=sphivals(jnphi);
end
end
end
% routine to specify and plot grid boundary
sibound

close all

%order of colors for plotting lines
co='rbgmcyrb';
%x and y positions of starting points
sx=sr.*cos(sphi);
sy=sr.*sin(sphi);

for jview=1:1

%create axes for plotting
g=subplot(2,1,jview);
set(g,'xtick',[]);
set(g,'ytick',[]);
set(g,'ztick',[]);

for jl=1:length(sr)
    jl
    %use streamline function to generate and plot field lines
    % data must be on cartesian mesh
    H=streamline(Y,X,Z,By,Bx,Bz,sx(jl),sy(jl),sz(jl),[0.05 90000]);
    set(H,'color',co(mod(jl-1,length(co))+1));
    set(H,'linewidth',2);
    hold on
    % make arrays of x,y, and z coordinated for the line
    line(jl).x=get(H,'XData');
    line(jl).y=get(H,'YData');
    line(jl).z=get(H,'Zdata');
    
    % routine to compue the length of the line
    line_length
end
%plots (part of) the 3D boundary of HIT-SI 
grid3plot
axis image
alpha(0.5)
end
%plotting effects for HIT-SI boundary
map=colormap('copper');
colormap(map(end-10,:));
material metal;


% get intersections of field lines with specified planes
phi_planes=[pi/4 5*pi/4];
subplot(2,1,2)
for jpl=1:2
    hold on
    phi_plane=phi_planes(jpl);
for jl=1:length(sr)
    % call routine to find plane crossings
    plane_cross
    if jpl==2
        rc=-rc;
    end
    if (~isempty(rc))
    plot(rc,zc,strcat(co(mod(jl-1,length(co))+1),'.'))
    hold on
end
end
end
sibound
axis image