From 402e6080bda7cf7b5030f97a1ef9bf7cc3ef3ce7 Mon Sep 17 00:00:00 2001
From: Jae-Joon Lee <lee.j.joon@gmail.com>
Date: Thu, 20 May 2021 12:17:43 +0900
Subject: [PATCH] [galactic_extinction] initial commit for lallerment2019

---
 galactic_extinction/README.md                 | 11 +++
 galactic_extinction/setup.py                  | 12 +++
 galactic_extinction/src/__init__.py           |  0
 .../src/galactic_extinction_lallerment2019.py | 87 +++++++++++++++++++
 4 files changed, 110 insertions(+)
 create mode 100644 galactic_extinction/README.md
 create mode 100644 galactic_extinction/setup.py
 create mode 100644 galactic_extinction/src/__init__.py
 create mode 100644 galactic_extinction/src/galactic_extinction_lallerment2019.py

diff --git a/galactic_extinction/README.md b/galactic_extinction/README.md
new file mode 100644
index 0000000..953f77f
--- /dev/null
+++ b/galactic_extinction/README.md
@@ -0,0 +1,11 @@
+This is small helper module for Galactic extinction data.
+For now, it includes helper for 3d cube data from Lallement+2019.
+
+# Install
+It has not been fully packaged yet. So, you need to manually copy the python file (it is a single file for now). And point it with a data path.
+
+The data can be downloaded from
+
+[lallerment2019_cube_converted.fits](http://210.110.233.67:9005/various-data-jjl/lallermant2019/lallerment2019_cube_converted.fits)
+
+
diff --git a/galactic_extinction/setup.py b/galactic_extinction/setup.py
new file mode 100644
index 0000000..9ed7386
--- /dev/null
+++ b/galactic_extinction/setup.py
@@ -0,0 +1,12 @@
+from distutils.core import setup
+
+setup(name='galactic_extinction',
+      version='0.1',
+      description='JJ\'s small python modules for astronomy',
+      author='Jae-Joon Lee',
+      author_email='lee.j.joon@gmail.com',
+      url='https://data.kasi.re.kr/gitlab/leejjoon/small-astronomy-modules',
+      packages=['galactic_extinction'],
+      package_dir={'galactic_extinction':'src'}
+     )
+
diff --git a/galactic_extinction/src/__init__.py b/galactic_extinction/src/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/galactic_extinction/src/galactic_extinction_lallerment2019.py b/galactic_extinction/src/galactic_extinction_lallerment2019.py
new file mode 100644
index 0000000..a245780
--- /dev/null
+++ b/galactic_extinction/src/galactic_extinction_lallerment2019.py
@@ -0,0 +1,87 @@
+import numpy as np
+import pandas as pd
+
+import astropy.io.fits as pyfits
+from scipy.interpolate import RegularGridInterpolator
+from astropy.coordinates import SkyCoord, Galactic
+
+
+def _load_data():
+    fn = "../STILISM_cube.fits.gz"
+    f = pyfits.open(fn)
+    d = f[1].data
+
+    nx = int((3000+3000)/5. + 1)
+    ny = int((3000+3000)/5. + 1)
+    nz = int((400 + 400)/5. + 1)
+    zmax = 400
+
+    X = d["X"].reshape((nz, ny, nx))
+    Y = d["Y"].reshape((nz, ny, nx))
+    Z = d["Z"].reshape((nz, ny, nx))
+
+    dA0_dD = d['dA0_dD'].reshape((nz, ny, nx))
+
+    return ((Z[:, 0, 0], Y[0, :, 0], X[0, 0, :]),
+            dA0_dD)
+
+
+def _save_converted_cube():
+    (z, y, x), dA0_dD = _load_data()
+    hdu = pyfits.PrimaryHDU(dA0_dD)
+    hdu_x = pyfits.ImageHDU(x, name="X")
+    hdu_y = pyfits.ImageHDU(y, name="y")
+    hdu_z = pyfits.ImageHDU(z, name="z")
+    outname = "lallerment2019_cube_converted.fits"
+    pyfits.HDUList([hdu, hdu_x, hdu_y, hdu_z]).writeto(outname, overwrite=True)
+
+
+def _load_converted_cube(converted_cube_name):
+    # outname = "lallerment2019_cube_converted.fits"
+    f = pyfits.open(converted_cube_name)
+    cube = f[0].data
+    x = f["X"].data
+    y = f["Y"].data
+    z = f["Z"].data
+
+    return (z, y, x), cube
+
+
+class ExtinctionCube(object):
+    def __init__(self, converted_cube_name=None):
+        if converted_cube_name is None:
+            converted_cube_name = "lallerment2019_cube_converted.fits"
+
+        (z, y, x), cube = _load_converted_cube(converted_cube_name)
+
+        self.zmax = np.max(z)
+        self.knots = RegularGridInterpolator((z, y, x), cube)
+
+    def get_a0_lb(self, l, b, unit="deg"):
+        # c = SkyCoord(30, 1, unit='deg', frame=Galactic)
+
+        c = SkyCoord(l, b, unit=unit, frame=Galactic)
+
+        xyz = c.cartesian.xyz.value
+
+        r_pc = np.arange(3000)
+        xyz_lb = xyz * r_pc[:, np.newaxis]
+        zmsk = np.abs(xyz_lb[:, -1]) < self.zmax
+
+        a0 = self.knots(xyz_lb[zmsk, ::-1])
+
+        return r_pc[zmsk], a0
+
+    def get_A0_lb(self, l, b, unit="deg"):
+        r_pc, a0 = self.get_a0_lb(l, b, unit=unit)
+        A0 = np.add.accumulate(a0)
+
+        return r_pc, A0
+
+
+if __name__ == "__main__":
+    l, b = 30, 1
+    ec = ExtinctionCube()
+    r_pc, A0 = ec.get_A0_lb(l, b)
+    print(r_pc, A0)
+
-- 
GitLab