<jupyter_start><jupyter_text>![astropy_banner.jpg](media/astropy_banner.jpg)  Sharing code is healthy for the community and the science it produces- a community-developed core library for professional astronomical research- combines many functionalities from a variety of astronomy packages and languages- Goals: Usability, Interoperability and Collaboration between packages Affiliated Packages- astronomy related python packages- adhere to interface standards of astropy but are not (yet) part of the astropy core package- currently 40 affiliated packages covering many areas of astronomical research Since astropy Version 3.0.0 the package only supports Python 3! Is now included by default in Anaconda  Queries, Coordinate Systems, Time and Units Observation planning with astropy  Observability of Fomalhaut from the VLT: initial setup<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import astropy_mpl_style
plt.style.use(astropy_mpl_style)
import astropy.units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation, AltAz

# Lets observe the star Fomalhaut with the ESO VLT - 8m Telescope in Chile
# Load the position of Fomalhaut from the Simbad database

fomalhaut = SkyCoord.from_name('Fomalhaut')

print(fomalhaut)

# Load the position of the Observatory. Physical units should be assigned via the
# units function
paranal = EarthLocation(lat=-24.62*u.deg, lon=-70.40*u.deg, height=2635*u.m) 

print(paranal)  # The coordinates are stored as geocentric (position relative to
# earth centre-of-mass) as a default<jupyter_output><empty_output><jupyter_text>Calculate Fomalhaut and solar altitude/azimuth<jupyter_code># We want to observe Fomalhaut next week. We will determine the position
# in the sky as seen from Paranal in a 24 hour window centred on local midnight Oct 15.
# (Oct 15 *starts* at 00:00:00)

midnight = Time('2019-10-15 00:00:00')
delta_midnight = np.linspace(-12, 12, 1000)*u.hour
times_Oct14_to_15 = midnight + delta_midnight
frame_Oct14_to_15 = AltAz(obstime=times_Oct14_to_15, location=paranal)

# Now we transform the Fomalhaut object to the Altitute/Azimuth coordinate system

fomalhaut_altazs_Oct14_to_15 = fomalhaut.transform_to(frame_Oct14_to_15)

#print(fomalhaut_altazs_Oct14_to_15)


# We also check the position of the sun in the sky over the same time range
from astropy.coordinates import get_sun

sunaltazs_Oct14_to_15 = get_sun(times_Oct14_to_15).transform_to(frame_Oct14_to_15)<jupyter_output><empty_output><jupyter_text>Determining the night-time observability of Fomalhaut<jupyter_code># Plot the sun altitude
plt.plot(delta_midnight, sunaltazs_Oct14_to_15.alt, color='r', label='Sun')
# Plot Fomalhaut's alt/az - use a colour map to represent azimuth
plt.scatter(delta_midnight, fomalhaut_altazs_Oct14_to_15.alt, 
            c=fomalhaut_altazs_Oct14_to_15.az, label='Fomalhaut', lw=0, s=8,
            cmap='viridis')
# Now plot the range when the sun is below the horizon, and at least 18 degrees below 
# the horizon - this shows the range of twilight (-0 to -18 deg) and night (< -18 deg)
plt.fill_between(delta_midnight.to('hr').value, 0, 90,
                 sunaltazs_Oct14_to_15.alt < -0*u.deg, color='0.7', zorder=0)
plt.fill_between(delta_midnight.to('hr').value, 0, 90,
                 sunaltazs_Oct14_to_15.alt < -18*u.deg, color='0.4', zorder=0)
plt.colorbar().set_label('Azimuth [deg]')
plt.legend(loc='upper left')
plt.xlim(-12, 12)
plt.xticks(np.arange(13)*2 -12)
plt.ylim(0, 90)
plt.xlabel('Hours from UT Midnight')
plt.ylabel('Altitude [deg]')
plt.show()<jupyter_output><empty_output><jupyter_text>Fitting models to data Health warning: This is a very basic intro - to do this properly and not make catastrophic errors in your interpretation and/or analysis, you need to do a statistics course! (e.g. for much more detail, see Statistical Methods course in Block 3)  Simple 1d model<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting

# Generate data with some random noise
x = np.linspace(-5., 5., 200)
y = 3 * np.exp(-0.5 * (x - 1.3)**2 / 0.8**2) # Gaussian function
y += np.random.normal(0., 0.2, x.shape) # add Normally distributed scatter

plt.plot(x, y, 'ko')
plt.show()
# Fit the data using a Gaussian
g_init = models.Gaussian1D(amplitude=1., mean=0, stddev=1.)
fit_g = fitting.LevMarLSQFitter() # Initialises fitting method as Levenberg-Marquardt
# least squares fitting


g = fit_g(g_init, x, y)  # Fits the data, note that fit_g is effectively an alias
# for the L-M fitting function, so must use the same arguments in addition to any
# pre-specified in the initial assignment of the function to fit_g above
# (none are included in this case)

plt.plot(x, y, 'ko')
plt.plot(x, g(x), "r-")
plt.show()<jupyter_output><empty_output><jupyter_text>2d model<jupyter_code>import warnings
import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting

# Generate fake data
np.random.seed(0)
y, x = np.mgrid[:128, :128]
z = 2. * x ** 2 - 0.5 * x ** 2 + 1.5 * x * y - 1. # Polynomial function
z += np.random.normal(0., 0.1, z.shape) * 50000.
# Fit the data using astropy.modeling
p_init = models.Polynomial2D(degree=2)
fit_p = fitting.LevMarLSQFitter()

with warnings.catch_warnings():
    # Ignore model linearity warning from the fitter
    warnings.simplefilter('ignore')
    p = fit_p(p_init, x, y, z)
    

# Plot the data with the best-fit model
plt.figure(figsize=(8, 2.5))
plt.subplot(1, 3, 1)
plt.imshow(z, origin='lower', interpolation='nearest', vmin=-1e4, vmax=5e4)
plt.title("Data")
plt.subplot(1, 3, 2)
plt.imshow(p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Model")
# We can also plot data - model to look for any systematic deviations of the data
# from the model that might suggest we need a better model or choice of parameters
plt.subplot(1, 3, 3)
plt.imshow(z - p(x, y), origin='lower', interpolation='nearest', vmin=-1e4,
           vmax=5e4)
plt.title("Residual")
plt.show()<jupyter_output><empty_output><jupyter_text>Some affiliated packages  astroquery<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from astroquery.vizier import Vizier

# Cepheids Period-Luminosity data from Bhardwaj et al. 2017
catalog = Vizier.get_catalogs('J/A+A/605/A100')
print(catalog)

period = np.array(catalog[0]['Period'])
log_period = np.log10(period)
k_mag = np.array(catalog[0]['__Ksmag_'])
k_mag_err = np.array(catalog[0]['e__Ksmag_'])

plt.errorbar(log_period, k_mag, k_mag_err, fmt='k.')
plt.xlabel(r'$\log_{10}$(Period [days])')
plt.ylabel('Ks')
plt.show()
# Lets now fit a simple model to the data
model = models.Linear1D()
fitter = fitting.LinearLSQFitter()
best_fit = fitter(model, log_period, k_mag, weights=1.0/k_mag_err**2)

plt.errorbar(log_period,k_mag,k_mag_err,fmt='k.')
plt.plot(log_period, best_fit(log_period), color='g', linewidth=3)
plt.xlabel(r'$\log_{10}$(Period [days])')
plt.ylabel('Ks')
plt.show()<jupyter_output><empty_output><jupyter_text>photutils<jupyter_code>from photutils import datasets
from photutils import aperture_photometry
from photutils import SkyCircularAperture
from astropy import units as u
from astropy.coordinates import SkyCoord
import astropy.io.fits as fits

hdu = fits.open("./data/spitzer.fits") 
catalog = datasets.load_spitzer_catalog()

plt.imshow(hdu[0].data, vmin=-1, vmax=30, origin="lower", cmap=plt.cm.gist_heat)
plt.show()

positions = SkyCoord(catalog['l'], catalog['b'], frame='galactic')
apertures = SkyCircularAperture(positions, r=4.8 * u.arcsec) 
phot_table = aperture_photometry(hdu[0], apertures) 

# conversion to flux per pixel with pixel scale of 1.2 arcsec/pixel
factor = (1.2 * u.arcsec) ** 2 / u.pixel   
converted_aperture_sum = (phot_table['aperture_sum'] *
                          factor).to(u.mJy / u.pixel)

# loading the catalog measurements
fluxes_catalog = catalog['f4_5']

import matplotlib.pyplot as plt
plt.scatter(fluxes_catalog, converted_aperture_sum.value)
plt.xlabel('Spitzer catalog PSF-fit fluxes ')
plt.ylabel('Aperture photometry fluxes')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10923450789793439, 'fraction_numerical': 0.031105710814094776, 'mean_word_length': 5.2403335860500375, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '576737', 'license': 'MIT', 'n_tokens_mistral': 3098, 'n_tokens_neox': 2830, 'n_words': 829, 'path': '6_astropy.ipynb', 'repo_name': 'philuttley/basic_linux_and_coding', 'size': '8230'}