Quiero hacer una proyección de púlsares de milisegundos en el plano galáctico, muy parecida a esta de Sala et al. 2004 :
He probado varios métodos y no he llegado a ninguna parte. Este es mi código actual, junto con lo que está produciendo:
import math as m
import numpy as np
import csv
import matplotlib.pyplot as plt
import pandas as pd
import astropy.coordinates as coord
import astropy.units as u
from astropy.io import ascii
from astropy.coordinates import SkyCoord
data = pd.read_csv('galacticwperiod.csv')
xarr = np.array(data.iloc[:,0])
yarr = np.array(data.iloc[:,1])
eq = SkyCoord(xarr[:], yarr[:], unit=u.deg)
gal = eq.galactic
#print(xarr)
#xarr = np.array(df.iloc[:,0])
#yarr = np.array(df.iloc[:,1])
#zarr = np.array(df.iloc[:,2])
#ra = coord.Angle(xarr[:], unit=u.hour)
#ra.degree
#ra = ra.wrap_at(180*u.degree)
#dec = coord.Angle(yarr[:], unit=u.deg)
#print(ra)
plt.figure(figsize=(6,5))
fig = plt.figure()
ax = fig.add_subplot(111, projection="aitoff")
plt.plot(gal.l.wrap_at(180*u.deg), gal.b.wrap_at(180*u.deg), linestyle='None')
ax.scatter(gal.l, gal.b, linestyle='None')
#ax.set_facecolor('xkcd:battleship grey')
#fig.patch.set_facecolor('white')
#ax.tick_params(axis='both', which='major', labelsize=10)
#ax.grid(color='b', linestyle='solid')
fig.show()
#plt.savefig('millisecondcoloraitoff.png', dpi=600)
Aquí hay algunas líneas del archivo de entrada 'galacticwperiod.csv':
Gl,Gb
111.383,-54.849
305.898,-44.888
305.913,-44.877
Esta es la imagen que produce:
Estoy casi seguro de que esto está mal porque no están distribuidos a lo largo del plano galáctico, como deberían estarlo. Los datos que estoy usando son del catálogo ATNF Pulsar .
Estos son los sitios que ya he mirado, como referencia:
https://stackoverflow.com/questions/33105898/astropy-matplotlib-and-plot-galactic-coordinates https://astropy4scipy2014.readthedocs.io/_static/notebooks/06_Celestial_Coordinates_solutions.html
Cualquier ayuda con esto sería muy apreciada.
Este código lee las coordenadas como ecuatoriales (ra, dec) y las transforma en galácticas (l, b):
eq = SkyCoord(xarr[:], yarr[:], unit=u.deg)
gal = eq.galactic
El contenido de 'galacticwperiod.csv' ya está en coordenadas galácticas y no debe transformarse. Algo como esto puede dar mejores resultados:
gal = SkyCoord(xarr[:], yarr[:], frame='galactic', unit=u.deg)
El otro problema es que las proyecciones geográficas de pyplot parecen esperar ángulos en radianes. Este código:
plt.subplot(111, projection='aitoff')
plt.grid(True)
plt.scatter(gal.l.wrap_at('180d').radian, gal.b.radian)
produce esta trama:
Mirando la trama que tienes, noto que hay una concentración de puntos en los "polos". La concentración debería estar a lo largo del ecuador, y particularmente hacia el centro de la galaxia.
Creo que ha cambiado sus ejes, por lo que los puntos a lo largo del ecuador galáctico se dibujan a lo largo del "meridiano principal" vertical.
from astropy.io import fits
from astropy import wcs
import matplotlib.pyplot as plt
x = fits.open('f160w_noopt_drz.fits')['SCI']
y = wcs.WCS(x.header)
plt.subplots(figsize=(8, 8), subplot_kw={'projection':y})
plt.imshow(x.data, cmap='viridis', origin='lower', vmin=0, vmax=10)
plt.show()
james k
María
Kornpob Bhirombhakdi