notebook.community

Edit and run

prolegomena





In [1]:



    


# --+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
# I need the time axis labels in italian,
# it's usually better to do this at first
from locale import setlocale, LC_TIME ; setlocale(LC_TIME, 'italian')

# astro's computations will be done using ephem 
import ephem as x

# not every important city is includes in ephem city database
from ephem.cities import _city_data, city
_city_data.update({
    'Monza':('45.577721','9.300896',160.0),
    'Cabiate':('45.679350','9.165230',255.0),
    'Meda':('45.661400','9.155734',223.0),
     })

# we need to fiddle a bit with the labeling of the time axis
from matplotlib.ticker import MultipleLocator
from matplotlib.dates import MonthLocator, DateFormatter

# the preferred format for inline plots, alternatives are 'png' and 'svg'
%config InlineBackend.figure_format = 'png'

# the time unit is a day, here we define day, hour, etc
d, h, m, s = 1.0, 1.0/24, 1.0/24/60, 1.0/24/60/60

Moonlight

Initialization

We need an observer, or better two observers, the second one will be used behind the scenes to do some computation.





In [2]:



    


my_city = 'Cabiate'
me  = city(my_city)
me2 = city(my_city)

Finally, we need a moon.





In [3]:



    


moon = x.Moon()

"Vertical" velocity of the moon

For a given date, i.e., moment in time,

  1. i compute the position of the moon in the sky of the observer (me2) a little before and a little later,
  2. i extract the height of the moon (moon.alt) in the sky of the observer in these two instants and
  3. i approximate the vertical velocity with the usual finite differences formula.




In [4]:



    


def velocity(date):
    me2.date=date-0.1*s
    moon.compute(me2) ; alt0 = moon.alt
    me2.date=date+0.1*s
    moon.compute(me2) ; alt1 = moon.alt
    return alt1/0.2/s - alt0/0.2/s

Day-by-day maximum heigth of the moon

Given that the maximum heigth is reached when the moon is close to transit, after a bit of initialization,

  1. the position of the moon is computed at the beginning of the day
  2. the transit time is extracted
  3. test if the moon transits during the day
    1. bracketing the transit time, use the default solver to find the time for which the vertical velocity is zero,
    2. save in a container the couple of values time, altitude of the moon (nb, the plotting routines expect that time is in a particular format, hence the .datetime() call)
  4. the date of the observer is incremented by one day (me.date = me.date + d),
  5. if the new date is in 2014 break otherwise repeat the cycle.

The container (a Python list) is converted to an array, so that we can index it with a [i,j] notation.





In [5]:



    


me.date = "2013/01/01 00:00:01"
next_year = x.date("2014/01/01 00:00:01")
container = []

while 1:
    moon.compute(me)
    mtt = moon.transit_time
    if mtt:
        time = x.Date(x.newton(velocity,mtt-0.2*h,mtt+0.2*h))
        container.append( (time.datetime(), 180*moon.alt/x.pi))
    me.date = me.date + d
    if me.date>next_year: break

container=array(container)

The phases of the moon

This is ad-hoc for 2013, because i know that the first phase is a last quarter, anyway... the dates of new, 1st q., full and last q. moons are stored in a list, together with an integer 0,..,3 denoting the phase. At the end, the list of moons is converted to an array.





In [6]:



    


date = x.Date('2013/01/01') ; moons = []

while 1:
    
    date = x.next_last_quarter_moon(date)
    if date>next_year: break
    moons.append((date.datetime(),3))
    
    date = x.next_new_moon(date)
    if date>next_year: break
    moons.append((date.datetime(),0))
    
    date = x.next_first_quarter_moon(date)
    if date>next_year: break
    moons.append((date.datetime(),1))
    
    date = x.next_full_moon(date)
    if date>next_year: break
    moons.append((date.datetime(),2))

moons = array(moons)

We want to superimpose the envelope curve with different symbols for each phase of the moon. For now





In [7]:



    


i = 0 ; delendo = []

for m in moons:
    tm, phase = m
    
    while container[i,0]<tm:
        i = i+1

    t0, alt0 = container[i-1,:]; t0 = x.Date(t0)
    t1, alt1 = container[ i ,:]; t1 = x.Date(t1)
    altm = alt0 + (x.Date(tm)-t0)*(alt1-alt0)/(t1-t0)
    delendo.append((tm,altm,phase))
    
moons = array(delendo) ; del delendo

When applying a test on some elements of an array, we obtain an array of boolean values (i.e., true vs false) of the same length of the tested sequence.





In [8]:



    


new  = moons[:,2]    == 0
full = moons[:,2]    == 2
qrtr = moons[:,2] %2 == 1

The important fact is that we can use these boolean arrays to index an array... e.g., we print the new array and then we use it to print the first 3 new moons.





In [9]:



    


print new
print moons[new][:3]



    


















    






[False  True False False False  True False False False  True False False
 False  True False False False  True False False False  True False False
 False  True False False False  True False False False  True False False
 False  True False False False  True False False False  True False False
 False]
[[datetime.datetime(2013, 1, 11, 19, 43, 37, 595710) 26.133027939214347 0]
 [datetime.datetime(2013, 2, 10, 7, 20, 6, 621791) 34.04232785394464 0]
 [datetime.datetime(2013, 3, 11, 19, 51, 0, 142128) 44.13846824549499 0]]

Plotting

Prepare the axis on which to plot

  • We want to plot 13 cycles of an almost periodic function, it's better to have an elongated x-axis.
  • The margins around the graph are by default a fixed ratio of the figure size and the horizontal margins are hence too large, fix this.
  • The range of ordinates is ~ 20 to 65 degrees so limits on the y axis would be too tight, so we specify a range 0-90 degrees and, while at it, we give a precise range to the x axis.
  • By default, we have a label every two months, we want a label for every month and the year in every label.
  • By defalt, no minor ticks... we want a minor tick for every day.

Show our current results.





In [10]:



    


# size of the figure in inches
figsize(10,2.5)

# in this figure there is one x-y graph
subplot(111)

# save the current figure and the current x-y graph
f = gcf() ; ax = gca() ; close()

# the margins of the x-y are fractions of the figure,
# having increased the figsize we narrow the margins
f.subplots_adjust(left=0.06, right=0.96, bottom=0.08, top=0.92)

# adjust the limits of the axes
ax.axis(ymin=0, ymax=90,
        xmin=x.Date("2013-01-01").datetime(),
        xmax=x.Date("2014/01/01").datetime())

ax.xaxis.set_minor_locator(MultipleLocator(1))
ax.xaxis.set_major_locator(MonthLocator(bymonth=range(13)))
ax.xaxis.set_major_formatter(DateFormatter("%b %Y"))
ax.set_axisbelow(True)
ax.grid(c='#B0B0B0', linestyle='-', linewidth=0.1)

show_all = True

if show_all: f.show()



    


















    






/usr/lib/pymodules/python2.7/matplotlib/figure.py:371: UserWarning: matplotlib is currently using a non-GUI backend, so cannot show the figure
  "matplotlib is currently using a non-GUI backend, "

Enhancing the axis

  • The minor ticks on x are too large.
  • The x labels almost run one into the other, choose a smaller font and, while we are at it, change also the size of the y labels.
  • I prefer y labels in vertical.




In [12]:



    


ax.tick_params(axis='x', which='minor', length=4, width=0.2, color='#909090')
[tik.set_size('xx-small') for tik in ax.get_xticklabels()]
[tik.set_size('xx-small') for tik in ax.get_yticklabels()]
[tik.set_rotation(90.0)   for tik in ax.get_yticklabels()]

if show_all: f.show()

Titles for the x-y graph and the individual axes

I feel no need for a specific label for the whole x-axis, month names are good enough to infere the meaning of the abscissae.





In [13]:



    


# the title and a label for the ordinates
ax.set_title(
    'La Luna vista da %s (lat = %s, lon = %s, altezza = %5.1f mslm).'%
    (me.name,x.degrees(me.lat),x.degrees(me.lon), me.elevation),
    size='small')
ax.set_ylabel('Inviluppo dell\'altezza lunare, in gradi', size='x-small')

if show_all: f.show()

Plot the envelope of the moon altitude over the horizon

We have a function that is defined only in a discrete set of points and we plot it as a continuous function with a very light color, then we superimpose small black dots in the position where we have found a local maximum.





In [14]:



    


# plot the envelope of moon altitudes two times, first with a continuous line
# & later with the smallest of the available dot typess
ax.plot(container[:,0],container[:,1],'-',color='lightgray', linewidth=0.3)
ax.plot(container[:,0],container[:,1],',',color="black")

if show_all: f.show()

Plot the New Moons, etc

First i plot differently colored circles, with a slight transparency so that the underlying envelope curve is still visible and later i superimpose in the center of each circle a small black dot to help the eye to exactly position each moon.





In [15]:



    


# plot the postion of a) new moons, b) full moons and c) both first quarter and last quarter
# of the moon

# first time, with circles of different colors
ax.plot(moons[new, 0],moons[new, 1],'o',markersize=5, color='#202090',   alpha=0.75)
ax.plot(moons[full,0],moons[full,1],'o',markersize=5, color="#ffff00",   alpha=0.75)
ax.plot(moons[qrtr,0],moons[qrtr,1],'o',markersize=5, color="lightgray", alpha=0.75)

# second time, with small black dots
ax.plot(moons[:,0],moons[:,1],'.k',markersize=2)

if show_all: f.show()

What's up with these random circles?

I know that the small circles represent the approximate maximum altitude of the moon when the moon enters a new phase and i think this is almost evident, but an annotation is easy to place.





In [16]:



    


xpos = x.Date('2013/07/1').datetime()

ax.annotate('''I cerchi sul grafico indicano la posizione delle lune piene (cerchi gialli),
delle lune nuove (cerchi blu) e dei quarti di luna (cerchi grigi).''', (xpos,12),
    size='x-small', ha='center', va='center',
    bbox=dict(boxstyle="round,pad=0.5",fc='white',ec='lightgray'))

if show_all: f.show()

A gratuitous comment...

It's in Italian, means that these two full moons are very high in the sky and, consequently, we will have a sort of "white nights" in the 3 days atound the full moon date.





In [17]:



    


ax.annotate('',(moons[-2,0],moons[-2,1]), xytext=(10./12.,0.85), textcoords='axes fraction', ha='right',
    arrowprops=dict(width=0.1,headwidth=4.0,shrink=0.15,color='gray'))
ax.annotate(u'''Le due ultime lune piene di autunno saliranno molto alte nel
cielo e le tre notti della luna piena saranno molto luminose.''',
    (moons[-6,0],moons[-6,1]), xytext=(10./12.,0.85), textcoords='axes fraction', ha='right',
    name='sans', size='xx-small', color='#404090',
    arrowprops=dict(width=0.1,headwidth=4.0,shrink=0.25,color='gray',),
    bbox={"boxstyle":"round,pad=0.5","fc":"#feffff","alpha":0.9,'ec':'#b0b0ff'})

f.show()

This fine plot deserves to be saved

And it's just a single line of code!





In [18]:



    


f.savefig('Moon_at_%s.pdf'%(me.name,))



    











In [18]:



    






    











In [18]:



    






    











In [ ]:

Content source: boffi/boffi.github.io

Similar notebooks:

notebook.community | gallery | about