Commit c221ca6d authored by Akke Viitanen's avatar Akke Viitanen
Browse files

Use 'trapz' instead of 'trapezoid' consistently

parent 9860b0a2
Loading
Loading
Loading
Loading
+48 −0
Original line number Diff line number Diff line
#!/usr/bin/env python3
# encoding: utf-8
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: <++>

"""
<+A python3 script template. Remove this line and add your description.+>
"""


import argparse
import glob
import math
import os
import random
import re
import subprocess
import sys
import time

import astropy as ap
import astropy.coordinates as c
import astropy.units as u
import fitsio
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import pyvo

tap_service = pyvo.dal.TAPService("https://datalab.noirlab.edu/tap")
query = f"""
SELECT TOP 100 * FROM lsst_sim.simdr2 WHERE
(0 <= pmode AND pmode <= 1) AND
(4 <= label AND label <= 6)
"""
print(query)

print("Getting the results")
tap_results = tap_service.search(query)

print("Converting to table")
table = tap_results.to_table()

print("Writing to file")
table.write("stars.fits")
+82 −0
Original line number Diff line number Diff line
#!/usr/bin/env python3
# encoding: utf-8
# Author: Akke Viitanen
# Email: akke.viitanen@helsinki.fi
# Date: 2024-06-14 12:21:18

"""
Get cepheid lightcurves
"""

import glob

from astropy.table import Table
import fitsio
import matplotlib.pyplot as plt
import numpy as np


def parse_filename_cepheid(filename):

    """
    Return the parameters from the filename
    """

    parts = filename.replace("/", "_").split("_")
    parts = parts[-10:]
    return {
        "z":      float(parts[0][1:]),
        "y":      float(parts[1][1:]),
        "logte":  np.log10(float(parts[3])),
        "logl":   float(parts[4]),
        "mass":   float(parts[5]),
        "pmode":  {"F": 0, "FO": 1}[parts[8]],
    }


def get_lightcurve_template_cepheid(source):

    """
    Return the closest lightcurve template according the given parameters
    """

    columns = ("z", "y", "logte", "logl", "mass")
    pmode = source["pmode"]
    dist = []

    vec1 = np.array([source[k] for k in columns])

    filenames = glob.glob(f"CEP_lsst_quasar_project/*/*_%s_lsst.csv" % (["F", "FO"][pmode]))
    for filename in filenames:
        parameters = parse_filename_cepheid(filename)
        vec2 = np.array([parameters[k] for k in columns])
        delta = vec1 - vec2
        dist.append(np.sum(delta ** 2))

    idx = np.argmin(dist)
    table = Table.read(filenames[idx])
    return table


import util
def get_lightcurve_normalized_cepheid(source, band):

    # Find the best lightcurve template
    template = get_lightcurve_template_cepheid(source)
    select = (0 <= template["phase"]) * (template["phase"] <= 1)
    template = template[select]

    # The template spells 'y' with a capital 'Y'
    b = band.replace("lsst-", "")
    key = "Y" if b == 'y' else b
    time1 = template["phase"] * template["P"]
    flux1 = util.mag_to_flux(10 ** (-0.4 * template[f"{key}_lsst"]))

    # Normalize the period to the period of the source
    period = source["period%d" % source["pmode"]]
    time2 = time1 / template["P"] * period

    # Normalize the flux to the mean flux of the source
    flux2 = flux1 / flux1.mean() * util.mag_to_flux(source[f"{b}mag"])

    return time2, flux2
+217 B

File added.

No diff preview for this file type.

+217 B

File added.

No diff preview for this file type.

+217 B

File added.

No diff preview for this file type.

Loading