Search code examples
pythonastronomyastropy

Astropy Units equivalencies - interferometry baselines

I am trying to create a astropy.units equivalency to convert between different units of UV-coordinates when working with interferometric astronomical data. The most common way of storing the coordinates is in seconds, however I usually convert directly to lambdas (dependent on the rest-wavelength/frequency). I want to be able to go between : (nano) seconds - (kilo) lambda - meters. The input that is needed for the conversion is the rest-frequency of the individual observations.

Initial description of the method:

What I have come up with so far is the following.

import astropy.units as un
import astropy.constants as co

restfreq_hz = 203e9 #203 Ghz
lambdas = un.def_unit('lambdas', format={'format' : r'\lambda'})
klambdas = un.def_unit('kilolambdas', format={'format' : r'k\lambda'})

# equivalency (from_unit, to_unit, forward, backward)
lambdas_equivalencies = [
    (lambdas, un.s, lambda x: x/restfreq_hz, lambda x: x*restfreq_hz),
    (lambdas, un.m, lambda x: x/restfreq_hz * co.c.to(un.m/un.s).value, lambda x: x/co.c.to(un.m/un.s).value * restfreq_hz),
    (lambdas, un.ns, lambda x: x/restfreq_hz * 1e9, lambda x: x / 1e-9*restfreq_hz ),
    (lambdas, klambdas, lambda x: x*1e-3, lambda x: x*1e3),
    (klambdas, un.s, lambda x: 1e3*x/restfreq_hz, lambda x: 1e-3*x*restfreq_hz),
    (klambdas, un.m, lambda x: 1e3*x/restfreq_hz * co.c.to(un.m/un.s).value, lambda x: 1e-3*x/co.c.to(un.m/un.s).value * restfreq_hz),
    (klambdas, un.ns, lambda x: 1e3*x/restfreq_hz * 1e9, lambda x: 1e-3*x / 1e-9*restfreq_hz ),
    (un.m, un.s, lambda x: x/co.c.to(un.m/un.s).value, lambda x: x*co.c.to(un.m/un.s).value),
    (un.m, un.ns, lambda x: x/co.c.to(un.m/un.ns).value, lambda x: x*co.c.to(un.m/un.ns).value)
]

So as an example I can now do:

In [10]: (100.*klambdas).to(un.m ,equivalencies=lambdas_equivalencies)
Out[10]: <Quantity 147.68101379310343 m>

In [13]: (12 * un.m).to(lambdas, equivalencies=lambdas_equivalencies)
Out[13]: <Quantity 8125.621359026984 lambdas>

In [29]: (1000000*un.ns).to(lambdas, equivalencies=lambdas_equivalencies)
Out[29]: <Quantity 203000000.0 lambdas>

Is this the preferred/best way of doing this, or am I missing something?Any additional tweaks/tips are welcome!

Additional issue(s):

I want to incorporate this into an object. So that I define an array (object attribute) with say unit 'klambda'. Then I want to be able to convert that on-the-fly to 'lambda' or 'm'. Can I do that without redefining the array class?

Solution

  • What you have currently works but you can actually simplify it a little. In particular, if astropy.units already knows how to convert e.g. s to ns then you don't need to define both m to s and m to ns, it will be able to figure it out. To simplify further, you can define klambdas as being defined as a multiple of lambdas. This gives:

    lambdas = un.def_unit('lambdas', format={'format' : r'\lambda'})
klambdas = un.def_unit('kilolambdas', 1e3 * lambdas, format={'format' : r'k\lambda'})

# equivalency (from_unit, to_unit, forward, backward)
lambdas_equivalencies = [
    (lambdas, un.s, lambda x: x/restfreq_hz, lambda x: x*restfreq_hz),
    (lambdas, un.m, lambda x: x/restfreq_hz * co.c.to(un.m/un.s).value, lambda x: x/co.c.to(un.m/un.s).value * restfreq_hz),
    (un.m, un.s, lambda x: x/co.c.to(un.m/un.s).value, lambda x: x*co.c.to(un.m/un.s).value),
]

    In reality you should probably have a function for the equivalency that takes the frequency:

    def lambdas_equivalencies(restfreq_hz):
    eq = [
    (lambdas, un.s, lambda x: x/restfreq_hz, lambda x: x*restfreq_hz),
    (lambdas, un.m, lambda x: x/restfreq_hz * co.c.to(un.m/un.s).value, lambda x: x/co.c.to(un.m/un.s).value * restfreq_hz),
    (un.m, un.s, lambda x: x/co.c.to(un.m/un.s).value, lambda x: x*co.c.to(un.m/un.s).value),
    ]
    return eq

    then use it as 

    (100.*klambdas).to(un.m ,equivalencies=lambdas_equivalencies(restfreq_hz))

    and you should also be able to have restfreq_hz be instead a quantity that you can convert to Hz if needed inside the function::

    def lambdas_equivalencies(restfreq):
    restfreq_hz = restfreq.to(u.Hz, equivalencies=u.spectral())
    ...

    then you can even pass wavelengths, etc.

    For your second question, I think you may have to create a new quantity class that inherits from Quantity and simply overloads to.