% Animation that rotates the Earth by 360 degrees, then removes an octant 
% to reveal the inner structure of the Earth of the lithosphere, upper 
% mantle, lower mantle, outer core, and inner core 

clear
close all

vid = VideoWriter('earthCrossSection.mp4','MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
set(gca,'linewidth',7);
lp = [0.7 0.7 0.7]; 
LL = 1.25; 

[xx,yy,zz] = sphere(3000); % Earth with octant removed 
[xx2,yy2,zz2] = sphere(3000); % Full spherical Earth 
[xx3,yy3,zz3] = sphere(3000); % Octant of Earth 

% Remove octant 
for ii = 1:3001
    for jj = 1:3001
        if (xx(ii,jj)>=0) && (yy(ii,jj)>=0) && (zz(ii,jj)<=0)
            xx(ii,jj) = NaN;
            yy(ii,jj) = NaN;
            zz(ii,jj) = NaN;
        end
    end
end

% Octant only 
for ii = 1:3001
    for jj = 1:3001
        if (xx3(ii,jj)<0) || (yy3(ii,jj)<0) || (zz3(ii,jj)>0)
            xx3(ii,jj) = NaN;
            yy3(ii,jj) = NaN;
            zz3(ii,jj) = NaN;
        end
    end
end 

% Read in representation of different components of Earth from lithosphere 
% to inner core as a set of coloured concentric circles. Drawn by the 
% author 
earthSlice = imread('earthSlice.png');
XX = [1 1; -1 -1];
YY = [1 -1; 1 -1];
ZZ = [0 0; 0 0];

% Get rid of regions that are black due to original image having 
% transparent regions outside circular representations of the different 
% depth regions of the Earth. Set them to NaN 
iimax = size(earthSlice,1);
jjmax = size(earthSlice,2);
for ii = 1:iimax
    for jj = 1:jjmax
        if (earthSlice(ii,jj,1)==0) && ...
                (earthSlice(ii,jj,2)==0) && ...
                (earthSlice(ii,jj,3)==0)
            earthSlice(ii,jj,1) = NaN;
            earthSlice(ii,jj,2) = NaN;
            earthSlice(ii,jj,3) = NaN;
        end
    end
end
AlphaData = earthSlice(:,:,1) + earthSlice(:,:,2) + earthSlice(:,:,3);

% Read in surface representation of continents and oceans from an image
% that was produced to map onto a spherical surface. Original image
% 'earth.jpg' licenced under Creative Commons at 
% https://commons.wikimedia.org/wiki/File:PathfinderMap_hires_(4996917742).jpg
earth = imread('earthMap.png'); 

% Full earth, no octant removed 
% Rotate this through a day 
for rot = 0:pi/60:2*pi
    newplot
    hold off
    fullEarth = surf(xx2,yy2,-zz2,earth,'FaceColor','texturemap','EdgeColor','none');
    rotate(fullEarth,[0 0 1],rot*180/pi,[0,0,0]);
    axis equal
    axis off
    hold on
    view(110,23);
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]); 
    set(gca, 'Projection','perspective');

    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);

end 
pauseMov(vid,60); 

% Now move the octant (approximately covering longitude from Greenwich to
% India and from the equator to North pole in latitude) away from the
% original sphere 
for ii = 0:0.005:0.25 
    newplot
    hold off
    view(110,23);
    surf(xx3+ii,yy3+ii,-zz3+ii,earth,'FaceColor','texturemap','EdgeColor','none')
    axis equal
    axis off
    hold on
    surf(xx,yy,-zz,earth,'FaceColor','texturemap','EdgeColor','none');
    
    surf(XX,YY,-ZZ,earthSlice*0.98,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    surf(ZZ,XX,-YY,earthSlice*1.1,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    surf(YY,ZZ,-XX,earthSlice*1.25,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    
    view(110,23);
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]); 
    set(gca, 'Projection','perspective');

    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
end 

% Fade the same quadrant away to reveal inner structure 
for ii = 1:-1/30:0 
    newplot
    hold off
    view(125,31);
    surf(xx3+0.25,yy3+0.25,-zz3+0.25,earth,'FaceColor','texturemap','EdgeColor','none')
    alpha(ii)
    axis equal
    axis off
    hold on
    surf(xx,yy,-zz,earth,'FaceColor','texturemap','EdgeColor','none');
    
    surf(XX,YY,-ZZ,earthSlice*0.98,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    surf(ZZ,XX,-YY,earthSlice*1.1,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    surf(YY,ZZ,-XX,earthSlice*1.25,'FaceColor','texturemap','EdgeColor','none','Alphadata',AlphaData,'FaceAlpha','texturemap');
    
    view(110,23);
    xlim([-LL LL]);
    ylim([-LL LL]);
    zlim([-LL LL]); 
    set(gca, 'Projection','perspective');
    set(gcf, 'color', 'w'); % Black background

    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
    
end 

close(vid) 
