%% Plots the location
close all
clear all

kp=3;

mltv=0:0.25:24;
%mltvf=[18:0.05:19];
mltvf=0:0.05:24;

lshv=linspace(3,25,50);
lshv=lshv';

[mltm,lshm]=meshgrid(mltv,lshv);
ep=modmcle(mltm,lshm,kp,1);
[xmat,ymat]=pol2cart(mltm*pi/12,lshm);

[lshb,epb]=modmcilpp2(mltvf,kp);
[xmatpp,ymatpp]=pol2cart(mltvf*pi/12,lshb);


mltv=0:24;
lshv=3:12;
lshv=lshv';

[mltm2,lshm2]=meshgrid(mltv,lshv);
[vth,vrd]=modmcilvel2(lshm2,mltm2,kp,0);
[xmat2,ymat2]=pol2cart(mltm2*pi/12,lshm2);


vrdx=vrd.*cos(mltm2*pi/12);
vrdy=vrd.*sin(mltm2*pi/12);
vthx=vth.*cos(mltm2*pi/12+pi/2);
vthy=vth.*sin(mltm2*pi/12+pi/2);

velx=vrdx+vthx;
vely=vrdy+vthy;


timv=0:0.5:23.5;
[xve,yve]=pol2cart(timv*pi/12,ones(size(timv)));
timv2=-6:0.5:6;
[xve2,yve2]=pol2cart(timv2*pi/12,ones(size(timv2)));

figure(1)
contour(xmat,ymat,ep,50,'k')
hold on
fill(xve,yve,'g')
fill(xve2,yve2,'k')
hold off
grid
axis equal
axis([-12.5 12.5 -12.5 12.5])
xlabel('Earth radii','fontweight','bold','fontsize',16)
ylabel('Earth radii','fontweight','bold','fontsize',16)
set(gca,'fontweight','bold','fontsize',14)

figure(2)
quiver(xmat2,ymat2,velx,vely,1.75)
axis equal
hold on
fill(xve,yve,'g')
fill(xve2,yve2,'k')
hold off
grid
axis equal
axis([-12.5 12.5 -12.5 12.5])
xlabel('Earth radii','fontweight','bold','fontsize',16)
ylabel('Earth radii','fontweight','bold','fontsize',16)
set(gca,'fontweight','bold','fontsize',14)

figure(3)
quiver(xmat2,ymat2,velx,vely,1.75,'k')
axis equal
hold on
plot(xve,yve,'k')
fill(xve2,yve2,'k')
hold off
grid
axis equal
axis([-12.5 12.5 -12.5 12.5])
xlabel('Earth radii','fontweight','bold','fontsize',16)
ylabel('Earth radii','fontweight','bold','fontsize',16)
set(gca,'fontweight','bold','fontsize',14)
text(-0.2,1.5,'6','fontweight','bold','fontsize',14)
text(-0.5,-1.5,'18','fontweight','bold','fontsize',14)
text(1.25,0,'0','fontweight','bold','fontsize',14)
text(-2.5,0,'12','fontweight','bold','fontsize',14)

figure(4)
quiver(xmat2,ymat2,velx,vely,1.75,'k')
axis equal
hold on
plot(xve,yve,'k')
fill(xve2,yve2,'k')
plot(xmatpp,ymatpp,'k','linewidth',2)

hold off
grid
axis equal
axis([-12.5 12.5 -12.5 12.5])
xlabel('Earth radii','fontweight','bold','fontsize',16)
ylabel('Earth radii','fontweight','bold','fontsize',16)
set(gca,'fontweight','bold','fontsize',14)
text(-0.2,1.5,'6','fontweight','bold','fontsize',14)
text(-0.5,-1.5,'18','fontweight','bold','fontsize',14)
text(1.25,0,'0','fontweight','bold','fontsize',14)
text(-2.5,0,'12','fontweight','bold','fontsize',14)