Skip to content
Snippets Groups Projects
Commit e673dc31 authored by Harm Munk's avatar Harm Munk
Browse files

first commit after some major revisions

parent 4acc78c0
Branches
No related tags found
No related merge requests found
#!/usr/bin/env python3
from typing import Any
import astropy.units as u
from astropy.coordinates import SkyCoord
import numpy as np
import time
import os
from tqdm import tqdm
from numpy import sin, cos
import matplotlib.pyplot as plt
import zmq
import sys
import warnings
import tomllib
import pathvalidate
warnings.filterwarnings("ignore")
def process_command_line() -> dict[str, Any] | None:
if len(sys.argv) > 1:
with open(sys.argv[1], "rb") as f:
return tomllib.load(f)
def main():
params: dict[str, Any] = \
{'object': {
'name': 'm31'
},
'observation': {
'integration_time': 60,
'offset': 0.25,
'nop_longitude': 25,
'nop_latitude': 7
}
}
params |= process_command_line()
# print(params)
INTEGRATION_TIME = params["observation"]["integration_time"]
OFFSET = params["observation"]["offset"]
N_ALONG = params["observation"]["nop_longitude"]
N_ORTHO = params["observation"]["nop_latitude"]
ANGLE = params ["observation"]["rotation_angle"]
try:
from telescope import Telescope
have_telescope = True
except ImportError:
have_telescope = False
print("Telescope class not imported")
if have_telescope:
try:
dt = Telescope(consoleHost="pc2019", setmode="J2000")
except OSError:
print("dt object not created.")
have_telescope = False
if not have_telescope:
print("Telescope not found!")
context = zmq.Context()
hi_viewer = context.socket(zmq.PUB)
hi_viewer.bind("tcp://127.0.0.1:5556")
obj = SkyCoord.from_name(params["object"]["name"])
center_ra = round(obj.ra.deg / OFFSET) * OFFSET * u.deg
center_dec = round(obj.dec.deg / OFFSET) * OFFSET * u.deg
center = SkyCoord(ra=center_ra, dec=center_dec)
ortho_offsets = (np.arange(N_ORTHO) - N_ORTHO // 2) * OFFSET
along_offsets = (np.arange(N_ALONG) - N_ALONG // 2) * OFFSET
angle = ANGLE * u.deg
offsets_grid = np.array(np.meshgrid(along_offsets, ortho_offsets))
rotation_matrix = np.array([[cos(angle), sin(angle)], [-sin(angle), cos(angle)]])
offsets = rotation_matrix @ offsets_grid.reshape(2, -1)
pointings = center.spherical_offsets_by(offsets[0] * u.deg, offsets[1] * u.deg)
indices = np.indices((N_ORTHO, N_ALONG)).reshape(2, -1).T
cwd = os.path.dirname(os.path.realpath(__file__))
fig, ax = plt.subplots()
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_aspect(1)
radeclist = []
for pointing, idx in tqdm(zip(pointings, indices), leave=False, total=len(pointings)):
# filename = cwd + "/" + pathvalidate.sanitize_filename(params["object"]["name"]) + "_" +\
# str(idx[0])[0:2] + "_" + str(idx[1])[0:2] + ".ecsv"
filename = f"{cwd}/{pathvalidate.sanitize_filename(params["object"]["name"])}_{idx[0]:02d}_{idx[1]:02d}.ecsv"
radeclist.append(
(pointing.ra.deg - center.ra.deg, pointing.dec.deg - center.dec.deg)
)
ax.plot([radeclist[-1][0]], [radeclist[-1][1]], ".")
if os.path.isfile(filename):
continue
if have_telescope:
dt.setRaDec(pointing)
time.sleep(1)
dt.waitUntilThere()
hi_viewer.send_string("reset")
hi_viewer.send_string(f"integration {INTEGRATION_TIME}")
for i in tqdm(range(INTEGRATION_TIME), leave=False):
plt.pause(1)
hi_viewer.send_string(f"save {filename}")
if __name__ == "__main__":
main()
#!/usr/bin/env python3
import astropy.units as u
from astropy.coordinates import SkyCoord
import numpy as np
import time
import os
from tqdm import tqdm
from numpy import sin, cos
import matplotlib.pyplot as plt
import zmq
import warnings
warnings.filterwarnings("ignore")
INTEGRATION_TIME = 60
OFFSET = 0.25
N_ALONG = 25
N_ORTHO = 7
try:
from telescope import Telescope
have_telescope = True
except ImportError:
have_telescope = False
if have_telescope:
try:
dt = Telescope(consoleHost="console", setmode="J2000")
except OSError:
have_telescope = False
if not have_telescope:
print("Telescope not found!")
context = zmq.Context()
hi_viewer = context.socket(zmq.PUB)
hi_viewer.bind("tcp://127.0.0.1:5556")
m31 = SkyCoord.from_name("M31")
center_ra = round(m31.ra.deg / OFFSET) * OFFSET * u.deg
center_dec = round(m31.dec.deg / OFFSET) * OFFSET * u.deg
center = SkyCoord(ra=center_ra, dec=center_dec)
ortho_offsets = (np.arange(N_ORTHO) - N_ORTHO // 2) * OFFSET
along_offsets = (np.arange(N_ALONG) - N_ALONG // 2) * OFFSET
angle = -52 * u.deg
offsets_grid = np.array(np.meshgrid(along_offsets, ortho_offsets))
rotation_matrix = np.array([[cos(angle), sin(angle)], [-sin(angle), cos(angle)]])
offsets = rotation_matrix @ offsets_grid.reshape(2, -1)
pointings = center.spherical_offsets_by(offsets[0] * u.deg, offsets[1] * u.deg)
indices = np.indices((N_ORTHO, N_ALONG)).reshape(2, -1).T
cwd = os.path.dirname(os.path.realpath(__file__))
fig, ax = plt.subplots()
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_aspect(1)
radeclist = []
for pointing, idx in tqdm(zip(pointings, indices), leave=False, total=len(pointings)):
filename = f"{cwd}/m31_{idx[0]:02d}_{idx[1]:02d}.ecsv"
radeclist.append(
(pointing.ra.deg - center.ra.deg, pointing.dec.deg - center.dec.deg)
)
ax.plot([radeclist[-1][0]], [radeclist[-1][1]], ".")
if os.path.isfile(filename):
continue
if have_telescope:
dt.setRaDec(pointing)
time.sleep(1)
dt.waitUntilThere()
hi_viewer.send_string("reset")
hi_viewer.send_string(f"integration {INTEGRATION_TIME}")
for i in tqdm(range(INTEGRATION_TIME), leave=False):
plt.pause(1)
hi_viewer.send_string(f"save {filename}")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment