% this plot routine will plot the desity contours for Titan

clf

load bcompsw -ascii
bx = bcompsw(:,1);
by = bcompsw(:,2);
bz = bcompsw(:,3);
b = bcompsw(:,4);

clear bx
clear by
clear bz
clear bcompsw

% convert the one dimensional array into a three dimensional array
ix = 0;
iy = 0;
iz = 0;

% Constants
rm = 339400000

% convert the one dimensional array into a three dimensional array
ix = 0;
iy = 0;
iz = 0;


% number of cells for full half a planet Titan run

%nx = 66+1;
%ny = 153+1;
%nz = 91+1;

nx = 110 + 1;
ny = 115 + 1;
nz = 85 + 1;


% the length of the grid for RUN92 and Run59
xlen = 3.0e9/rm
ylen = 3.0e9/rm
zlen = 2.55e9/rm

% Cell size 
dx = xlen/(nx-1)
dy = ylen/(ny-1)
dz = zlen/(nz-1)

iskipx = 1
iskipy = 1
iskipz = 1


%dx = dx/rm
%dy = dy/rm
%dz = dz/rm

bmax = 0.0;
bmin = 100;

bmag = zeros(nx,ny,nz);

% set up the data into a 3-D array
% Note:  reducevolume could be used to do this
for k = 1:iskipz:nz
   iz = iz+1
   iy = 0;
   for j = 1:iskipy:ny
      iy = iy+1;
      ix = 0;
      for i = 1:iskipx:nx
      ix = ix+1;
         index = i + (j-1)*nx + (k-1)*nx*ny;
         if (b(index) > bmax) 
            bmax = b(index);
         end
         if (b(index) < bmin) 
            bmin = b(index);
         end
         bmag(ix,iy,iz) = b(index) / 1.0e-5;
      end
   end
end

rm2 = 2*rm;

% Create a 3-d plaid mesh for the slice routine

[y,x,z] = meshgrid(-1*ylen/2:dy:ylen/2,-1*xlen/2:dx:xlen/2,-1*zlen/2:dz:zlen/2);

ylb = -1*xlen/2;
yub = xlen/2;
xlb = -1*ylen/2;
xub = ylen/2;
zlb = -1*zlen/2;
zub = zlen/2;

slice(y,x,z,bmag,ylb,0,zlb);
set(gca,'YDir','rev')
shading interp;
colorbar

xlabel('Y Rt')
ylabel('X Rt')
zlabel('Z Rt')

%caxis([0 3])
%invrtgray = flipud(gray(60))
%colorbar('vert')
%colormap(invrtgray)


%    p = patch(isosurface(y,x,z,bmag,6.0));
%    p2 = patch(isosurface(x,y,z,dens,dmax));
%    p3= patch(isosurface(x,y,z,dens,0.1e-4));
%    p4 = patch(isosurface(x,y,z,dens,0.05e-4));
         
%isonormals(y,x,z,bmag,p)
%isonormals(y,x,z,dens,p2)
%set(p,'FaceColor','blue','EdgeColor','none','AmbientStrength',0.5);
%set(p2,'FaceColor','cyan','EdgeColor','none','AmbientStrength',0.5);
%set(p3,'FaceColor','cyan','EdgeColor','none','AmbientStrength',0.5);
%set(p4,'FaceColor','blue','EdgeColor','none','AmbientStrength',0.5);

%h = legend('2.e-4','0.9e-4','0.1e-4','0.05e-4',2)

hold on

%sphere

daspect([1 1 1])
% view(3)
% camlight
% camlight(180,30) 
%lighting phong

% lighting gouraud
 
stop


light('Position',[0 0 -1],'Style','infinite');



% set up the data into a 2-D array

for k = 1:iskipz:nz
      for j = 1:iskipy:ny       
        Bmagxplan(j,k) = Bmagplot(1,j,k);
   end
end

% Set up the mesh axis & convert to units of rm
for k = 1:iskipz:nz
   z(k) = (k-1)*dz;
end

for j = 1:iskipy:ny
   y(j) = (j-1)*dy;
end
%[z2d,y2d] = meshgrid(z,y);
[y2d,z2d] = ndgrid( 1:iskipy:ny, 1:iskipz:nz);

numcont = 30
ylabel('Distance Along Fieldline (Rm)');
xlabel('Polar Direction (Rm)');
zlabel('Flow Direction from the Sun (Rm)');

%contour(y2d,z2d,Bmagxplan,6);
%pcolor(z2d,y2d,Bmagxplan);
% invert the gray scale map by using flipud

invrtgray = flipud(gray(numcont))
contourf(y2d,z2d,Bmagxplan,numcont);
xlabel('Polar Direction (Rm)');
ylabel('Flow Direction from the Sun (Rm)');
colorbar('vert')
colormap(invrtgray)