calculations.py 5.05 KB
Newer Older
1
2
3
4
5
6
7
8
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Created by schaffer at 8/22/19

"""

9
import numpy as np
10
import sunpy.sun.constants as sun_const
11
from astropy import units as u
12
from astropy.coordinates import SkyCoord, Angle
13
from astropy.wcs import WCS
14
15
16
17
from sunpy.coordinates import Helioprojective, frames, sun


def get_solar_radius(date):
Carl Schaffer's avatar
Carl Schaffer committed
18
    """Get apparent solar radius for a given date
19
20
21
22
23
24
25
26
27
28
29

    Args:
        date: observation date string

    Returns:
        r: apparent solar radius in arcsec
    """
    r = sun.angular_radius(date).to(u.arcsec)
    return r.value


Carl Schaffer's avatar
Carl Schaffer committed
30
31
def calc_heliocentric_angle(tx, ty, date="", rsun_ang=945.57):
    """Calculate heliocentric angle from helioprojective cartesian coordinates theta_x and theta_y
32
33
34
35
        When the date is given SunPy transformation to heliocentric-radial (cylindrical representation) is made
        to get rho radius (in AU by default) and then mu/theta using formula Eq.8 from Thompson 2006 paper.
        If date is not given, use apparent solar radius rsun and formulas provided by Markus Schmassmann.
        Difference between two methods is approx. 20 arcsec near the limb (theta > 85 deg).
Carl Schaffer's avatar
Carl Schaffer committed
36

37
38
39
40
41
    Args:
        tx: helioprojective westaward angle in arcsec
        ty: helioprojective northward angle in arcsec
        date: observation date string
        rsun_ang: apparent solar radius in arcsec for selected date
Carl Schaffer's avatar
Carl Schaffer committed
42

43
    Returns:
Taras Yakobchuk's avatar
Taras Yakobchuk committed
44
        theta, mu: heliocentric angle in degrees and its cosine
45
    """
Carl Schaffer's avatar
Carl Schaffer committed
46
47
48
49
50
51
52
53
    if date != "":
        hp_crd = SkyCoord(
            tx * u.arcsec,
            ty * u.arcsec,
            frame=frames.Helioprojective,
            observer="earth",
            obstime=date,
        )
54
        hc_crd = hp_crd.transform_to(frames.Heliocentric)
Carl Schaffer's avatar
Carl Schaffer committed
55
        hc_crd.representation_type = "cylindrical"
56
57
58
59
60
61
62
63
64
65
66
67
        hc_rho = hc_crd.rho.to(u.meter).value
        R_SUN = sun_const.radius.value
        theta = np.arcsin(hc_rho / R_SUN)
        mu = np.cos(theta)
        theta = (theta * u.radian).to(u.degree).value
    else:
        hp_r = (np.sqrt(tx ** 2 + ty ** 2) * u.arcsec).to(u.radian).value
        r_sun = (rsun_ang * u.arcsec).to(u.radian).value
        theta = np.arcsin(np.sin(hp_r) / np.sin(r_sun)) - hp_r
        mu = np.cos(theta)
        theta = (theta * u.radian).to(u.degree).value
    return theta, mu
68

69

70
71
72
73
74
def to_heliographic(coords, date):
    """
    Transform a coordinate tuple from Helioprojective Cartesian to Heliographic coordinates

    """
75
    tx, ty = coords
Carl Schaffer's avatar
Carl Schaffer committed
76
77
78
    coords_in = Helioprojective(
        tx * u.arcsec, ty * u.arcsec, obstime=date, observer="Earth"
    )
79
    output = coords_in.transform_to(frames.HeliographicStonyhurst)
80
81
    lon, lat = output.lon.degree, output.lat.degree
    return lon, lat
82
83


84
85
def calc_wcs_bbox(header):
    """Calculate WCS bounding box of the image
86
87
88
89
90
91
92
93
94

    Args:
      header: FITS header

    Returns:
      NumPy array of dimensionless min and max values of x and y
      with shape [[xmin, ymin], [xmax, ymax]]

    """
95
96
97
    wcs = WCS(header, naxis=2)
    naxis1 = wcs.pixel_shape[0]
    naxis2 = wcs.pixel_shape[1]
Carl Schaffer's avatar
Carl Schaffer committed
98
99
100
101
102
103
    rect_pixel = np.array(
        [[1, 1], [1, naxis2], [naxis1, naxis2], [naxis1, 1]], dtype=np.float64
    )
    rect_deg = Angle(wcs.pixel_to_world_values(rect_pixel) * u.degree).wrap_at(
        180 * u.degree
    )
104
105
    rect_arcsec = rect_deg.to(u.arcsec).value

Carl Schaffer's avatar
Carl Schaffer committed
106
107
108
109
    return [
        [min(rect_arcsec[:, 0]), min(rect_arcsec[:, 1])],
        [max(rect_arcsec[:, 0]), max(rect_arcsec[:, 1])],
    ]
Carl Schaffer's avatar
Carl Schaffer committed
110
111


112
def rotate_around_first_coord(X, Y, angle_rad):
113
114
115
116
117
    """
    Rotate a set of coordinates by a given angle (in radians). The rotation is performed around (X[0], Y[0])
    Args:
        X: numpy array of x values
        Y: numpy array of y values
118
        angle_rad: rotational angle in radians
119
120
121
122

    Returns:
        coords: numpy array of shape (len(X),2) with rotation results
    """
Carl Schaffer's avatar
Carl Schaffer committed
123
124
125
    x_cent = X - X[0]
    y_cent = Y - Y[0]

126
    c, s = np.cos(angle_rad), np.sin(angle_rad)
Carl Schaffer's avatar
Carl Schaffer committed
127
128
129
130
131
    R = np.array(((c, -s), (s, c)))
    x_rot, y_rot = np.dot(R, [x_cent, y_cent])
    x_rot += X[0]
    y_rot += Y[0]

132
133
    return x_rot, y_rot

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

def rotate_around_origin(xy, angle_rad):
    """
    Rotate a point clockwise around the origin (0, 0).
    Args:
        xy: tuple of floats (x, y)
        angle_rad: angle of rotation in radians

    Returns:
        (x_rotated, y_rotated)
    """
    x, y = xy
    xx = x * np.cos(angle_rad) + y * np.sin(angle_rad)
    yy = -x * np.sin(angle_rad) + y * np.cos(angle_rad)

149
    return xx, yy
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167


def estimate_telescope_drift(hours_since_center):
    """
    Gives an estimate of the drift of the telescope centering information. The parametrization of the degradation is
    based on studies performed by Olivier Grassin and backed up by a validation by Carl Schaffer based on ~500 cross
    correlations between GREGOR data and HMI. Further information can be found in the kis_tools repository under
    GRIS coordinate study
    Args:
        hours_since_center: hours since the centering procedure was done

    Returns:
        degradation: float mean drift in arcseconds

    """
    offset, slope = 23.4, 21.7
    degradation = offset + slope * hours_since_center
    return degradation