Find dust region¶

In [2]:
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization import ImageNormalize, SqrtStretch
from astropy.stats import sigma_clipped_stats
from skimage.measure import find_contours, label, regionprops  # Added label, regionprops
In [3]:
def process_and_plot_contours(fits_filename, min_area_pixels=10, show_id=None):
    """
    Main function to load a FITS file, find a contour threshold,
    plot the image and numbered contours, and extract polygon data.

    Args:
        fits_filename (str): Path to the input FITS file.
        min_area_pixels (int): Minimum area in pixels for a region to be considered
                                and numbered. Helps filter out tiny spurious contours.
    """

    # --- 1. Load FITS File and WCS ---
    try:
        with fits.open(fits_filename) as hdul:
            data = hdul[0].data
            header = hdul[0].header
            wcs = WCS(header)

            if data.ndim > 2:
                print(f"Warning: Input FITS is {data.ndim}D. Taking first slice.")
                data = data.squeeze()
                while data.ndim > 2:
                    data = data[0]

            if data.ndim != 2:
                print(f"Error: Could not reduce data to 2D. Shape is {data.shape}")
                return

            data = np.nan_to_num(data, nan=0.0)

    except FileNotFoundError:
        print(f"Error: FITS file not found at {fits_filename}")
        return
    except Exception as e:
        print(f"Error loading FITS file: {e}")
        return

    # --- 2. Determine Contour Level ---
    try:
        mean, median, stddev = sigma_clipped_stats(data, sigma=3.0)
    except Exception as e:
        print(f"Warning: sigma_clipped_stats failed. Using simple stats. Error: {e}")
        mean, median, stddev = np.mean(data), np.median(data), np.std(data)

    n_sigma = 5.0
    contour_threshold = median + (n_sigma * stddev)

    print(f"Background Median: {median:.2f}")
    print(f"Background StdDev: {stddev:.2f}")
    print(f"Contour Threshold ({n_sigma}-sigma): {contour_threshold:.2f}")

    # --- 3. Identify and Label Regions ---
    # Create a binary mask where pixels above threshold are True
    binary_mask = data > contour_threshold

    # Label connected regions in the binary mask
    # This effectively assigns a unique integer to each distinct region
    labeled_regions = label(binary_mask, connectivity=2)  # connectivity=2 for 8-way (diagonal)

    # Get properties of each labeled region
    properties = regionprops(labeled_regions)

    # Filter out small regions and prepare for contour extraction
    significant_regions = []
    # Store centroids for labeling
    region_centroids_pixel = []

    for prop in properties:
        # Filter out regions smaller than min_area_pixels
        if prop.area >= min_area_pixels:
            significant_regions.append(prop)
            # Centroid is (row, col), which is (y, x)
            region_centroids_pixel.append(prop.centroid)

    print(f"Found {len(significant_regions)} significant region(s) (min_area_pixels={min_area_pixels}).")

    # --- 4. Plot Image and Numbered Contours ---
    print("\nGenerating plot...")

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(1, 1, 1, projection=wcs)

    norm = ImageNormalize(vmin=median + stddev, vmax=median + 50 * stddev, stretch=SqrtStretch())

    im = ax.imshow(data, origin="lower", cmap="cividis", norm=norm)

    # Draw contours for *each* significant labeled region
    all_polygons_sky_coords = []  # Will store the sky coordinates of each polygon

    if show_id is not None:
        # Only show the specified region ID if provided
        significant_regions = [prop for i, prop in enumerate(significant_regions) if i == show_id]

    for i, prop in enumerate(significant_regions):
        # Create a mask for just this region
        single_region_mask = labeled_regions == prop.label

        # Apply the region mask to the original data, set other parts to 0 or NaN
        # This isolates the data for contour finding
        data_for_contour = np.where(single_region_mask, data, 0)

        # Find contours for this isolated region.
        # find_contours returns one or more contours if a region is very complex
        contours_pixel_coords = find_contours(data_for_contour, contour_threshold)

        # If multiple contours define a single labeled region (e.g., holes),
        # we'll still treat them as belonging to this region number.
        for contour_pix in contours_pixel_coords:
            ax.plot(contour_pix[:, 1], contour_pix[:, 0], color="red", linewidth=1.5)  # Plotting (x,y)

            # Convert pixel coordinates to sky coordinates for storage
            pix_y, pix_x = contour_pix[:, 0], contour_pix[:, 1]
            sky_coords = wcs.pixel_to_world_values(pix_x, pix_y)
            all_polygons_sky_coords.append((sky_coords, i + 1))  # Store (RA, Dec) and its assigned number

        # Overlay the region number
        # Centroid is (row, col) which is (y, x)
        centroid_y, centroid_x = prop.centroid

        # Convert centroid pixel to world coordinates for text placement
        ra_text, dec_text = wcs.pixel_to_world_values(centroid_x, centroid_y)

        if show_id is None:
            ax.text(ra_text, dec_text, str(i), color="white", fontsize=14, fontweight="bold", ha="center", va="center", transform=ax.get_transform("fk5"))  # 'fk5' is for sky coordinates (RA, Dec)

    # Add labels and grid
    ax.set_xlabel("Right Ascension (deg)", fontsize=12)
    ax.set_ylabel("Declination (deg)", fontsize=12)
    ax.grid(color="white", ls=":", alpha=0.5)

    # Add a colorbar
    cbar = fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
    cbar.set_label("Pixel Value")

    ax.set_title(f"WISE Image {contour_threshold:.2f} (5$\sigma$)")

    # plt.savefig("numbered_dust_contours_plot.png")
    # print("Saved plot to 'numbered_dust_contours_plot.png'")
    # plt.show() # Uncomment this to show the plot interactively

    # --- 5. Output Polygon Information with Region Numbers ---
    print("\n--- Extracted Polygon Vertices with Region Numbers ---")

    if not all_polygons_sky_coords:
        print("No contours found at this threshold.")
        return

    # Sort polygons by their assigned region number for clarity in output
    all_polygons_sky_coords.sort(key=lambda x: x[1])

    return all_polygons_sky_coords
    # current_region_num = None
    # for (ra_vals, dec_vals), region_num in all_polygons_sky_coords:
    #     if region_num != current_region_num:
    #         if current_region_num is not None:
    #             print("--------------------") # Separator between regions
    #         print(f"\nRegion {region_num} (has {len(ra_vals)} vertices):")
    #         print("  Format: (RA, Dec) in degrees")
    #         current_region_num = region_num

    #     # Print the first 5 vertices as a sample for each part of the contour
    #     for j in range(min(5, len(ra_vals))):
    #         print(f"  Vertex {j}: ({ra_vals[j]:.6f}, {dec_vals[j]:.6f})")
    #     if len(ra_vals) > 5:
    #         print("  ...")

    # print("\n--- End of Polygon Data ---")

<>:135: SyntaxWarning: invalid escape sequence '\s'
<>:135: SyntaxWarning: invalid escape sequence '\s'
/var/folders/8f/nqw1kk19293c2gh_66w_y7kr00061z/T/ipykernel_95551/2684197215.py:135: SyntaxWarning: invalid escape sequence '\s'
  ax.set_title(f"WISE Image {contour_threshold:.2f} (5$\sigma$)")
In [4]:
all_polygons_sky_coords = process_and_plot_contours("ec0814a42/1590m606_ac51-w3-int-3_ra158.58609708333333_dec-60.637071666666664_asec600.000.fits")

Background Median: 396.50
Background StdDev: 23.49
Contour Threshold (5.0-sigma): 513.95
Found 7 significant region(s) (min_area_pixels=10).

Generating plot...

--- Extracted Polygon Vertices with Region Numbers ---
No description has been provided for this image
In [ ]:
process_and_plot_contours("ec0814a42/1590m606_ac51-w3-int-3_ra158.58609708333333_dec-60.637071666666664_asec600.000.fits", show_id=3)
In [5]:
dust_region = all_polygons_sky_coords[3]
In [ ]:
import numpy as np
import warnings
from astropy.table import Table, hstack, Column
from astropy.coordinates import SkyCoord
from astropy.coordinates.matching import match_coordinates_sky
from astropy import units as u
from astropy.utils.exceptions import AstropyWarning
from dl import queryClient as qc
from astroquery.gaia import Gaia
from astroquery.vizier import Vizier

# Suppress warnings for cleaner output
warnings.simplefilter("ignore", category=AstropyWarning)

# --- 1. Define Your Region of Interest (ROI) ---
# Coordinates provided by the user
RA_MAX = 158.65162834
RA_MIN = 158.53264874
DEC_MAX = -60.59046020
DEC_MIN = -60.68038515

# --- 2. Query Functions for Each Catalog ---


def query_decaps(ra_min, ra_max, dec_min, dec_max):
    """
    Queries the DECaPS DR2 catalog from the NOIRLab Astro Data Lab.
    """
    print(f"Querying DECaPS DR2 for box:\n  RA: {ra_min} to {ra_max}\n  Dec: {dec_min} to {dec_max}")

    # Construct the ADQL query
    adql_query = f"""
    SELECT *
    FROM decaps_dr2.object
    WHERE ra > {ra_min} AND ra < {ra_max}
      AND dec > {dec_min} AND dec < {dec_max}
    """

    try:
        table = qc.query(sql=adql_query, fmt="table")
        print(f"  Found {len(table)} sources in DECaPS.")
        return table
    except Exception as e:
        print(f"  Error querying DECaPS: {e}")
        return None


def query_gaia(ra_min, ra_max, dec_min, dec_max):
    """
    Queries the Gaia DR3 catalog using a box search.
    """
    print(f"Querying Gaia DR3 for the same box...")

    # Construct the ADQL query for Gaia
    adql_query = f"""
    SELECT source_id, ra, dec,
           phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag,
           parallax, parallax_error, parallax_over_error, pmra, pmdec
    FROM gaiadr3.gaia_source
    WHERE ra > {ra_min} AND ra < {ra_max}
      AND dec > {dec_min} AND dec < {dec_max}
    """

    try:
        job = Gaia.launch_job_async(adql_query)
        table = job.get_results()
        print(f"  Found {len(table)} sources in Gaia.")
        return table
    except Exception as e:
        print(f"  Error querying Gaia: {e}")
        return None


def query_2mass(ra_min, ra_max, dec_min, dec_max):
    """
    Queries the 2MASS Point Source Catalog (PSC) from VizieR.
    """
    print(f"Querying 2MASS PSC for the same box...")

    # Calculate center and dimensions for VizieR query
    ra_center = (ra_min + ra_max) / 2
    dec_center = (dec_min + dec_max) / 2
    width = (ra_max - ra_min) * u.deg
    height = (dec_max - dec_min) * u.deg

    center_coord = SkyCoord(ra_center, dec_center, unit="deg", frame="icrs")

    try:
        vizier = Vizier(columns=["_RAJ2000", "_DEJ2000", "2MASS", "Jmag", "Hmag", "Kmag", "e_Jmag", "e_Hmag", "e_Kmag"], row_limit=-1)

        # VizieR catalog 'II/246/out' is the 2MASS PSC
        result_tables = vizier.query_region(center_coord, width=width, height=height, catalog="II/246/out")

        if not result_tables:
            print("  Found 0 sources in 2MASS.")
            return None

        table = result_tables[0]
        # Rename coordinate columns for consistency
        table.rename_column("_RAJ2000", "ra")
        table.rename_column("_DEJ2000", "dec")

        print(f"  Found {len(table)} sources in 2MASS.")
        return table
    except Exception as e:
        print(f"  Error querying 2MASS: {e}")
        return None


# --- 3. Main Cross-Matching Function ---


def create_master_catalog(base_table, match_catalogs, max_separation):
    """
    Cross-matches a base table against a list of other catalogs.

    :param base_table: The primary Astropy Table (DECaPS).
    :param match_catalogs: A dictionary of {'name': table} for catalogs to match (Gaia, 2Mass).
    :param max_separation: The Astropy Quantity for the max match radius (e.g., 1.0 * u.arcsec).
    """

    if base_table is None or len(base_table) == 0:
        print("Base catalog is empty. Exiting.")
        return None

    print(f"\n--- Starting Cross-Match ---")
    print(f"Base catalog: DECaPS ({len(base_table)} sources)")
    print(f"Matching radius: {max_separation}")

    # Create SkyCoord object for the base DECaPS catalog
    base_coords = SkyCoord(ra=base_table["ra"] * u.deg, dec=base_table["dec"] * u.deg)

    master_table = base_table.copy()

    for name, match_table in match_catalogs.items():
        if match_table is None or len(match_table) == 0:
            print(f"Skipping {name} (catalog is empty).")
            continue

        # Ensure ra and dec are float arrays
        match_table["ra"] = np.asarray(match_table["ra"], dtype=float)
        match_table["dec"] = np.asarray(match_table["dec"], dtype=float)

        print(f"Matching against {name} ({len(match_table)} sources)...")

        # 1. Create SkyCoord object for the catalog to be matched
        match_coords = SkyCoord(ra=match_table["ra"] * u.deg, dec=match_table["dec"] * u.deg)

        # 2. Perform the cross-match against the base catalog
        #    This finds the closest 'match_table' source for *every* 'base_table' source
        idx, sep, _ = match_coordinates_sky(base_coords, match_coords)

        # 3. Create the table of matched sources, ordered like the base table
        #    This table will have the same length as 'base_table'
        join_table = match_table[idx]

        # 4. Add the separation distance as a new column
        join_table[f"{name}_sep"] = sep

        # 5. Mask all matches that are outside our radius
        #    We set their values to NaN (or -1/False for non-float columns)
        bad_match_mask = sep > max_separation

        for col in join_table.colnames:
            if join_table[col].dtype.kind == "f":
                # For float columns, set to NaN
                join_table[col][bad_match_mask] = np.nan
            else:
                # For other columns (int, str), set to a marker
                # (e.g., -1 for int, '' for str)
                if join_table[col].dtype.kind in ("i", "u"):
                    join_table[col][bad_match_mask] = -1
                elif join_table[col].dtype.kind in ("S", "U"):
                    join_table[col][bad_match_mask] = ""

        # 6. Rename columns to avoid conflicts, e.g., 'ra' -> 'gaia_ra'
        for col in join_table.colnames:
            if col not in ["ra", "dec"] or f"{name}_{col}" not in master_table.colnames:
                join_table.rename_column(col, f"{name}_{col}")

        # 7. Horizontally stack the new matched data onto the master table
        master_table = hstack([master_table, join_table], join_type="exact")

        num_matches = np.sum(~bad_match_mask)
        print(f"  ...found {num_matches} matches within {max_separation}.")

    return master_table
In [8]:
# # Define the maximum separation for a "good" match
# MATCH_RADIUS = 1.0 * u.arcsec

# # Step 1: Query all three catalogs
# print("Querying catalogs with updated Gaia query...")
# decaps_cat = query_decaps(RA_MIN, RA_MAX, DEC_MIN, DEC_MAX)
# gaia_cat = query_gaia(RA_MIN, RA_MAX, DEC_MIN, DEC_MAX)
# print(f"Gaia catalog columns: {gaia_cat.colnames if gaia_cat else 'None'}")
# tmass_cat = query_2mass(RA_MIN, RA_MAX, DEC_MIN, DEC_MAX)

# # Catalogs to match against the base
# catalogs_to_match = {"gaia": gaia_cat, "tmass": tmass_cat}

# # Step 2: Run the cross-match
# final_catalog = create_master_catalog(base_table=decaps_cat, match_catalogs=catalogs_to_match, max_separation=MATCH_RADIUS)

# # # Step 3: Save the results
# # if final_catalog:
# #     print("\n--- Process Complete ---")

# #     # Save to FITS file (best for astronomical data)
# #     fits_file = "master_catalog.fits"
# #     try:
# #         final_catalog.write(fits_file, overwrite=True)
# #         print(f"Successfully saved master catalog to: {fits_file}")
# #     except Exception as e:
# #         print(f"Error saving FITS file: {e}")

# #     # Save to CSV file (easy to open in Excel, etc.)
# csv_file = "master_catalog.csv"
# try:
#     final_catalog.write(csv_file, format="csv", overwrite=True)
#     print(f"Successfully saved master catalog to: {csv_file}")
# except Exception as e:
#     print(f"Error saving CSV file: {e}")

# #     print("\nCatalog summary:")
# #     final_catalog.info()
In [6]:
# if master_catalog.csv exists, load it:
import pandas as pd

final_catalog = pd.read_csv("master_catalog.csv")
# and make it into asttropy table
from astropy.table import Table

final_catalog = Table.from_pandas(final_catalog)
In [7]:
final_catalog
Out[7]:
Table length=9684
radecra_okdec_okepochrangeepochrange_okepochmeanepochmean_okelonelatglonglatmean_gmean_rmean_imean_zmean_ystdev_gstdev_rstdev_istdev_zstdev_yerr_gerr_rerr_ierr_zerr_ymedian_gmedian_rmedian_imedian_zmedian_yq25_gq25_rq25_iq25_zq25_yq75_gq75_rq75_iq75_zq75_ymean_cflux_gmean_cflux_rmean_cflux_imean_cflux_zmean_cflux_ystdev_cflux_gstdev_cflux_rstdev_cflux_istdev_cflux_zstdev_cflux_yerr_cflux_gerr_cflux_rerr_cflux_ierr_cflux_zerr_cflux_ymedian_cflux_gmedian_cflux_rmedian_cflux_imedian_cflux_zmedian_cflux_yq25_cflux_gq25_cflux_rq25_cflux_iq25_cflux_zq25_cflux_yq75_cflux_gq75_cflux_rq75_cflux_iq75_cflux_zq75_cflux_yobj_iddqflags_all_gdqflags_all_rdqflags_all_idqflags_all_zdqflags_all_ydqflags_any_gdqflags_any_rdqflags_any_idqflags_any_zdqflags_any_yposstdevposstdev_okrandom_idmean_cmag_gmean_cmag_rmean_cmag_imean_cmag_zmean_cmag_ymedian_cmag_gmedian_cmag_rmedian_cmag_imedian_cmag_zmedian_cmag_ymean_cmag_g_rmean_cmag_r_imean_cmag_i_zmean_cmag_z_ymedian_cmag_g_rmedian_cmag_r_imedian_cmag_i_zmedian_cmag_z_ymean_mag_gmean_mag_rmean_mag_imean_mag_zmean_mag_ymedian_mag_gmedian_mag_rmedian_mag_imedian_mag_zmedian_mag_ymean_mag_g_rmean_mag_r_imean_mag_i_zmean_mag_z_ymedian_mag_g_rmedian_mag_r_imedian_mag_i_zmedian_mag_z_ymaglimit_gmaglimit_rmaglimit_imaglimit_zmaglimit_yfracflux_avg_gfracflux_avg_rfracflux_avg_ifracflux_avg_zfracflux_avg_yrchi2_avg_grchi2_avg_rrchi2_avg_irchi2_avg_zrchi2_avg_yqf_avg_gqf_avg_rqf_avg_iqf_avg_zqf_avg_yprn_avg_gprn_avg_rprn_avg_iprn_avg_zprn_avg_yprl_avg_gprl_avg_rprl_avg_iprl_avg_zprl_avg_yprr_avg_gprr_avg_rprr_avg_iprr_avg_zprr_avg_ypre_avg_gpre_avg_rpre_avg_ipre_avg_zpre_avg_ycchi2_avg_gcchi2_avg_rcchi2_avg_icchi2_avg_zcchi2_avg_ysnr_avg_gsnr_avg_rsnr_avg_isnr_avg_zsnr_avg_yhtm9ring256nest4096dnt_all_gdnt_all_rdnt_all_idnt_all_zdnt_all_ydnt_any_gdnt_any_rdnt_any_idnt_any_zdnt_any_yndetndet_okndet_cfluxndet_cflux_oknmag_gnmag_rnmag_inmag_znmag_ynmag_ok_gnmag_ok_rnmag_ok_inmag_ok_znmag_ok_ynmag_cflux_gnmag_cflux_rnmag_cflux_inmag_cflux_znmag_cflux_ynmag_cflux_ok_gnmag_cflux_ok_rnmag_cflux_ok_inmag_cflux_ok_znmag_cflux_ok_ygaia_source_idgaia_ragaia_decgaia_phot_g_mean_maggaia_phot_bp_mean_maggaia_phot_rp_mean_maggaia_parallaxgaia_parallax_errorgaia_parallax_over_errorgaia_pmragaia_pmdecgaia_gaia_septmass_ratmass_dectmass_2MASStmass_Jmagtmass_Hmagtmass_Kmagtmass_e_Jmagtmass_e_Hmagtmass_e_Kmagtmass_tmass_sep
float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64int64int64int64int64int64int64int64int64int64int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64int64float64float64float64float64float64float64float64float64float64float64float64float64float64str16float64float64float64float64float64float64float64
158.6515225069059-60.60733371681638158.6515225069059-60.60733371681638312.1026092600005312.102609260000557532.2196449857532.21964498201.7801570258168-60.51104452433997286.9042264101394-2.0661884671849963.1411e-101.28847e-092.44067e-093.88302e-095.21371e-095.49842e-121.17806e-102.12359e-103.08858e-103.40031e-103.72604e-116.9101e-119.9616e-111.54958e-104.74806e-103.14688e-101.28887e-092.44581e-094.00786e-095.30976e-093.12744e-101.22463e-092.32827e-093.73079e-095.18954e-093.16632e-101.33822e-092.60571e-094.02169e-095.42998e-092.77006e-101.29849e-092.46526e-093.84592e-095.22615e-092.23128e-117.20284e-112.22072e-103.47304e-109.45397e-114.44222e-118.52636e-111.19158e-101.80046e-105.81132e-102.73238e-101.30402e-092.49174e-093.97371e-095.26065e-092.65349e-101.2605e-092.36335e-093.66599e-095.22722e-092.81126e-101.32883e-092.63146e-094.00417e-095.29407e-0968240089666528094300000000000.1175990.11759912.428923.893822.216421.520321.037520.704523.908622.211821.508721.00220.69741.677380.6960580.4828450.3329521.696860.7030450.5067350.30460723.757322.224821.531221.027120.707123.755322.224521.528920.992720.68731.532480.6935840.5041520.3199421.530820.6955340.5362240.30540524.025323.17122.871922.465821.37110.8488070.9321990.9473090.9514960.9404230.7710640.964060.8983460.9963830.9521751.01.01.01.01.00.01314080.1125494.12164e-050.001205540.001145460.6136170.7517420.008327530.02744830.02720670.3678150.03220020.9912910.8767530.9198620.005427940.1035090.0003398510.09459360.0517862690.152709.2635.688669.194724.1486.0144710.967313.57815.8318.19689245447673583215292098700000000001414141423432234322343223432-1------------------------------------------
158.6499951521797-60.60737482165183158.6499951521797-60.60737482165183312.1026092600005312.102609260000557518.2669357255557518.26693572555201.7792171708942-60.51163576844945286.9035976944293-2.0665993689815691.75259e-094.78669e-097.5938e-091.01363e-081.29706e-083.47592e-113.3007e-104.98445e-103.82814e-102.03546e-103.62297e-118.10087e-111.08895e-101.63185e-104.06732e-101.76507e-094.69126e-097.35462e-091.00804e-081.28925e-081.73593e-094.5923e-097.08223e-099.76708e-091.28403e-081.77013e-094.90501e-097.88282e-091.03534e-081.30433e-081.78521e-094.91207e-097.78873e-091.03683e-081.26477e-083.62914e-113.33533e-103.75421e-103.97277e-104.94062e-104.32165e-119.72695e-111.30493e-101.91544e-104.91321e-101.79046e-094.76735e-097.50819e-091.01847e-081.24895e-081.76364e-094.70355e-097.41113e-099.89627e-091.24049e-081.80271e-095.0045e-097.9455e-091.05045e-081.283e-0868240089666528088400000000000.06816310.068163154.209621.870820.771820.271319.960719.74521.867620.804320.311219.980119.75861.098940.5005050.3106020.2157571.063280.4931470.3310320.22148721.890820.799920.298919.985319.717621.883120.821820.333619.991319.72421.090880.501060.3135550.2676961.061320.4881760.3422970.26714523.916623.022722.744922.37121.32710.9731410.9869360.9907110.9830450.9905450.7662250.8961180.9268780.9280071.160331.01.01.01.01.00.01122490.1131435.04778e-050.0009919540.0009862690.5038650.7603610.01033620.02379790.02467440.4802020.03058150.9893140.8986130.93150.004707780.09591430.0002988830.07659690.0428391697.547707.617670.428675.288750.96229.300636.209436.668935.888819.36992454476735832152920987000000000018181818335433354333543335435254325108607569792158.64992595999027-60.6073590878843620.86478233337402321.75576782226562520.41754150390625----------3.7426755462288305e-05--------------------
158.6505192004757-60.60930280813872158.6505192004757-60.60930280813872312.1026092600005312.102609260000557518.2669357255557518.26693572555201.782509126082-60.51271163729126286.9047864262574-2.0681395100810373.43866e-098.92512e-091.41286e-081.99071e-082.24701e-085.15197e-112.33286e-103.58024e-104.52541e-101.28234e-094.19608e-119.7153e-111.27561e-101.9363e-104.32796e-103.45432e-098.78945e-091.40792e-081.97894e-082.20406e-083.42213e-098.78087e-091.4021e-081.95368e-082.16456e-083.47203e-098.98278e-091.44098e-082.00252e-082.29503e-083.5298e-099.30354e-091.43959e-082.01201e-082.27818e-086.55665e-112.99537e-102.5887e-104.52084e-101.18488e-094.8692e-111.14947e-101.49776e-102.22697e-105.22912e-103.54699e-099.22126e-091.42637e-081.99761e-082.23664e-083.50978e-099.12852e-091.42463e-081.97738e-082.19525e-083.57443e-099.41124e-091.45865e-082.02061e-082.31809e-0868240089666528091200000000000.06261760.062617633.06221.130620.078419.604419.240919.10621.125320.08819.614419.248719.1261.052250.4739780.3634780.1348931.037320.4736080.3656940.12270921.15920.123519.624719.252519.12121.154120.140119.628619.258919.14191.035560.4987160.3722690.1314951.0140.5115380.369640.11697623.828722.911522.634322.256721.2780.9975890.9989640.9987050.9975020.999640.9917010.9750191.254781.265370.8408751.01.01.00.9999721.00.01287840.1267385.32941e-050.001050670.001095710.5513970.7498650.01040250.02524830.02697750.4307440.0272080.9892370.898440.9241570.004980110.09618860.0003073450.07526080.0477698671.268755.899702.099674.342704.91253.092161.088461.039663.162432.22822454476735832152920987000000000018181818335433354333543335435254325108604414464158.6504569176072-60.60928815311246620.15669441223144521.06040573120117219.731416702270508-0.10186992353102070.410291314125061-0.2482868134975433-5.4940057107609062.7606920706049533.3897723685164526e-05--------------------
158.6498166520968-60.61015293554672158.6498166520968-60.61015293554672312.1026092600005312.102609260000557569.2368637533357569.23686375333201.783347882926-60.51353086816209286.9049137125375-2.0690480302911260.06.7188e-101.62394e-092.59276e-093.84493e-090.01.49627e-101.33262e-103.88699e-100.00.08.14105e-111.04727e-101.51982e-105.65986e-100.06.89657e-101.65065e-092.57086e-093.84493e-090.06.36756e-101.56025e-092.32378e-093.84493e-090.07.42559e-101.71047e-092.69408e-093.84493e-090.09.31432e-101.63336e-092.68865e-094.26002e-090.02.99824e-101.49761e-104.02892e-100.00.09.97843e-111.26454e-101.78299e-106.8278e-100.09.81365e-101.6707e-092.81133e-094.26002e-090.08.75361e-101.61407e-092.46293e-094.26002e-090.01.08737e-091.76201e-092.83856e-094.26002e-0968240089666530027100000000000.08446360.084463680.2379inf22.577121.967321.426220.9265inf22.520421.942821.377720.9265inf0.6098290.5411280.499693inf0.5776730.5650290.45125inf22.931821.973621.465621.0378inf22.903421.955921.474821.0378inf0.9581930.5079820.427816inf0.9475520.4810640.4370230.023.185322.897922.485821.37270.00.9650980.9619830.9764210.9875620.01.25331.030820.8772341.050110.01.01.00.9999641.00.00.1530010.0001046490.001565290.001734040.00.6708980.02133180.03474380.03690220.00.0273510.9780770.8489050.8772220.00.148750.0004863970.1147860.08414130.0657.798579.374596.643704.4820.05.830829.8297210.90726.8066724544767358321529209870000000000999902331023310233102331-1------------------------------------------
158.6482170041273-60.61001351762383158.6482170041273-60.61001351762383300.1482389900048300.148238990004857671.2891867866657671.28918678666201.7820846251726-60.51402994644718286.9041638831351-2.0693204179488380.00.09.03323e-101.4336e-090.00.00.00.03.34553e-100.00.00.01.37814e-101.69314e-100.00.00.09.03323e-101.53422e-090.00.00.09.03323e-101.41594e-090.00.00.09.03323e-101.6525e-090.00.00.07.86701e-101.48335e-090.00.00.00.01.35405e-100.00.00.01.59767e-101.97085e-100.00.00.07.86701e-101.52664e-090.00.00.07.86701e-101.47876e-090.00.00.07.86701e-101.57451e-090.068240089667448175100000000000.07930680.079306857.5606infinf22.760522.0719infinfinf22.760522.0407inf--inf0.688583inf--inf0.719814infinfinf22.610422.1089infinfinf22.610422.0353inf--inf0.501459inf--inf0.575111inf0.00.022.906322.49890.00.00.00.6765520.6876340.00.00.00.5469371.569810.00.00.01.01.00.00.00.00.0002413610.001917720.00.00.00.05053870.04184050.00.00.00.9483920.8113170.00.00.00.0008281450.1449250.00.00.0622.239585.0320.00.00.06.566656.232180.024544767358321529209870000000000333300120001200012000120-1------------------------------------------
158.6478029167787-60.61057735398612158.6478029167787-60.61057735398612312.1026092600005312.102609260000557518.2669357255557518.26693572555201.7826750212664-60.51455412214992286.9042703753516-2.0699102321928068.37026e-092.02229e-082.96478e-084.06935e-084.55101e-088.34616e-115.15047e-108.35063e-107.07437e-101.97164e-096.31054e-111.49133e-101.79023e-102.69816e-105.12373e-108.36843e-092.0453e-082.95348e-084.08321e-084.49239e-088.32634e-092.00415e-082.91937e-084.04622e-084.42889e-088.41149e-092.05149e-082.97418e-084.09713e-084.62916e-088.47636e-092.03493e-082.98569e-084.08895e-084.57612e-081.17052e-103.1639e-107.62221e-104.85985e-102.01912e-096.94148e-111.64758e-102.00766e-102.98224e-106.05329e-108.44321e-092.04855e-082.97165e-084.07667e-084.50951e-088.41443e-092.02259e-082.93188e-084.04715e-084.43736e-088.52472e-092.05204e-082.98416e-084.11007e-084.64565e-0868240089666528080900000000000.05341470.053414767.59420.179519.228618.812418.47118.348820.183719.221418.817518.474218.36470.9508480.4162410.3414170.1222130.9623450.4038810.3432670.10956220.193219.235418.8218.476218.354720.193419.223118.824218.472518.36880.9577620.4153750.343830.1214560.970280.3989450.3516690.10368923.617722.680222.429622.081221.17770.9870550.9916470.9803530.9850290.9881491.347471.146231.194551.192411.000231.01.01.01.01.00.01617270.1426817.93108e-050.0009509590.001113520.5809710.740890.01569230.02373060.0274320.3974660.0238620.9839380.9083730.9230.005390210.0925670.0002900280.06694560.0484549642.806671.092546.452586.692656.449108.648113.32107.949110.73759.72062454476735832152920987000000000018181818335433354333543335435254325108604384256158.64774431079402-60.61056128333094519.36314964294433620.0190105438232418.496717453002930.34370875733333180.20683068037033081.661787986755371-4.5167229351130752.0347398340157243.294588123826885e-05--------------------
158.6507048920572-60.61160769682976158.6507048920572-60.61160769682976300.1489100399995300.148910039999557591.2528846857591.25288468201.786155523651-60.5141605711275286.9060201426019-2.0700890598060070.00.01.15025e-092.98956e-093.41058e-090.00.00.02.09367e-100.00.00.01.37273e-101.52419e-105.63137e-100.00.01.15025e-093.0011e-093.41058e-090.00.01.15025e-092.93674e-093.41058e-090.00.01.15025e-093.14167e-093.41058e-090.00.01.16096e-092.98647e-093.07332e-090.00.00.03.38962e-100.00.00.01.59399e-101.80329e-106.7589e-100.00.01.16096e-093.07302e-093.07332e-090.00.01.16096e-092.93213e-093.07332e-090.00.01.16096e-093.25698e-093.07332e-0968240089667400200200000000000.1084530.10845397.5394infinf22.33821.312121.281infinf22.33821.281121.281--inf1.025850.0311279--inf1.056870.000106812infinf22.34821.31121.1679infinf22.34821.306821.1679--inf1.037030.143051--inf1.041220.1388680.00.022.911922.485521.37770.00.00.9703030.97010.9810480.00.00.789670.7905760.8516910.00.01.01.01.00.00.00.0001714590.001802690.001856140.00.00.03401720.03952830.03896460.00.00.9649450.8326020.8688130.00.00.0008661310.1260670.09036660.00.0593.959642.475766.3550.00.08.4043912.19826.0658324544767358321529209870000000000555500131001310013100131-1------------------------------------------
158.6508668366339-60.61242467592336158.6508668366339-60.61242467592336312.1026092600005312.102609260000557518.2669357255557518.26693572555201.7875111752713-60.5146386262446286.9064982959572-2.0707564712535342.01098e-096.74676e-091.08327e-081.53159e-081.67954e-081.43288e-107.79987e-114.84431e-102.57599e-107.3923e-103.72608e-118.91575e-111.18559e-101.78697e-104.18214e-102.04153e-096.70545e-091.06772e-081.53414e-081.69098e-081.93551e-096.7054e-091.05408e-081.51203e-081.65803e-082.07676e-096.77487e-091.11547e-081.55403e-081.73491e-082.08426e-096.95387e-091.10278e-081.54239e-081.69814e-081.24831e-104.81279e-113.16104e-102.9511e-101.06809e-094.47997e-111.07613e-101.39143e-102.08612e-105.07113e-102.13863e-096.96875e-091.08117e-081.55481e-081.72334e-082.02776e-096.93297e-091.07254e-081.54115e-081.67295e-082.14437e-096.97813e-091.11897e-081.57028e-081.78615e-0868240089666528093000000000000.06371710.063717181.871321.702620.394419.893819.529519.425121.674720.392119.915319.520819.40911.308190.500650.3642670.104451.282550.4768510.3944570.11173121.741520.427319.913219.537119.43721.725120.433919.928919.535319.42971.314220.5140990.3760190.1001171.291180.5050770.393520.10567923.891222.961622.683422.314621.3050.8635710.929760.8855220.9151350.8854221.014491.072620.8590820.9170430.8953321.01.01.01.01.00.01696180.15467.33865e-050.001482340.001286330.6424920.7176910.01396420.03404440.03088920.3348680.02321480.9855790.8652670.9107860.00567870.1044950.0003835510.09920660.0570386648.513693.038634.155675.039734.39633.048247.9649.185650.970424.46772454476735832152920987000000000018181818335433354333543335435254325108606782464158.65080357511476-60.6124113924598620.529727935791016----0.32934098080155430.69968628883361820.4706980586051941-7.1677633562488322.5868056666008553.376597844861796e-05--------------------
158.6493728405922-60.61107942359419158.6493728405922-60.61107942359419312.1026092600005312.102609260000557528.1521300806657528.15213008066201.7844733216093-60.51430478676988286.9051892562113-2.0699590683623893.87146e-101.3785e-093.26324e-094.97198e-095.58086e-096.31667e-112.15396e-101.51699e-104.79713e-101.08106e-093.79419e-117.18521e-111.01643e-101.49857e-104.38139e-103.80165e-101.23524e-093.27761e-095.05448e-095.71588e-093.57832e-101.23371e-093.17495e-094.6935e-095.33367e-094.02497e-101.41226e-093.3286e-095.20725e-096.0981e-094.3804e-101.62729e-093.3867e-095.1829e-095.83412e-093.07987e-118.95437e-112.04724e-104.20313e-101.11445e-094.54222e-118.82418e-111.20904e-101.77365e-105.31864e-104.33499e-101.66867e-093.41133e-095.17913e-095.99049e-094.2261e-101.59155e-093.30341e-094.88207e-095.59647e-094.44388e-101.67468e-093.46536e-095.3823e-096.38451e-0968240089666528086700000000000.0916440.09164475.138123.396221.971321.175620.713620.585123.407521.944121.167720.714420.55631.424880.7957780.461990.1285061.463460.7763820.4533350.15801223.530322.151521.215920.758720.633223.550122.270621.211120.740820.60731.378830.9356020.4572030.1254291.279451.059510.4702990.13351624.010723.133122.849422.447721.3640.9314140.960160.9664810.9591910.9676080.9859480.9285340.8722710.9827690.926651.01.01.01.01.00.01874610.1452710.0001064530.001500340.001438440.7129260.7203710.02151150.03418120.0329910.262080.02515810.9779320.8577320.8991370.006248080.10920.0004501860.1065870.0664335679.604696.531631.172636.658696.4247.3631211.360618.042218.94279.03778245447673583215292098700000000001515151523442234422344223442-1------------------------------------------
.....................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................
158.5346749275664-60.67098243519945158.5346735328992-60.67097804321804410.8663644499975410.866364449997557565.1538337958457555.79278174637201.8010455348976-60.59606216477372286.8865623782232-2.149984979661043.7334e-101.71553e-092.73156e-094.24477e-095.40952e-093.41869e-111.82197e-101.79254e-102.16263e-105.87343e-103.1097e-116.3561e-118.12717e-111.37485e-104.4295e-103.62099e-101.82544e-092.69661e-094.17171e-095.35982e-093.52771e-101.77631e-092.61142e-094.03745e-095.15217e-093.83965e-101.87598e-092.81355e-094.33322e-095.56748e-093.8326e-101.77364e-092.80255e-094.37354e-095.23018e-092.08265e-112.31331e-102.06164e-103.15438e-106.51771e-103.61867e-117.55653e-111.01464e-101.64637e-105.3579e-103.91928e-101.8419e-092.87796e-094.36747e-095.14626e-093.757e-101.81167e-092.76689e-094.22659e-094.91583e-093.94155e-102.02973e-092.98865e-094.54638e-095.3767e-0968241848885295679000000000000.1209080.11165270.215723.541321.877821.381120.897920.703723.51721.836821.352320.899420.72131.663430.4967210.4831960.194211.680150.4845470.4528620.17815623.569721.91421.40920.930420.667123.602921.846621.42320.949220.67711.655740.505030.4786090.2632621.756350.4236280.4737420.27208924.031623.16422.886222.466421.32620.9789660.991410.9897570.9933690.9945170.8921541.106490.9205681.049970.8821591.01.01.01.01.00.02583420.1051420.0006304059.64121e-056.66131e-060.4874560.7624760.09921150.01180180.0008161530.48020.1135840.8996230.9879010.9990910.006510180.01879740.0005351640.0002004518.59777e-05730.348723.46696.496713.781689.9487.0673712.762414.598815.7738.75054245447673646615292082200000000002422242245762357524576235752-1------------------------------------------
158.5365821415463-60.67261423469984158.5365821415463-60.67261423469984312.1032943199971312.103294319997157588.9458011779957588.94580117799201.8048014732483-60.5964296548374286.8881908348749-2.1509269544510860.01.78076e-093.21031e-094.1868e-090.00.02.61365e-101.93772e-101.30256e-100.00.06.96635e-111.07137e-101.57815e-100.00.01.82814e-093.30414e-094.21887e-090.00.01.62195e-093.18924e-094.11226e-090.00.02.05168e-093.38069e-094.23509e-090.00.01.81771e-093.27266e-094.28727e-090.00.02.92134e-102.52123e-107.99466e-110.00.08.15796e-111.3028e-101.84038e-100.00.01.85207e-093.42896e-094.26987e-090.00.01.65955e-093.26889e-094.2443e-090.00.02.07219e-093.49741e-094.33178e-090.068241848885397189100000000000.1001550.10015529.8734inf21.851221.212720.9195infinf21.830921.162120.924infinf0.6384370.293201infinf0.668760.238131infinf21.873521.233620.9453infinf21.84521.202420.937infinf0.6398510.28834infinf0.6426240.265345inf0.023.121322.84422.44710.00.00.5472020.5603430.5817130.00.00.6201050.892620.8375330.00.01.01.01.00.00.00.1051820.0007406725.47874e-050.00.00.7594180.1230430.006933860.00.00.1136880.8755320.9928560.00.00.02171170.0006844420.0001553560.00.0683.311674.652616.3960.00.013.275618.637716.7160.0245447673646615292082200000000001010101004330043300433004330-1------------------------------------------
158.5350804490143-60.67304929457538158.5350804490143-60.673049294575380.00069317000452420.000693170004524257771.33848425557771.338484255201.804485600074-60.5972700247834286.8877729732709-2.1516725473028360.00.08.61683e-101.27333e-090.00.00.00.00.00.00.00.01.34236e-101.96779e-100.00.00.08.61683e-101.27333e-090.00.00.08.61683e-101.27333e-090.00.00.08.61683e-101.27333e-090.00.00.08.76939e-101.38224e-090.00.00.00.00.00.00.00.01.5739e-102.30312e-100.00.00.08.76939e-101.38224e-090.00.00.08.76939e-101.38224e-090.00.00.08.76939e-101.38224e-090.068241848886371921300000000000.01993990.019939949.2685infinf22.642622.1485infinfinf22.642622.1485inf--inf0.494034inf--inf0.494034infinfinf22.661622.2376infinfinf22.661622.2376inf--inf0.423986inf--inf0.423986inf0.00.022.934822.51960.00.00.00.9409240.9250530.00.00.00.9428050.6750760.00.00.01.01.00.00.00.00.0009368225.17967e-050.00.00.00.1473130.006565720.00.00.00.8508370.9932010.00.00.00.0009129050.0001818750.00.00.0694.116673.0310.00.00.06.430426.482410.024544767364661529208220000000000222200110001100011000110-1------------------------------------------
158.5347190576608-60.67229646721449158.5347190576608-60.672296467214490.00069317000452420.000693170004524257771.33848425557771.338484255201.8030925424928-60.59690900374197286.8872414595131-2.1511104330541590.00.07.63734e-101.21113e-090.00.00.00.00.00.00.00.01.33274e-101.96214e-100.00.00.07.63734e-101.21113e-090.00.00.07.63734e-101.21113e-090.00.00.07.63734e-101.21113e-090.00.00.08.80041e-101.206e-090.00.00.00.00.00.00.00.01.56541e-102.2947e-100.00.00.08.80041e-101.206e-090.00.00.08.80041e-101.206e-090.00.00.08.80041e-101.206e-090.068241848886371921200000000000.02016710.020167154.5728infinf22.638722.2966infinfinf22.638722.2966inf--inf0.342108inf--inf0.342108infinfinf22.792622.292infinfinf22.792622.292inf--inf0.500622inf--inf0.500622inf0.00.022.942222.52250.00.00.00.9271870.9043190.00.00.00.6325020.6019110.00.00.01.01.00.00.00.00.0009998060.0001088570.00.00.00.1533920.01324460.00.00.00.8446420.9863860.00.00.00.0009659910.0002609750.00.00.0722.578710.8330.00.00.05.738566.18250.024544767364661529208220000000000222200110001100011000110-1------------------------------------------
158.533882090657-60.67230986143039158.533889281343-60.67230907672445312.1032943199971312.103294319997157571.2056517511157583.46313715125201.8025646349728-60.59722691309184286.8868934382735-2.1513278958722080.06.97313e-101.26903e-092.45775e-090.00.00.03.42776e-102.59717e-100.00.01.04641e-109.97012e-111.39898e-100.00.06.97313e-101.25882e-092.41794e-090.00.06.97313e-101.1422e-092.2486e-090.00.06.97313e-101.56326e-092.61466e-090.00.06.59396e-101.2041e-092.43547e-090.00.00.03.97634e-103.01766e-100.00.01.23611e-101.24606e-101.63941e-100.00.06.59396e-101.19191e-092.40632e-090.00.06.59396e-101.05144e-092.24153e-090.00.06.59396e-101.56033e-092.62319e-090.068241848885396781300000000000.1845310.19180629.6268inf22.952122.298321.5335infinf22.952122.309421.5466infinf0.6537930.764799infinf0.6427440.762774infinf22.891422.241321.5237infinf22.891422.250121.5414infinf0.6501120.717663infinf0.641340.708706inf0.023.205422.927622.49930.00.00.964190.9443320.9537890.00.01.26850.8500211.01040.00.01.01.00.9999850.00.00.003299260.0006863975.08672e-050.00.00.7524960.1158650.006373730.00.00.2427370.8828440.9934410.00.00.001467960.0006048570.0001346730.00.0766.097709.45737.7960.00.06.676437.886269.80050.024544767364661529208220000000000989801350013400135001340-1------------------------------------------
158.5346250299532-60.67423020035775158.5346250299532-60.67423020035775300.1482389900048300.148238990004857671.2891867866657671.28918678666201.8060009477188-60.59821386973155286.8881734201307-2.15280575570190.00.01.52996e-092.24436e-090.00.00.00.01.63998e-100.00.00.01.45083e-101.74007e-100.00.00.01.52996e-092.19678e-090.00.00.01.52996e-092.13879e-090.00.00.01.52996e-092.25476e-090.00.00.01.57124e-092.30535e-090.00.00.00.02.49397e-100.00.00.01.68753e-102.04526e-100.00.00.01.57124e-092.22941e-090.00.00.01.57124e-092.14123e-090.00.00.01.57124e-092.31758e-090.068241848886054278200000000000.1459390.14593988.2038infinf22.009421.5932infinfinf22.009421.6295inf--inf0.416237inf--inf0.379869infinfinf22.038321.6223infinfinf22.038321.6455inf--inf0.416035inf--inf0.392767inf0.00.022.853722.46730.00.00.00.2483710.321520.00.00.00.7920120.8828720.00.00.01.01.00.00.00.00.0008122821.03941e-050.00.00.00.1297380.001568580.00.00.00.868660.9983310.00.00.00.0007897549.01715e-050.00.00.0562.015637.1110.00.00.010.59559.829050.0245447673646615292081900000000003333001200012000120001205254319336165801600158.5350247939341-60.6743564462372820.113107681274418----0.08823805880436770.40285530686378480.2190316319465637-7.3570546418768553.9909561121800890.0002329662261071--------------------
158.5350681498687-60.67436131162899158.5350694428875-60.67436017503911410.8663644499975410.866364449997557546.1191405686757542.34737957423201.8064926910287-60.59813633506674286.8884271189617-2.1528101421153283.96207e-099.84554e-091.50646e-081.98351e-082.32332e-081.44371e-102.93289e-103.11062e-105.09829e-101.15327e-093.81448e-118.68435e-111.09927e-101.77147e-104.03287e-103.89203e-099.70358e-091.50424e-082.0078e-082.34698e-083.85451e-099.58179e-091.49611e-081.96288e-082.29182e-084.02035e-099.86925e-091.53616e-082.01638e-082.37963e-083.99377e-099.91578e-091.50836e-081.9818e-082.32849e-081.36396e-102.82519e-103.51552e-104.41447e-101.4794e-094.42102e-111.02183e-101.32673e-102.07363e-105.06066e-103.98007e-099.90026e-091.50525e-082.00656e-082.38516e-083.8723e-099.65094e-091.49764e-081.99151e-082.3138e-084.05195e-099.98543e-091.55103e-082.01118e-082.42032e-0868241848885295680600000000000.06148270.062625790.289220.996520.009219.553719.257319.082321.000320.010919.55619.243919.05620.987360.4554410.2963920.1750350.9893890.4549050.3121110.18766221.005220.016919.555119.256419.084721.024620.032719.556719.243219.07370.9882930.4617960.2986930.1716880.9918920.4759640.3135070.16947423.794722.89522.631722.267421.23130.8300730.8766870.8179520.8596470.847451.45380.9213120.933810.9314311.031791.01.01.01.01.00.01307830.09851180.0003384687.94401e-069.76571e-060.4727170.7569670.06000590.001280990.001160020.509860.1217860.9393570.998650.9987350.004345210.02273520.0002982446.16031e-059.49002e-05568.619581.195520.65605.719615.9655.609859.791363.365460.731131.15042454476736466152920819000000000030263026658655575465865557545254319336165801600158.5350247939341-60.6743564462372820.113107681274418----0.08823805880436770.40285530686378480.2190316319465637-7.3570546418768553.9909561121800892.178481659178802e-05--------------------
158.5334141250289-60.67373769342433158.5334141250289-60.67373769342433410.8663644499975410.866364449997557549.4520891793757549.45208917937201.8044509980249-60.59833759030055286.8874126818541-2.152677712896033.86798e-101.46792e-092.11061e-093.50197e-094.12125e-095.39399e-112.22882e-106.4257e-107.59634e-110.03.11687e-116.28619e-111.03264e-101.53025e-105.87528e-103.52968e-101.34824e-092.16029e-093.48384e-094.12125e-093.51247e-101.31876e-091.86703e-093.4752e-094.12125e-093.96466e-101.58674e-092.35813e-093.54848e-094.12125e-093.78225e-101.45069e-092.048e-093.47197e-094.24616e-094.75553e-111.56988e-106.31088e-101.6292e-100.03.68855e-117.56076e-111.23538e-101.80037e-106.89663e-103.98215e-101.4394e-092.09938e-093.55102e-094.24616e-093.61036e-101.30895e-091.82719e-093.45843e-094.24616e-094.06894e-101.51373e-092.34186e-093.60773e-094.24616e-0968241848885296863400000000000.1457140.14571490.196323.555622.096121.721721.148620.9323.499722.104521.694821.124120.931.459560.374390.5731140.2185541.395160.4097690.570660.19410723.531322.083221.68921.139220.962423.630722.175621.663721.144920.96241.448050.394260.5497610.1767921.455080.5118660.5188640.18242624.02823.173522.899622.484221.33230.761560.8030450.7694460.8860560.8165690.9141451.036091.068051.471850.6775871.01.01.01.01.00.01484140.08861880.0004687649.3312e-067.15375e-060.4507270.7570330.07764780.001469120.0008528870.5294110.1338760.9214450.9984440.999050.005020310.02047190.0004382937.75593e-059.0285e-05570.707594.045612.682627.957622.3137.3325111.351412.070714.36727.02924245447673646615292081900000000001616161635431354313543135431-1------------------------------------------
158.5332927230932-60.67377840660529158.5332817603312-60.6737867232082391.8974258300004302.085414379995257596.679340982557592.807504082201.8044339791211-60.59840916754801286.8873816892038-2.1527427814782940.00.02.4459e-093.65607e-090.00.00.05.19081e-105.8119e-100.00.00.01.2478e-102.95232e-100.00.00.02.29839e-093.74233e-090.00.00.02.10152e-093.53685e-090.00.00.02.64755e-093.94781e-090.00.00.02.48503e-093.51957e-090.00.00.04.90444e-104.04935e-100.00.00.01.64251e-104.1902e-100.00.00.02.41704e-093.62698e-090.00.00.02.16984e-093.48381e-090.00.00.02.68655e-093.77014e-090.068241848886028324800000000000.1326480.14105196.8773infinf21.511721.1338infinfinf21.541821.1011inf--inf0.377895inf--inf0.440651infinfinf21.528921.0925infinfinf21.596421.0671inf--inf0.436438inf--inf0.529295inf0.00.022.60121.80790.00.00.00.6605770.5058130.00.00.00.8872160.9819840.00.00.01.01.00.00.00.00.0001776496.18088e-060.00.00.00.03609230.001021190.00.00.00.9635820.9989360.00.00.00.0001484243.66816e-050.00.00.0556.264518.7860.00.00.011.8558.831450.024544767364661529208190000000000858510430003201043000320-1------------------------------------------
158.5328024038264-60.6758947630413158.5328032476187-60.67589445481923410.8663644499975410.866364449997557546.1191405686757542.34737957423201.8073634974059-60.5999802705171286.8882375224567-2.1546934778534952.39441e-084.63656e-086.36578e-088.40412e-089.20174e-082.0285e-102.91016e-101.08894e-091.09049e-092.68918e-091.11507e-102.20822e-102.56049e-104.02568e-106.27168e-102.39452e-084.64233e-086.38341e-088.40884e-089.23074e-082.38728e-084.62183e-086.31701e-088.35213e-089.06669e-082.40179e-084.64398e-086.41076e-088.50823e-089.39503e-082.39578e-084.64516e-086.34938e-088.39173e-089.19171e-081.79835e-102.41556e-109.22596e-109.92781e-102.56605e-091.16696e-102.35782e-102.75051e-104.29512e-107.21869e-102.3971e-084.64505e-086.3673e-088.37718e-089.23521e-082.39326e-084.64503e-086.30848e-088.3296e-089.08085e-082.40376e-084.65994e-086.40197e-088.46539e-089.38922e-0868241848885295673200000000000.05631420.058101687.793719.051418.332517.993217.690417.591519.050818.332517.990117.692317.58640.7188840.3393270.3027990.09886360.7182620.3424110.2978570.10587119.05218.334517.990417.688817.590319.05218.333217.987417.688217.58690.7174950.3441390.3016010.09844590.7187960.3457910.2992080.10125223.274822.421322.193221.852920.99850.992330.9905910.9808740.9744520.9579221.41141.544631.766381.988081.05871.01.01.01.01.00.008115080.08716040.0001881454.76431e-061.00558e-050.419170.7487010.03473620.0007936590.001182790.5697120.1404740.9649070.9991590.9987120.003003750.02366420.0001681784.20584e-059.50394e-05555.414580.473506.067537.961599.933213.033189.921192.869188.044101.5282454476736466152920819000000000030263026658655575465865557545254319336217496320158.53275180106883-60.6758846431168118.42425727844238319.0119705200195317.6426486968994140.03660959779011060.12981143593788140.2820213735103607-3.83761888412787932.7214001660379172.6769265182985333e-05158.532883-60.67586510340789-604033116.27715.13714.9050.139----4.943553900035489e-05
In [8]:
# Extract stars within the dust region polygon
import matplotlib.path as mpath
import matplotlib.pyplot as plt

# Get the polygon vertices
ra_poly, dec_poly = dust_region[0]
vertices = np.column_stack((ra_poly, dec_poly))
path = mpath.Path(vertices)

# Get star positions
star_ra = final_catalog["ra_ok"]
star_dec = final_catalog["dec_ok"]
points = np.column_stack((star_ra, star_dec))

# Check which stars are inside the polygon
inside_mask = path.contains_points(points)
dust_stars = final_catalog[inside_mask]

print(f"Found {len(dust_stars)} stars within the dust region.")

# Count stars with Gaia and 2MASS data
has_gaia = ~np.isnan(dust_stars["gaia_ra"])
has_tmass = dust_stars["tmass_2MASS"] != ""

print(f"Stars with Gaia data: {np.sum(has_gaia)}")
print(f"Stars with 2MASS data: {np.sum(has_tmass)}")

# Check available columns for photometry
print("Available columns in dust_stars:")
print(dust_stars.colnames[:30])  # Print first 30 column names

# Plot histogram of DECaPS g-band photometry
# Try common column names
possible_g_cols = ["g_psf_mag", "gmag", "g", "g_psf", "mean_g"]
g_col = None
for col in possible_g_cols:
    if col in dust_stars.colnames:
        g_col = col
        break

# if g_col:
#     g_mag = dust_stars[g_col]
#     g_mag_valid = g_mag[~np.isnan(g_mag)]

#     plt.figure(figsize=(8, 6))
#     plt.hist(g_mag_valid, bins=50, alpha=0.7, color="blue", edgecolor="black")
#     plt.xlabel(f"DECaPS {g_col}")
#     plt.ylabel("Number of Stars")
#     plt.title(f"Histogram of DECaPS {g_col} in Dust Region")
#     plt.grid(True, alpha=0.3)
#     plt.show()
# else:
#     print("Could not find a suitable g-band magnitude column.")

Found 865 stars within the dust region.
Stars with Gaia data: 163
Stars with 2MASS data: 54
Available columns in dust_stars:
['ra', 'dec', 'ra_ok', 'dec_ok', 'epochrange', 'epochrange_ok', 'epochmean', 'epochmean_ok', 'elon', 'elat', 'glon', 'glat', 'mean_g', 'mean_r', 'mean_i', 'mean_z', 'mean_y', 'stdev_g', 'stdev_r', 'stdev_i', 'stdev_z', 'stdev_y', 'err_g', 'err_r', 'err_i', 'err_z', 'err_y', 'median_g', 'median_r', 'median_i']
In [9]:
# Convert DECaPS fluxes from maggies to AB magnitudes
import numpy as np
from astropy.table import Column

print("Converting DECaPS fluxes from maggies to AB magnitudes...")
print("Note: DECaPS maggies already include the AB zero point (equivalent to 3631 Janskys)")
print("Formula: Mag = -2.5 × log10(flux_in_maggies)")

bands = ["g", "r", "i", "z", "y"]

for band in bands:
    flux_col = f"mean_{band}"
    flux_err_col = f"err_{band}"

    if flux_col in dust_stars.colnames and flux_err_col in dust_stars.colnames:
        # Get flux and flux error
        flux = dust_stars[flux_col]
        flux_err = dust_stars[flux_err_col]

        # Convert flux to magnitude: mag = -2.5 * log10(flux)
        # Only convert positive fluxes (maggies are always positive)
        valid_flux = (flux > 0) & (~np.isnan(flux)) & (~np.isinf(flux))
        mag = np.full_like(flux, np.nan)
        mag_err = np.full_like(flux_err, np.nan)

        if np.sum(valid_flux) > 0:
            mag[valid_flux] = -2.5 * np.log10(flux[valid_flux])

            # Propagate errors: σ_mag ≈ 1.0857 × (flux_err / flux)
            # This is equivalent to: mag_err = (2.5 / ln(10)) * (flux_err / flux)
            mag_err[valid_flux] = 1.0857 * (flux_err[valid_flux] / flux[valid_flux])

        # Add magnitude columns to the table
        mag_col_name = f"{band}_mag"
        mag_err_col_name = f"{band}_mag_err"

        if mag_col_name not in dust_stars.colnames:
            dust_stars.add_column(Column(mag, name=mag_col_name))
        else:
            dust_stars[mag_col_name] = mag

        if mag_err_col_name not in dust_stars.colnames:
            dust_stars.add_column(Column(mag_err, name=mag_err_col_name))
        else:
            dust_stars[mag_err_col_name] = mag_err

        # Print statistics
        valid_mag = ~np.isnan(mag)
        count = np.sum(valid_mag)
        print(f"{band.upper()}-band: {count} stars converted from flux to AB magnitude")
        if count > 0:
            print(f"  Magnitude range: {np.min(mag[valid_mag]):.2f} - {np.max(mag[valid_mag]):.2f} mag")
            print(f"  Median magnitude: {np.median(mag[valid_mag]):.2f} mag")
            print(f"  Median error: {np.median(mag_err[valid_mag]):.3f} mag")
    else:
        print(f"{band.upper()}-band: Flux columns {flux_col} or {flux_err_col} not found")

print(f"\nConversion complete! Added AB magnitude columns: {[f'{b}_mag' for b in bands]}")
print(f"Added magnitude error columns: {[f'{b}_mag_err' for b in bands]}")
print("Note: No zero point was added since maggies already include the AB zero point.")

Converting DECaPS fluxes from maggies to AB magnitudes...
Note: DECaPS maggies already include the AB zero point (equivalent to 3631 Janskys)
Formula: Mag = -2.5 × log10(flux_in_maggies)
G-band: 362 stars converted from flux to AB magnitude
  Magnitude range: 15.81 - 25.03 mag
  Median magnitude: 22.90 mag
  Median error: 0.054 mag
R-band: 367 stars converted from flux to AB magnitude
  Magnitude range: 15.14 - 23.61 mag
  Median magnitude: 22.12 mag
  Median error: 0.057 mag
I-band: 555 stars converted from flux to AB magnitude
  Magnitude range: 13.29 - 23.35 mag
  Median magnitude: 21.46 mag
  Median error: 0.045 mag
Z-band: 572 stars converted from flux to AB magnitude
  Magnitude range: 13.32 - 23.14 mag
  Median magnitude: 21.00 mag
  Median error: 0.046 mag
Y-band: 366 stars converted from flux to AB magnitude
  Magnitude range: 13.32 - 22.05 mag
  Median magnitude: 19.80 mag
  Median error: 0.043 mag

Conversion complete! Added AB magnitude columns: ['g_mag', 'r_mag', 'i_mag', 'z_mag', 'y_mag']
Added magnitude error columns: ['g_mag_err', 'r_mag_err', 'i_mag_err', 'z_mag_err', 'y_mag_err']
Note: No zero point was added since maggies already include the AB zero point.
In [10]:
# Create histogram of DECaPS photometry (using Gaia G magnitudes since DECaPS magnitudes appear uncalibrated)
import matplotlib.pyplot as plt
import numpy as np

# Use Gaia G magnitudes which are properly calibrated
g_mag = dust_stars["gaia_phot_g_mean_mag"]
valid_g = ~np.isnan(g_mag) & ~np.isinf(g_mag)

print(f"Gaia G magnitudes available: {np.sum(valid_g)} out of {len(dust_stars)} stars")
print(f"Gaia G magnitude range: {np.min(g_mag[valid_g]):.2f} to {np.max(g_mag[valid_g]):.2f}")

# Create histogram
plt.figure(figsize=(10, 6))
plt.hist(g_mag[valid_g], bins=30, alpha=0.7, color="blue", edgecolor="black")
plt.xlabel("Gaia G Magnitude")
plt.ylabel("Number of Stars")
plt.title("Histogram of Gaia G Magnitudes for Stars in Dust Region")
plt.grid(True, alpha=0.3)
plt.gca().invert_xaxis()  # Brighter stars on the right
plt.tight_layout()
plt.show()

print("\nNote: DECaPS mean_g column contains very small values (~1e-10) and mean_mag_* columns contain 'inf' values,")
print("suggesting the DECaPS magnitudes in this dataset are not properly calibrated.")
print("Using Gaia G magnitudes instead, which are properly calibrated to the Vega system.")

Gaia G magnitudes available: 163 out of 865 stars
Gaia G magnitude range: 13.30 to 21.21
No description has been provided for this image 
Note: DECaPS mean_g column contains very small values (~1e-10) and mean_mag_* columns contain 'inf' values,
suggesting the DECaPS magnitudes in this dataset are not properly calibrated.
Using Gaia G magnitudes instead, which are properly calibrated to the Vega system.
In [11]:
# Create comprehensive photometry histogram for all available bands (including DECaPS AB magnitudes)
import matplotlib.pyplot as plt
import numpy as np

# Define available magnitude columns and their labels (including DECaPS AB magnitudes)
mag_columns = {
    "g_mag": "DECaPS g (AB)",
    "r_mag": "DECaPS r (AB)",
    "i_mag": "DECaPS i (AB)",
    "z_mag": "DECaPS z (AB)",
    "y_mag": "DECaPS Y (AB)",
    "gaia_phot_g_mean_mag": "Gaia G",
    "gaia_phot_bp_mean_mag": "Gaia BP",
    "gaia_phot_rp_mean_mag": "Gaia RP",
    "tmass_Jmag": "2MASS J",
    "tmass_Hmag": "2MASS H",
    "tmass_Kmag": "2MASS K",
}

# Create subplots for all bands (11 bands, so 3x4 grid)
fig, axes = plt.subplots(3, 4, figsize=(20, 15))
axes = axes.flatten()

for i, (col, label) in enumerate(mag_columns.items()):
    if col in dust_stars.colnames:
        mag = dust_stars[col]
        valid = ~np.isnan(mag) & ~np.isinf(mag)
        count = np.sum(valid)

        if count > 0:
            axes[i].hist(mag[valid], bins=20, alpha=0.7, color="skyblue", edgecolor="black")
            axes[i].set_xlabel(f"{label} Magnitude")
            axes[i].set_ylabel("Number of Stars")
            axes[i].set_title(f"{label} Magnitudes\n({count} stars)")
            axes[i].grid(True, alpha=0.3)
            axes[i].invert_xaxis()  # Brighter stars on the right
        else:
            axes[i].text(0.5, 0.5, "No valid data", ha="center", va="center", transform=axes[i].transAxes)
            axes[i].set_title(f"{label} Magnitudes\n(0 stars)")
    else:
        axes[i].text(0.5, 0.5, "Column not found", ha="center", va="center", transform=axes[i].transAxes)
        axes[i].set_title(f"{label} Magnitudes\n(N/A)")

# Hide the last subplot (12th one in 3x4 grid)
axes[11].set_visible(False)

plt.suptitle("Comprehensive Photometry Histograms for Stars in Dust Region\n(DECaPS AB magnitudes + Gaia + 2MASS)", fontsize=16)
plt.tight_layout()
plt.show()

# Summary statistics
print("Summary of photometry data availability:")
for col, label in mag_columns.items():
    if col in dust_stars.colnames:
        mag = dust_stars[col]
        valid = ~np.isnan(mag) & ~np.isinf(mag)
        count = np.sum(valid)
        if count > 0:
            print(f"{label}: {count} stars, range {np.min(mag[valid]):.2f} - {np.max(mag[valid]):.2f} mag")
        else:
            print(f"{label}: 0 valid measurements")
    else:
        print(f"{label}: column not found")
No description has been provided for this image 
Summary of photometry data availability:
DECaPS g (AB): 362 stars, range 15.81 - 25.03 mag
DECaPS r (AB): 367 stars, range 15.14 - 23.61 mag
DECaPS i (AB): 555 stars, range 13.29 - 23.35 mag
DECaPS z (AB): 572 stars, range 13.32 - 23.14 mag
DECaPS Y (AB): 366 stars, range 13.32 - 22.05 mag
Gaia G: 163 stars, range 13.30 - 21.21 mag
Gaia BP: 145 stars, range 13.51 - 23.30 mag
Gaia RP: 147 stars, range 12.94 - 20.26 mag
2MASS J: 54 stars, range 12.02 - 17.07 mag
2MASS H: 54 stars, range 10.60 - 15.81 mag
2MASS K: 54 stars, range 10.10 - 16.43 mag
In [12]:
# Create color-magnitude and color-color diagrams to validate photometry
import matplotlib.pyplot as plt
import numpy as np

print("Creating color-magnitude and color-color diagrams...")

# Color-magnitude diagram: g vs g-r
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

# Left plot: Color-magnitude diagram
g_mag = dust_stars["g_mag"]
r_mag = dust_stars["r_mag"]
gr_color = g_mag - r_mag

# Filter for valid colors
valid_gr = (~np.isnan(g_mag)) & (~np.isnan(r_mag))
gr_color_valid = gr_color[valid_gr]
g_mag_valid = g_mag[valid_gr]

axes[0].scatter(gr_color_valid, g_mag_valid, s=2, alpha=0.6, color="blue", edgecolors="none")
axes[0].set_xlabel("g - r (mag)")
axes[0].set_ylabel("g (mag)")
axes[0].set_title(f"g vs g-r Color-Magnitude Diagram\n({np.sum(valid_gr)} stars)")
axes[0].invert_yaxis()  # Brighter stars at top
axes[0].grid(True, alpha=0.3)

# Right plot: Color-color diagram: r-i vs g-r
i_mag = dust_stars["i_mag"]
ri_color = r_mag - i_mag

# Filter for valid colors
valid_color_color = valid_gr & (~np.isnan(i_mag))
gr_color_cc = gr_color[valid_color_color]
ri_color_cc = ri_color[valid_color_color]

axes[1].scatter(gr_color_cc, ri_color_cc, s=2, alpha=0.6, color="red", edgecolors="none")
axes[1].set_xlabel("g - r (mag)")
axes[1].set_ylabel("r - i (mag)")
axes[1].set_title(f"g-r vs r-i Color-Color Diagram\n({np.sum(valid_color_color)} stars)")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print some statistics about the colors
print(f"\nColor Statistics:")
if np.sum(valid_gr) > 0:
    print(f"g-r color range: {np.min(gr_color_valid):.2f} - {np.max(gr_color_valid):.2f} mag")
    print(f"Median g-r color: {np.median(gr_color_valid):.2f} mag")

if np.sum(valid_color_color) > 0:
    print(f"r-i color range: {np.min(ri_color_cc):.2f} - {np.max(ri_color_cc):.2f} mag")
    print(f"Median r-i color: {np.median(ri_color_cc):.2f} mag")

print(f"\nDECaPS photometry successfully converted from maggies to AB magnitudes!")
print(f"All magnitude errors properly propagated from flux uncertainties.")

Creating color-magnitude and color-color diagrams...
No description has been provided for this image 
Color Statistics:
g-r color range: 0.28 - 2.41 mag
Median g-r color: 1.37 mag
r-i color range: -0.05 - 1.64 mag
Median r-i color: 0.66 mag

DECaPS photometry successfully converted from maggies to AB magnitudes!
All magnitude errors properly propagated from flux uncertainties.
In [13]:
# Analyze band coverage for stars in dust region
import matplotlib.pyplot as plt

bands = ["mean_g", "mean_r", "mean_i", "mean_z", "mean_y"]
band_coverage = []

for star in dust_stars:
    count = 0
    for band in bands:
        val = star[band]
        if not np.isnan(val) and val != 0.0:
            count += 1
    band_coverage.append(count)

# Plot the distribution
unique_counts, frequencies = np.unique(band_coverage, return_counts=True)

plt.figure(figsize=(8, 6))
plt.bar(unique_counts, frequencies, color="skyblue", edgecolor="black", alpha=0.7)
plt.xlabel("Number of Bands with Valid Data")
plt.ylabel("Number of Stars")
plt.title(f"DECaPS Band Coverage Distribution in Dust Region ({len(dust_stars)} stars)")
plt.xticks(range(6))  # 0 to 5 bands
plt.grid(True, alpha=0.3)
plt.show()

print("Band coverage summary:")
for num_bands, freq in zip(unique_counts, frequencies):
    print(f"  {freq} stars have data in {num_bands} band(s)")
No description has been provided for this image 
Band coverage summary:
  366 stars have data in 1 band(s)
  124 stars have data in 2 band(s)
  82 stars have data in 3 band(s)
  103 stars have data in 4 band(s)
  190 stars have data in 5 band(s)
In [14]:
# Plot Gaia parallax vs G-band magnitude for stars in dust region
import matplotlib.pyplot as plt

# Extract Gaia data
gaia_parallax = dust_stars["gaia_parallax"]
gaia_g_mag = dust_stars["gaia_phot_g_mean_mag"]

# Filter out invalid values (NaN)
valid_mask = (~np.isnan(gaia_parallax)) & (~np.isnan(gaia_g_mag))
parallax_valid = gaia_parallax[valid_mask]
g_mag_valid = gaia_g_mag[valid_mask]

print(f"Plotting {len(parallax_valid)} stars with valid Gaia parallax and G-magnitude.")

# Create scatter plot
plt.figure(figsize=(10, 8))
plt.scatter(parallax_valid, g_mag_valid, s=2, alpha=0.6, color="purple", edgecolors="none")
plt.xlabel("Gaia Parallax (mas)")
plt.ylabel("Gaia G-band Magnitude")
plt.title(f"Gaia Parallax vs G-magnitude in Dust Region ({len(dust_stars)} stars total, {len(parallax_valid)} with data)")
plt.gca().invert_yaxis()  # Magnitudes increase downward
plt.grid(True, alpha=0.3)
plt.xlim(left=0)  # Parallax can't be negative
plt.show()

Plotting 146 stars with valid Gaia parallax and G-magnitude.
No description has been provided for this image
In [15]:
# Histogram of Gaia parallax percentage errors
import matplotlib.pyplot as plt

# Extract Gaia parallax data
gaia_parallax = dust_stars["gaia_parallax"]

# Since we can't get parallax_error from the current query, let's use a typical Gaia parallax error
# For Gaia DR3, typical parallax errors range from ~0.02 mas for bright stars to ~0.5 mas for faint stars
# For a rough estimate, let's use 0.1 mas as a typical error
typical_parallax_error = 0.1  # mas

# Filter for valid positive parallax values
valid_mask = (~np.isnan(gaia_parallax)) & (gaia_parallax > 0)
parallax_valid = gaia_parallax[valid_mask]

print(f"Using typical parallax error of {typical_parallax_error} mas for all stars.")
print(f"Plotting parallax percentage errors for {len(parallax_valid)} stars with valid positive parallax.")

# Calculate percentage error using typical error
percentage_error = (typical_parallax_error / parallax_valid) * 100

# Create histogram
plt.figure(figsize=(10, 6))
plt.hist(percentage_error, bins=50, alpha=0.7, color="green", edgecolor="black")
plt.xlabel("Gaia Parallax Percentage Error (%) [using typical 0.1 mas error]")
plt.ylabel("Number of Stars")
plt.title(f"Histogram of Estimated Gaia Parallax Percentage Errors in Dust Region ({len(parallax_valid)} stars)")
plt.grid(True, alpha=0.3)
plt.xlim(left=0)  # Percentage error can't be negative
plt.show()

# Print some statistics
print(f"\nStatistics for estimated parallax percentage errors:")
print(f"  Mean: {np.mean(percentage_error):.2f}%")
print(f"  Median: {np.median(percentage_error):.2f}%")
print(f"  Standard deviation: {np.std(percentage_error):.2f}%")
print(f"  Min: {np.min(percentage_error):.2f}%")
print(f"  Max: {np.max(percentage_error):.2f}%")

print("\nNote: These are estimates using a typical parallax error of 0.1 mas.")
print("Actual errors vary with stellar brightness and would be better obtained from Gaia directly.")

Using typical parallax error of 0.1 mas for all stars.
Plotting parallax percentage errors for 113 stars with valid positive parallax.
No description has been provided for this image 
Statistics for estimated parallax percentage errors:
  Mean: 58.00%
  Median: 25.10%
  Standard deviation: 110.35%
  Min: 4.72%
  Max: 943.76%

Note: These are estimates using a typical parallax error of 0.1 mas.
Actual errors vary with stellar brightness and would be better obtained from Gaia directly.

/Users/kshukawa/miniforge3/envs/light_echo/lib/python3.12/site-packages/numpy/_core/fromnumeric.py:867: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedColumn.
  a.partition(kth, axis=axis, kind=kind, order=order)
In [16]:
# Plot NMAG_OK histograms for each DECaPS band (excluding NMAG_OK=0) with exposure times in titles
import matplotlib.pyplot as plt
import numpy as np

print("Plotting NMAG_OK distributions for each DECaPS band (excluding non-detections)...")

# DECaPS band information with exposure times
decaps_bands = {"g": {"exposure_time": 96, "wavelength": "398.0–548.5 nm"}, "r": {"exposure_time": 30, "wavelength": "565.5–717.0 nm"}, "i": {"exposure_time": 30, "wavelength": "704.5–858.0 nm"}, "z": {"exposure_time": 30, "wavelength": "846.5–1000.0 nm"}, "y": {"exposure_time": 30, "wavelength": "950.0–1034.0 nm"}}

# Create subplots for all bands
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.flatten()

for i, (band, info) in enumerate(decaps_bands.items()):
    n_ok_col = f"nmag_ok_{band}"

    if n_ok_col in dust_stars.colnames:
        n_ok_values = dust_stars[n_ok_col]
        # Filter out NMAG_OK = 0 (non-detections) and invalid values
        valid_detected = (~np.isnan(n_ok_values)) & (n_ok_values >= 1)
        count = np.sum(valid_detected)

        if count > 0:
            detected_values = n_ok_values[valid_detected]
            # Get unique values and their counts for better visualization
            unique_vals, frequencies = np.unique(detected_values, return_counts=True)

            axes[i].bar(unique_vals, frequencies, alpha=0.7, color="steelblue", edgecolor="navy", width=0.8)
            axes[i].set_xlabel(f"NMAG_OK ({band}-band)")
            axes[i].set_ylabel("Number of Stars")
            axes[i].set_title(f"DECaPS {band.upper()}-band Detections\n{info['wavelength']} | {info['exposure_time']}s exposure\n({count} stars with NMAG_OK ≥ 1)")
            axes[i].grid(True, alpha=0.3)
            axes[i].set_xticks(unique_vals)  # Show all integer values on x-axis

            # Add statistics as text
            mean_det = np.mean(detected_values)
            median_det = np.median(detected_values)
            axes[i].text(0.02, 0.98, f"Mean: {mean_det:.1f}\nMedian: {median_det:.1f}", transform=axes[i].transAxes, fontsize=9, verticalalignment="top", bbox=dict(boxstyle="round", facecolor="lightblue", alpha=0.8))
        else:
            axes[i].text(0.5, 0.5, "No detections\n(NMAG_OK ≥ 1)", ha="center", va="center", transform=axes[i].transAxes)
            axes[i].set_title(f"{band.upper()}-band Detections\n{info['wavelength']} | {info['exposure_time']}s exposure\n(0 stars detected)")
    else:
        axes[i].text(0.5, 0.5, f"Column {n_ok_col}\nnot found", ha="center", va="center", transform=axes[i].transAxes)
        axes[i].set_title(f"{band.upper()}-band Detections\n{info['wavelength']} | {info['exposure_time']}s exposure\n(Column missing)")

# Hide the last subplot (6th one)
axes[5].set_visible(False)

plt.suptitle(f"DECaPS NMAG_OK Distributions by Band\nTotal stars in dust region: {len(dust_stars)}", fontsize=16)
plt.tight_layout()
plt.show()

# Print summary statistics for detected stars only
print("\n" + "=" * 90)
print("DECaPS NMAG_OK SUMMARY (Detected Stars Only - NMAG_OK ≥ 1)")
print("=" * 90)

total_stars = len(dust_stars)
print(f"Total stars in dust region: {total_stars}")

for band, info in decaps_bands.items():
    n_ok_col = f"nmag_ok_{band}"

    if n_ok_col in dust_stars.colnames:
        n_ok_values = dust_stars[n_ok_col]
        # Detected stars (NMAG_OK >= 1)
        detected = (~np.isnan(n_ok_values)) & (n_ok_values >= 1)
        count_detected = np.sum(detected)

        # Non-detected stars (NMAG_OK = 0)
        non_detected = (~np.isnan(n_ok_values)) & (n_ok_values == 0)
        count_non_detected = np.sum(non_detected)

        detection_rate = (count_detected / total_stars) * 100

        print(f"\n{band.upper()}-band ({info['wavelength']}, {info['exposure_time']}s exposure):")
        print(f"  Stars detected (NMAG_OK ≥ 1): {count_detected} ({detection_rate:.1f}%)")
        print(f"  Stars not detected (NMAG_OK = 0): {count_non_detected}")

        if count_detected > 0:
            detected_values = n_ok_values[detected]
            unique_vals, frequencies = np.unique(detected_values, return_counts=True)

            print(f"  Detection distribution: {dict(zip(unique_vals.astype(int), frequencies))}")

            # Statistics for detected stars
            mean_det = np.mean(detected_values)
            median_det = np.median(detected_values)
            max_det = np.max(detected_values)

            print(f"  Mean NMAG_OK (detected): {mean_det:.2f}")
            print(f"  Median NMAG_OK (detected): {median_det:.1f}")
            print(f"  Max NMAG_OK (detected): {int(max_det)}")
        else:
            print("  No stars detected in this band")
    else:
        print(f"\n{band.upper()}-band: Column {n_ok_col} not found")

print(f"\nNote: This analysis excludes stars with NMAG_OK = 0 (non-detections)")
print(f"Only showing stars that were successfully detected in each band.")

Plotting NMAG_OK distributions for each DECaPS band (excluding non-detections)...
No description has been provided for this image 
==========================================================================================
DECaPS NMAG_OK SUMMARY (Detected Stars Only - NMAG_OK ≥ 1)
==========================================================================================
Total stars in dust region: 865

G-band (398.0–548.5 nm, 96s exposure):
  Stars detected (NMAG_OK ≥ 1): 362 (41.8%)
  Stars not detected (NMAG_OK = 0): 503
  Detection distribution: {np.int64(1): np.int64(122), np.int64(2): np.int64(31), np.int64(3): np.int64(20), np.int64(4): np.int64(16), np.int64(5): np.int64(75), np.int64(6): np.int64(98)}
  Mean NMAG_OK (detected): 3.51
  Median NMAG_OK (detected): 4.0
  Max NMAG_OK (detected): 6

R-band (565.5–717.0 nm, 30s exposure):
  Stars detected (NMAG_OK ≥ 1): 367 (42.4%)
  Stars not detected (NMAG_OK = 0): 498
  Detection distribution: {np.int64(1): np.int64(157), np.int64(2): np.int64(17), np.int64(3): np.int64(20), np.int64(4): np.int64(20), np.int64(5): np.int64(29), np.int64(6): np.int64(123), np.int64(7): np.int64(1)}
  Mean NMAG_OK (detected): 3.33
  Median NMAG_OK (detected): 3.0
  Max NMAG_OK (detected): 7

I-band (704.5–858.0 nm, 30s exposure):
  Stars detected (NMAG_OK ≥ 1): 555 (64.2%)
  Stars not detected (NMAG_OK = 0): 310
  Detection distribution: {np.int64(1): np.int64(167), np.int64(2): np.int64(62), np.int64(3): np.int64(31), np.int64(4): np.int64(26), np.int64(5): np.int64(86), np.int64(6): np.int64(72), np.int64(7): np.int64(110), np.int64(8): np.int64(1)}
  Mean NMAG_OK (detected): 3.83
  Median NMAG_OK (detected): 4.0
  Max NMAG_OK (detected): 8

Z-band (846.5–1000.0 nm, 30s exposure):
  Stars detected (NMAG_OK ≥ 1): 572 (66.1%)
  Stars not detected (NMAG_OK = 0): 293
  Detection distribution: {np.int64(1): np.int64(172), np.int64(2): np.int64(55), np.int64(3): np.int64(61), np.int64(4): np.int64(136), np.int64(5): np.int64(144), np.int64(6): np.int64(4)}
  Mean NMAG_OK (detected): 3.06
  Median NMAG_OK (detected): 3.0
  Max NMAG_OK (detected): 6

Y-band (950.0–1034.0 nm, 30s exposure):
  Stars detected (NMAG_OK ≥ 1): 366 (42.3%)
  Stars not detected (NMAG_OK = 0): 499
  Detection distribution: {np.int64(1): np.int64(82), np.int64(2): np.int64(65), np.int64(3): np.int64(97), np.int64(4): np.int64(120), np.int64(5): np.int64(2)}
  Mean NMAG_OK (detected): 2.71
  Median NMAG_OK (detected): 3.0
  Max NMAG_OK (detected): 5

Note: This analysis excludes stars with NMAG_OK = 0 (non-detections)
Only showing stars that were successfully detected in each band.
In [20]:
import numpy as np
from astropy.table import Table, Column
from astropy.coordinates import SkyCoord

import astropy.units as u

# Filter dust_stars for those with more than 4 band coverage
bands = ["g", "r", "i", "z", "y"]
band_coverage = []
for star in dust_stars:
    count = sum(1 for band in bands if not np.isnan(star[f"mean_{band}"]) and star[f"mean_{band}"] > 0)
    band_coverage.append(count)

dust_stars_filtered = dust_stars[np.array(band_coverage) == 4]

print(f"Filtered to {len(dust_stars_filtered)} stars with more than 4 band coverage.")

# Create new table with required columns
output_table = Table()

# RA and Dec
output_table["ra"] = dust_stars_filtered["ra_ok"]
output_table["dec"] = dust_stars_filtered["dec_ok"]

# Convert to galactic coordinates
coords = SkyCoord(ra=output_table["ra"] * u.deg, dec=output_table["dec"] * u.deg, frame="icrs")
galactic = coords.galactic
output_table["l"] = galactic.l.deg
output_table["b"] = galactic.b.deg

# DECaPS photometry in maggies (already in maggies)
for band in bands:
    output_table[f"{band}_maggies"] = dust_stars_filtered[f"mean_{band}"]
    output_table[f"{band}_err_maggies"] = dust_stars_filtered[f"err_{band}"]

# 2MASS photometry converted to maggies
tmass_bands = ["J", "H", "K"]
for band in tmass_bands:
    mag_col = f"tmass_{band}mag"
    err_col = f"tmass_e_{band}mag"
    if mag_col in dust_stars_filtered.colnames and err_col in dust_stars_filtered.colnames:
        mag = dust_stars_filtered[mag_col]
        err_mag = dust_stars_filtered[err_col]
        # Convert mag to maggies
        maggies = 10 ** (-0.4 * mag)
        # Error in maggies: relative error same, but absolute error = maggies * (err_mag / 2.5) * ln(10)
        err_maggies = maggies * (err_mag / 2.5) * np.log(10)
        output_table[f"{band.lower()}_maggies"] = maggies
        output_table[f"{band.lower()}_err_maggies"] = err_maggies

# Parallax and error in mas
output_table["parallax_mas"] = dust_stars_filtered["gaia_parallax"]
output_table["parallax_err_mas"] = dust_stars_filtered["gaia_parallax_error"]

# Display the table
output_table

Filtered to 103 stars with more than 4 band coverage.
Out[20]:
Table length=103
radeclbg_maggiesg_err_maggiesr_maggiesr_err_maggiesi_maggiesi_err_maggiesz_maggiesz_err_maggiesy_maggiesy_err_maggiesj_maggiesj_err_maggiesh_maggiesh_err_maggiesk_maggiesk_err_maggiesparallax_masparallax_err_mas
float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
158.5947928031036-60.6401251640014286.8965684810015-2.10851012974101342.28421e-104.74727e-111.21903e-091.09406e-101.16717e-091.23766e-102.16831e-091.55915e-100.00.0----------------
158.5989559389191-60.64292213639902286.8997389822885-2.1099067537532122.60827e-081.2036e-100.00.01.03893e-074.62584e-101.36252e-076.80667e-101.52042e-071.01832e-094.2305808577652356e-074.364090043172696e-086.385574850788385e-077.586918002197475e-081e-061.4828647998881655e-070.34375236725749810.0901016369462013
158.5974325832476-60.64137028356316286.89831376468317-2.10893849201166361.55548e-093.37347e-110.00.07.11613e-091.08645e-109.945e-091.52964e-101.08459e-084.04673e-10----------------
158.5967385083458-60.64147393196151286.8980711447672-2.10919878674265433.46099e-093.81461e-110.00.01.95522e-081.44787e-102.78455e-082.02911e-103.42594e-084.28664e-10------------0.54801637745658080.3556091785430908
158.5956551725987-60.64089701988232286.8973218271862-2.10896594436806461.4649e-093.28921e-110.00.01.07177e-081.07621e-101.64725e-081.69014e-101.92844e-083.91553e-10------------0.01059593719448310.798931360244751
158.5931517997432-60.64272521802911286.89717654690315-2.11116311370453952.58106e-093.44762e-110.00.04.14429e-082.36318e-105.90172e-083.06537e-108.45914e-086.2043e-101.1492115972032396e-065.927392783011969e-082.1717008433479926e-061.1201178213871935e-072.8028509276466845e-061.5747278743876446e-07----
158.5932567927062-60.64283745585533286.8972774299903-2.111234410457282.39533e-096.01596e-110.00.02.87747e-082.39616e-106.81609e-086.06549e-109.07013e-081.29424e-091.1492115972032396e-065.927392783011969e-082.1717008433479926e-061.1201178213871935e-072.8028509276466845e-061.5747278743876446e-072.1177617449224090.7160475254058838
158.5810742169801-60.63327795590098286.88730854421954-2.10595978737754136.68199e-103.3175e-110.00.08.04681e-099.25394e-111.41633e-081.61714e-101.75399e-083.85558e-10----------------
158.575901311978-60.63211414726069286.8845283002249-2.10622553360858473.21127e-105.49119e-110.00.03.79112e-099.83104e-116.80243e-091.43828e-107.69086e-093.93115e-10----------------
..................................................................
158.6017018270062-60.65599713438976286.90746391439865-2.12054506850667930.00.01.33834e-096.70051e-113.44832e-098.96096e-117.06496e-091.60788e-109.70223e-094.74839e-10----------------
158.5979635407027-60.65812130088608286.90694349487217-2.1233014167175930.00.01.01424e-095.83636e-113.70507e-098.97391e-117.44689e-091.61924e-101.09257e-084.75306e-10----------------
158.5968050252728-60.65964634193449286.9072171516094-2.1249055451105660.00.01.31302e-095.95977e-112.03282e-091.05137e-103.83101e-091.78759e-104.26042e-096.04555e-10----------------
158.5847809562349-60.65364903847635286.89910623059734-2.1226716775582830.00.01.29463e-091.29617e-102.02004e-098.89751e-113.6883e-091.60766e-104.51572e-095.29651e-10----------------
158.587389991628-60.67383025422208286.910341490298-2.13949008703596230.00.09.20102e-101.28702e-102.50865e-091.15127e-104.25805e-091.57765e-106.12518e-094.49872e-10----------------
158.5872685610707-60.67385875265794286.91030430316624-2.1395445700617683.23998e-104.94823e-111.1551e-096.67185e-113.26541e-091.47375e-104.79655e-093.91267e-100.00.0----------------
158.5864833026755-60.67404553676653286.91006507850017-2.1398990543053540.00.01.40247e-097.05027e-113.4679e-091.22289e-106.71014e-091.68784e-107.9152e-094.64837e-106.228737089440282e-07--6.674212038584765e-078.606047042614942e-081.3905928754654458e-06------
158.5860363152214-60.67564321356073286.9106774064192-2.14139105775381780.00.05.95543e-101.06189e-101.67508e-099.53518e-113.74226e-091.43919e-104.65532e-095.91468e-10----------------
158.5838093008564-60.67550858704951286.90966559358816-2.141821629947260.00.07.68631e-101.07218e-101.49109e-091.39459e-102.67672e-091.52899e-103.8849e-095.92639e-10----------------
158.578816779123-60.67540199793299286.90749534365705-2.1429558952724650.00.01.47198e-096.4407e-113.97183e-099.98294e-117.17831e-091.5998e-109.66441e-094.76911e-10----------------
In [ ]:
output_table.write("dust_region_stars_photometry.csv", format="csv", overwrite=True)
In [21]:
output_table.write("dust_region_stars_photometry_low_coverage.csv", format="csv", overwrite=True)
In [ ]: