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

# In[2]:


get_ipython().run_line_magic('matplotlib', 'inline')


# PyGSM quickstart
# ------------------
# 
# Firstly, we import the `pygsm` module:

# In[5]:


from pygsm import GlobalSkyModel


# We can very quickly get started with generating a sky map and plotting it:

# In[6]:


gsm = GlobalSkyModel()
gsm.generate(500)
gsm.view(logged=True)


# We might then want to write this to a FITS file (in healpix format):

# In[7]:


gsm.write_fits("gsm_500mhz.fits")


# There are also a few options we can play around with. For example:
# 
# * let's change to use a different interpolation method (pchip) 
# * instead of using the Haslam 408 MHz map for structure, let's use WMAP 23 GHz
# * instead of using MHz to define the units, let's use GHz

# In[11]:


gsm = GlobalSkyModel(freq_unit='GHz', interpolation='pchip', basemap='wmap')
gsm.generate(30) # Generate at 30 GHz
gsm.view(logged=True)


# You can set things on the fly, too, once you've got a `gsm` object:

# In[13]:


gsm.set_basemap('5deg')
gsm.view(logged=True)


# Just in case you're interested in spectral cubes, you can even pass a frequency range:

# In[19]:


gsm = GlobalSkyModel()
freqs = np.linspace(10, 1000, 100)
map_cube = gsm.generate(freqs)

print map_cube.shape
plt.loglog(map_cube[:,0])
plt.loglog(map_cube[:,12345])
plt.xlabel("Frequency [MHz]")
plt.ylabel("Temperature [K]")


# Be careful with this one though, as it can eat up a lot of memory!

# ## Generate observed sky for a given lat, long
# 
# A common task is to find out what the sky looks like at a given lat, long and date. `PyGSM` gives a quick method to do this:

# In[49]:


from pygsm import GSMObserver
from datetime import datetime

# Setup observatory location - in this case, Parkes Australia
(latitude, longitude, elevation) = ('-32.998370', '148.263659', 100)
ov = GSMObserver()
ov.lon = longitude
ov.lat = latitude
ov.elev = elevation
ov.date = datetime(2000, 1, 1, 23, 0)


# Now generate a sky model and view an all-sky orthographic plot

# In[53]:


ov.generate(1400)
d = ov.view(logged=True)


# And we can view galactic Mollweide projection too, with below-horizon data blanked out:

# In[55]:


d = ov.view_observed_gsm(logged=True)