Commit 3f8ec44f authored by Carl Schaffer's avatar Carl Schaffer
Browse files

Implementing lucia's calculation for the lit orientation

parent 8630231c
......@@ -251,12 +251,23 @@ class GrisFitsFile(FitsFile):
Returns: zenit float in degrees
Raises Value error if elevation is not contained in the header
"""
zenit = 90 - self.elevation
return zenit
@property
def elevation(self):
"""
Get elevation angle
Returns: in degrees
Raises Value error if elevation is not contained in the header
"""
elev = self._get_value(['ELEV_ANG', 'ELEVATIO'])
if elev == '':
raise ValueError(f'No elevation tracked for {self.__repr__()}')
zenit = 90 - elev
return zenit
return elev
@property
def slit_orientation(self):
......@@ -267,39 +278,20 @@ class GrisFitsFile(FitsFile):
Returns:
slitorie: float (degrees)
"""
# Pull values from header
# p0 angle, rotates over the course of the year
# Pull Values from header
p0_angle = self.p0_angle
# Azimut and zenit at time of observation start
rotangle = self.ROTANGLE if 'ROTANGLE' in self.header else None
azimut = self._get_value(['AZIMUT'])
zenit = self.zenit
# Slitorientation from header
slitorie = self.SLITORIE
# Calculate declination of disk center from time
gregor_observer = Observer(location=gregor_coords)
date_tz = pytz.utc.localize(self.obs_time)
time = gregor_observer.datetime_to_astropy_time(date_tz)
sun = get_sun(time)
dec = sun.dec
# Calculate parallactic angle as described in the GREGOR derotator manual
# http://www.leibniz-kis.de/fileadmin/user_upload/observatorien/gre_docs/man/GRE-KIS-MAN-0008_v001_derotator.pdf
parallactic_angle = np.degrees(
np.arcsin(np.sin(np.radians(azimut)) * np.cos(gregor_coords.lat.to(u.rad)) / np.cos(dec.to(u.rad))))
self.parallactic_angle = parallactic_angle
const_term = -130 # determined by trial and error to match results from cross correlation with HMI
alpha = -1 * parallactic_angle.to(u.degree).value + azimut - zenit - p0_angle + const_term
if self.derotator_in:
rotangle = self.ROTANGLE
alpha = -2 * rotangle + (const_term / 2)
angle = slitorie - alpha
elev = self.elevation
# Calculate Angel
angle = calculate_slit_angle(
rotangle=rotangle,
azimuth=azimut,
elevation=elev,
p0_angle=p0_angle,
date=self.obs_time
)
# *******************************************************************************************
# **************WARNING!*********************************************************************
......
......@@ -15,11 +15,17 @@ from itertools import chain
from os.path import basename
from pathlib import Path
import astropy.units as u
import numpy as np
import pandas as pd
import pytz
from astroplan import Observer
from astropy.coordinates import SkyCoord, GCRS, get_sun
from astropy.io.fits.hdu.hdulist import fitsopen as fitsopen
from sunpy.coordinates import Helioprojective
from kis_tools.util.calculations import rotate_around_origin
from kis_tools.util.constants import gregor_coords
from kis_tools.util.util import gris_run_number
......@@ -177,6 +183,7 @@ def extract_position_kws_from_manolo(raw_gris_file):
return pos
def telescope_to_hpc(slitposx=None, slitposy=None, p0_angle=None, xcenter=None, ycenter=None):
"""
Transform GRIS header coordinate data to HPC coordinates
......@@ -238,3 +245,57 @@ def calculate_slit_translation(slit_number: int, stepsize_perp: int, slit_orient
dy = slit_number * stepsize_perp * np.cos(np.radians(slit_orientation))
dxdy = np.array((dx, dy))
return dxdy
# Define a utility function to calculate the parallactic angle from a datetime and azimut
def get_parallactic_angle(datetime, az_deg, hpc_x=None, hpc_y=None):
"""
Calculates parallactic angle for GREGOR based on a time and helioprojective coordinates
parallactic angle defined in GREGOR derotator manual:
http://www.leibniz-kis.de/fileadmin/user_upload/observatorien/gre_docs/man/GRE-KIS-MAN-0008_v001_derotator.pdf
"""
gregor_observer = Observer(location=gregor_coords)
date_tz = pytz.utc.localize(datetime)
time = gregor_observer.datetime_to_astropy_time(date_tz)
# Build Sky Coord object, default to sun center at time
if hpc_x and hpc_y:
sun_coords = SkyCoord(hpc_x * u.arcsec, hpc_y * u.arcsec, frame=Helioprojective, observer='earth', obstime=time)
sun_coords.representation_type = 'spherical'
sun_coords = sun_coords.transform_to(GCRS)
else:
sun_coords = get_sun(time)
dec = sun_coords.dec.to(u.rad)
lat = gregor_coords.lat.to(u.rad)
angle = np.degrees(np.arcsin(np.sin(np.radians(az_deg)) * np.cos(lat.value) / np.cos(dec.value)))
return angle
def calculate_slit_angle(rotangle=None, azimuth=None, elevation=None, p0_angle=None, date=None,
**kwargs):
"""
Calculate the slit orientation as described in
https://wolke7.leibniz-kis.de/f/1477681
Args:
rotangle: header value for derotator setting
azimuth: float degrees, azimuth in gregor configuration (south to west; S=0° W=90°)
elevation: float, degrees elevation/altitude angle
p0_angle: float, degrees p0 angle from gris ehader
date: datetime, time of observation
**kwargs: no effect
Returns: slit orientation_angle, degrees: clockwise angle between solar equator and slit
"""
# print(rotangle, parallactic_angle, azimuth, elevation, p0_angle, date)
parallactic_angle = get_parallactic_angle(date, azimuth)
const_offset = 50.1 if date < datetime.datetime(2019, 12, 31) else 45.67 # changes in optical setup
if rotangle:
alpha = 2 * rotangle - parallactic_angle - azimuth - elevation + p0_angle + 2 * const_offset
else:
alpha = -parallactic_angle - azimuth - elevation + p0_angle + 2 * const_offset
return alpha
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment