8

Looking at the Sky

In this chapter we give a brief introduction to astronomy using some of the functions available in Mathematica. In the process we will also learn to model the orbits of eclipsing binary stars by creating and analyzing light curves, very useful when searching for extrasolar planets.

8.1 A Short Astronomical Walk

Contemplating the sky on a starry night, away from light pollution, is one of the most fascinating shows that we can enjoy. It’s not strange then that for millennia humans have looked up at the sky and realized that the sun and the moon follow regular cycles. The majority of ancient cultures came up with mythological or religious explanations about celestial objects but they also used those objects for practical purposes such as the creation of calendars to schedule harvests. For example, the Egyptians associated the appearance of Sirius, the brightest star in the norther hemisphere (now we know it’s a binary system), with the annual flooding of the Nile. Despite the reduced number of celestial bodies that can be seen with the naked eye, it’s surprising how much knowledge of the sky that astronomers already had before the invention of the telescope.

The arrival of the telescope about 400 years ago (1609) heralded a revolution: the number of known stars went from thousands to millions, the moon, far from being a perfect sphere, was found to be covered with mountains and craters, the wandering stars (the planets) were more than mere bright dots, they actually had satellites, and Saturn had rings (it’s not the only planet). Later on we learned to analyze the composition of the light coming from space. Many of the objects that looked like blurry spots under the telescope became galaxies with myriads of stars. It was discovered that those galaxies were moving apart from each other: the universe was expanding, an expansion that started around 14 billion years ago. It turned out that the visible stars and galaxies contain less than 5% of the mass of the Universe. The rest is dark energy and matter, we know it exists (at least the dark matter) but we don’t really know what it is. The idea of the existence of dark energy (“a riddle wrapped inside a mystery”) started with an amazing discovery (1998): The expansion of the universe is accelerating, a discovery that was worth a Nobel Prize in 2011. And the discovery race continues, promising new surprises: the important thing is still the journey.

For astronomical calculations, Mathematica has several functions: CometData, ExoplanetData, GalaxyData, MinorPlanetData, MoonPosition Sunset, NebulaData, PlanetData, PulsarData, SiderealTime, StarClusterData, StarData, SunPosition, SupernovaData, Sunrise and many more. These functions replace and extend the AstronomicalData function available until Mathematica 9 although you can still use it with versions 10 and 11.

The syntax of all these functions is similar to that of the other computable functions. Example: Function ["name", "property"] where "name" or entity (Entity) is an object and " property" is a characteristic or parameter of the object.

Dates and hours are of crucial importance in Astronomy. The existence throughout history, even nowadays, of different calendars makes it difficult to compare dates of astronomical and historical events. To avoid this problem we use Julian dates. A Julian date (not to be confused with the Julian calendar: Julian year) indicates the days that have gone by since January 1, 4713 B.C. at 12:00 PM UT (November 24, 4714 BC, in the proleptic Gregorian calendar). The name of the Julian date comes from Julius Scaliger (not Julius Caesar) who proposed it in 1583 after realizing that such a day was the least common multiple of 3 calendar cycles used by the Julian calendar. Using 12:00 PM as the starting time instead of 12:00 AM has the advantage that astronomical observations, usually done at night, are recorded on the same day. Starting with Mathematica 10.2 the function JulianDate is available. The function converts dates from our calendar, the Gregorian one, to Julian dates. It takes into account the modifications to the Julian calendar done by Gregorio XIII (1582) that improved the handling of leap years, but in which the year 0 doesn’t exist.

We can stop the dynamic updating by unchecking it in the Evaluation menu: Evaluation ►Dynamic Updating Enabled.

Star positions in the sky don’t follow an annual cycle exactly. They are subject to perturbations due to their own motions and other celestial movements such as the precession of equinoxes. Because of this, it’s necessary to regularly set up dates to use when referring to celestial coordinates. These standard reference dates are known as epochs. The last one was set on January 1, 2000 at 11:58:55.816 AM. When referring to this epoch we use J2000.0. StarData and other astronomical functions refer to the position of stars using this latest epoch without taking into account either the precession of equinoxes or the stars’ own motions. For practical purposes these adjustments are irrelevant but we need to consider them when making historical projections.

8.2 Stargazing

 
Clear ["Global` *"]

An incredible experience is to look at the sky away from light-polluted cities, trying to imagine the hidden celestial mechanism generating the incredible view. We will soon realize the greatness of those night owl geniuses (Hipparchus, Ptolemy, Copernicus, Brahe, Kepler and others) that established the foundations of modern astronomy without any other means than their plain sight and basic instruments such as the astrolabe to measure the stars’ positions. Everything changed in 1609 when Galileo used the telescope for astronomical purposes for the first time. Another equally important invention was the pendulum clock that enabled the measurement of time with great precision. Centuries later, telescopes mounted on satellites and atomic clocks heralded a new revolution.

8.2.1 Naked-Eye Astronomy

When looking at the sky, almost all of the objects that can be seen with the naked eye are stars. Throughout history, different cultures imagined that some stars resembled figures, named constellations. The same phenomenon occurred in civilizations that were not connected in any way. For example: The seven brightest stars from the Ursa Major were interpreted the same way in very different places. The inhabitants of the British islands said that those stars were the legendary King Arthur’s cart. For the Germans, they represented a wagon pulled by three horses. The Greeks came up with a more imaginative story: In a Greek legend, the god Zeus and the mortal Callisto had a son named Arcas. Hera, Zeus’ jealous wife, turned Callisto into a female bear. Arcas, unaware that the bear was actually her mother almost killed her. Afterward, Zeus also turned Arcas into another female bear. Callisto is Ursa Major and Arcas is Ursa Minor. Greek mythology is full of such stories to explain almost all the stars that appear near each other in the sky. Many of those myths have actually given names to the constellations. Over the centuries, many constellations were added until they eventually covered the entire visible sky from any place on Earth. There are 88 constellations in total. In 1930, the International Astronomical Union established the limits of each of them. Obviously, they are imaginary lines covering the entire firmament, including both hemispheres. As in maps, we don’t need to keep all the constellations’ details. They and the names of their most important stars are represented in a planisphere, a celestial map that shows us the sky each night for a given date and time.

  • The function below returns the sky chart for our location at the time given in the input. We can get Mathematica to create it for us by typing “Sky chart at 22 GMT” in free-form input, right-clicking on the sky chart and selecting Paste input for ► Subpod content.

    WolframAlpha["Sky chart at 22 GMT", {{"SkyMap:PlanetData", 1}, "Content"}]

    images

  • With the instruments and methods available to astronomers nowadays, the information that we can obtain about stars is very extensive. In the next example we show the properties most commonly used when describing stars:

    properties = {"AbsoluteMagnitude", "ApparentMagnitude",
       "RightAscension", "Declination", "Altitude", "Azimuth",
       "Constellation", "DistanceFromSun", "SpectralClass"};

  • Let’s find those properties for Polaris, the North Star:

    StarData["Polaris", properties, "PropertyAssociation"]

    images

The apparent magnitude gives how bright an astronomical object appears to an observer on regardless of its intrinsic brightness (http://scienceworld.wolfram.com/astronomy/ApparentMagnitude.html).

The smaller the magnitude, the bigger its luminosity. For the most luminous stars we use negative values. Most of the people looking at the sky with the naked eye and with little light pollution just see stars with apparent magnitude of less than 5.

  • Here are the 10 closest stars to Earth, many of them not visible to the naked eye:

    StarData[EntityClass ["Star", "StarNearest10"]]

    images

    Naturally, the closest one is the Sun. Usually for solar system bodies distances are shown by default in astronomical units (au), the average distance between the Earth and the Sun. The Earth rotates around the Sun following an ellipse, so its distance changes, although as its eccentricity is very small, in many calculations we can assume that it follows a circle. Its average value is 1 au.

  • The command below returns the name of the star along with its distance to earth in (ly), except for the Sun, and its apparent magnitude.

    {#, StarData[#, "DistanceFromEarth"],
       StarData[#, "ApparentMagnitude"]} & /@
    StarData[EntityClass["Star", "StarNearest10"]]

    images

    We can immediately notice that most of them are not visible since their apparent magnitude is bigger than 5. With the exception of the Sun, the star with the biggest apparent magnitude is Sirius, a star that we can see during winter nights (in the Northern hemisphere) next to the Orion constellation. The closest one is Proxima Centauri that is only visible in the Southern hemisphere (Proxima Centauri rotates around the double system Centaurus A and B in a 500,000-year period. By chance, it’s currently in its orbital position closest to the solar system).

  • The number of stars visible without the use of technology is very small compared to the number of stars in our galaxy. The command below generates a plot representing the naked-eye stars and their distances to Earth in light years. (This type of query may take a long time to execute).

    ListPlot[
      [Sort [StarData [StarData [EntityClass ["Star", "NakedEyeStar"], "Entities"],
         "DistanceFromSun"]], AxesLabel →{"stars", "ly"}]

    images

    From the graph above we can see that almost all of the visible stars are less than 1,500 ly away. Our galaxy, the Milky Way, has a width of approximately 100,000 ly. Therefore, with our bare eyes we only see a very small fraction of the stars in our own galaxy, most of them the ones closest to Earth. The previous figure refers to the number of potentially visible stars without light pollution or moonlight. In practice, the actual number is 2 or 3 thousand, since we only see part of the sky depending on our position.

Besides stars with apparent magnitude 5 or less, the only other celestial bodies visible using only our eyesight are: the Moon, the five nearest planets and, sporadically, some comets. In the Northern hemisphere, the only object visible that doesn’t belong to our galaxy is the Andromeda galaxy. In short, we see an insignificant fraction of our own galaxy, which in turn is just one among billions of galaxies.

We’ve seen that among the properties of StarData and other astronomical functions are: “RightAscension”, “Declination”, “Altitude” and ”Azimuth”. These properties are commonly used to indicate the position of a celestial body. In equatorial coordinates, “RightAscension” and “Declination” are similar to latitude and longitude but they refer to the celestial sphere so they are independent of the observer (for further details see: The Celestial Sphere, http://demonstrations.wolfram.com/TheCelestialSphere, by Jeff Bryant). Another type of coordinates widely used are the alt-azimuthal (alt/az) that depend on the observer: The “Altitude”(alt) is the height of the star over the horizon. It goes from 0º to 90º and has a positive sign for stars located above the horizon and a negative sign for the ones located below it. The ”Azimuth” (az) is the arc in the horizon measured counterclockwise from the South point until the object’s vertical. Its value ranges from 0º to 360º.

  • The next example shows Sirius’ displacement on the first of each month at 20:00 h using alt/az coordinates.

    sirius = Table[QuantityMagnitude[
          AstronomicalData["Sirius", {#, {2017, i, 1}, {20, 00}}]] & /@
        {"Azimuth", "Altitude"}, {i, 1, 12}]
     
    {{138.2, 42.32}, {179.3, 53.3}, {217.9, 44.60}, {241.2, 22.46},
     {253.4, -3.222}, {262.1, -31.27}, {270.6, -59.0}, {329.6, -86.2},
     {87.8, -62.7}, {96.8, -35.01}, {105.3, -6.85}, {116.7, 19.11}}
     
    Graphics [{Orange, Point[sirius]}, Frame → True,
     FrameLabel→{"azimuth", "altitude"}, AspectRatio → 1/ 2]

    images

    As we can see, there are negative values that correspond to the period of the year in which Sirius falls below the horizon and as a consequence is not visible. As a matter of fact, the Egyptians considered that the appearance of the star marked the beginning of a new year.

  • In the previous example we didn’t specify the observer’s location. Remember that in cases like that the function uses the location given by FindGeoLocation.

    FindGeoLocation[]
    GeoPosition[{40.42, -3.68}]

  • However, it’s important to keep in mind that the function actually gives us the position of the server that we are using to connect to the Internet based on the IP address. Sometimes it may be far away from our real position. This would be the case if, for example, we were using a phone to surf the net. We can check it by showing the obtained position in a map:


    GeoGraphics[
     GeoVisibleRegionBoundary[FindGeoLocation []], GeoZoomLevel -> 12]

    images

Naked-eye stars display an amazing regularity. It’s possible to forecast the location of each star without prior knowledge of celestial mechanics. If we observe a star at a certain hour, for example at 24:00 h, ignoring the displacement due to the inclination of the Earth’s axis with respect to the ecliptic, the next day it will be on the same spot four minutes earlier, that is at 23:56 h, completing a cycle in a year. The explanation that the Greek civilization and others provided for this fact was to assume that all the stars (which they called fixed, to distinguish them from the planets or wandering stars that did not exhibit such behavior) were glued to a dome that rotated with a daily cycle of 23 h and 56 min, the sidereal day.

images

1 sidereal days

  • We can use UnitConvert to find out the actual duration. When executing the command above, the predictive interface will do the conversion automatically.

    UnitConvert[1 sidereal day, MixedRadix ["Hours", "Minutes", "Seconds"]]
     
     23 h 56 min 4.09054 s

However, over very long periods we can notice discrepancies that at first sight may seem insignificant during a person’s lifetime but become quite obvious after several generations. Hipparchus of Nice (c. 190 BC – c. 120 BC) compared the stellar maps available at the time (celestial cartography is over 2,500 years old!) and realized that there was a relative movement of the stars with respect to the ecliptic. This movement is nowadays known as the precession of the equinoxes. This is probably his most famous discovery. The displacement happens when the Earth’s axis, moving along a circumference with respect to the ecliptic, rotates with a period of 25,771 years.

An excellent interactive demonstration (Figure 8.1) showing this effect can be downloaded from http://demonstrations.wolfram.com/PrecessionOfTheEquinoxes.

images

Figure 8.1     Precession of the Equinoxes demonstration.

As a result, the North celestial pole keeps on moving. Today it is close to the North Star but 4,800 years ago it was pointing toward Thuban (Alpha Draconis). William Shakespeare did not realize that, when in his play Julius Caesar stated: “But I am constant as the northern star, Of whose true-fix’d and resting quality There is no fellow in the firmament”.

In reality, all the stars seen from the Earth have a slow displacement motion as a result of the precession and their own orbits in the galaxy. This last type of movement is known as proper motion.

  • We can find the proper motion of the North Star as follows:

    images

     11.75 mas/yr

    “mas/yr” means milliarc seconds per year. It’s a value that cannot be observed during a lifetime. We need to keep in mind that a telescope on the ground can rarely see details smaller than 1 arc second. However, during long periods even the shapes of the constellations change.

Astrology (a pseudoscience) attributed predictive powers to the constellations. For example, they were supposed to determine the future of people born under them. The firmament was divided into 12 signs corresponding approximately to the number of lunar cycles in a year. Each division was assigned a symbol called a zodiac sign. The names of the zodiac are associated with the first constellations named by the Greeks during the 5th century BC, even though the Babylonians were the first ones to mentioned them 4,000 years ago. The starting point was the constellation pointing toward Aries (the moment when the Sun moves from the south celestial hemisphere to the north one coinciding with the spring equinox) in the 5th century BC. Nowadays, we still use the same division although the constellations are now in a different location to where they were 2,500 years ago. The real astronomical zodiac dates correspond to the constellation located behind the solar disk during the spring equinox, in the direction opposite to the Earth’s direction (the Sun actually passes through 13 constellations and not 12, the 13th one is Ophiuchus). The zodiac signs must be adjusted by a month to adapt to the current astronomical reality, so chances are you may have to revise your zodiac sign (although it wouldn’t be very useful).

  • The natural language input below returns the location of the Ursa Minor stars in the year 100,000 compared to their current position (dashed line). The rest of the stars experience similar changes.

    images

8.2.2 Wandering Planets

Since antiquity, people have realized that there are a small number of stars exhibiting a behavior different from the one of the fixed stars and their annual cycles. These stars were once called wandering stars. Nowadays we know that they correspond to the 5 visible planets: Mercury, Venus, Mars, Jupiter and Saturn (Uranus was not added to the list of planets until 1783 given how difficult it is to observe it with the naked eye).

To get information about planets we use the function PlanetData but before we learn more about this function, let’s see an example first available in the Wolfram Cloud to showcase the use of Dataset and CloudGet. Remember that a dataset has a hierarchical structure and that this general structure is quite common in practice. In this case, the dataset (http://wolfr.am/7FxLgPm5) contains information about the mass and radius of each planet and its corresponding moons.

  • Let’s import the dataset:

    planets = CloudGet["http://wolfr.am/7FxLgPm5"];

  • Here we show the moons of Mars:

    planets["Mars", "Moons"]

     

    Mass

    Radius

    Phobos

    1.072 ×1016 kg

    11.1 km

    Deimos

    1.5 ×1015 kg

    6.2 km

  • The next command displays the number of moons for each planet:


    planets[All, "Moons", Length]

    Mercury

    0

    Venus

    0

    Earth

    1

    Mars

    2

    Jupiter

    63

    Saturn

    61

    Uranus

    27

    Neptune

    13

  • Another way to obtain the Uranus moons is:

    Length /@ PlanetData[PlanetData[], "Satellites", "EntityAssociation"]

    images

  • Note that there is a discrepancy in the number of Neptunian satellites. In fact, PlanetData uses curated information but in this case we are using http://wolfr.am/7FxLgPm5 as an example to show how to import a dataset. For further information we can check Wikipedia:

    TextSentences[WikipediaData ["Moons of Neptune"]][[ ;; 2]]
     
    {Neptune has 14 known moons, which are
      named for minor water deities in Greek mythology.,
     By far the largest of them is Triton, discovered by William
      Lassell on October 10, 1846, just 17 days after the
      discovery of Neptune itself; over a century passed before
      the discovery of the second natural satellite, Nereid.}

As we mentioned previously, if we periodically observe a planet with the naked eye at the same time and from the same place and we record (sometimes we can take pictures, as we saw with the example of the Sun) its location with respect to the horizon, the resulting figure is called an analemma. In Chapter 5 we built the one corresponding to the Sun. Here we are going to do the same for Mars and Venus.

  • We need first to define our position (latitude and longitude).

    zone = "Location" -> GeoPosition[{40.96, -5.66}];

  • Next, we display the Mars and Venus positions at 20:00 h each 5 days for the first 90 days of 2017.

    analemmamars = Table[{
        QuantityMagnitude [PlanetData ["Mars", EntityProperty["Planet", "Azimuth",
            {"Date" → DateObject[DateList[{2017, 1, i, 20}]], zone}]],
          "AngularDegrees"], QuantityMagnitude[
          PlanetData ["Mars", EntityProperty["Planet", "Altitude",
            {"Date" → DateObject[DateList[{2017, 1, i, 20}]], zone}]],
          "AngularDegrees"]}, {i, 1, 90, 5}];
    analemmavenus = Table[{
        QuantityMagnitude[
         PlanetData ["Venus", EntityProperty["Planet", "Azimuth",
           {"Date" → DateObject[DateList[{2017, 1, i, 20}]], zone}]],
         "AngularDegrees"], QuantityMagnitude [
         PlanetData ["Venus", EntityProperty["Planet", "Altitude",
           {"Date" → DateObject[DateList[{2017, 1, i, 20}]], zone}]],
         "AngularDegrees"]}, {i, 1, 90, 5}];

  • With GraphicsRow we can display both graphs in the same row. The graph for Mars has been modified using AspectRatio to make its visualization easier.

    GraphicsRow[{
     Graphics [{Orange, Point[analemmamars]},
      Frame→ True, FrameLabel →{"azimuth", "altitude"},
      AspectRatio→ 1 /2, PlotLabel → Style["Mars Analemma", Bold]],
     Graphics [{Orange, Point[analemmavenus]}, Frame → True,
      FrameLabel→{"azimuth", "altitude"},
      PlotLabel→ Style["Venus Analemma", Bold]]}, Frame → All]

    images

    The points with negative altitude correspond to those times when the planet is below the horizon and therefore invisible. We cannot see the planet when it’s behind the Sun either.

We can also see that in the case of Venus, the planet is only visible a few degrees over the horizon. This is because Venus is an interior planet and as such, when observed from the Earth, it will never be very high above the horizon. This means that we will never be able to see it next to the zenith since the sunlight would blind us. Obviously, this doesn’t happen with exterior planets such as Mars.

  • The graph below shows the orbits of the terrestrial planets or inner planets: Mercury, Venus, the Earth and Mars, with the Sun (not in scale) in the center.

    images

    InputForm[%]
    EntityClass["Planet", "InnerPlanet"]
    terrestialPlanets=
      PlanetData[EntityClass ["Planet", "InnerPlanet"], "OrbitPath"];
    Graphics3D
     {{Yellow, Sphere[{0, 0, 0}, 0.1]}, terrestialPlanets}, Boxed → False]

    images

  • The next four functions, based on examples included in the PlanetData help page, enable us to calculate and display graphically the distance from Earth and the apparent magnitude of the visible planets over a year:

    visibleplanets[planet _, date _, property _] :=
     PlanetData [planet, EntityProperty ["Planet", property, {"Date" → date}]]
     
    ts[planet _, d1 _, d2_, property _] := TimeSeries[
      {#, visibleplanets[planet, #, property]} & /@ DateRange[d1, d2, "Week"]]
     
    DateListPlot[
     AssociationMap[ts [#, {2017, 1, 1}, {2017, 12, 31}, "DistanceFromEarth"] &,
      {"Mercury", "Venus", "Mars", "Jupiter", "Saturn"}],
     PlotLegends→ "Expressions", FrameLabel → Automatic]

    images

    DateListPlot [
     AssociationMap[ts [#, {2017, 1, 1}, {2017, 12, 31}, "ApparentMagnitude"] &,
      {"Mercury", "Venus", "Mars", "Jupiter", "Saturn"}],
     PlotLegends→ "Expressions", FrameLabel → Automatic]

    images

Even though naked-eye observation and the use of angle-measuring instruments (sextant, astrolabe, etc.) advanced the science of Astronomy, the arrival of the telescope and other instruments such as precision clocks was no less important. They opened new possibilities for exploring the firmament. Later on, in the 20th century, radio telescopes and spacial astronomy would arrive taking the science to a whole new level.

8.2.3. Dwarf Planets

We’ve seen that PlanetData doesn’t include Pluto, demoted to the “Dwarf Planet” category in a controversial meeting of the International Astronomical Union on August 24, 2006. This action was motivated mainly by the discovery of planets beyond Pluto, in an area known as the Kuiper Belt that probably includes thousands of planetoids. Pluto is considered to be part of this belt. All the planets that didn’t fit the new definition were categorized as dwarf planets. There are reasons to justify that Pluto doesn’t belong to the same category as the classical planets but the dwarf label doesn’t seem the most adequate one since there are dwarf planets that are most likely bigger than the classical planet Mercury.

  • The celestial bodies that orbit around the Sun but are not planets have been included in MinorPlanetData.

    images

  • Next, we use a couple of functions to know how long it would take Pluto to complete one orbit starting in 1970-01-01 and how far away in astronomical units (au) it would be from the Earth each year.

    MinorPlanetData["Pluto", "OrbitPeriod"]
    247.92065 Julian years

    DateListPlot[
    {{#, 01, 01}, MinorPlanetData["Pluto", EntityProperty["MinorPlanet",
    "DistanceFromEarth", {"Date" → #}]]} & /@
    Range[1970, 1970 + 248, 10]]

    images

    Notice that during our lifetimes Pluto will be further away each year. In July 2015 a probe, New Horizons, visited it for the first time.

  • One of the most interesting characteristics of Pluto and other celestial bodies in the Kuiper Belt is that their orbits are normally slanted with respect to the ecliptic plane, next to the ones from the classical planets, as we can see with the following command:

    images

    In the coming years the number of newly discovered dwarf planets will most likely increase substantially.

  • We dynamically simulate the movement of the dwarf planets starting on January 1, 2017 using MinorPlanetData. Notice the use of Tooltip to see the name when placing the mouse over a planet. The function is slow.

    minorplanets = MinorPlanetData[EntityClass["MinorPlanet", "Plutoid"]];
    Manipulate[
     Graphics3D[{(Tooltip[Sphere[#[[2]], 1], #[[1]]]) & /@
      Transpose[{minorplanets,
       QuantityMagnitude[
            MinorPlanetData[EntityClass ["MinorPlanet", "Plutoid"],
       EntityProperty ["MinorPlanet", "HelioCoordinates",{"Date" ->
        DateObject[{2017, 1, date}]}]]]}], ColorData[1, 1],
       MinorPlanetData[#, "OrbitPath"] & /@ minorplanets},
    PlotLabel -> DateString[DateList[{2017, 1, date}]]],
     {date, 1, 300 × 365.25}, SaveDefinitions → True]

    images

8.2.4 Exoplanets

To obtain information about planets outside of the solar system we can use the function ExoplanetData.

  • According to a paper published in Nature on August 25, 2016, there’s an Earth-sized planet orbiting our nearest neighboring star other than the sun, Proxima Centauri.

    ExoplanetData[EntityClass ["Exoplanet",  {EntityProperty ["Exoplanet",
    "DistanceFromEarth"]  -> TakeSmallest[1]}]]

    images

  • Next all the available information using ExoplanetData regarding this planet as of 2016-12-04 is shown:

    images

Although the planet orbits (0.04 au) closer to its star than Mercury does to the Sun, the star itself is far fainter than the Sun. The surface temperature would allow the presence of liquid water. However, the conditions on the surface may be strongly affected by the ultraviolet and X-ray flares from the star far more intense than the Earth experiences from the Sun.

8.2.5 Galaxies and Nebulae

Nowadays, we know that the Sun is just one of the stars in a galaxy known as the Milky Way. However, the notion of galaxy was not introduced until the 1920s. Before then, the difference between galaxies and nebulae was not clear. Seen through the telescope both appeared like blurry objects with shapes resembling clouds. That’s why the term nebulae was used to refer to both. Galaxies were considered a type of nebula and our galaxy was the entire Universe.

  • Let’s see the definition of galaxy:

    images

    1

    noun

    a splendid assemblage (especially of famous people)

    2

    noun

    tufted evergreen perennial herb having spikes of tiny white flowers and glossy green round to heart-shaped leaves that become coppery to maroon or purplish in fall

    3

    noun

    (astronomy) a collection of star systems; any of the billions of systems each having many stars and nebulae and dust

Based on the 3rd definition, nebulae are just part of galaxies. With telescopes we can only see the ones located in our own galaxy. The generic name nebulae actually includes two types of very different structures: the planetary nebulae (“PlanetaryNebula”) generated by the accumulation of dust left after a supernova explosion, and nebulae (“Nebula”) star nurseries, such as the Orion Nebula. As we can see the name is confusing as planetary nebulae have nothing to do with the formation of planets. Other very important galaxy elements are clusters or globular clusters (groups of stars). Information about galaxies, clusters and nebulae can be found using GalaxyData, StarClusterData and NebulaData.

  • A very interesting nebula is Orion, visible specially during Winter in the northern hemisphere. It is an enormous gas cloud where new stars are constantly being born. The command below displays some of its properties (“//Short” should be removed to see the complete output).

    DeleteMissing[
      Transpose [{NebulaData["Properties"], NebulaData["OrionNebula", #] & /@
        NebulaData["Properties"]}], 1, 1] // Short

    images

  • We can find out the diameter of our galaxy. To express it in terms of light-years we use UnitConvert.

    {GalaxyData ["MilkyWay", "Diameter"],
     UnitConvert [GalaxyData ["MilkyWay", "Diameter"], "LightYears"]}
     
    {9.5 × 1017 km, 1.0 × 105 ly}

  • We’ve seen that big distances, such as the ones between galaxies, are usually measured in parsec (pc), the equivalent of the distance corresponding to one arc second. The closest galaxies to our own are shown below, note the use of Dataset. They all form a group of galaxies known as the local group. Many of them are quite small.

    images

  • We select the galaxies closer than 100 kpc:

    lggaxies[Select[# < 100 kpc &]]

    Boötes Dwarf Galaxy

    60.3591 kpc

    Canis Major Dwarf Galaxy

    7.65979 kpc

    Draco Dwarf Galaxy

    81.8065 kpc

    Large Magellanic Cloud

    49.9418 kpc

    Milky Way

    7.61113 kpc

    Sagittarius Dwarf Elliptical Galaxy

    19.9154 kpc

    Sculptor Dwarf Galaxy

    79.049 kpc

    Sextans Dwarf Galaxy

    85.7896 kpc

    Small Magellanic Cloud

    60.5736 kpc

    Ursa Minor Dwarf Galaxy

    65.8742 kpc

    Willman I

    36.767 kpc

  • To check which galaxies are visible to the naked eye (Apparent Magnitude < 5) we sort the apparent magnitudes of the local group members (remember that the smaller the magnitude the higher the luminosity).

    images

    Large Magellanic Cloud

    0.9

    M31

    3.5

    Sagittarius Dwarf Elliptical Galaxy

    4.5

    Small Magellanic Cloud

    2.2

As mentioned earlier, with our own eyes it is very difficult to see magnitudes greater than 5. This means that in the northern hemisphere the only visible galaxy is M31 (Andromeda), the rest of the galaxies with high luminosity (and correspondingly low apparent magnitude) are only visible from the southern hemisphere.

8.3 Application: Determining the Color of the Stars

Clear["Global` *"]

The French philosopher Auguste Comte (1798–1857) believed that the composition of stars was going to be beyond human knowledge. In his Cours de Philosophie he used to say that stars would only be known as specks of light in the sky since they were very far away. However, the analysis of stellar light has enabled us to know reasonably well stars’ composition, temperature and evolution.

A star is what in Physics is known as a black body and therefore has a light spectrum whose characteristics are related to its surface temperature following Planck’s law.

We are going to build a dynamic representation of Planck’s law that shows the relationship between color and temperature.

A more complete model by Jeff Bryant can be found in:

"Blackbody Spectrum" (http://demonstrations.wolfram.com/BlackbodySpectrum).

8.4 The Measurement of Distances across the Universe

Clear ["Global` *"]

Historically, the measurement of distances across the Universe was one of the fundamental problems in Astronomy. The first step toward solving it was to measure the Earth’s size but before even that, someone had to realize that our planet is a sphere. It is believed that it was the Pythagoreans, mysterious people that didn’t leave any written works behind, who for the first time, circa 430 B.C., came up with the idea that the Earth was round, based probably on two observations: a) Sailors traveling from Greece to Northern Africa had noticed that several constellations, such as Ursa Major, would appear higher in the horizon in Egypt than in Athens for a given date and time (although there are other geometric shapes that may explain this fact, the sphere is the one that provides the simplest explanation), b) The way the Moon hides during lunar eclipses can be easily explained if one assumes that it’s a shadow projected by a sphere, the Earth, upon the Moon.

Once it was clear that the Earth was round, the next step was its measurement. It was Eratosthenes, who eventually became the head of the library of Alexandria, the first person to measure the Earth’s size. The method he used is legendary and well-known among Astronomy enthusiasts, but it’s worth describing it in the next few lines. He started measuring the distance between Alexandria and Syene (nowadays Aswan), that were located approximately on the same meridian (they really differ by 3º). He also assumed that the Sun was distant enough that its rays would hit the Earth in a parallel manner. He knew that during the summer solstice in Syene the light would illuminate the bottom of the wells in the city. At that moment the sunbeams would hit the ground perpendicularly while at the same time in Alexandria they would form an angle of approximately 1/50th of the length of a circumference (number deduced by the projection of the shadows). He then used a basic trigonometric identity to figure out that the distance between Alexandria and Syene was 1/50th of the Earth’s circumference. Since that distance was 5,000 stadia, the Earth’s circumference was 250,000 stadia. There’s no agreement about the exact conversion rate to meters. Depending on the reference, if we take a value ranging from 185 m to 157.2 m, the error is between 17% and 1% of the actual size, quite an extraordinary approximation in any case. In 1492, Christopher Columbus assumed a considerably smaller figure while trying to convince Queen Isabella of Spain to finance his trip against the opinion of her advisors, who knew Eratosthenes’ estimate. He persisted in his error and ended up discovering America.

Shortly before Eratosthenes’ death, Hipparchus of Nice was born (c. 190 BC). He’s not only famous due to his discovery of the precession of equinoxes, as we have already mentioned, but also because of his quite accurate measurement of the distance between the Earth and the Moon using the shadow projected by the former onto the later during Moon eclipses. He also initiated an empirical study of the apparent magnitude of the stars and their positions. For that purpose, he invented the ecliptic coordinates (different from the elliptic ones) and used the astrolabe (invention attributed to him). Many centuries later (from 1989 to 1993) a satellite, carrying its name, repeated his work with technology from the late 20th century. Its successor: Gaia (http://www.esa.int/Our_Activities/Space_Science/Gaia_overview), currently in mission, will measure the position of the objects in our galaxy with extreme precision.

We can compute the distance to the nearest stars using geometric methods (Figure 8.2). We measure the position of the star seen from a certain location and 6 months later we will observe its apparent displacement. This phenomenon is known as parallax, and its explanation is similar to what happens when we place our thumb between and object and our eyes. When closing either eye, we will notice the apparent movement of the object. In the case of a stellar parallax, each eye corresponds to the vertices of the base of an isosceles triangle created when joining the Earth position in two vertices A and B, separated by 6 months along the orbit of the star, with vertex C, the position of the star, considered fixed. The parallax is half the angle ACB.

images

Figure 8.2     Calculating the distance to the nearest stars.

The distance, d, from the Sun to the star with parallax, p, in radians, can be computed as follows:

d=1AUtan p={given that the angle is very small,tan p≃p}≃1pAU.

A big leap forward in the measurement of astronomical distances happened when Henrietta Swan Leavitt (1868–1921) noticed that certain stars, now known as Cepheid variables, exhibited regular changes in their luminosity from which their absolute magnitude could be calculated.

The relationship period/luminosity for the Cepheid variables has been revised several times since Henrietta Leavitt’s original measurements. Currently, the following function is usually applied:

M = -2.78 log (P) - 1.35

where M is the absolute magnitude of the star (the luminosity of the star if it was 1 pc away) and P is the average period in days.

Applying a similar criterion to the one used to determine the distance between a bulb with known brightness and the brightness experienced by an observer, we can deduct the following relationship between M, the apparent magnitude m (magnitude seen from the Earth, calculated as the average between the maximum and minimum observed magnitudes) and the distance D (in pc):

m - M = 5 log (D/10)= 5 log(D) - 5

This method can be applied to stars in our galaxy with the help of big ground telescopes or space ones. We can also use it to measure the distances to our nearest galaxies assuming we can find Cepheid variable stars in them.

Example: In Astroex (http://www.astroex.org), part of the European Southern Observatory (ESO), we can find a set of basic astronomy exercises using real data. In one of them we are given the results of the calculations done by Hubble in several Cepheid variables located in the M100 galaxy. According to the data, one star has a period of 53.5 days and an apparent magnitude of 24.5. Stars of such magnitude can only be observed with large telescopes.

It’s one of the most distant objects ever measured using this method.

For remote galaxies the method used is based on measuring the apparent luminosity of a type Ia supernovae and relating it to the real magnitude (that can be deduced using theory).

Another key parameter related to the previous one is the redshift experienced by light. This is an example of the Doppler effect caused by galaxy movements, mostly associated with the expansion of the universe.

Astronomer Edwin Powell Hubble (1889–1953) (the famous telescope was named after him) realized that galaxies frequently presented redshift. The amount changed depending on the galaxy. Later on, this astronomical phenomenon became known as Hubble’s law.

cz = H0D

with:

z (redshift) = (λob - λe)/λe, where λob corresponds to the observed wavelength and λe to the emitted one; z is therefore a dimensionless number.

c is the speed of light.

D is the actual distance to the galaxy (in Mpc).

H0 is the Hubble constant at the moment of observation (during our lifetime and for the duration of our civilization, this value will be practically constant).

Once z is known, we can calculate D using basic algebra:

D=czH0

Hubble’s law is considered a fundamental relation between recessional velocity and distance. However, the relation between recessional velocity and redshift depends on the cosmological model adopted, and is not established except for small redshifts.

A modified version of Hubble’s law and other properties of the universe are now included in UniverseModelData.

We can see that there are some discrepancies. This is because galaxies are also affected by the gravitational forces of their galactic neighbors, e.g., Andromeda influencing the Milky Way. In fact, Andromeda and our galaxy are approaching each other and will collide in approximately 4 billion years from now.

The possibilities of the astronomical functions in Mathematica don’t end here. Feel free to continue exploring.

8.5 Application: Binary Systems and the Search for Extrasolar Planets

Quit[]

In this section we will show how we can use Mathematica for other Astronomy related matters. We will describe in particular several features developed by the author for the study of binary systems and, by extension, the search for extrasolar planets.

Since 1995, when the first extrasolar planet was found, more than 400 have been detected. NASA’s Kepler spacecraft (http://www.nasa.gov/kepler) launched on March 6, 2009, was specifically designed to search for extrasolar planets. Probably, by the time you read these lines, their total number will have increased significantly.

To find exoplanets we mainly use two methods:

1) Radial velocity: A star with a planet will move in its own small orbit in response to the planet’s gravity leading to variations in the radial velocity of the star with respect to Earth. These movements are only detectable for giant planets next to stars.

2) Transit photometry: When a planet crosses in front of its parent star’s disk the observed brightness of the star drops by a small amount. This last method is becoming the most effective one. The principle is similar to the one used when studying binary star systems.

More than half of the stars that we see are actually systems of two or more stars. Normally it’s not possible to distinguish between different stars in multiple systems (not to confuse them with naked-eye binaries that are usually stars that seem close to each other when seen from Earth but are actually located in different systems altogether). The detection is usually done with the transit photometry method. This approach is very effective when dealing with eclipsing binaries, systems of two stars with coplanar orbits when observed from Earth.

Next we are going to describe how to build this kind of system using Mathematica and study its brightness. In the example that follows, we use as reference Eclipsing Binary Stars by Dan Bruton: http://www.physics.sfasu.edu/astro/ebstar/ebstar.html.

Figure 8.3 represents a binary system (although we discuss the case of two stars, the method is similar if we consider a star and a planet). The coordinates f the centers are (x1, y1, z1) for star 1 and (x1, y1, z1) for star 2 (assuming homogeneous spheres, the centers of mass coincide with the geometric centers). We consider {0, 0, 0} as the coordinates for the center of mass of the system. These coordinates have been chosen so that an observer, located on Earth, is on the axis OZ. That is: the one that goes from the center of mass to Earth. For the OY axis, we take the perpendicular to OZ in the star 1 plane when it forms an angle θ = 0. This happens when the star is closest to the observer.

images

Figure 8.3     A binary system.

With:

i = orbital inclination: angle, in radians, between the line of vision and the orbital plane axis. R = distance between the centers of both stars. It’s expressed in arbitrary units, usually the average distance between the stars. In the case of a constant orbit R = 1.

m1, m2 = stars’ masses (we can choose arbitrary units but it’s important that they keep the equivalence q = m2/m1).

r1, r2 = radii. For convenient purposes we choose r1 > r2. In practice, we can choose stars 1 and 2 so that the condition is always met. It’s convenient to express the radii in fractions of the major semi-axis. In the case of a circular orbit, obviously the major semi-axis is greater than or equal to the minor one and identical to the radius R.

θ = Phase, corresponding to the angle between the star and OZ, using θ = 0, the angle where the center of star 2 is the closest to the observer. In a binary system with circular orbit, θ is a variable parameter (it’s time dependent) while the previous parameters are constant (although that is not always the case, for example: There might be some intrinsic changes to the brightness or even the mass could change as in the case of a mass transfer between stars. However in the example that follows we’ll use the simplest case, that happens to be the most common one). θ is known as the azimuthal angle and follows θ=2π(t−t0)P, with P being the orbital period.

In a binary system, the coordinates {x1, y1, z1} and {x2, y2, z2} for the location of the stars are given by:

x1=−x1+(1/q),y1=−y1+(1/q),z1=−z1+(1/q);x2=x1+q,y2=y1+q,z2=z1+q;

The initial moment θ = 0 is the time when both stars are aligned and can be seen from the OZ axis (where the observer is located), that is when the interposition of the stars is at its maximum. In a real system, all the parameters would be fixed and only θ would change as a function of time.

8.6 Light Curves

8.6.1 Light Curves in Binary Systems

Quit[]

The transit photometry method discussed in the previous section is based on the study of the changes in brightness when a star or a planet interposes itself between another star and the observer. As mentioned before, since stars are not visually distinguishable, even less a star and a planet, what we will see is a variation in the apparent brightness.

The method described here, with the arrival of astronomical CCD cameras, is within the reach of Astronomy fans. There have been instances where fans have even been able to detect exoplanets using this technique. This is a great example of how science can be done without the need for major resources.

The luminosity l of a star is defined as the amount of energy emitted by the star per unit of time. The flow F is the energy emitted per unit of surface and per unit of time, that is F1=l14πR12 and F2=l24πR22 where l1, l2 = luminosities (in arbitrary units, as long as they keep the ratio l2/l1).

Then an observer will see the system with a luminosity given by:

l = K (F1 A1 + F2 A2)

where A1 and A2 represent each of the areas of the stars’ disks as seen by the observer and K that can be determined from the observer’s detector area, on the OZ axis, and the distance between the Earth and the binary system.

To find these areas we need to know the apparent distance ρ as seen by the observer (located on the OZ axis). That distance can be calculated as follows:

ρ=(x1−x2)2+(y1−y2)2

  • The next function computes Ρ(R, i, θ, q):

    rho[R _, i _, θ_, q_] = Module[{x, y, x1, x2, y1, y2}, x = R Sin[θ];
       y = R Cos[i] Cos[θ];
       {x1,y1}={−x1+(1/q),−y1+(1/q)}{x2,y2}={x1+q,y1+q};
       Sqrt[(x2 - x1) ^2 +(y2 - y1) ^2]];
    rho[2, 2 Pi 85 /360, θ, 0.1/ 0.6]
    √(0.0303845 Cos[θ]2 + 4. Sin[θ]2)

  • With the example data (R = 2, i = 2 Pi 85 /360, m1 = 0.6, m2 = 0.1) the apparent distance Ρ changes according to the graph below. Based on our problem formulation, the lowest apparent size corresponds to the maximum interposition that happens when θ = 0, Pi, 2 Pi,..., n Pi.

    Plot [rho[2, 2 Pi 85/360, θ, 0.1/0.6], {θ, 0, 4 Pi}, AxesLabel → (θ, ρ}]

    images

To calculate A1 and A2 it is necessary to take into account the current stage in the cycle.

  • The two figures below show the star with the bigger radius (r1) in red, and the one with the smaller radius (r2) in blue.

    Graphics [{Red, Disk[{0, 0}, 2], Blue, Disk[{4, 0}, 1]}]

    images

Without eclipse.- Corresponds to the period in which all the light emitted by both stars in the direction of the observer reaches him. Therefore, this is when the luminosity is at its maximum. In this stage, the following relationship holds: ρ > r1 + r2 with:

A1=π r12 and A2=π r22

Partial eclipse.- This happens when one star partially covers the other one. During this stage, r1 + r2 > ρ > r1 - r2.

  • Figure 8.4 was created using the following function:

    Graphics [{{Red, Circle[{0, 0}, 1]}, {Blue, Circle[{1.15, 0}, 0.5,
    {-2 Pi/3, 2 Pi/3}]}, {Dashed, Blue, Circle[{1.15, 0}, 0.5]},
    Line[{{0, 0}, {0.9, 0.425}, {0.9, -0.425}, {0, 0}}], Line[{{0.9,
    0.425}, {1.15, 0}, {0.9, -0.425}}]}]

    We added legends to the output of the function using Mathematica’s 2D Drawing Tools: images in Windows or images in OS X.

    images

    Figure 8.4     Partial eclipse diagram.

To calculate the amount of light that arrives from the eclipsed star, we have to subtract the eclipsed area from the total area. This means we have to compute ΔA1 and ΔA2 corresponding to the star 1 and star 2 segments respectively using the following formulas:

ΔA1=12R12(φ1−Sin[φ1]); ΔA2=12R22(φ2−Sin[φ2])

The following expressions enable us to compute ϕ1 and ϕ2:

R22=R12+ρ2−2R1 ρ Cos[φ12],R12=R22+ρ2−2R2 ρ Cos[φ22]

  • The output below represents the different stages of a partial eclipse:

    {Graphics [{Blue, Disk[{2, 0}, 1], Red, Disk[{0, 0}, 2]}],
    Graphics [{Blue, Disk[{-2, 0}, 1], Red, Disk[{0, 0}, 2]}],
    Graphics [{Red, Disk[{0, 0}, 2], Blue, Disk[{-2, 0}, 1]}],
    Graphics [{Red, Disk[{0, 0}, 2], Blue, Disk[{2.5, 0}, 1]}]}

    images

The two figures on the left correspond to the scenario where star 1 (red) is the closest one, that is z1 > z2. In this case the light that reaches the observer will be:

A1=πr12 and A2=πr22−ΔA1−ΔA2A1=πr12 and A2=πr22−ΔA1−ΔA2

The rest of the figures corresponds to the case where star 1 (red) is the most distant one, that is z1 < z2 then:

A1=πr12−ΔA1−ΔA2 and A2=πr22A1=πr12−ΔA1−ΔA2 and A2=πr22

Total or annular.- Corresponds to the period when a star is shaded by the other one completely. In this stage: ρ ≤ r1 - r2 with: A1=πr12 and A2 = 0 for z1>z2 and A1=πr12−πr22 and A2=πr22

All these conditions are summarized in the table below:

images

Now we are going to include all the criteria described earlier in the function area[R, i, m2/m1, r1, r2, l1, l2, θ], where: m2 > m1 and r1 > r2. The function returns the area that the observer will measure.

Without eclipse. When ρ≥ r1 + r2:

   area[R_, i_, q_, r1_, r2_, l1_, l2_, θ_] ≔
    Module[{A1, A2}, {A1, A2}={Pi r1^2, Pi r2^2};
     l1 A1/r1^2 + l2 A2/r2^2]/; rho[R, i, θ, q] ≥ r1 +r2

Partial eclipe. When r1 + r2 > ρ > r1 - r2:

  • We first calculate ΔA1 and ΔA2

    ΔA1[R_, i_, q_, r1_, r2_, θ_] ≔ Module[{ϕ1, ρ, phi, phi1},
       phil1=Last[Solve[r22 == r12+ρ2−2 r1 ρ Cos[phi2],phi]/.
            ρ -> rho[R, i, 0, q]][[1, 2]];
       12r12(φ1−Sin[φ1]) /. φ1→phil];
     
    ΔA2[R_, i_, q_, r1_, r2_, 0_] ≔ Module [{φ2, ρ, phi, phi1},
       phi1=Last[Solve[r12 == r22+ρ2−2 r2 ρ Cos[phi2],phi]/.
           ρ -> rho[R, i, 0, q]][[1, 2]];
       12r22(φ2−Sin[φ2])/. φ2→phi1];

  • In this type of eclipse, for z1 > z2: A1=πr12 and A2=πr22−ΔA1−ΔA2 while for z1<z2:A1=πr12−ΔA1−ΔA2 and A2=πr22.


    area[R_, i_, q_, r1_, r2_, l1_, l2_, 0] :=
      Module[{z1,z2,a1,a2,A1,A2},z1=−R Cos[θ]Sin[i]1+1q;z2=R Cos[θ] Sin[i]1+q;
      a1 = ΔA1[R, i, q, r1, r2, 0];
      a2 = ΔA2[R, i, q, r1, r2, 0];
       {A1, A2}=
        If[z1 > z2, {Pi r1^2, Pi r2^2 - a1 -a2}, {Pi r1^2 - a1 -a2, Pi r2^2}];
       l1 A1/r1^2 + 12 A2/r2^2] /; r1 + r2 > rho[R, i, 0, q] > r1 - r2

Total and / or annular eclipse.

  • When ρ ≤ r1 - r2 one of the stars will be directly in front of the other one, interposing it. Then a total eclipse will happen if z1 > z2 with A1=πr12 and A2 = 0, or an annular one will occur if z1 < z2 with A1' = A1 - A2 and A2' = A2.

    area[R_, i_, q_, r1_, r2_, l1_, l2_, θ_] :=
     Module[{z1,z2,A1,A2},z1=−R Cos[θ]Sin[i]1+1q;z2=R Cos[θ] Sin[i]1+q;
      {A1, A2} = If[z1 > z2, {Pi r1^2, θ, 0}, {Pi r1^2 - Pi r2^2, Pi r2^2}];
      l1 A1/r1^2 + l2 A2/r2^2] /; rho[R, i, θ, q] ≤ r1 - r2

Let’s simulate the ideal case of an eclipsing binary. We can distinguish the repeating phases of an eclipse: one corresponds to the period without eclipse (t1 to t2), in t2 the eclipse starts and reaches its plenitude between t3 and t4, from t4 to t5 we are in the partial eclipse phase, from t6 to t9 is the transit phase, and when t9 = t1 the cycle starts again.

  • We use Plot[area[2, 90 Degree, 0.1/0.6, 0.6, 0.3, 0.6, 0.2, 2 Pi t { t, 0, 12}, AxesLabel -> {"t", "I"}, AxesOrigin -> {0, 0}] and include the legends using the Drawing Tools option in the Graphics menu: Graphics ►Drawing Tools.

images

Figure 8.5     Light curve.

Example: With R = 2, i = 85º, m1 = 0.6, m2 = 0.1, R1 = 0.4, R2 = 0.3, l1 = 0.6, l2 = 0.2, we’d like to show I/Imax.

  • We first calculate Imax.

    max = NMaximize[area[2, 2 Pi 85 / 360, 0.1/ 0.6, 0.4, 0.3, 0.6, 0.2, θ],
    {θ, 0, 2Pi}][[1]] // Quiet
    2.51327

  • Then we display I/Imax.

    Plot [area [2, 85 Degree, 0.1/0.6, 0.4, 0.3, 0.6, 0.2, 2 Pi θ] /max,
    {θ, 0, 1.5}, AxesLabel →{"Phase", "I/Imax"},
    AxesOrigin→{0, 0}] // Quiet

    images

8.6.2 Application: Using Light Curves to Determine the Period-Luminosity Relation

Besides eclipsing variables or planets interposition, periodic changes in luminosity can happen for other reasons. Probably the best known example is the change in luminosity in variable stars, such as Cepheids. The usual method for building an experimental light curve consists of measuring the luminosity at different times t and deduce the period from them.

  • The following expression simulates magnitude as a function of time (in practice, since we often don’t know this function we try to deduce it from the experimental data). The chosen function is very simple and we only use it for didactic purposes).

    mag[t _] = 5 + 0.14 Cos[0.3 t + 0.1];

  • Let’s assumed that we have already taken measurements at different times ti.

    data = Table[{t, Round[mag[t], 0.01]}, {t, 0, 200, 5}];

  • We represent graphically the previous data. If we only had the previous data points it would be very difficult to figure out the underlying pattern (remove Joined → True), joining the dots makes it easier to observe the cyclical nature of the data. It’s normal for many intermediate points to be missing making the analysis even more challenging.

    ListPlot[data, Joined → True]

    images

The analysis method frequently used is to represent the previous phenomenon in a phase diagram.

To calculate the phase we need to apply this equation:

ϕ=Decimal part of[t−tθP]

where:

ϕ represents the phase in cycles.

t0 (normally called epoch) represents the origin of the time variable expressed as a Julian Date (DJ).

t is the moment when the measurement is taken (normally we will measure the apparent magnitude, m).

P is the period, the time it takes for the cycle to repeat itself.

An equivalent expression to the previous one and the one that we will be using is:

φ=Mod [t−tθP]

  • With the function below Mathematica returns the phase given the initial data points and the period.

    phi[list _, p_] ≔ Module[{data}, data = list;
      Transpose [{Mod[data[[All, 1]] /p, 1], data[[All, 2]]}]];

In this example we are assuming that the period is known, but in practice that’s what we are trying to find out. We can check that by using an approximate value of the period the representation is more diffuse.

  • With real observations, the period has to be determined and that is not always easy. The next function will help us to compute it. We keep in mind that for period-magnitude curves, the order of representation in the axis OY is inverted. That is, the magnitudes are sorted in descending order not in the usual ascending one: -3 >-2 > 0 > 1 > 2.

    curveL[data_, per_, kmax_] ≔
     Module[{mag, a, b, k, i}, mag = Round[Transpose[data][[2]], 0.1];
      Manipulate[ListPlot[Join[phi[data, per + k]/. {a_, b_} -> {-(1 - a), -b},
         phi[data, per + k] /. {a_, b_} -> {-(1- a), -b}], AxesLabel ->
         {"Period", "Magnitude"}, AxesOrigin -> {-1, -Max[mag]- 0.1},
        Ticks -> {Automatic, Table[{i, -i}, {i, -Max[mag]- 0.1,
            -Min[mag] + 0.1, (Max[mag] - Min[mag]) / 10}]}],
       {{k, 0, "Deviation from the period (in days)"}, 0, kmax}]]

  • Let’s assume a complete cycle (2Π). We use the previous function to simulate the experimental measures and include a random component. We make the assumption that the period is 20 days (the duration of a complete orbit).

    simulation = Table[
        {θ, area[2, 2 Pi 85 /360, 0.1 /0.6, 0.4, 0.3, 0.6, 0.2, 2 Pi θ/20] +
          RandomReal[NormalDistribution[0, 0.005]]},
        {θ, 0, 60, 0.05}]; // Quiet

  • We represent the data in a phase diagram:

    simulation1 = Sort[phi[simulation, 20]];
    ListPlot [simulation1, Joined → True]

    images

The previous luminosity data have to be multiplied by a parameter k to convert them into a magnitude measure. Remember that higher luminosity numerically implies lower magnitude. For that reason, to represent the data in a phase/luminosity diagram we assume k = -1.

  • If you execute the following function and use the slider to choose a 5-day period, you’ll see that the curve is very similar to the previous figure.


    curveL [Map[Times [{1, -1}, #] &, simulation], 15, 20]

    images

8.6.3 Application: Finding the Light Curve of a Known Eclipsing Binary

Let’s apply the newly created function to the luminosity data from NSV 03199 obtained by Garcia-Melendo, E., Henden, A., 1998, IBVS, No. 4546. The data are in the file: NSV03199.xls in our data directory. It has been downloaded from: http://www.astrogea.org/VARIABLE/variables.htm.

  • Let’s import the luminosity measures:

    vsdata =
    Drop [Import [FileNameJoin [{NotebookDirectory[], "Data", "NSV03199.xls"}],
    {"XLS", "Data", 1}], 1];

  • The data have Julian dates so we use as reference JHD=2450510.46542.

    vsdata1 = vsdata /. {a_, b_} -> {a - 2450510.46542, b};

  • We cannot really see the cycle pattern:

    ListPlot [vsdata1]

    images

However, when using a period-magnitude curve, the periodic behavior corresponding to an eclipsing binary can be seen clearly.

  • We find that the optimum period is 1.046400 days:

    curveL[vsdata1, 1.046400, 0.5]

    images

All the previous functions along with other examples in the book have been put together in a package that can be downloaded from diarium.usal.es/guillermo.

8.7. Additional Resources

Astronomy demonstrations:

http://demonstrations.wolfram.com/topic.html?topic=Astronomy

Eric Weisstein’s World of Astronomy: http://scienceworld.wolfram.com/astronomy

To follow objects with variable luminosity (variable stars, supernovae and many more), the best place is the website of the American Association of Variable Star Observers (AAVSO): http://www.aavso.org

Author’s website astronomy section (in Spanish):

http://diarium.usal.es/guillermo/astronomia