MATLAB: Fixing axes for a rotating globe under an orbit

animationmodelingrotations

Hi,
Please forgive me if this is a dumb question, I am very new to MATLAB and am trying to learn! I'm writing a program to plot the time-dependent orbital position of a satellite in the geocentric frame of reference. I'm trying to get the globe model of the Earth to rotate underneath the orbit in a realistic way. Then, I want to output the frames into a video. The issue I'm having is that the rotation is moving the axes wildly around instead of rotating the globe in place about the z-axis (or perhaps it is the view that is being moved? I'm really not sure).
How do I change this so it will simply rotate in place? As I said, I'm trying to learn so please explain where I'm going wrong!
To summarize: I'm trying to plot an orbit around the Earth while simultaneously having the Earth rotate on its z-axis. Then, I want to export this to video. However, the way I'm doing this is doesn't seem to rotate the Earth in place.
NOTE: The function "trajectory" just calculates the time-dependent position of an orbiting object. It outputs a N x 3 matrix of positions.
function GeocentricTrack(R0, V0, dt, step,mu)
      load topo
      if nargin == 4
          mu  = 398600;        % [km^3/s^2] Standard Gravitational Parameter
      end
      [Ri,Vi] = trajectory(R0,V0,dt,step,mu);
      %%Instantiate Movie Maker
      str = input('Would you like a movie?\n','s');
      if strcmpi(str,'y') || strcmpi(str,'yes')
          movName = input('Please name the movie\n','s');
          con     = strcat('~/Documents/matlab/movie/',movName);
          vid = VideoWriter(con,'MPEG-4');
          open(vid);
          grs80 = referenceEllipsoid('grs80','km');
          figure('Renderer','opengl')
          ax = axesm('globe','Geoid',grs80,'Grid','off', ...
              'GLineWidth',1,'GLineStyle','-',...
              'Gcolor',[0.9 0.9 0.1],'Galtitude',100);
          ax.Position = [0 0 1 1];
          axis equal off
          h1 = geoshow(topo,topolegend,'DisplayType','texturemap');
          demcmap(topo)
          land = shaperead('landareas','UseGeoCoords',true);
          h2 = plotm([land.Lat],[land.Lon],'Color','black');
          rivers = shaperead('worldrivers','UseGeoCoords',true);
          h3 = plotm([rivers.Lat],[rivers.Lon],'Color','blue');
          n          = dt/step;
          rate       = 360/86164;       % Earth's rotation rate [deg/sec]
          rotperstep = n*rate;
          totalrot   = rotperstep*step; % Total degrees of rotation
          rotation   = 360:-rotperstep:totalrot;
          direction  = [0, 0, 1];
          set(gca,'units','pixels','position', [0,0,560,420])
          rotate(h1,direction,rotation(1))
          rotate(h2,direction,rotation(1))
          rotate(h3,direction,rotation(1))
          view(0,23.5)
          drawnow
          scatter3(Ri(1,1),Ri(1,2),Ri(1,3),'.r');
          frame = getframe(gca);
          writeVideo(vid,frame);
          hold on
          for i = 2:length(Ri)
              rotate(h1,direction,rotation(i))
              rotate(h2,direction,rotation(i))
              rotate(h3,direction,rotation(i))
              view(0,23.5)
              drawnow
              scatter3(Ri(i,1),Ri(i,2),Ri(i,3),'.r');
              frame = getframe(gca);
              writeVideo(vid,frame);
          end
          close(vid);
      else
          grs80 = referenceEllipsoid('grs80','km');
          figure('Renderer','opengl')
          ax = axesm('globe','Geoid',grs80,'Grid','off', ...
              'GLineWidth',1,'GLineStyle','-',...
              'Gcolor',[0.9 0.9 0.1],'Galtitude',100);
          ax.Position = [0 0 1 1];
          axis equal off
          view(3)
          geoshow(topo,topolegend,'DisplayType','texturemap')
          demcmap(topo)
          land = shaperead('landareas','UseGeoCoords',true);
          plotm([land.Lat],[land.Lon],'Color','black')
          rivers = shaperead('worldrivers','UseGeoCoords',true);
          plotm([rivers.Lat],[rivers.Lon],'Color','blue')
          for i = 1:length(Ri)
              scatter3(Ri(i,1),Ri(i,2),Ri(i,3),'.r');
          end
      end
  end

Best Answer

I can’t run you code because I don’t have the Mapping Toolbox.
I offer this code snippet as inspiration:
[X,Y,Z] = sphere;
x = linspace(-1, 1, 75);
y = -x.^2;
figure(1)
for k2 = 1:2
    for k1 = 1:length(x)
        surf(X+y(k1), Y+sin(2*pi*x(k1)), Z+cos(2*pi*x(k1)))
        axis([-5  5    -5  5    -5  5])
        axis square
        view([-20  20])
        refreshdata
        drawnow
    end
end
You probably need to define the reference with an axis function call in every iteration.
Related Question