%	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);
%               r = sqrt((i-21)^2 + (j-21)^2);
%               icy3(i,j,k) = icy(i+30,j+30,k+30);

%                if (r < 15)
%                  icy3(i,j,k) = 0.0;
%                end

%	      	m=m+1;
%	      end
%	   end
       %	end

	x=[-10:.5:10];
	y=[-10:.5:10];
	z=[-10:.5:10];
    
	x2=[-34.75:.5:34.75];
	y2=[-34.75:.5:34.75];
	z2=[-9.75:.5:9.75];
    
f=figure;
set(f, 'Color', 'black')

lighting phong
set(gcf,'Renderer','zbuffer')

	
%h=slice(y,x,z,icy2,[],[],0);

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

%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)

%pii=patch(isosurface(x,y,z,icy3,20)); 
%isonormals(x,y,z,icy3,pii);
%set(pii,'FaceColor','magenta','EdgeColor','none','FaceAlpha',0.3,...
%          'BackFaceLighting','lit','AmbientStrength',.6);

hold on
%        clear icy
%	clear icy2
%        clear icy3

%	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

%	load OH2 -ascii;
%       	temp1 = OH2(:,4);
%	clear OH2
%	m=1
%	for k = 1:101;
%	   for j = 1:101;
%	      for i = 1:101;
%	      	oha(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;
%	      	oh(i,j,k) = oha(i+30,j+30,k+30);
%	      	m=m+1;
%	      end
%	   end
%	end
%clear oha
%p2=patch(isosurface(x,y,z,oh,8)); 
%isonormals(x,y,z,oh,p2);
%set(p2,'FaceColor','yello','EdgeColor','none','FaceAlpha',0.5,...
%             'BackFaceLighting','lit','AmbientStrength',.6);

%	clear oh
%hold on

	load titann -ascii;
       	temp1 = titann(:,4);
	clear titann
	m=1
	for k = 1:40;
	   for j = 1:140;
	      for i = 1:140;
	      	tita(i,j,k) = temp1(m);
	      	m=m+1;
	      end
	   end
	end
	clear temp1

p2=patch(isosurface(x2,y2,z2,tita,2)); 
isonormals(x2,y2,z2,tita,p2);
set(p2,'FaceColor','yello','EdgeColor','none','FaceAlpha',0.5,...
             'BackFaceLighting','lit','AmbientStrength',.6);

	clear tita
hold on

	load testith3 -ascii;
	wion = testith3(:,4);
	clear testith3
	m = 1
	for k = 1:21;
	   for i = 1:101;
	      for j = 101:-1:1;
		if wion(m) < 500 ;
			wion(m) = 0;
		end
	      	ion(i,j,k) = wion(m);
	      	m=m+1;
	      end
	   end
	end

	ion2=smooth3(ion);
% 	clear ion

	x=[-25:.5:25];
	y=[-25:.5:25];
	z=[-5:.5:5];
p=patch(isosurface(x,y,z,ion2,4000)); 
isonormals(x,y,z,ion2,p);
set(p,'FaceColor','green','EdgeColor','none','FaceAlpha',0.5,...
             'BackFaceLighting','lit','AmbientStrength',.6);

hold on
clear ion
clear wion
clear ion2


% Plot RINGS

for ir = 1:99
   radius = 1.112 + (ir-1) * (1.92-1.112)/99;
   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:50
   radius = 2.02 + (ir-1) * (2.27-2.02)/51;
   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',3)

   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',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([-45 45 -45 45 -45 45])

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])
light('Position',[30 0 -14.4])
lighting phong




% Mimas

 %rgb=imread('europa.jpg');
 %x64=double(rgb)/255;
 %[xtet,ytet,ztet] = sphere(25);
  %   xtet = xtet*529/60268*10 + 1.534;
  %   ytet = ytet*529/60268*10 - 2.692;
  %   ztet = ztet*529/60268*10 + 0.066;

%surface(xtet,ytet,ztet,flipdim(x64,1),...
%        'FaceColor','texturemap',...
%        'EdgeColor','none',...
 %       'CDataMapping','direct','AmbientStrength',0.2)

% Enceladus

% [xtet,ytet,ztet] = sphere(25);
 %    xtet = xtet*529/60268*10 + 3.793;
 %    ytet = ytet*529/60268*10 + 1.165;
 %    ztet = ztet*529/60268*10 - 0.000;

%surface(xtet,ytet,ztet,flipdim(x64,1),...
%        'FaceColor','texturemap',...
 %       'EdgeColor','none',...
 %       'CDataMapping','direct','AmbientStrength',0.2)

% Tethys

 %[xtet,ytet,ztet] = sphere(25);
 %    xtet = xtet*529/60268*10 - 3.813;
 %    ytet = ytet*529/60268*10 - 3.059;
 %    ztet = ztet*529/60268*10 + 0.092;

%surface(xtet,ytet,ztet,flipdim(x64,1),...
 %       'FaceColor','texturemap',...
 %       'EdgeColor','none',...
 %       'CDataMapping','direct','AmbientStrength',0.2)

% Dione

 %[xtet,ytet,ztet] = sphere(25);
 %    xtet = xtet*529/60268*10 + 2.350;
 %    ytet = ytet*529/60268*10 - 5.808;
 %    ztet = ztet*529/60268*10 - 0.003;

%surface(xtet,ytet,ztet,flipdim(x64,1),...
 %       'FaceColor','texturemap',...
 %       'EdgeColor','none',...
 %      'CDataMapping','direct','AmbientStrength',0.2)

% Rea

 %[xtet,ytet,ztet] = sphere(25);
  %   xtet = xtet*529/60268*10 + 8.007;
  %   ytet = ytet*529/60268*10 + 3.495;
  %   ztet = ztet*529/60268*10 - 0.047;

%surface(xtet,ytet,ztet,flipdim(x64,1),...
 %       'FaceColor','texturemap',...
 %       'EdgeColor','none',...
 %       'CDataMapping','direct','AmbientStrength',0.2)

 rgb=imread('titanmap.gif');
 x64=double(rgb)/255;
 [xtet,ytet,ztet] = sphere(25);
     xtet = xtet*2575/60268*1 + 18.656;
     ytet = ytet*2575/60268*1 - 7.196;
     ztet = ztet*2575/60268*1 - 0.130;

surface(xtet,ytet,ztet,flipdim(x64,1),...
        'FaceColor','texturemap',...
        'EdgeColor','none',...
        'CDataMapping','direct','AmbientStrength',0.2)



camproj perspective 
camva(20)
camlight

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

load oTAext -ascii
xt = double(oTAext(:,1));
yt = double(oTAext(:,2));
zt = double(oTAext(:,3));
t = double(oTAext(:,4));

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

axis vis3d
axis off
stop



%c = [1 1 1]


%for i=1:30
%	s1 = rand(1) * 50 - 25;
%        s2 = rand(1) * 50 - 25;
%        s3 = rand(1) * 50 - 25;
%        sh = rand(1) * 25;
%        ps = fix(rand(1)*3+1);
%        ub = 45;
%	lb = -45; 
%	scatter3(s1,s2,ub,ps,c,'.');
%	scatter3(lb,s2,s3,ps,c,'.');
%	scatter3(s3,s2,lb,ps,c,'.');
%	scatter3(s3,ub,s2,ps,c,'.');
%	scatter3(s2,lb,s1,ps,c,'.');
%	scatter3(ub,s1,s3,ps,c,'.');
%end

%        scatter3(45,0,21.6,5,c,'*');

hold on

%orbit camera

%campos([40,0,15])

%for i=1:36
%    camorbit(10,0,'data',[0 0 1])
%    drawnow
%    F(i) = getframe(gcf);
%end

%campos([15,0,5.625])

%for j=1:36
%    camorbit(10,0,'data',[0 0 1])
%    drawnow
%    F(i+j) = getframe(gcf);
%end

%disp('Please wait: saving as AVI');
%movie2avi(F,'SOI-flyby','fps',10,'compression','none');






 %stop
tiny=1e-10
tinyx = min(abs(xt))/100
tinyy = min(abs(yt))/100
tinyz = min(0.005)/100
j=2;

campos([xt(1),yt(1),zt(1)])
camtarget([xt(1)+min(xt)/100,yt(1)+min(yt)/100,zt(1)+tinyz])
drawnow
 F(1) = getframe(gcf);
for i=2:length(xt)-1
i
  campos([xt(i),yt(i),zt(i)])
%%  camtarget([xt(i+5)+min(xt)/100,yt(i),zt(i)])
%%     r = sqrt((xt(i+2)-xt(i))^2+(yt(i+2)-yt(i))^2+(zt(i+2)-zt(i))^2);
%%     dx = (xt(i+2)-xt(i))/r*1e-13;   %/abs(xt(i+1)-xt(i)+tiny)%*tinyx
%%     dy = (yt(i+2)-yt(i))/r*1e-13;   %/abs(yt(i+1)-yt(i)+tiny)%*tinyy
%%     dz = (zt(i+2)-zt(i))/r*1e-13;   %/abs(zt(i+1)-zt(i)+tiny)%*tinyz

     xv = (xt(i+1)-xt(i-1))/(t(i+1)-t(i-1))
     yv = (yt(i+1)-yt(i-1))/(t(i+1)-t(i-1))
     zv = (zt(i+1)-zt(i-1))/(t(i+1)-t(i-1))
 
     dx = xv*1e-5;   %/abs(xt(i+1)-xt(i)+tiny)%*tinyx
     dy = yv*1e-5;   %/abs(yt(i+1)-yt(i)+tiny)%*tinyy
     dz = zv*1e-5;   %/abs(zt(i+1)-zt(i)+tiny)%*tinyz

     camtarget([xt(i)+dx,yt(i)+dy,zt(i)+dz])
     if(i > 175)
            camtarget([xt(i)-dx,yt(i)-dy,zt(i)-dz])
     end

  %camlight('headlight')
  drawnow
  F(j) = getframe(gcf);
j = j+1;
end
disp('Please wait: saving as AVI');
movie2avi(F,'SOI-flythrough1d','fps',10,'compression','none');

movie(F)

stop
