Star_catalog_lv7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For stars dimmer\n",
    "\n",
    "```sql\n",
    "with x as\n",
    "(\n",
    "\tSELECT G.source_id, G.ra, G.dec, G.parallax, G.parallax_over_error, G.pmra, G.pmdec, G.phot_g_mean_mag, G.bp_rp, G.astrometric_params_solved, G.ruwe, G.phot_g_mean_mag - (0.01426 * POWER(G.bp_rp, 3) - 0.2156 * POWER(G.bp_rp, 2) + 0.01424 * POWER(G.bp_rp, 1) - 0.02704) as v_mag \n",
    "    FROM gaiadr3.gaia_source AS G\n",
    "\tWHERE G.phot_g_mean_mag IS NOT NULL\n",
    ")\n",
    "SELECT * \n",
    "FROM x\n",
    "WHERE (x.bp_rp IS NOT NULL AND x.v_mag > 15.5 AND x.v_mag <= 16.75) OR (x.bp_rp IS NULL AND x.phot_g_mean_mag > 15.5  AND x.phot_g_mean_mag <= 16.75)\n",
    "```\n",
    "\n",
    "correspond to file named `1733602856427O-result.fits`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "57534724\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 57534724/57534724 [36:05<00:00, 26565.60it/s]\n"
     ]
    }
   ],
   "source": [
    "import pathlib\n",
    "import struct\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from astropy.table import Table\n",
    "from py.geodesic import radec2xyz, GeodesicGrid\n",
    "import tqdm\n",
    "from py.gaia import gbprp_to_bv\n",
    "import trimesh\n",
    "from astropy.time import Time\n",
    "\n",
    "\n",
    "gaia_base_path = pathlib.Path(\"./gaia_query_results\")\n",
    "starcatalog_base_path = pathlib.Path(\"./star_catalogs\")\n",
    "\n",
    "# =============================================================================\n",
    "# Process the Gaia DR3 catalog\n",
    "# =============================================================================\n",
    "gaia_t = Table.read(gaia_base_path / \"1733602856427O-result.fits\", format=\"fits\")\n",
    "good_astrometry_idx = (gaia_t[\"astrometric_params_solved\"] == 31) | (\n",
    "    gaia_t[\"astrometric_params_solved\"] == 95\n",
    ") & (gaia_t[\"ruwe\"] < 1.4)\n",
    "df_gaia = gaia_t.to_pandas()\n",
    "b_v = np.zeros(len(df_gaia)) * np.nan\n",
    "\n",
    "# read synthetic photometry\n",
    "synth_phot_df = pd.read_csv(gaia_base_path / \"Gaia_XP_JKC.csv\")\n",
    "# try to get synthetic photometry\n",
    "matched_source_id, idx1, idx2 = np.intersect1d(\n",
    "    df_gaia[\"source_id\"].values, synth_phot_df[\"source_id\"].values, return_indices=True\n",
    ")\n",
    "# if we have synthetic photometry, use it\n",
    "b_v[idx1] = (synth_phot_df.loc[idx2, \"Jkc_mag_B\"].values - synth_phot_df.loc[idx2, \"Jkc_mag_V\"].values)\n",
    "# baseline B-V calculation, inferior to the synthetic photometry\n",
    "_, b_v_est = gbprp_to_bv(df_gaia[\"phot_g_mean_mag\"], df_gaia[\"bp_rp\"], red_correction=True)\n",
    "b_v = np.where(np.isnan(b_v), b_v_est, b_v)\n",
    "df_gaia.fillna({\"b_v\": b_v_est}, inplace=True)\n",
    "df_gaia[\"b_v\"] = b_v\n",
    "\n",
    "# make df_gaia\n",
    "df = pd.DataFrame(\n",
    "    index=range(len(df_gaia)),\n",
    "    data={\n",
    "        \"source_id\": df_gaia[\"source_id\"].values,\n",
    "        \"ra\": df_gaia[\"ra\"].values,\n",
    "        \"dec\": df_gaia[\"dec\"].values,\n",
    "        \"epoch\": np.zeros(len(df_gaia), dtype=float)\n",
    "        + 2016.0,  # we have propagated to J2000.0\n",
    "        \"parallax\": df_gaia[\"parallax\"].values,\n",
    "        \"parallax_error\": df_gaia[\"parallax\"].values\n",
    "        / df_gaia[\"parallax_over_error\"].values,\n",
    "        \"pmra\": df_gaia[\"pmra\"].values,\n",
    "        \"pmdec\": df_gaia[\"pmdec\"].values,\n",
    "        \"b_v\": df_gaia[\"b_v\"].values,\n",
    "        \"vmag\": df_gaia[\"v_mag\"].values,\n",
    "    },\n",
    ")\n",
    "\n",
    "# fill bad data with 0\n",
    "df = df.fillna(0)\n",
    "assert df[\"source_id\"].dtype == np.int64\n",
    "neg_parallax_idx = df[\"parallax\"] < 0\n",
    "df.loc[neg_parallax_idx, \"parallax\"] = 0\n",
    "df.loc[neg_parallax_idx, \"parallax_error\"] = 0\n",
    "\n",
    "min_vmag = [15.5]\n",
    "max_vmag = [16.75]\n",
    "major_version = [0]\n",
    "minor_version = [3]\n",
    "data_type = [1]\n",
    "# mag_ranges = [12800]\n",
    "# mag_steps = [256]\n",
    "levels = [7]\n",
    "suffixes = [None]\n",
    "catalog_epoch_jd_fromJ2000 = Time(2016, format=\"jyear\").jd\n",
    "# =============================================================================\n",
    "# Write to file\n",
    "# =============================================================================\n",
    "\n",
    "for lv, datatype, majver, minver, min_v, max_v, suff in zip(\n",
    "    levels,\n",
    "    data_type,\n",
    "    major_version,\n",
    "    minor_version,\n",
    "    min_vmag,\n",
    "    max_vmag,\n",
    "    suffixes\n",
    "):\n",
    "    if suff is not None:\n",
    "        suff = f\"_{suff}\"\n",
    "    else:\n",
    "        suff = \"\"\n",
    "    df6 = df[(df[\"vmag\"] > min_v) & (df[\"vmag\"] <= max_v)].reset_index(drop=True)\n",
    "    print(len(df6))\n",
    "\n",
    "    # =============================================================================\n",
    "    # Deal with geodesic grid\n",
    "    # =============================================================================\n",
    "    grid = GeodesicGrid(level=lv)\n",
    "\n",
    "    # setup ray-tracing\n",
    "    xyz = radec2xyz(df6[\"ra\"], df6[\"dec\"]).T\n",
    "    vertices = np.vstack(grid.vertices[-1])\n",
    "    faces = np.vstack(grid.faces[-1])\n",
    "    mesh = trimesh.Trimesh(vertices=vertices, faces=faces)\n",
    "    ray_origins = np.tile([[0, 0, 0]], (len(df6),1))\n",
    "    ray_directions = np.array(xyz)\n",
    "    intersected_faces, ray_indices = mesh.ray.intersects_id(\n",
    "        ray_origins=ray_origins,\n",
    "        ray_directions=ray_directions,\n",
    "        multiple_hits=False,\n",
    "    )\n",
    "    df6.loc[ray_indices, f\"zone{lv}\"] = intersected_faces\n",
    "    # find which number ray_indices are missing\n",
    "    missing = np.setdiff1d(np.arange(len(df6)), ray_indices)\n",
    "    # if there are missing values, we need to do a more expensive search with algebra\n",
    "    zones = grid.search_zone(xyz[missing])\n",
    "    df6.loc[missing, f\"zone{lv}\"] = zones\n",
    "\n",
    "    # sort by zone, within each zone sort by vmag\n",
    "    df6 = df6.sort_values([f\"zone{lv}\", \"vmag\"]).reset_index(drop=True)\n",
    "    df6[\"pmra_wo_cosdec\"] = df6[\"pmra\"] / np.cos(np.radians(df6[\"dec\"]))\n",
    "    \n",
    "    f = open(starcatalog_base_path / f\"./stars_{lv}_{datatype}v{majver}_{minver}{suff}.cat\", \"w+b\")\n",
    "    f.write(b\"\\n\\x04_\\x83\")  # Magic Number\n",
    "    f.write(np.int32(datatype).tobytes())  # Data Type\n",
    "    f.write(np.int32(majver).tobytes())  # Major Version\n",
    "    f.write(np.int32(minver).tobytes())  # Minor Version\n",
    "    f.write(np.int32(lv).tobytes())  # Level\n",
    "    f.write(np.int32(min_v * 1000).tobytes())  # Magnitude Minimum\n",
    "    f.write(np.float32(catalog_epoch_jd_fromJ2000).tobytes())  # Catalog Epoch\n",
    "    n_zones = 20 * 4**lv + 1  # plus 1 global zone\n",
    "\n",
    "    # count number of stars in each zone in df6\n",
    "    zone_info = df6[f\"zone{lv}\"].value_counts().sort_index()\n",
    "    for z in range(n_zones):\n",
    "        f.write(struct.pack(\"i\", zone_info.get(z, 0)))\n",
    "\n",
    "    max_records = sum(zone_info)\n",
    "\n",
    "    if datatype == 1:\n",
    "        # round to nearest integer\n",
    "        df6[\"ra_3600000\"] = (df6[\"ra\"] * 3_600_000).round().astype(int)\n",
    "        df6[\"dec_3600000\"] = (df6[\"dec\"] * 3_600_000).round().astype(int)\n",
    "        df6[\"pmra_wo_cosdec_1000\"] = (df6[\"pmra_wo_cosdec\"] * 1000).round().astype(int)\n",
    "        df6[\"pmdec_1000\"] = (df6[\"pmdec\"] * 1000).round().astype(int)\n",
    "        df6[\"b_v_1000\"] = (df6[\"b_v\"] * 1000).round().astype(int)\n",
    "        df6[\"vmag_1000\"] = (df6[\"vmag\"] * 1000).round().astype(int)\n",
    "        df6[\"parallax_100\"] = (df6[\"parallax\"] * 100).round().astype(int)\n",
    "        df6[\"parallax_error_100\"] = (df6[\"parallax_error\"] * 100).round().astype(int)\n",
    "        for i in tqdm.tqdm(range(max_records)):\n",
    "            # star_header = (\n",
    "            #     \"gaia_id\",\n",
    "            #     \"x0\",\n",
    "            #     \"x1\",\n",
    "            #     \"dx0\",\n",
    "            #     \"dx1\",\n",
    "            #     \"b_v\",\n",
    "            #     \"mag\",\n",
    "            #     \"plx\",\n",
    "            #     \"plx_err\"\n",
    "            # )\n",
    "            f.write(\n",
    "                struct.pack(\n",
    "                    \"qiiiihhHH\",\n",
    "                    df6.loc[i, \"source_id\"],\n",
    "                    df6.loc[i, \"ra_3600000\"],\n",
    "                    df6.loc[i, \"dec_3600000\"],\n",
    "                    df6.loc[i, \"pmra_wo_cosdec_1000\"],\n",
    "                    df6.loc[i, \"pmdec_1000\"],\n",
    "                    df6.loc[i, \"b_v_1000\"],\n",
    "                    df6.loc[i, \"vmag_1000\"],\n",
    "                    df6.loc[i, \"parallax_100\"],  # switch form 50 to 100 beacuse proxima centauri is not here\n",
    "                    df6.loc[i, \"parallax_error_100\"],\n",
    "                )\n",
    "            )\n",
    "    else:\n",
    "        raise ValueError(\"Data Type must be 0 or 1\")\n",
    "    f.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}