	ion = zeros(140,140,40);
	icy = zeros(140,140,40);
	rings = zeros(41,41,41);

	load testit4-icy -ascii;
	temp1 = testit4_icy(:,5);
	clear testith4_icy

	m=1
	for k = 1:101;
	   for j = 1:101;
	      for i = 1:101;
	      	icy(i,j,k) = temp1(m);
	      	m=m+1;
	      end
	   end
	end
	clear temp1

	for k = 1:41;
	   for j = 1:41;
	      for i = 1:41;
	      	icy2(i,j,k) = icy(i+30,j+30,k+30);
	      	m=m+1;
	      end
	   end
	end

	x=[-10:.5:10];
	y=[-10:.5:10];
	z=[-10:.5:10];

f=figure;
set(f, 'Color', 'black')
	
h=slice(y,x,z,icy2,[],[],0);

set(h,'EdgeColor','none','FaceColor','interp',...
    'FaceAlpha',0.4,'AmbientStrength',.7)
%     alpha(0.4)

	axis([-10 10 -10 10 -10 10])

	xlabel('X (Saturn Radii)')
	ylabel('Y (Saturn Radii)')
	zlabel('Z (Saturn Radii)')

title('H_2O Number Density Isosurface for Saturn''s Rings and Icy Satellites')

pi=patch(isosurface(x,y,z,icy2,65)); 
isonormals(x,y,z,icy2,pi);
set(pi,'FaceColor','magenta','EdgeColor','none','FaceAlpha',0.3,...
          'BackFaceLighting','lit','AmbientStrength',.6);
%     alpha(0.3)

hold on
	clear icy

	load testit4-rings -ascii;
       	temp1 = testit4_rings(:,5);
	clear testith4_rings

	m=1
	for k = 1:101;
	   for j = 1:101;
	      for i = 1:101;
	      	ring(i,j,k) = temp1(m);
	      	m=m+1;
	      end
	   end
	end
	clear temp1
	for k = 1:41;
	   for j = 1:41;
	      for i = 1:41;
	      	rings(i,j,k) = ring(i+30,j+30,k+30);
	      	m=m+1;
	      end
	   end
	end

p2=patch(isosurface(x,y,z,rings,30)); 
isonormals(x,y,z,rings,p2);
set(p2,'FaceColor','red','EdgeColor','none','FaceAlpha',0.5,...
             'BackFaceLighting','lit','AmbientStrength',.6);

	clear rings
hold on

for ir = 1:49
   radius = 1.112 + (ir-1) * (1.92-1.112)/50;
   for t=1:101
	thet = 2*3.141592654*(t-1)/(100);
        xr(t) = radius * cos(thet);
        yr(t) = radius * sin(thet);
        zr(t) = 0.0;
   end
   hra = plot3(xr,yr,zr);
   set(hra,'color',[0.3 0.3 0.3],'linewidth',1)
end

for ir = 1:20
   radius = 2.02 + (ir-1) * (2.27-2.02)/21;
   for t=1:101
	thet = 2*3.141592654*(t-1)/(100);
        xr(t) = radius * cos(thet);
        yr(t) = radius * sin(thet);
        zr(t) = 0.0;
   end
   hrb = plot3(xr,yr,zr);
   set(hrb,'color',[0.3 0.3 0.3],'linewidth',1)
end

   radius = 2.3264;
   for t=1:101
	thet = 2*3.141592654*(t-1)/(100);
        xr(t) = radius * cos(thet);
        yr(t) = radius * sin(thet);
        zr(t) = 0.0;
   end
   hrf = plot3(xr,yr,zr);
   set(hrf,'color',[0.2 0.2 0.2],'linewidth',2)

   radius = 2.751;
   for t=1:101
	thet = 2*3.141592654*(t-1)/(100);
        xr(t) = radius * cos(thet);
        yr(t) = radius * sin(thet);
        zr(t) = 0.0;
   end
   hrg = plot3(xr,yr,zr);
   set(hrg,'color',[0.3 0.3 0.3],'linewidth',2)


        load operi74 -ascii
xt = operi74(:,1);
yt = operi74(:,2);
zt = operi74(:,3);

H = plot3(xt,yt,zt);
	set(H,'Color',[1 0.5 0,],'LineWidth',3)



 rgb=imread('sat0fds1a.jpg');

   x64=double(rgb)/255;

%h=findobj('Type','surface');rot = flipdim(x64,1);
%%   temp = rot(:,1:1080,:);
%%   rot(:,1:1080,:)=[];
%%   rot = [rot temp];
%   set(h,'CData',rot,'FaceColor','texturemap','AmbientStrength',0.7)
%%   daspect([1 1 1]);
%hold on

[x2,y2,z2]=sphere(20);

surface(x2,y2,z2,flipdim(x64,1),...
        'FaceColor','texturemap',...
        'EdgeColor','none',...
        'CDataMapping','direct','AmbientStrength',0.2)

axis('square')
set(gca,'Color',[0 0 0]);
set(gca,'XColor','white',...
        'YColor','white',...
        'ZColor','white');
set(gca,'FontSize',20);

light('Position',[2.5 0 -1.2])
lighting phong
%camlight

print -djpeg80 trajper74.jpg

stop
