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

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

"""

import math

11
import numpy as np
12
from astropy import units as u
13
from astropy.wcs import WCS
14
15
from sunpy.coordinates import Helioprojective, frames

16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

def calc_theta(x, y, r_ref=945.57):
    """Calculate the value of the heliocentric angle theta

    Args:
      x: helio-projective cartesian x posn in arcsec
      y: helio-projective cartesian y posn in arcsec
      r_ref: apparent angular disk radiusin arcsun

    Returns:
      theta: float

    """
    r_obs = math.sqrt(math.pow(x, 2) + math.pow(y, 2))
    radius_ratio = min([r_obs / r_ref, 1])
    theta_value = math.asin(radius_ratio) / math.pi * 180
    return theta_value


def calc_mu(x, y, r_ref=945.57):
    """Calculate the value of mu (cos of heliocentric angle)

    Args:
      x: helio-projective cartesian x posn in arcsec
      y: helio-projective cartesian y posn in arcsec
      r_ref: apparent angular disk radiusin arcsun

    Returns:
      mu: float, range 1 to 0; 1 means disk center

    """
    mu = math.cos(calc_theta(x, y, r_ref))
    return mu
49
50
51
52
53
54
55
56


def to_heliographic(coords, date):
    """
    Transform a coordinate tuple from Helioprojective Cartesian to Heliographic coordinates

    """
    x, y = coords
57
    input = Helioprojective(x * u.arcsec, y * u.arcsec, obstime=date, observer='Earth')
58
59
60
    output = input.transform_to(frames.HeliographicStonyhurst)
    lon, lat = output.lon.degree, output.lat.degree
    return lon, lat
61
62
63


def calc_wcs_footprint(header):
64
    """Calculate the WCS footprint of the image
65
66
67
68
69
70
71
72
73
74
75

    Args:
      header: FITS header

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

    """
    wcs = WCS(header)
    pixel_shape = wcs.pixel_shape
76
    n_x = pixel_shape[0]
77
    n_y = pixel_shape[1]
78
    naxis = len(pixel_shape)
79

80
    slit_rect = np.array([[0, n_x - 1] + [0] * (naxis - 2), [0, n_y - 1] + [0] * (naxis - 2)])
81
82
83
84
85
86
87
88
    slit_world = np.array(wcs.pixel_to_world_values(slit_rect))
    slit_deg = u.degree * slit_world
    slit_arcsec = slit_deg.to(u.arcsec)

    return np.array([[min(slit_arcsec[0, 0], slit_arcsec[1, 0]).value,
                      min(slit_arcsec[0, 1], slit_arcsec[1, 1]).value],
                     [max(slit_arcsec[0, 0], slit_arcsec[1, 0]).value,
                      max(slit_arcsec[0, 1], slit_arcsec[1, 1]).value]])