<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Using pyND.gbt.GBTspec
This code reads individual datasets from a GBTIDL-style FITS file.
Step1:  Loading GBT spectra from GBTIDL ASCII output
Step2:  Metadata...
Step3:  Loading GBT spectra from GBTIDL FITS format
Loading from the list of objects
Step4:  Loading with an object name
Step5:  Resample the results to a coarser velocity grid
This is a flux-conserving process.
Step6:  Compare the two results
In this example, the results are slightly different, as the ASCII data are saved in the OPTICAL-LSR frame, while the GBTIDL data are saved using the RADI-LSR, the radio astronomical definition of the LSR.
Step7:  The following is not yet working...
Change OPTICAL to RADIO
N.B. This approach doesn't seem to shift the spectra by enough to be the source of the difference...
Step8:  Change RADIO to OPTICAL
Step9:  Change RADIO to OPTICAL with built-in function
<ASSISTANT_TASK:>
Python Code:
# Typical imports here
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
import astropy.constants as c
%load_ext autoreload
%autoreload 2
from pyND.gbt import GBTspec
Explanation: Using pyND.gbt.GBTspec
This code reads individual datasets from a GBTIDL-style FITS file.
End of explanation
input_filename = 'data/HS0033+4300_GBT.dat'
x = GBTspec.from_ascii(input_filename)
x.plotspectrum()
x.velocity[0:5]
x.Tb[0:5]
Explanation: Loading GBT spectra from GBTIDL ASCII output
End of explanation
x.meta.keys()
x.meta['object'],x.meta['RA'],x.meta['DEC']
Explanation: Metadata...
End of explanation
input_filename = 'data/GBTdata.fits'
y = GBTspec.from_GBTIDLindex(input_filename)
y.plotspectrum()
Explanation: Loading GBT spectra from GBTIDL FITS format
Loading from the list of objects
End of explanation
input_filename = 'data/GBTdata.fits'
object_name = 'HS0033+4300'
z = GBTspec.from_GBTIDL(input_filename,object_name)
z.plotspectrum()
Explanation: Loading with an object name:
End of explanation
z_new = z.copy()
new_velocity = np.arange(-400,100,10.)
z_new.resample(new_velocity,masked=True)
plt.figure(figsize=(8,5))
plt.plot(z.velocity,z.Tb,drawstyle='steps-mid',label='Original')
plt.plot(z_new.velocity,z_new.Tb,drawstyle='steps-mid',label='Resampled',lw=4)
plt.xlim(-100,50)
plt.legend(loc='upper left')
plt.xlim(-200,50)
plt.xlabel('Velocity')
plt.ylabel('Tb [K]');
Explanation: Resample the results to a coarser velocity grid
This is a flux-conserving process.
End of explanation
plt.plot(z.velocity,z.Tb,drawstyle='steps-mid',label='RADIO')
plt.plot(x.velocity,x.Tb,drawstyle='steps-mid',label='OPTICAL')
plt.xlim(-100,50)
plt.legend(loc='upper left')
Explanation: Compare the two results
In this example, the results are slightly different, as the ASCII data are saved in the OPTICAL-LSR frame, while the GBTIDL data are saved using the RADI-LSR, the radio astronomical definition of the LSR.
End of explanation
light_speed = np.float64(c.c.to('m/s').value)
nu0 = np.float64(1420405800.0000000000000000000)
# First calculate frequency from optical:
nu = (nu0/(1+(x.velocity)*1000./light_speed))
# Calculate radio definition
vrad =light_speed*((nu0-nu)/nu0)/1000.
plt.figure(figsize=(8,5))
plt.plot(y.velocity,y.Tb,drawstyle='steps-mid',label='RADIO')
plt.plot(vrad,x.Tb,drawstyle='steps-mid',label='OPTICAL-->RADIO',zorder=0,linewidth=3)
plt.xlim(-100,50)
plt.legend(loc='upper left')
Explanation: The following is not yet working...
Change OPTICAL to RADIO
N.B. This approach doesn't seem to shift the spectra by enough to be the source of the difference...
End of explanation
light_speed = (c.c.to('m/s').value)
nu0 = (1420405800.0000000000000000000)
# Frequency from radio:
nu = nu0*(1-(x.velocity)*1000./light_speed)
# Calculate radio definition
vopt = (light_speed/1000.)*((nu0-nu)/nu)
plt.plot(vopt,y.Tb,drawstyle='steps-mid',label='RADIO-->OPTICAL')
plt.plot(x.velocity,x.Tb,drawstyle='steps-mid',label='OPTICAL',zorder=0,linewidth=3)
plt.xlim(-100,50)
plt.legend(loc='upper left')
Explanation: Change RADIO to OPTICAL
End of explanation
input_filename = 'data/AMIGA-GBT.fits'
object_name = 'RBS2055'
y = GBTspec.from_GBTIDL(input_filename,object_name)
input_filename = 'data/RBS2055_GBT.dat'
x = GBTspec.from_ascii(input_filename)
x.change_veldef()
plt.plot(y.velocity,y.Tb,drawstyle='steps-mid',label='RADIO-->OPTICAL')
plt.plot(x.velocity,x.Tb,drawstyle='steps-mid',label='OPTICAL',zorder=0,linewidth=3)
plt.xlim(-100,50)
plt.legend(loc='upper left')
Explanation: Change RADIO to OPTICAL with built-in function
End of explanation
<END_TASK>