Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Nishikawa/kisa_io/#71 #95

Merged
merged 9 commits into from
Apr 22, 2021
Merged

Nishikawa/kisa_io/#71 #95

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
9615af5
#67#71#72 resolve all those issues
KaoruNishikawa Apr 21, 2021
2eac70f
#67 support numpy datetime format
KaoruNishikawa Apr 21, 2021
f34f08e
#67 prepare for future implementation
KaoruNishikawa Apr 21, 2021
68c3e01
#71 add dependency
KaoruNishikawa Apr 21, 2021
a1ca802
#71 fix typo
KaoruNishikawa Apr 22, 2021
d0c985f
#71 fix typo
KaoruNishikawa Apr 22, 2021
328c3b4
#71 change reference of constants
KaoruNishikawa Apr 22, 2021
2760715
#71 change reference of constants
KaoruNishikawa Apr 22, 2021
31b7c5a
#71 fix invalid condition
KaoruNishikawa Apr 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 56 additions & 58 deletions nasco_analysis/doppler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,12 @@
from functools import wraps
from datetime import datetime

from astropy.coordinates import EarthLocation
from astropy.coordinates import SkyCoord, CartesianDifferential, LSR
from astropy.time import Time
import astropy.constants as const
from astropy.time import Time
import astropy.units as u
import numpy as np
import n_const.constants as n2const


class Doppler(object):
Expand All @@ -24,12 +24,12 @@ class Doppler(object):
rest_freq : astropy.units.quantity.Quantity
Rest frequency of the line spectrum to be evaluated.
species : str
Name of the observed species; only CO isotopes (12CO_10,
13CO_10, C18O_10, 12CO_21, 13CO_21, C18O_21) are supported.
Name of the observed species; only CO isotopes (J10_12CO,
J10_13CO, J10_C18O, J21_12CO, J21_13CO, J21_C18O) are supported.
LO1st_freq : astropy.units.quantity.Quantity
Frequency of 1st local oscilator.
Frequency of 1st local oscillator.
LO1st_factor : int or float
Factor of frequency multiplier for 1st local oscilator.
Factor of frequency multiplier for 1st local oscillator.
LO2nd_freq : astropy.units.quantity.Quantity
Frequency of 2nd local oscilator.
obstime : float, datetime.datetime or astropy.units.quantity.Quantity
Expand Down Expand Up @@ -72,25 +72,18 @@ class Doppler(object):
<Quantity -36.03406182 km / s>

""" # noqa: E501

__VARS = [
'spectrometer', 'rest_freq', 'LO1st_freq', 'LO1st_factor',
'LO2nd_freq', 'obstime', 'ra', 'dec', 'species'
"spectrometer",
"rest_freq",
"LO1st_freq",
"LO1st_factor",
"LO2nd_freq",
"obstime",
"ra",
"dec",
"species",
]
# observatory location
LOC_NANTEN2 = EarthLocation(
lon=-67.70308139 * u.deg,
lat=-22.96995611 * u.deg,
height=4863.85 * u.m,
)
SPECIES_to_REST_FREQ = {
'12co_10': 115.271202 * u.GHz,
'13co_10': 110.201353 * u.GHz,
'c18o_10': 109.782173 * u.GHz,
'12co_21': 230.538000 * u.GHz,
'13co_21': 220.398681 * u.GHz,
'c18o_21': 219.560354 * u.GHz,
# reference: Naomasa Nakai et al. 2009, ISBN978-4-535-60766-8
}

def __init__(self, args=None, **kwargs):
# set instance variables
Expand All @@ -117,20 +110,20 @@ def set_args(self, args=None, **kwargs):
if key in self.__VARS:
setattr(self, key, val)
else:
raise NameError(f'invalid argument : {key}')
raise NameError(f"invalid argument : {key}")
return

def __species2freq(func):
@wraps(func)
def translate(inst, *args):
if hasattr(inst, 'species'):
inst.rest_freq = inst.SPECIES_to_REST_FREQ[inst.species.lower()]
if hasattr(inst, "species"):
inst.rest_freq = n2const.REST_FREQ[inst.species.lower()]
return func(inst, *args)

return translate

def __check_args(self, args):
"""

Check if the parameters needed to run function are declared.

Parameters
Expand All @@ -141,7 +134,6 @@ def __check_args(self, args):
------
AttributeError
When undeclared parameter exists.

"""
undeclared = []
for arg in args:
Expand All @@ -150,9 +142,7 @@ def __check_args(self, args):
except AttributeError:
undeclared.append(arg)
if undeclared:
raise AttributeError(
f'Parameter {undeclared} is not given. Use `set_args`'
)
raise AttributeError(f"Parameter {undeclared} is not given. Use `set_args`")
return

@__species2freq
Expand Down Expand Up @@ -193,18 +183,22 @@ def heterodyne(self):
""" # noqa: E501
# check if necessary parameters have already been declared
__vars = [
'spectrometer', 'rest_freq', 'LO1st_freq', 'LO1st_factor', 'LO2nd_freq'
"spectrometer",
"rest_freq",
"LO1st_freq",
"LO1st_factor",
"LO2nd_freq",
]
self.__check_args(__vars)
# set spectrometer constants
if self.spectrometer.lower() == 'xffts':
self.ch_num = 32768
self.band_width = 2 * u.GHz
elif self.spectrometer.lower() == 'ac240':
self.ch_num = 16384
self.band_width = 1 * u.GHz
if self.spectrometer.lower() == "xffts":
self.ch_num = n2const.XFFTS.ch_num
self.band_width = n2const.XFFTS.bandwidth
elif self.spectrometer.lower() == "ac240":
self.ch_num = n2const.AC240.ch_num
self.band_width = n2const.AC240.bandwidth
else:
raise(ValueError(f"Invalid spectrometer name : {self.spectrometer}"))
raise (ValueError(f"Invalid spectrometer name : {self.spectrometer}"))
# calculate frequency shift
freq_after1st = self.rest_freq - self.LO1st_freq * self.LO1st_factor
self.sideband_factor = np.sign(freq_after1st)
Expand Down Expand Up @@ -250,13 +244,14 @@ def ch_speed(self):
spec_freq = self.heterodyne()
freq_resolution = self.band_width / self.ch_num
speed_resolution = -1 * self.band_width / self.rest_freq / self.ch_num * const.c
v_0GHz = -1 * speed_resolution * spec_freq \
/ freq_resolution * self.sideband_factor
v_0GHz = (
-1 * speed_resolution * spec_freq / freq_resolution * self.sideband_factor
)
v_base = np.arange(self.ch_num) * speed_resolution * self.sideband_factor
v_apparent = v_base + v_0GHz
v_obs = self.v_obs()
v_lsr = v_apparent + v_obs
return v_lsr.to(u.km/u.s)
v_lsr = v_apparent + np.atleast_2d(v_obs).T
return v_lsr.to(u.km / u.s)

def v_obs(self):
"""Observer's verocity relative to LSR.
Expand Down Expand Up @@ -294,7 +289,7 @@ def v_obs(self):

""" # noqa: E501
# check if necessary parameters have already been declared
__vars = ['obstime', 'ra', 'dec']
__vars = ["obstime", "ra", "dec"]
self.__check_args(__vars)

# galactic rotation parameter (solar motion)
Expand All @@ -303,40 +298,43 @@ def v_obs(self):
ra=18 * 15 * u.deg,
dec=30 * u.deg,
frame="fk4",
equinox=Time('B1900'),
equinox=Time("B1900"),
).galactic
U = v_sun * np.cos(dir_sun.b) * np.cos(dir_sun.l)
V = v_sun * np.cos(dir_sun.b) * np.sin(dir_sun.l)
W = v_sun * np.sin(dir_sun.b)
v_bary = CartesianDifferential(U, V, W)

# set target position
if isinstance(self.obstime, datetime):
obstime_repr = list(np.array(self.obstime))[0]
if isinstance(obstime_repr, (datetime, np.datetime64)):
self.obstime = Time(self.obstime)
if not isinstance(self.obstime, Time): # to interpret both int and float
self.obstime = Time(self.obstime, format='unix')
if not isinstance(obstime_repr, Time): # to catch both int and float
self.obstime = Time(self.obstime, format="unix")
target = SkyCoord(
ra=self.ra,
dec=self.dec,
frame="fk5",
obstime=self.obstime,
location=self.LOC_NANTEN2,
location=n2const.LOC_NANTEN2,
)

# calculate V_obs
v_obs = SkyCoord(
self.LOC_NANTEN2.get_gcrs(self.obstime)
).transform_to(
LSR(v_bary=v_bary)
).velocity
v_obs = (
SkyCoord(n2const.LOC_NANTEN2.get_gcrs(self.obstime))
.transform_to(LSR(v_bary=v_bary))
.velocity
)

# normal component of V_obs to stellar body
v_correction = v_obs.d_x * np.cos(target.icrs.dec) * np.cos(target.icrs.ra) + \
v_obs.d_y * np.cos(target.icrs.dec) * np.sin(target.icrs.ra) + \
v_obs.d_z * np.sin(target.icrs.dec)
v_correction = (
v_obs.d_x * np.cos(target.icrs.dec) * np.cos(target.icrs.ra)
+ v_obs.d_y * np.cos(target.icrs.dec) * np.sin(target.icrs.ra)
+ v_obs.d_z * np.sin(target.icrs.dec)
)

return v_correction


if __name__ == '__main__':
if __name__ == "__main__":
pass
Loading