{"id":2655,"date":"2014-05-14T08:51:01","date_gmt":"2014-05-14T13:51:01","guid":{"rendered":"https:\/\/blogs.mathworks.com\/community\/?p=2655"},"modified":"2014-07-29T12:42:22","modified_gmt":"2014-07-29T17:42:22","slug":"whats-that-in-the-sky-build-a-matlab-planetarium","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/community\/2014\/05\/14\/whats-that-in-the-sky-build-a-matlab-planetarium\/","title":{"rendered":"What&#8217;s that in the sky? Build a MATLAB Planetarium"},"content":{"rendered":"<div class=\"content\"><!--introduction--><p>I look up at the sky just after sunset and I see an especially bright star. It's probably a planet. But which one?<\/p><p>This question gives me a good opportunity to play around with MATLAB. Let's do a visualization that shows where the planets are relative to the earth and the sun. In the process, we'll use JSON services, the File Exchange, MATLAB graphics, and 3-D vector mathematics cribbed from Wikipedia.<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/explanation1.png\" alt=\"\"> <\/p><p>Here is the basic grade-school illustration of the solar system, the one that shows the planets rolling around the sun like peas on a plate. For simplicity, we're just showing the sun, the earth, the moon, Venus, and Mars.<\/p><p>But we never see anything like this with our own eyes. Instead, we see bright spots on a dark background somewhere \"up there.\" So let's simplify our problem to determining what direction each naked-eye planet is in. This leads to an image like this.<\/p><p><img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/explanation3.png\" alt=\"\"> <\/p><p>Our goal is to make an accurate up-to-date version of this diagram. Specifically, relative to the sun, where should we look to find the moon and the naked-eye planets (Mercury, Venus, Mars, Jupiter, and Saturn)? To do this, we need to solve a few problems.<\/p><div><ol><li>Find the planets<\/li><li>Find the unit vector pointing from earth to each planet<\/li><li>Squash all these vectors onto a single plane<\/li><li>Visualize the resulting disk<\/li><\/ol><\/div><!--\/introduction--><h3>Contents<\/h3><div><ul><li><a href=\"#efb551df-af2d-4c78-a73b-8f4806ab5a2c\">Where Are the Planets?<\/a><\/li><li><a href=\"#242720c5-da18-40ca-aea0-145e7d041532\">Parse the data<\/a><\/li><li><a href=\"#69db2a03-0af9-4e72-a4ce-17e6a7258579\">Aerospace Toolbox Ephemeris<\/a><\/li><li><a href=\"#7e098dad-e06c-4f38-867e-048051a55cd7\">Build the Solar System Structure<\/a><\/li><li><a href=\"#983013aa-7264-436a-babb-d7995a335893\">Plot the planets<\/a><\/li><li><a href=\"#691012af-279b-45e4-9237-e7a4fff3b0a4\">Plot the result<\/a><\/li><\/ul><\/div><h4>Where Are the Planets?<a name=\"efb551df-af2d-4c78-a73b-8f4806ab5a2c\"><\/a><\/h4><p>First of all, where are the planets? There's a free JSON service for just about everything these days. I found planetary data hosted on Davy Wybiral's site here:<\/p><p><a href=\"http:\/\/davywybiral.blogspot.com\/2011\/11\/planetary-states-api.html\">http:\/\/davywybiral.blogspot.com\/2011\/11\/planetary-states-api.html<\/a><\/p><pre class=\"codeinput\">url = <span class=\"string\">'http:\/\/www.astro-phys.com\/api\/de406\/states?bodies=sun,moon,mercury,venus,earth,mars,jupiter,saturn'<\/span>;\r\njson = urlread(url);\r\n<\/pre><h4>Parse the data<a name=\"242720c5-da18-40ca-aea0-145e7d041532\"><\/a><\/h4><p>Now we can use <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/42236-parse-json-text\">Joe Hicklin's JSON parser<\/a> from the File Exchange. It returns a well-behaved MATLAB structure.<\/p><pre class=\"codeinput\">data = JSON.parse(json)\r\n<\/pre><pre class=\"codeoutput\">\r\ndata = \r\n\r\n       date: 2.4568e+06\r\n    results: [1x1 struct]\r\n       unit: 'km'\r\n\r\n<\/pre><p>The payload is in the \"results\" field. Each entry has three position components and three velocity components.<\/p><pre class=\"codeinput\">data.results\r\n<\/pre><pre class=\"codeoutput\">\r\nans = \r\n\r\n    mercury: {{1x3 cell}  {1x3 cell}}\r\n        sun: {{1x3 cell}  {1x3 cell}}\r\n       moon: {{1x3 cell}  {1x3 cell}}\r\n    jupiter: {{1x3 cell}  {1x3 cell}}\r\n       mars: {{1x3 cell}  {1x3 cell}}\r\n      earth: {{1x3 cell}  {1x3 cell}}\r\n      venus: {{1x3 cell}  {1x3 cell}}\r\n     saturn: {{1x3 cell}  {1x3 cell}}\r\n\r\n<\/pre><p>The distances are in kilometers, and I'm not even sure how this representation is oriented relative to the surrounding galaxy. But it doesn't really matter, because all I care about is the relative positions of the bodies in question.<\/p><h4>Aerospace Toolbox Ephemeris<a name=\"69db2a03-0af9-4e72-a4ce-17e6a7258579\"><\/a><\/h4><p>Side note: We could also have used the <a href=\"https:\/\/www.mathworks.com\/help\/aerotbx\/ug\/planetephemeris.html\">Aerospace Toolbox<\/a> to get the same information.<\/p><p><tt>[pos,vel] = planetEphemeris(juliandate(now),'Sun','Earth')<\/tt><\/p><h4>Build the Solar System Structure<a name=\"7e098dad-e06c-4f38-867e-048051a55cd7\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"comment\">% List of bodies we care about<\/span>\r\nssList = {<span class=\"string\">'sun'<\/span>,<span class=\"string\">'moon'<\/span>,<span class=\"string\">'mercury'<\/span>,<span class=\"string\">'venus'<\/span>,<span class=\"string\">'earth'<\/span>,<span class=\"string\">'mars'<\/span>,<span class=\"string\">'jupiter'<\/span>,<span class=\"string\">'saturn'<\/span>};\r\n\r\nss = [];\r\n<span class=\"keyword\">for<\/span> i = 1:length(ssList)\r\n    ssObjectName = ssList{i};\r\n    ss(i).name = ssObjectName;\r\n    ssData = data.results.(ssObjectName);\r\n    ss(i).position = [ssData{1}{:}];\r\n    ss(i).velocity = [ssData{2}{:}];\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4>Plot the planets<a name=\"983013aa-7264-436a-babb-d7995a335893\"><\/a><\/h4><pre class=\"codeinput\"><span class=\"comment\">% Plot in astronomical units<\/span>\r\nau = 149597871;\r\nk = 5;\r\n\r\ncla\r\n<span class=\"keyword\">for<\/span> i = 1:length(ss)\r\n    p = ss(i).position\/au;\r\n\r\n    line(p(1),p(2),p(3), <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'Marker'<\/span>,<span class=\"string\">'.'<\/span>,<span class=\"string\">'MarkerSize'<\/span>,24)\r\n\r\n    text(p(1),p(2),p(3),[<span class=\"string\">'   '<\/span> ss(i).name]);\r\n<span class=\"keyword\">end<\/span>\r\n\r\nview(3)\r\ngrid <span class=\"string\">on<\/span>\r\naxis <span class=\"string\">equal<\/span>\r\n<\/pre>\r\n\r\n<img decoding=\"async\" vspace=\"5\" hspace=\"5\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/calculations_01.png\" alt=\"\"> \r\n\r\n<p>This is accurate, but not yet very helpful. Let's now calculate the geocentric position vectors of each planet. To do this, we'll put the earth at the center of the system. Copernicus won't mind, because A) he's dead, and B) we admit this reference frame orbits the sun.<\/p><p>We're also going to use another file from the File Exchange. Georg Stillfried's <a href=\"https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/25372-marrow3-m-easy-to-use-3d-arrow\">mArrow3<\/a> will help us make nice 3-D arrows in space.<\/p><pre class=\"codeinput\">clf\r\npEarth = ss(5).position;\r\n<span class=\"keyword\">for<\/span> i = 1:length(ss)\r\n    <span class=\"comment\">% pe = position relative to earth<\/span>\r\n    <span class=\"comment\">% (i.e. a vector pointing from earth to body X)<\/span>\r\n    pe = ss(i).position - pEarth;\r\n    <span class=\"comment\">% pne = normalized position relative to earth<\/span>\r\n    pne = pe\/norm(pe);\r\n    ss(i).pne = pne;\r\n\r\n    mArrow3([0 0 0],pne, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'stemWidth'<\/span>,0.01,<span class=\"string\">'FaceColor'<\/span>,[1 0 0]);\r\n\r\n    text(pne(1),pne(2),pne(3),ss(i).name, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'HorizontalAlignment'<\/span>,<span class=\"string\">'center'<\/span>);\r\n    hold <span class=\"string\">on<\/span>\r\n<span class=\"keyword\">end<\/span>\r\nlight\r\nhold <span class=\"string\">off<\/span>\r\naxis <span class=\"string\">equal<\/span>\r\naxis <span class=\"string\">off<\/span>\r\naxis([-1.2 1.2 -0.8 1.1 -0.2 0.8])\r\n<\/pre>\r\n\r\n<img decoding=\"async\" loading=\"lazy\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/calculations_021.png\" alt=\"calculations_02\" width=\"468\" height=\"423\" class=\"alignnone size-full wp-image-2701\" \/>\r\n\r\n<p>These are unit vectors pointing out from the center of the earth towards each of the other objects. It's a little hard to see, but these vectors are all lying in approximately the same plane.<\/p><p>If we change our view point to look at the system edge-on, we can see the objects are not quite co-planar. For simplicity, let's squash them all into the same plane. For this, we'll use the plane defined by the earth's velocity vector crossed with its position relative to the sun. This defines \"north\" for the solar system.<\/p><pre class=\"codeinput\">pEarth = ss(5).position;\r\npSun = ss(1).position;\r\nvEarth = ss(5).velocity;\r\n\r\nearthPlaneNormal = cross(vEarth,pSun - pEarth);\r\n\r\n<span class=\"comment\">% Normalize this vector<\/span>\r\nearthPlaneNormalUnit = earthPlaneNormal\/norm(earthPlaneNormal);\r\nmArrow3([0 0 0],earthPlaneNormalUnit, <span class=\"keyword\">...<\/span>\r\n    <span class=\"string\">'stemWidth'<\/span>,0.01,<span class=\"string\">'FaceColor'<\/span>,[0 0 0]);\r\nview(-45,15)\r\naxis([-1.2 1.2 -0.8 1.1 -0.2 0.8])\r\n<\/pre>\r\n\r\n\r\n<img decoding=\"async\" loading=\"lazy\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/calculations_031.png\" alt=\"calculations_03\" width=\"428\" height=\"384\" class=\"alignnone size-full wp-image-2702\" \/>\r\n\r\n<p>Now we project the vectors onto the plane defined by earth's motion around the sun. I learned what I needed from Wikipedia here: <a href=\"http:\/\/en.wikipedia.org\/wiki\/Vector_projection\">Vector Projection<\/a>.<\/p><p>Since we are working with the normal, we are technically doing a <a href=\"http:\/\/en.wikipedia.org\/wiki\/Vector_projection#Vector_rejection_2\">vector rejection<\/a>. Using Wikipedia's notation, this is given by<\/p><p>$$ \\mathbf{a_2} = \\mathbf{a} - \\frac{\\mathbf{a} \\cdot \\mathbf{b}}{\\mathbf{b} \\cdot \\mathbf{b}} \\mathbf{b} $$<\/p><pre class=\"codeinput\">hold <span class=\"string\">on<\/span>\r\n<span class=\"keyword\">for<\/span> i = 1:length(ss)\r\n    pne = ss(i).pne;\r\n    pneProj = pne - dot(pne,earthPlaneNormalUnit)\/dot(earthPlaneNormalUnit,earthPlaneNormalUnit)*earthPlaneNormalUnit;\r\n    ss(i).pneProj = pneProj;\r\n\r\n    mArrow3([0 0 0],pneProj, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'stemWidth'<\/span>,0.01,<span class=\"string\">'FaceColor'<\/span>,[0 0 1]);\r\n<span class=\"keyword\">end<\/span>\r\nhold <span class=\"string\">off<\/span>\r\naxis <span class=\"string\">equal<\/span>\r\n<\/pre>\r\n\r\n<img decoding=\"async\" loading=\"lazy\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/calculations_042.png\" alt=\"calculations_04\" width=\"520\" height=\"453\" class=\"alignnone size-full wp-image-2705\" \/>\r\n\r\n<p>We're close to the end now. Let's just calculate the angle between the sun and each element. Then we'll place the sun at the 12:00 position of our planar visualization and all the other planets will fall into place around it.<\/p><p>Calculate the angle between the sun and each of the bodies. Again, from the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Vector_projection#Definitions_in_terms_of_a_and_b\">Wikipedia article<\/a>, we have<\/p><p>$$ cos \\theta = \\frac{\\mathbf{a} \\cdot \\mathbf{b}}{|\\mathbf{a}||\\mathbf{b}|} $$<\/p><pre class=\"codeinput\">sun = ss(1).pneProj;\r\nss(1).theta = 0;\r\n\r\n<span class=\"keyword\">for<\/span> i = 1:length(ss)\r\n    pneProj = ss(i).pneProj;\r\n    cosTheta = dot(sun,pneProj)\/(norm(sun)*norm(pneProj));\r\n    theta = acos(cosTheta);\r\n\r\n    <span class=\"comment\">% The earth-plane normal vector sticks out of the plane. We can use it<\/span>\r\n    <span class=\"comment\">% to determine the orientation of theta<\/span>\r\n\r\n    x = cross(sun,pneProj);\r\n    theta = theta*sign(dot(earthPlaneNormalUnit,x));\r\n    ss(i).theta = theta;\r\n<span class=\"keyword\">end<\/span>\r\n<\/pre><h4>Plot the result<a name=\"691012af-279b-45e4-9237-e7a4fff3b0a4\"><\/a><\/h4><pre class=\"codeinput\">cla\r\n\r\nk1 = 1.5;\r\nk2 = 1.2;\r\n<span class=\"keyword\">for<\/span> i = 1:length(ss)\r\n    beta = ss(i).theta + pi\/2;\r\n    x = cos(beta);\r\n    y = sin(beta);\r\n    mArrow3([0 0 0],[x y 0], <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'stemWidth'<\/span>,0.01, <span class=\"keyword\">...<\/span>\r\n        <span class=\"string\">'FaceColor'<\/span>,[0 0 1]);\r\n    line([0 k1*x],[0 k1*y],<span class=\"string\">'Color'<\/span>,0.8*[1 1 1])\r\n    text(k1*x,k1*y,ss(i).name,<span class=\"string\">'HorizontalAlignment'<\/span>,<span class=\"string\">'center'<\/span>)\r\n<span class=\"keyword\">end<\/span>\r\nt = linspace(0,2*pi,100);\r\nline(k2*cos(t),k2*sin(t),<span class=\"string\">'Color'<\/span>,0.8*[1 1 1])\r\nline(0,0,1,<span class=\"string\">'Marker'<\/span>,<span class=\"string\">'.'<\/span>,<span class=\"string\">'MarkerSize'<\/span>,40,<span class=\"string\">'Color'<\/span>,[0 0 1])\r\n\r\naxis <span class=\"string\">equal<\/span>\r\naxis(2*[-1 1 -1 1])\r\n<\/pre>\r\n\r\n<img decoding=\"async\" loading=\"lazy\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/calculations_051.png\" alt=\"calculations_05\" width=\"480\" height=\"472\" class=\"alignnone size-full wp-image-2704\" \/>\r\n\r\n<p>And there you have it: an accurate map of where the planets are in the sky for today's date. In this orientation, planets \"following\" the sun through the sky are on the left side of the circle. So Jupiter will be high in the sky as the sun sets.<\/p><p>And that is a very satisfying answer to my question, by way of vector math, JSON feeds, MATLAB graphics, and the File Exchange.<\/p><script language=\"JavaScript\"> <!-- \r\n    function grabCode_cada7fe7de3f4d3b9bb4ad223d735d1e() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='cada7fe7de3f4d3b9bb4ad223d735d1e ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' cada7fe7de3f4d3b9bb4ad223d735d1e';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2014 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }   \r\n     --> <\/script><p style=\"text-align: right; font-size: xx-small; font-weight:lighter;   font-style: italic; color: gray\"><br><a href=\"javascript:grabCode_cada7fe7de3f4d3b9bb4ad223d735d1e()\"><span style=\"font-size: x-small;        font-style: italic;\">Get \r\n      the MATLAB code <noscript>(requires JavaScript)<\/noscript><\/span><\/a><br><br>\r\n      Published with MATLAB&reg; R2014a<br><\/p><\/div><!--\r\ncada7fe7de3f4d3b9bb4ad223d735d1e ##### SOURCE BEGIN #####\r\n%% Sky Clock calculations\r\n%\r\n% I look up at the sky just after sunset and I see an especially bright\r\n% star. It's probably a planet. But which one? \r\n%\r\n% This question gives me a good opportunity to play around with MATLAB.\r\n% Let's do a visualization that shows where the planets are relative to the\r\n% earth and the sun. In the process, we'll use JSON services, the File\r\n% Exchange, MATLAB graphics, and 3-D vector mathematics cribbed from\r\n% Wikipedia.\r\n%\r\n% <<explanation1.png>>\r\n%\r\n% Here is the basic grade-school illustration of the solar system, the one\r\n% that shows the planets rolling around the sun like peas on a plate. For\r\n% simplicity, we're just showing the sun, the earth, the moon, Venus, and\r\n% Mars. \r\n%\r\n% But we never see anything like this with our own eyes. Instead, we see\r\n% bright spots on a dark background somewhere \"up there.\" So let's simplify\r\n% our problem to determining what direction each naked-eye planet is in.\r\n% This leads to an image like this.\r\n%\r\n% <<explanation3.png>>\r\n%\r\n% Our goal is to make an accurate up-to-date version of this diagram.\r\n% Specifically, relative to the sun, where should we look to find the moon\r\n% and the naked-eye planets (Mercury, Venus, Mars, Jupiter, and Saturn)? To\r\n% do this, we need to solve a few problems.\r\n%\r\n% # Find the planets\r\n% # Find the unit vector pointing from earth to each planet\r\n% # Squash all these vectors onto a single plane\r\n% # Visualize the resulting disk\r\n\r\n%% Where Are the Planets?\r\n% First of all, where are the planets? There's a free JSON service for just\r\n% about everything these days. I found planetary data hosted on Davy\r\n% Wybiral's site here:\r\n%\r\n% <http:\/\/davywybiral.blogspot.com\/2011\/11\/planetary-states-api.html>\r\n\r\nurl = 'http:\/\/www.astro-phys.com\/api\/de406\/states?bodies=sun,moon,mercury,venus,earth,mars,jupiter,saturn';\r\njson = urlread(url);\r\n\r\n%% Parse the data\r\n% Now we can use\r\n% <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/42236-parse-json-text\r\n% Joe Hicklin's JSON parser> from the File Exchange. It returns a\r\n% well-behaved MATLAB structure.\r\n\r\ndata = JSON.parse(json)\r\n\r\n%%\r\n% The payload is in the \"results\" field. Each entry has three position\r\n% components and three velocity components. \r\n\r\ndata.results\r\n\r\n%%\r\n% The distances are in kilometers, and I'm not even sure how this\r\n% representation is oriented relative to the surrounding galaxy. But it\r\n% doesn't really matter, because all I care about is the relative positions\r\n% of the bodies in question.\r\n\r\n%% Aerospace Toolbox Ephemeris\r\n% Side note: We could also have used the Aerospace Toolbox to get the same\r\n% information.\r\n%\r\n% |[pos,vel] = planetEphemeris(juliandate(now),'Sun','Earth')|\r\n\r\n%% Build the Solar System Structure\r\n\r\n% List of bodies we care about\r\nssList = {'sun','moon','mercury','venus','earth','mars','jupiter','saturn'};\r\n\r\nss = [];\r\nfor i = 1:length(ssList)\r\n    ssObjectName = ssList{i};\r\n    ss(i).name = ssObjectName;\r\n    ssData = data.results.(ssObjectName);\r\n    ss(i).position = [ssData{1}{:}];\r\n    ss(i).velocity = [ssData{2}{:}];\r\nend\r\n\r\n%% Plot the planets\r\n\r\n% Plot in astronomical units\r\nau = 149597871;\r\n[x,y,z] = sphere;\r\nk = 5;\r\n \r\ncla\r\nfor i = 1:length(ss)\r\n    p = ss(i).position\/au;\r\n\r\n    line(p(1),p(2),p(3), ...\r\n        'Marker','.','MarkerSize',24)\r\n \r\n    text(p(1),p(2),p(3),['   ' ss(i).name]);\r\nend\r\n \r\nview(3)\r\ngrid on\r\nshading flat\r\naxis equal\r\nlight\r\n\r\n%%\r\n% This is accurate, but not yet very helpful. Let's now calculate the\r\n% geocentric position vectors of each planet. To do this, we'll put the\r\n% earth at the center of the system. Copernicus won't mind, because A) he's\r\n% dead, and B) we admit this reference frame orbits the sun.\r\n%\r\n% We're also going to use another file from the File Exchange. Georg\r\n% Stillfried's <https:\/\/www.mathworks.com\/matlabcentral\/fileexchange\/25372-marrow3-m-easy-to-use-3d-arrow\r\n% mArrow3> will help us make nice 3-D arrows in space.\r\n\r\nclf\r\npEarth = ss(5).position;\r\nfor i = 1:length(ss)\r\n    % pe = position relative to earth\r\n    % (i.e. a vector pointing from earth to body X)\r\n    pe = ss(i).position - pEarth;\r\n    % pne = normalized position relative to earth\r\n    pne = pe\/norm(pe);\r\n    ss(i).pne = pne;\r\n    \r\n    mArrow3([0 0 0],pne, ...\r\n        'stemWidth',0.01,'FaceColor',[1 0 0]);\r\n    \r\n    text(pne(1),pne(2),pne(3),ss(i).name, ...\r\n        'HorizontalAlignment','center');\r\n    hold on\r\nend\r\nlight\r\nhold off\r\naxis equal\r\naxis off\r\naxis([-1.2 1.2 -0.8 1.1 -0.2 0.8])\r\n\r\n%%\r\n% These are unit vectors pointing out from the center of the earth towards\r\n% each of the other objects. It's a little hard to see, but these vectors\r\n% are all lying in approximately the same plane.\r\n%\r\n% If we change our view point to look at the system edge-on, we can see the\r\n% objects are not quite co-planar. For simplicity, let's squash them all\r\n% into the same plane. For this, we'll use the plane defined by the earth's\r\n% velocity vector crossed with its position relative to the sun. This defines\r\n% \"north\" for the solar system.\r\n\r\npEarth = ss(5).position;\r\npSun = ss(1).position;\r\nvEarth = ss(5).velocity;\r\n\r\nearthPlaneNormal = cross(vEarth,pSun - pEarth);\r\n\r\n% Normalize this vector\r\nearthPlaneNormalUnit = earthPlaneNormal\/norm(earthPlaneNormal);\r\nmArrow3([0 0 0],earthPlaneNormalUnit, ...\r\n    'stemWidth',0.01,'FaceColor',[0 0 0]);\r\nview(-45,15)\r\naxis([-1.2 1.2 -0.8 1.1 -0.2 0.8])\r\n\r\n%%\r\n% Now we project the vectors onto the plane defined by earth's motion\r\n% around the sun. I learned what I needed from Wikipedia here:\r\n% <http:\/\/en.wikipedia.org\/wiki\/Vector_projection Vector Projection>.\r\n%\r\n% Since we are working with the normal, we are technically doing a\r\n% <http:\/\/en.wikipedia.org\/wiki\/Vector_projection#Vector_rejection_2 vector\r\n% rejection>. Using Wikipedia's notation, this is given by\r\n%\r\n% $$ \\mathbf{a_2} = \\mathbf{a} - \\frac{\\mathbf{a} \\cdot \\mathbf{b}}{\\mathbf{b} \\cdot \\mathbf{b}} \\mathbf{b} $$\r\n\r\nhold on\r\nfor i = 1:length(ss)\r\n    pne = ss(i).pne;\r\n    pneProj = pne - dot(pne,earthPlaneNormalUnit)\/dot(earthPlaneNormalUnit,earthPlaneNormalUnit)*earthPlaneNormalUnit;\r\n    ss(i).pneProj = pneProj;\r\n    \r\n    mArrow3([0 0 0],pneProj, ...\r\n        'stemWidth',0.01,'FaceColor',[0 0 1]);\r\n%     text(pneProj(1),pneProj(2),pneProj(3),ss(i).name);\r\nend\r\nhold off\r\naxis equal\r\n\r\n%%\r\n% We're close to the end now. Let's just calculate the angle between the\r\n% sun and each element. Then we'll place the sun at the 12:00 position of\r\n% our planar visualization and all the other planets will fall into place\r\n% around it.\r\n% \r\n% Calculate the angle between the sun and each of the bodies. Again, from\r\n% the\r\n% <http:\/\/en.wikipedia.org\/wiki\/Vector_projection#Definitions_in_terms_of_a_and_b\r\n% Wikipedia article>, we have\r\n%\r\n% $$ cos \\theta = \\frac{\\mathbf{a} \\cdot \\mathbf{b}}{|\\mathbf{a}||\\mathbf{b}|} $$\r\n\r\nsun = ss(1).pneProj;\r\nss(1).theta = 0;\r\n\r\nfor i = 1:length(ss)\r\n    pneProj = ss(i).pneProj;\r\n    cosTheta = dot(sun,pneProj)\/(norm(sun)*norm(pneProj));\r\n    theta = acos(cosTheta);\r\n    \r\n    % The earth-plane normal vector sticks out of the plane. We can use it\r\n    % to determine the orientation of theta\r\n    \r\n    x = cross(sun,pneProj);\r\n    theta = theta*sign(dot(earthPlaneNormalUnit,x));\r\n    ss(i).theta = theta;    \r\nend\r\n\r\n%% Plot the result\r\n\r\ncla\r\n\r\nk1 = 1.5;\r\nk2 = 1.2;\r\nfor i = 1:length(ss)\r\n    beta = ss(i).theta + pi\/2;\r\n    x = cos(beta);\r\n    y = sin(beta);\r\n    mArrow3([0 0 0],[x y 0], ...\r\n        'stemWidth',0.01, ...\r\n        'FaceColor',[0 0 1]);\r\n    line([0 k1*x],[0 k1*y],'Color',0.8*[1 1 1])\r\n    text(k1*x,k1*y,ss(i).name,'HorizontalAlignment','center')\r\nend\r\nt = linspace(0,2*pi,100);\r\nline(k2*cos(t),k2*sin(t),'Color',0.8*[1 1 1])\r\nline(0,0,1,'Marker','.','MarkerSize',40,'Color',[0 0 1])\r\n\r\naxis equal\r\naxis(2*[-1 1 -1 1])\r\n\r\n%%\r\n% And there you have it: an accurate map of where the planets are in the\r\n% sky for today's date. In this orientation, planets \"following\" the sun\r\n% through the sky are on the left side of the circle. So Jupiter will be\r\n% high in the sky as the sun sets.\r\n%\r\n% And that is a very satisfying answer to my question, by way of vector\r\n% math, JSON feeds, MATLAB graphics, and the the File Exchange.\r\n\r\n##### SOURCE END ##### cada7fe7de3f4d3b9bb4ad223d735d1e\r\n-->","protected":false},"excerpt":{"rendered":"<div class=\"overview-image\"><img decoding=\"async\"  class=\"img-responsive\" src=\"https:\/\/blogs.mathworks.com\/community\/files\/explanation1.png\" onError=\"this.style.display ='none';\" \/><\/div><!--introduction--><p>I look up at the sky just after sunset and I see an especially bright star. It's probably a planet. But which one?... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/community\/2014\/05\/14\/whats-that-in-the-sky-build-a-matlab-planetarium\/\">read more >><\/a><\/p>","protected":false},"author":69,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2655"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/users\/69"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/comments?post=2655"}],"version-history":[{"count":20,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2655\/revisions"}],"predecessor-version":[{"id":2789,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/posts\/2655\/revisions\/2789"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/media?parent=2655"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/categories?post=2655"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/community\/wp-json\/wp\/v2\/tags?post=2655"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}