#!/usr/bin/env python
# coding: utf-8

# # Units and Quantities
# 
# ## Objectives
# 
# - Use units
# - Create functions that accept quantities as arguments
# - Create new units

# ## Basics
# 
# How do we define a Quantity and which parts does it have?

# In[ ]:


import astropy.units as u


# In[ ]:


# Define a quantity length
length = 26.2 * u.meter
# print it
print(length) # length is a quantity


# In[ ]:


# Type of quantity
type(length)


# In[ ]:


# Type of unit
type(u.meter)


# In[ ]:


# Quantity
length


# In[ ]:


# value
length.value


# In[ ]:


# unit
length.unit


# In[ ]:


# information
length.info


# Quantities can be converted to other units systems or factors by using `to()`

# In[ ]:


# Convert it to: km, lyr
print(length.to(u.km))
print(length.to(u.lightyear))


# We can do arithmetic operations when the quantities have the compatible units:

# In[ ]:


# arithmetic with distances
distance_start = 10 * u.mm
distance_end = 23 * u.km
length = distance_end - distance_start
print(length)


# Quantities can also be combined, for example to measure speed

# In[ ]:


# calculate a speed
time = 15 * u.minute
speed = length / time
print(speed)


# In[ ]:


# decompose it
print(speed.decompose())
print(speed.si)


# 
# <section class="challenge panel panel-success">
# <div class="panel-heading">
# <h2><span class="fa fa-pencil"></span> Unit conversions</h2>
# </div>
# 
# 
# <div class="panel-body">
# 
# <ol>
# <li>Convert the speed in imperial units (miles/hour) using:
#     <code>from astropy.units import imperial</code></li>
# <li>Calculate whether a pint is more than half litre.
#     <em>You can compare quantities as comparing variables.</em>
#     Something strange? Check what deffinition of <a href="https://en.wikipedia.org/wiki/Pint">pint</a> astropy is using.</li>
# <li>Does units work with areas? calculate the area of a rectangle of 3 km of side and 5 meter of width. Show them in $m^2$ and convert them to yards$^2$</li>
# </ol>
# 
# </div>
# 
# </section>
# 

# In[ ]:


#1
from astropy.units import imperial
print(speed.to(imperial.mile/u.hour))

#2
print(imperial.pint > 0.5 * u.l)
print(imperial.pint.to(u.ml))
 # A liquid pint in US is 473 ml; in UK is 568 ml

#3
rectangle_area = 3 * u.km * 5 * u.m
print(rectangle_area)
print(rectangle_area.decompose())
print(rectangle_area.to(imperial.yard ** 2))


# ## Composed units
# 
# Many units are compositions of others, for example, one could create new combinationes for ease of use:

# In[ ]:


# create a composite unit
cms = u.cm / u.s
speed.to(cms)


# In[ ]:


# and in the imperial system
mph = imperial.mile / u.hour
speed.to(mph)


# and others are already a composition:

# In[ ]:


# what can be converted from s-1?
(u.s ** -1).compose()


# In[ ]:


# or Jules?
(u.joule).compose()


# In[ ]:


# Unity of R
(13.605692 * u.eV).to(u.Ry)


# Sometime we get *no units* quantitites

# In[ ]:


# no units
nounits = 20. * u.cm / (1. * u.m)
nounits


# What happen if we add a number to this?

# In[ ]:


# arithmetic with no units
nounits + 3


# In[ ]:


# final value of a no unit quantity
nounits.decompose() # It's a unitless quantity


# ## Equivalencies
# 
# Some conversions are not done by a conversion factor as between miles and kilometers, for example converting between wavelength and frequency.

# In[ ]:


# converting spectral quantities
(656.281 * u.nm).to(u.Hz) # Fails because they are not compatible


# In[ ]:


# but doing it right
(656.281 * u.nm).to(u.Hz, equivalencies=u.spectral())


# Other built-in equivalencies are: 
#  - `parallax()`
#  - Doppler (`dopplr_radio`, `doppler_optical`, `doppler_relativistic`)
#  - spectral flux density
#  - brigthness temperature
#  - temperature energy
#  - and you can [build your own](http://astropy.readthedocs.org/en/stable/units/equivalencies.html#writing-new-equivalencies)

# In[ ]:


# finding the equivalencies
u.Hz.find_equivalent_units()


# In[ ]:


# but also using other systems
u.Hz.find_equivalent_units(equivalencies=u.spectral())


# ## Printing the quantities

# In[ ]:


# Printing values with different formats
print("{0.value:0.03f} {0.unit:FITS}".format(speed))
print("{0.value:0.03f} {0.unit:latex_inline}".format(speed))


# ## Arrays
# 
# Quantities can also be applied to arrays

# In[ ]:


# different ways of defining a quantity for a single value
length = 44 * u.m
time = u.Quantity(23, u.s)
speed = length / time
speed


# In[ ]:


# now with lists
length_list = [1, 2, 3] * u.m

# and arrays
import numpy as np
time_array = np.array([1, 2, 3]) * u.s

# and its arithmetics
length_list / time_array


# In[ ]:


# angles are smart!
angle = u.Quantity(np.arange(180), u.deg)
print(angle[[0, -1]])
print(np.sin(angle[[0, -1]]))


# ## Plotting quantities
# 
# To work nicely with matplotlib we need to do as follows:

# In[ ]:


# allowing for plotting
from astropy.visualization import quantity_support
quantity_support()

# loading matplotlib
get_ipython().run_line_magic('matplotlib', 'inline')
from matplotlib import pyplot as plt


# In[ ]:


# Ploting the previous array
plt.plot(angle, np.sin(angle))


# ## Creating functions with quantities as units
# 
# We want to have functions that contain the information of the untis, and with them we can be sure that we will be always have the *right* result.

# In[ ]:


# Create a function for the Kinetic energy
@u.quantity_input(mass=u.kg, speed=u.m/u.s)
def kinetic(mass, speed):
    return (mass * speed ** 2 / 2.).to(u.joule)


# In[ ]:


# run without units
kinetic(5, 10) # Fails! it doesn't specify the units.


# In[ ]:


# Run with units
kinetic(5 * u.kg, 100 * cms)


# 
# <section class="challenge panel panel-success">
# <div class="panel-heading">
# <h2><span class="fa fa-pencil"></span> Using `quantity_input`</h2>
# </div>
# 
# 
# <div class="panel-body">
# 
# <ol>
# <li>Create a function that calculates potential energy where $g$ defaults to Earth value, but could be used for different planets. Test it for any of the $g$ values for any other <a href="http://www.physicsclassroom.com/class/circles/Lesson-3/The-Value-of-g&quot;">planets</a>.</li>
# </ol>
# 
# </div>
# 
# </section>
# 

# In[ ]:


#4
@u.quantity_input(mass=u.kg, height=u.m, g=u.m / u.s ** 2)
def potential(mass, height, g=9.8 * u.m / u.s **2):
    return (0.5 * mass * g * height).to(u.joule)


# In[ ]:


# run it for some values
potential(5 * u.kg, 30 *u.cm)


# In[ ]:


# on Mars:
potential(5 * u.kg, 1 * u.m, g=3.75 * u.m/u.s**2)


# ## Create your own units
# 
# Some times we want to create our own units:

# In[ ]:


# Create units for a laugh scale
titter = u.def_unit('titter')
chuckle = u.def_unit('chuckle', 5 * titter)
laugh = u.def_unit('laugh', 4 * chuckle)
guffaw = u.def_unit('guffaw', 3 * laugh)
rofl = u.def_unit('rofl', 4 * guffaw)
death_by_laughing = u.def_unit('death_by_laughing', 10 * rofl)
print((1 * rofl).to(titter))


# 
# <section class="challenge panel panel-success">
# <div class="panel-heading">
# <h2><span class="fa fa-pencil"></span> Area with units</h2>
# </div>
# 
# 
# <div class="panel-body">
# 
# <ol>
# <li>Convert the area calculated before <code>rectangle_area</code> in <a href="https://en.wikipedia.org/wiki/Hectare">hectares</a> (1 hectare = 100 ares; 1 are = 100 $m^2$).</li>
# </ol>
# 
# </div>
# 
# </section>
# 

# In[ ]:


#5
ares = u.def_unit('ares', (10 * u.m)**2)
hectar = u.def_unit('hectares', 100 * ares)
print(rectangle_area.to(hectar))