• Author / Uploaded
• Felipe Ávila Vera

Citation preview

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Universidad De Concepci´ on Software y An´ alisis de Datos Astronomicos

Tarea y Laboratorio 5

Alumno: Felipe Avila

25 de abril de 2019

Software y analisis de datos astronomicos

1

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Laboratorio 1. Encuentren el articulo que define los archivos FITS. Cual es el primer autor del articulo? Cuantas citas tiene el articulo, y cuantas veces a sido “leido”? Que forma aproximada tienen los histogramas de citas y lecturas por a˜ no de este articulo? Creen que el numero total de citas a este articulo es razonable? Por que? Lean secci´ on 4 del articulo. 2. Usando la imagen “ch1 merged clean.fits”, averig¨ uen dentro de Python a que coordenadas (x,y) corresponden estas coordenadas absolutas: (R.A., Decl.) = 5:35:16.340, -5:23:20.07 (R.A., Decl.) = 5:35:14.134,-5:22:22.91 3. Usando la imagen OrionA all spire250 nh mask corr apex.fits y la lista de coordenadas de protoestrellas hops 1.5.reg, conviertan las coordenadas de protoestrellas a coordenadas (x,y) en la imagen. Den las coordenadas (x,y) de las fuentes 102, 099, 057, 235, 243. Incluyan los c´odigos de E2 y E3. 4. E4: Make a figure of thee N(H) map like in the example, with proper coordinates. Overplot contour levels 1E22 y 2E23 cm -2

5. Extra credit: can you figure out automatically which protostars do not fall within the OrionA all spire250 nh mask corr a map area? Where are these sources located? What cloud? (Hint: You can figure this out by looking in Stutz et al. 2013, for example). 1 2 3

#Pregunta 1 #El papper que define fits es "FITS: A flexible image transport system", escrito publicad en 1981 y con primer autor D.C Wells,ha sido citado 264 veces con su pick maximo el ao 2018 y haciendo una mirada a grueso modo podemos observar 2 gausianas con picks en 1995 y la otra tendria que ser dentro de los proximos 10 aoo,registra 5330 lecturas y podriamos decir que ha grandes rasgos las lecturas tienen una forma gausiana centrada en el 2009 con su maximo pick.Creo que es demciado poco el numero de citaciones para un articulo que define el formato mas utilizado por la astronomia moderna a la hora de trabajar los datos. Habla sobre las minimos keywords que tiene que tener una imagen tipo fits, para describir correctamente la imagen.

Software y analisis de datos astronomicos

2

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 1: Histograma de las citas al papper

Figura 2: Histograma de las lecturas del documento

Software y analisis de datos astronomicos

3

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 3: Imagen del papper original

1 2 3 4 5

6 7 8 9

10 11 12

13 14 15

#Pregunta 2 from astropy import units as u import numpy as np from astropy.coordinates import Angle, SkyCoord from regions import CircleSkyRegion, write_ds9,read_ds9, RectangleSkyRegion, PixCoord, CirclePixelRegion from astropy.io import fits from astropy.wcs import WCS header =fits.getheader(’ch1_merged_clean.fits’) c1 = SkyCoord(’5:35:16.340 -5:23:20.07’, unit=(u.hourangle, u.deg),frame=’fk5’) #ASI ESCRIBES EL COORDENADAS c2 = SkyCoord(’5:35:14.134 -5:22:22.91’, unit=(u.hourangle, u.deg),frame=’fk5’) w=WCS(header) #x=1819,y=2187 c1x, c1y=w.wcs_world2pix(c1.ra.deg,c1.dec.deg,1) #poner 1 es importante para no tner grandes errores c2x, c2y=w.wcs_world2pix(c2.ra.deg,c2.dec.deg,1) #hay que separar las coordenadas ra y dec print(’Coordenadas x,y de primera coordenada ’,c1x,c1y) print(’Coordenadas x,y de segunda coordenada ’,c2x,c2y)

Software y analisis de datos astronomicos

4

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 4: Cambio de coordenadas

1 2 3 4 5 6

7 8 9

#Pregunta 3 from astropy import units as u import numpy as np from astropy.coordinates import Angle, SkyCoord from astropy import wcs from regions import CircleSkyRegion, write_ds9,read_ds9, RectangleSkyRegion, PixCoord, CirclePixelRegion from astropy.io import fits from astropy.wcs import WCS import sys

10 11 12 13 14 15 16 17 18 19 20

filename = ’hops_1.5.reg’ regions = read_ds9(filename) w = WCS(’OrionA_all_spire250_nh_mask_corr_apex.fits’) #PARA VOLVER LA IMAGEN WCS print(’Las coordenadas de 102 sera =’,regions[74].to_pixel(w).center) print(’Las coordenadas de 099 sera =’,regions[71].to_pixel(w).center) print(’Las coordenadas de 057 sera =’,regions[39].to_pixel(w).center) print(’Las coordenadas de 235 sera =’,regions[182].to_pixel(w).center) print(’Las coordenadas de 243 sera =’,regions[191].to_pixel(w).center) #Den las coordenadas (x,y) de las #fuentes 102, 099, 057, 235, 243.

Software y analisis de datos astronomicos

5

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 5: Cambio de coordenadas para la region

1 2 3 4 5 6

#Pregunta 4 import matplotlib.pyplot as plt import numpy as np from astropy.wcs import WCS from astropy.io import fits from astropy.utils.data import get_pkg_data_filename

7 8 9 10 11 12 13 14 15 16 17 18 19 20

filename = ’OrionA_all_spire250_nh_mask_corr_apex.fits’ hdu = fits.open(filename)[0] wcs = WCS(hdu.header) H=hdu.data H=H.astype(’float’) H[H==0] =np.nan ax=plt.subplot(projection=wcs) #PUEDES VER LAS POSICIONES EN EL EJE USANDO WCS ax.imshow(H, vmin=0, vmax=1e+22, origin=’lower’) #vmax nivel de contraste ax.grid(color=’white’, ls=’solid’) ax.contour(H, levels=(1e+22,3e+22), colors=’black’, alpha=1) ax.set_xlabel(’RA’) ax.set_ylabel(’DEC’) plt.show()

Software y analisis de datos astronomicos

6

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 6: Imagen de la nebulosa con contornos

Tarea 1. E1: Usando simbad, averig¨ uen a que objetos corresponden estas coordenadas absolutas: a) (R.A., Decl.) = 5:35:16.340, -5:23:20.07 b) (R.A., Decl.) = 5:35:14.134,-5:22:22.91

Respuesta Al buscar las coordenadas en simbad, obtenemos una lista de objetos y dentro del top 20 para la primera coordenada el objeto mas referenciado es NAME Trapezium Cluster con 1412, para las coordenadas b) el objeto mas referenciado es NAME Orion-KL con 1897.(Las imagenes de los resultados de las busquedas a continuacion) a) http://simbad.u-strasbg.fr/simbad/sim-id?Ident=%40808485&Name=NAME%20Trapezium%20Cluster& submit=submit b) http://simbad.u-strasbg.fr/simbad/sim-id?Ident=%40808511&Name=NAME%20Orion-KL&submit= submit

Figura 7: a) coord 5:35:16.340 -5:23:20.07

Software y analisis de datos astronomicos

7

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 8: b) coord 5:35:14.134 -5:22:22.91 2. Write a brief paragraph on what each object is. Hint: the KL object is also called the “Becklin-Neugebauer Object”.

Respuesta a) Las primeras coordenadas, nos dan como resultado el C ¸ umulo trapecio”, que es un cumulo abierto ubicado en la nebulosa de Orion. Este es relativamente joven y posee un diametro de ∼ 30kpc. b) El segundo set es el que correspondiente a el Becklin-Neugebauer object. Fue descubierto en 1967 y se estima que puede ser una proto-estrella con una masa de ∼ 7M 3. If BN-KL is moving in the plane-of-the sky at about 30 km/s, what is it’s proper motion at the distance of Orion in units of milli-as/yr (milli-arcsec per year)? (Eventually we will want to compare the Gaia proper motion sensitivity to the stars in Orion.)

Respuesta 1

import numpy as np

2 3 4

def stoyr(x): #segundos a aos return x/(3600*24*365.25)

5 6

7 8 9 10 11 12 13

#v_t=4.74*u*d_ori ecuacion que relaciona distancia en pc, velocidad tangencial y movimiento propio #v_t/(4.74*d_ori) despejada v_t=30#km/s d_ori=412#pc u=v_t/(4.74*d_ori)#arcsec/s #pero la profesora pide en miliarcosegundos/yr print(’movimient propio u=’,stoyr(u*1000),’mas/yr’) #Resultado en pantalla "movimient propio u= 4.867900917725965e-07 mas/yr"

4. Encuentren el articulo reciente sobre la explosion en BNKL basado en datos de ALMA (Bally et al. 2017). ˆ ˆ Para buscar por primer autor en ADS, beben incluir un delante del apellido, e.g.: bally, j.Basado en este articulo, respondan a las siguientes preguntas.

Software y analisis de datos astronomicos

8

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Respuesta El articulo nombrado es el siguiente https://iopscience.iop.org/article/10.3847/1538-4357/aa5c8b/ pdf a) Cuando ocurri´ o la explosion? La explosion se data hace unos 500 a˜ nos atras. b) A que velocidad maxima se esta moviendo el material en la explosion? Se registraron jets de gas que viajan a velocidades maximas de 150

km s

debido a la explosion

c) Coincide el centro de la explosion con nuestras coordenadas de arriba? El centro de la explosion es en 05:35:14.11, 5: 22:18.7 coordenada cercana a la coordenadas b) de la pregunta 1. d ) La conversion mas fundamental entre velocidad, tiempo y distancia es f´acil de calcular. 1 km/s es igual a cuantos pc/Myr? Deben de redondear el resultado al numero entero mas cercano. Recuerden esta conversion. 1 2

#pregunta d import numpy as np

3 4 5 6 7 8 9

def kmtopc(x): return x/(3.086e+13) def stoMyr(x): return x/(1e6*3600*24*365.25) def kmstopcMyr(x): return stoMyr(x)/kmtopc(x) # se aplica stomyr/kmtopc y se obtiene el resultado

10 11

12

print(kmstopcMyr(1),’usualmente se aproxima a 1km/s a 1pc/Myr pero por aproximaciones dentro del calculo no obtengo una cifra mas precisa’) #resultado en consola "0.9778943899409335 usualmente se aproxima a 1km/s a 1pc/Myr pero por aproximaciones dentro del calculo no obtengo una cifra mas precisa"

e) Para la velocidad maxima, y una escala de tiempo de 500 a˜ nos, cuantos pc (y AU) puede cubrir el materia en la explosion? 1 2 3 4 5 6 7 8 9 10

11 12 13 14 15 16

#Pregunta e import numpy as np def kmtopc(x): return x/(3.086e+13) def yrtos(x): return x*(3600*24*365.25) def pctoUA(x): return x*206264.8062454 #Para la velocidad maxima, y una escala de tiempo de 500 a nos, cuantos pc (y AU) puede cubrir elmateria en la explosion? #con v_max=150km/s v_max=150#km/s t_yr=500#yr t_s=yrtos(t_yr) d_max=v_max*t_s #rapidez es v=d/t entonces d=v*t print(’la distancia maxima en pc sera’,kmtopc(d_max),’pc’,’y la distancia maxima en UA sera’,pctoUA(kmtopc(d_max)),’UA’)

Software y analisis de datos astronomicos

9

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

17

Universidad de Concepci´ on

#Resultado en consola sera "la distancia maxima en pc sera 0.0766953985742061 pc y la distancia maxima en UA sera 15819.561526822346 UA"

f ) Que linea estan observando Bally et al. con ALMA? Con que resoluci´on angular? Cual es la resolucion espacial a la distancia de Orion? Como se compara esta a la resolucion angular y espacial de los datos de 13CO de la semana pasada? Cual es su frecuencia en GHz de los datos de Bally et al.? Y rango de onda en mm? ”The OMC1 outflow was mapped with the Atacama Large Millimeter Array (ALMA) in the 1.3 mm atmospheric window (ALMA Band 6) using four 1.875 GHz wide bands centered at 216, 228, 231, and 233 GHz” Entonces observan en las bandas 216,228,231 y 233 GHz. Usando antenas alrededor de un diametro de 1km, ahora procedo a realizar los calculos solicitados y se comparan con los de la semana anterior, como comentario del programa. 1 2

#pregunta f import numpy as np

3 4

5 6

7 8 9 10 11 12 13 14 15 16

17 18

19

20 21 22 23 24 25

26 27 28 29 30 31 32 33

def AngRes(lamda,D):#Calcula la resolucion angular a partir de la distancia de apertura(D) y la longitud de onda(lamda), ambas en la misma unidad de radianes return (1.22 * lamda/D) #angulo en radianes def MaxDist(Size,angle):#Se le entrega una resolucion angular(angle) y el tamao de la fuente(Size), ofrece la mayor distancia a la que se puede resolver return Size/np.tan(angle) #metros def RadsToArcmin(x): #ingreso angulo en radianes return x*(60*180)/(np.pi) #retorna angulo en arcominutos def RadianToArcsec(x): #ingreso angulo en radianes return x*(3600*180)/(np.pi) #arcosegundo def ArcsecToRadian(x): #ingrso arco segundos return (x*np.pi)/(3600*180) # devuelve radianes def frectolamda(x): #frecuencia en 1/s=hz return 3e8/x #en metros def AperNeed(lamda,Res):#Apertura necesaria dada por una longitud de onda y una cierta resolucion return 1.22*lamda/Res #resultado en metros def resolucionespacial(d,angle):#Entrega tamao del objeto respecto a la resolucion angular entregada d es la distancia a el obj return d*np.tan(angle) #La distancia devuelta depende de las unidades usadas por d def PctoUA(x): #Transforma Pc a UAs return x/4.848136e-6 def wienT (lamda): #Ley de wien return 2.897e-3/lamda #Devuelve temperaturas en kelvin #d= DIAMETRO Y L= lamda #using four 1.875 GHz wide bands centered at 216, 228, 231, and233 GHz, distancia antenas=1000m d_ant=1000#m f_T=np.array([2.16,2.28,2.31,2.33])*(10**11)#hz l_T=frectolamda(f_T) #longitud de onda en m A_T=AngRes(l_T,d_ant) #en radianes d_ori=400#pc resolucionespacial_T=resolucionespacial(d_ori,A_T) #en pc e9=1e9 e3=1e3

34

Software y analisis de datos astronomicos

10

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

35 36 37 38

Universidad de Concepci´ on

print(f_T/e9,’Frecuencias en GHz’) print(l_T*e3,’longitudes de onda en mm’) print(A_T,’resoluciones angulares en rad’) print(resolucionespacial_T,’resoluciones espaciales en pc’)

39 40 41 42 43

44 45

46

47

48 49 50 51 52

53

54

55 56 57 58

59 60 61

62 63 64 65

#Pregunta semana pasada print(’Todo semana pasada’) resolucionangular_2D=ArcsecToRadian(8)#radianes f_2d=np.array([115.27,110.2])*(10**9) #Frecuencia Co12 (1-0) 115 Ghz,Frecuencia CO13 (1-0) 110 Ghz ambas en hz apminima_2d=AperNeed(frectolamda(f_2d),resolucionangular_2D) print(’Las aperturas necesarias serian’,apminima_2d,’en metros de 12CO y 13CO respectivamente’) # Los tamaos de las antenas de carma son de3.5-m, 6.1-m, and 10.4-m y NRO es de 45m, pero como se hace por medio de la interferometria #la resolucion aumentara, haciendo posible resolver estas longitudes de onda sin necesidad de utilizar un lente de 80m de diametro #Ejercicio 3 print(’Ejercicio 3’) d_e3=400 #pc resolucionespacial_2D=resolucionespacial(d_e3,resolucionangular_2D) print(resolucionangular_2D,’resolucion angular en radianes de pregunta semana pasada ’) print(’Resolucion espacial 2d’,resolucionespacial_2D,’en pc’,PctoUA(resolucionespacial_2D),’en UA’) resolucionespacial_3=resolucionespacial(d_e3,2*np.pi/(360)) #como 2pirad es 360 grados entonces 1grado=2pi/360 radianes #Resultados en consola [216. 228. 231. 233.] Frecuencias en GHz [1.38888889 1.31578947 1.2987013 1.28755365] longitudes de onda en mm [1.69444444e-06 1.60526316e-06 1.58441558e-06 1.57081545e-06] resoluciones angulares en rad [0.00067778 0.00064211 0.00063377 0.00062833] resoluciones espaciales en pc Todo semana pasada Las aperturas necesarias serian [81.86531522 85.63171403] en metros de 12CO y 13CO respectivamente Ejercicio 3 3.878509448876288e-05 resolucion angular en radianes de pregunta semana pasada Resolucion espacial 2d 0.01551403780328432 en pc 3200.0005369660257 en UA #Finalmente en general son similares en resolucion angular, resolucion espacial, aun que sus frecuencias tengan una diferencia de unos 100GHz

g) Podr´ a Gaia observar esta linea? Por que? Dificilmente puesto que gaia estudia el espectro visible y estas lineas estan en el espectro infrarojo. h) Hagan una region en DS9 indicando el centro de la explosion, y usen DS9 para comparar esta con la imagen de 13CO y N(H). Dada la velocidad de la explosion y la escala de tiempo (el tiempo pasado desde que ocurri´ o le explosion), puede llegar material de la explosion a 9’ al sur de la explosion? 05:35:14.11, 5: 22:18.7,tienes q scar el tama˜ no angular de cada pixel luego ver el centro en pixel y compararlo con el lugar dnd supuestamente tiene q estar es la mejor idea q tengo

Software y analisis de datos astronomicos

11

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

Figura 9: Region de 9’ en ambas imagenes

Figura 10: region de 9’ en 13CO

Figura 11: region de 9’ en n(h)

1 2 3

import numpy as np def radtoarmic(x): return x*60*180/(np.pi) #de radianes a arcominutos

Software y analisis de datos astronomicos

12

Tarea V

Facultad de Ciencias F´ısicas y Matem´ aticas. Departamento de Astronomia.

Universidad de Concepci´ on

4 5 6 7 8 9

d=400#distancia en pc a nube de orion dt=0.0766953985742061#pc angulo=np.arctan(dt/d) print(radtoarmic(angulo),’angulo maximo de recorrido en arcominutos’ ) #resultado en consola ’0.6591483888787851 angulo maximo de recorrido en arcominutos’, entonces no llegara a los 9’

Software y analisis de datos astronomicos

13

Tarea V