aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/utils/convert.py
blob: 02a4ffd908b735b251dde7986ff9c398d60b0e43 (plain)
1
2
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
34
35
36
37
38
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
# Copyright (c) 2016 Weitian LI <liweitianux@live.com>
# MIT license

"""
Utilities for conversion among common astronomical quantities.
"""

import numpy as np
import astropy.units as au
import numba


def Fnu_to_Tb(Fnu, omega, freq):
    """Convert flux density to brightness temperature, using the
    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.

    Parameters
    ----------
    Fnu : `~astropy.units.Quantity`
        Input flux density, e.g., `1.0*au.Jy`
    omega : `~astropy.units.Quantity`
        Source angular size/area, e.g., `1.0*au.sr`
    freq : `~astropy.units.Quantity`
        Frequency where the flux density measured, e.g., `1.0*au.MHz`

    Returns
    -------
    Tb : `~astropy.units.Quantity`
        Brightness temperature, with default unit `au.K`

    References
    ----------
    - Brightness and Flux
      http://www.cv.nrao.edu/course/astr534/Brightness.html
    - Wikipedia: Brightness Temperature
      https://en.wikipedia.org/wiki/Brightness_temperature
    - NJIT: Physics 728: Introduction to Radio Astronomy: Lecture #1
      https://web.njit.edu/~gary/728/Lecture1.html
    - Astropy: Equivalencies: Brightness Temperature / Flux Density
      http://docs.astropy.org/en/stable/units/equivalencies.html
    """
    equiv = au.brightness_temperature(omega, freq)
    Tb = Fnu.to(au.K, equivalencies=equiv)
    return Tb


def Sb_to_Tb(Sb, freq):
    """Convert surface brightness to brightness temperature, using the
    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.

    Parameters
    ----------
    Sb : `~astropy.units.Quantity`
        Input surface brightness, e.g., `1.0*au.Jy/au.sr`
    freq : `~astropy.units.Quantity`
        Frequency where the flux density measured, e.g., `1.0*au.MHz`

    Returns
    -------
    Tb : `~astropy.units.Quantity`
        Brightness temperature, with default unit `au.K`
    """
    omega = 1.0 * au.arcsec**2
    Fnu = (Sb * omega).to(au.Jy)  # [Jy]
    return Fnu_to_Tb(Fnu, omega, freq)


@numba.jit(nopython=True)
def Sb_to_Tb_fast(Sb, freq):
    """Convert surface brightness to brightness temperature, using the
    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.

    This function does the calculations explicitly, and does NOT rely
    on the `astropy.units`, therefore it is faster.  However, the input
    parameters must be in right units.

        Tb = Sb * c^2 / (2 * k_B * nu^2)

    where `SB` is the surface brightness density measured at a certain
    frequency (unit: [ Jy/rad^2 ] = [ erg/s/cm^2/Hz/rad^2 ]).

    Parameters
    ----------
    Sb : float
        Input surface brightness, unit [ Jy/deg^2 ]
    freq : float
        Frequency where the flux density measured, unit [ MHz ]

    Returns
    -------
    Tb : float
        Calculated brightness temperature, unit [ K ]
    """
    # NOTE: `radian` is dimensionless
    rad2_to_deg2 = np.rad2deg(1.0) * np.rad2deg(1.0)
    Sb_rad2 = Sb * rad2_to_deg2  # unit: [ Jy/rad^2 ] -> [ Jy ]
    c = 29979245800.0  # speed of light, [ cm/s ]
    k_B = 1.3806488e-16  # Boltzmann constant, [ erg/K ]
    coef = 1e-35  # take care the unit conversions
    Tb = coef * (Sb_rad2 * c*c) / (2 * k_B * freq*freq)  # unit: [ K ]
    return Tb


@numba.jit(nopython=True)
def Fnu_to_Tb_fast(Fnu, omega, freq):
    """Convert flux density to brightness temperature, using the
    Rayleigh-Jeans law, in the Rayleigh-Jeans limit.

    This function does the calculations explicitly, and does NOT rely
    on the `astropy.units`, therefore it is faster.  However, the input
    parameters must be in right units.

    Parameters
    ----------
    Fnu : float
        Input flux density, unit [ Jy ]
    omega : float
        Source angular size/area, unit [ deg^2 ]
    freq : float
        Frequency where the flux density measured, unit [ MHz ]

    Returns
    -------
    Tb : float
        Calculated brightness temperature, unit [ K ]
    """
    Sb = Fnu / omega
    return Sb_to_Tb_fast(Sb, freq)