Using EIT HEC information

gistfile1.txt

#
# Example EIT HEC to HGS conversion
#
import datetime

import numpy as np

import astropy.units as u
from astropy.coordinates import CartesianRepresentation, HeliocentricTrueEcliptic

from sunpy.coordinates import HeliographicStonyhurst, get_sunearth_distance
from sunpy.data.sample import EIT_195_IMAGE
from sunpy.time import parse_time
from sunpy.net import Fido
from sunpy.net import attrs as a
import sunpy.map

base_date = parse_time('2001-06-07')

n_step = 4

step_size = (1 * u.year).to(u.s).value / n_step
one_day = (1 * u.day).to(u.s).value

for i in range(0, n_step):

    # Search
    print(' ')
    print('Searching for example {:n} out of {:n}...'.format(i+1, n_step))
    start_date = base_date + datetime.timedelta(seconds=i*step_size)
    end_date = start_date + datetime.timedelta(seconds=one_day)
    q = Fido.search(a.Instrument('eit'), a.Time(start_date, end_date), a.Wavelength(195*u.AA))

    # Download one file
    print('Downloading...')
    f = Fido.fetch(q[0, 0])

    # Get the map
    m = sunpy.map.Map(f[0])

    hec_x = m.meta['hec_x'] * u.km
    hec_y = m.meta['hec_y'] * u.km
    hec_z = m.meta['hec_z'] * u.km

    cr = CartesianRepresentation(hec_x, hec_y, hec_z)

    hec_distance = np.sqrt(hec_x**2 + hec_y**2 + hec_z**2)

    hte = HeliocentricTrueEcliptic(cr)

    hgs_frame = HeliographicStonyhurst(obstime=m.date)

    hgs_coord = hte.transform_to(hgs_frame)

    # Image time
    print(' ')
    print('EIT 195 image time', m.date)

    #
    print(' ')
    print('[0] Coordinate Round trip')
    print('    (a) original Cartesian', cr)
    print('     ')
    print('    (b) distance', hec_distance)
    print('     ')
    print('    (c) heliocentric true ecliptic', hte)
    print('     ')
    print('    (d) transformed to HGS', hgs_coord)
    print('     ')
    print('    (e) back to heliocentric true ecliptic', hgs_coord.transform_to(HeliocentricTrueEcliptic))
    print(' ')

    # Should be around 0.99, since SOHO is about 1% closer to
    # the Sun than the Earth is
    print('[1] Ratio HEC distance / Sun-Earth distance at ', m.date)
    print('    should be around 0.99')
    gsed = (get_sunearth_distance(m.date)).to(u.km)
    print('    HEC distance / Sun-Earth distance', hec_distance/gsed)
    print('    HGS radius / Sun-Earth distance', hgs_coord.radius/gsed)
    print(' ')

    # Should be very very close to each other, since the two
    # co-ordinate systems represent the same co-ordinate in space
    print('[2] HEC distance and HGS Sun-Earth distance should be almost identical')
    print('   ', hec_distance, hgs_coord.radius)
    print('    HEC distance - HGS radius = ', hec_distance - hgs_coord.radius)
    print(' ')

    # Should be close to zero since SOHO sits close to the Earth/Sun line
    print('Longitude and Latitude from both sources should be close to each other')
    print('Longitudes should be close to zero since SOHO is close to Earth/Sun line')
    print('Latitudes should follow annual variation within the limits of the solar B0 angle range')
    print('[3] HGS Longitude', hgs_coord.lon.to(u.deg).value)
    print('    Longitude from FITS file', m.heliographic_longitude)
    print(' ')
    print('[4] HGS Latitude', hgs_coord.lat.to(u.deg).value)
    print('    latitude from FITS file', m.heliographic_latitude)