## MAXI All-Sky Map Generation from Event Data

This playbook describes how to generate all-sky X-ray maps by accumulating event-by-event photon data from MAXI.
The source dataset is [MAXI data for standard analysis](https://darts.isas.jaxa.jp/datasets/darts:maxi-stddata/).

For visualizing the pre-generated map products released by the MAXI team, see `maxi-allskymap-products`.

---

### Data Structure

Processed event files are organized by MJD under:

```
https://data.darts.isas.jaxa.jp/pub/maxi/mxdata/obs/MJD{outer}/MJD{mjd}/processed/{instrument}/
```

where `{outer}` is the MJD rounded down to the nearest 1000 (e.g., MJD 56001 → outer = MJD56000),
and `{instrument}` is `gsc_low` or `ssc_med`.

Each daily directory contains up to 768 event files (one per HEALPix tile).
If a tile has no events for that day, its file is absent.

---

### Energy Conversion

| Instrument | Conversion | Default band | RGB bands |
|---|---|---|---|
| GSC (`gsc_low`) | 1 keV = 20 PI | 2–20 keV | R: 2–4, G: 4–8, B: 8–16 keV |
| GSC (`gsc_med`) | 1 keV = 20 PI | 2–20 keV | R: 2–4, G: 4–8, B: 8–16 keV |
| SSC (`ssc_med`) | 1 PI = 0.00365 keV | 0.7–7 keV | R: 0.7–1, G: 1–2, B: 2–4 keV |

---

### Generating a Map with Python

#### Requirements

```
pip install astropy healpy numpy matplotlib requests
```

#### Helper functions

```python
import re
import numpy as np
import requests
import healpy as hp
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.time import Time
import astropy.units as u

BASE_URL = "https://data.darts.isas.jaxa.jp/pub/maxi/mxdata/obs"

GSC_PI_PER_KEV = 20.0
SSC_KEV_PER_PI = 0.00365

DEFAULT_BANDS = {
    "gsc_low": (2.0, 20.0),
    "ssc_med": (0.7, 7.0),
}
RGB_BANDS = {
    "gsc_low": [(2.0, 4.0), (4.0, 8.0), (8.0, 16.0)],
    "ssc_med": [(0.7, 1.0), (1.0, 2.0), (2.0, 4.0)],
}

def energy_to_pi(instrument, e_kev):
    if instrument == "gsc_low":
        return e_kev * GSC_PI_PER_KEV
    else:
        return e_kev / SSC_KEV_PER_PI

def dir_url(instrument, mjd):
    outer = (mjd // 1000) * 1000
    return f"{BASE_URL}/MJD{outer}/MJD{mjd}/processed/{instrument}/"

def list_evt_files(url):
    resp = requests.get(url, timeout=30)
    if resp.status_code != 200:
        return []
    return [url + name for name in re.findall(r'href="(mx_\w+\.evt\.gz)"', resp.text)]

def read_events(file_url, pi_min, pi_max, t_start=None, t_stop=None):
    with fits.open(file_url) as hdul:
        ev = hdul["EVENTS"].data
    mask = (ev["PI"] >= pi_min) & (ev["PI"] < pi_max)
    if t_start is not None:
        mask &= ev["TIME"] >= t_start
    if t_stop is not None:
        mask &= ev["TIME"] < t_stop
    return ev["RA"][mask], ev["DEC"][mask]

def events_to_hpmap(ra_deg, dec_deg, nside=256):
    coords = SkyCoord(ra=ra_deg * u.deg, dec=dec_deg * u.deg, frame="icrs").galactic
    ipix = hp.ang2pix(nside, np.pi / 2 - coords.b.rad, coords.l.rad)
    hpmap = np.zeros(hp.nside2npix(nside))
    np.add.at(hpmap, ipix, 1)
    return hpmap
```

#### Single-band map

```python
def make_allsky_map(instrument, mjd_start, mjd_stop,
                    e_min=None, e_max=None, nside=256,
                    t_start=None, t_stop=None):
    if e_min is None:
        e_min, e_max = DEFAULT_BANDS[instrument]
    pi_min = energy_to_pi(instrument, e_min)
    pi_max = energy_to_pi(instrument, e_max)

    ra_all, dec_all = [], []
    for mjd in range(mjd_start, mjd_stop + 1):
        for furl in list_evt_files(dir_url(instrument, mjd)):
            try:
                ra, dec = read_events(furl, pi_min, pi_max, t_start, t_stop)
                ra_all.append(ra)
                dec_all.append(dec)
            except Exception:
                pass

    if not ra_all:
        return np.zeros(hp.nside2npix(nside))
    return events_to_hpmap(np.concatenate(ra_all), np.concatenate(dec_all), nside)


# Example: GSC 2–20 keV, January 2012
mjd_start = int(Time("2012-01-01").mjd)
mjd_stop  = int(Time("2012-01-31").mjd)

hpmap = make_allsky_map("gsc_low", mjd_start, mjd_stop)

disp = np.where(hpmap > 0, hpmap, np.nan)
hp.mollview(disp, coord=["G"], cmap="afmhot", norm="log",
            title="MAXI/GSC All-Sky Map (2–20 keV, 2012-01)",
            unit="counts/pixel")
hp.graticule()
plt.savefig("maxi_gsc_allsky.png", dpi=150, bbox_inches="tight")
plt.close()
```

#### RGB three-color map

```python
def make_rgb_map(instrument, mjd_start, mjd_stop, nside=256):
    channels = []
    for e_min, e_max in RGB_BANDS[instrument]:
        hpmap = make_allsky_map(instrument, mjd_start, mjd_stop,
                                e_min=e_min, e_max=e_max, nside=nside)
        channels.append(hpmap)
    return channels

def plot_rgb(channels, title="MAXI All-Sky RGB Map"):
    def norm(arr):
        arr = np.where(arr > 0, arr, 0.0)
        vmax = np.percentile(arr[arr > 0], 99) if np.any(arr > 0) else 1.0
        return np.clip(arr / vmax, 0, 1)

    maps_2d = []
    for ch in channels:
        m2d = hp.mollview(norm(ch), coord=["G"], return_projected_map=True)
        plt.close()
        maps_2d.append(np.where(m2d.mask, 0.0, m2d.data))

    rgb = np.stack(maps_2d, axis=-1)
    fig, ax = plt.subplots(figsize=(12, 6))
    ax.imshow(rgb, origin="lower", aspect="auto")
    ax.set_title(title)
    plt.tight_layout()
    plt.savefig("maxi_allsky_rgb.png", dpi=150, bbox_inches="tight")
    plt.close()


# Example: SSC RGB, January 2012
channels = make_rgb_map("ssc_med", mjd_start, mjd_stop)
plot_rgb(channels, title="MAXI/SSC All-Sky RGB Map (R: 0.7–1, G: 1–2, B: 2–4 keV)")
```

---

### Acknowledgement

When using this data in publications, please include:

> "This research has made use of MAXI data provided by RIKEN, JAXA and the MAXI team."

and cite:

- Matsuoka, M. et al. (2009), PASJ, 61, 999. [doi:10.1093/pasj/61.5.999](https://doi.org/10.1093/pasj/61.5.999)
- Sugizaki, M. et al. (2011), PASJ, 63, S635. [doi:10.1093/pasj/63.sp3.S635](https://doi.org/10.1093/pasj/63.sp3.S635) — for GSC
- Tomida, H. et al. (2010), PASJ, 62, 1371. [doi:10.1093/pasj/62.6.1371](https://doi.org/10.1093/pasj/62.6.1371) — for SSC