''' A main script to plot the FOV map, trajectory, and E-t (JDC) and C-t (JNA).
This is originally from a copy of :mod:`apps120816_flyover.flyover`.
Related scripts are:
- :mod:`apps120718_jna_fov.traj_on_surface`.
- :mod:`apps120803_jdcfly.app08_jdc_jna_flyover`.
.. image:: ../../../src/scripts/apps120816_flyover/flyover_yaw_70steer.png
:width: 90%
'''
import datetime
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
#import pyana.pep.juice_spice1205 as jspice
import pyana.juice.jspice as js
js.init()
import pyana.pep.pep_attitude
import pyana.jna.frame0 as frame
import pyana.jna.fov0 as jnafov
import pyana.jna.frame0 as jnaframe
import pyana.util.intersection
def pla2geo(vecs_pla):
vecs_geo = np.ma.zeros(vecs_pla.shape)
vecs_geo[0] = vecs_pla[1]
vecs_geo[1] = -vecs_pla[0]
vecs_geo[2] = vecs_pla[2]
return vecs_geo
def geo2pla(vecs_geo):
vecs_pla = np.ma.zeros(vecs_geo.shape)
vecs_pla[0] = -vecs_geo[1]
vecs_pla[1] = vecs_geo[0]
vecs_pla[2] = vecs_geo[2]
return vecs_pla
def vec2lonlat(vecs_geo):
lons = np.arctan2(vecs_geo[1], vecs_geo[0]) * 180. / np.pi
# lons is -180 to 180. Lons should range -90 to 270.
lons = np.array(lons)
lons = np.where(lons < -90, lons + 360, lons)
r = np.sqrt((vecs_geo ** 2).sum(0))
lats = np.arcsin(vecs_geo[2] / r) * 180. / np.pi
return lons, lats
def segmentize_lons(lons):
''' Return index boundary that separates the segment
Segment is separated by a boundary that has difference of
longitudal angle more than 180 degrees.
For example,
>>> lons = np.array([-175, -178, -179, 179, 177, 173, 178, 179, -178, -175])
>>> print(lons.shape)
(10,)
will be separated into three segments, namely first three negatives,
following five positives, and the last two negatives.
Returned is boundary index. In this example, [0, 3, 8, 10] will be returned.
The number of the segment is the length of the returned minus 1.
>>> idx = segmentize_lons(lons)
>>> print(idx)
[0, 3, 8, 10]
You can use the returned as
>>> nseg = len(idx) - 1
>>> segmented_lons = []
>>> for i in range(nseg):
... segmented_lons.append(lons[idx[i]:idx[i+1]])
>>> print(segmented_lons)
[array([-175, -178, -179]), array([179, 177, 173, 178, 179]), array([-178, -175])]
Obviously, the array as
>>> lons = np.arange(1000)
will return the boundary as [0, 1000], while it exceeds 360 degrees.
>>> print(segmentize_lons(lons))
[0, 1000]
'''
nlons = len(lons)
diff = np.abs(lons[:-1] - lons[1:])
segb = (np.where(diff > 180)[0] + 1).tolist()
segb.insert(0, 0)
segb.append(nlons)
return segb
def main(type="nadir", get_yaw=None, output=None):
'''Main script
get_yaw should be a function that eats the latitude_degrees
returning the steering angle for the specific time.
'''
t0 = datetime.datetime(2033, 4, 4, 18, 50)
### Plot on map
fig = plt.figure(figsize=(24, 10))
ax = fig.add_subplot(111)
### Note that the longitude here is
m = Basemap(projection='cyl', lon_0=90, resolution='c', ax=ax)
m.drawparallels(np.arange(-90, 120, 30), labels=[True, False, False, False])
m.drawmeridians(np.arange(-180, 210, 30), labels=[False, False, False, True])
m.warpimage('../../pyana/pep/images/ganymede.jpg')
### Map
import pyana.pep.mhddata
fmap = pyana.pep.mhddata.PrecipitationFlux(type=type)
lonlist = fmap.reshape(fmap.lonlist()) # In plasma frame
lonlist = lonlist - 90 # In geo-frame
latlist = fmap.reshape(fmap.latlist())
flxlist = fmap.reshape(fmap.flist())
lflxlist = np.ma.log10(flxlist)
x, y = m(lonlist, latlist)
# LOG scale
pc = m.contourf(x, y, lflxlist, np.linspace(7.9, 8.5, 7))
cb = fig.colorbar(pc)
cb.set_label('Log of flux [/cm2 s]')
from matplotlib import cm
pc2 = m.contour(x, y, lflxlist, np.linspace(7.5, 7.9, 5), cmap=cm.Blues)
plt.clabel(pc2)
# LINEAR scale
# pc = m.contourf(x, y, flxlist, [1e8, 1.25e8, 1.5e8, 1.75e8, 2e8, 3e8])
# fig.colorbar(pc)
### Open-close line
import pyana.pep.mhddata
lval = pyana.pep.mhddata.LValue1201()
ocb = lval.opencloseboundaries()
lon = np.arange(-90, 270.001, 2)
x, y = m(lon, ocb[0])
m.plot(x, y, 'y', lw=3)
x, y = m(lon, ocb[1])
m.plot(x, y, 'y', lw=3)
### Track.
for tlen, style in zip((1080, 3600), ('w-', 'w--')):
track_time = t0 + datetime.timedelta(seconds=10) * np.arange(tlen) # 3 or 10 hours
poslist = np.array([js.get_position(t, origin='GANYMEDE', frame='IAU_GANYMEDE') for t in track_time])
poslist = poslist.T
lon, lat = vec2lonlat(poslist)
segidx = segmentize_lons(lon)
seglon = [lon[segidx[i]:segidx[i+1]] for i in range(len(segidx) - 1)]
seglat = [lat[segidx[i]:segidx[i+1]] for i in range(len(segidx) - 1)]
for slon, slat in zip(seglon, seglat):
print(slon, slat)
x, y = m(slon, slat)
m.plot(x, y, style)
### Track every 10 minutes
tick_time = t0 + datetime.timedelta(minutes=10) * np.arange(60)
poslist = np.array([js.get_position(t, origin='GANYMEDE', frame='IAU_GANYMEDE') for t in tick_time])
poslist = poslist.T
lon, lat = vec2lonlat(poslist)
x, y = m(lon, lat)
m.plot(x, y, 'wo')
### JNA FoV on map
fovtlist = t0 + datetime.timedelta(minutes=2) * np.arange(40)# + datetime.timedelta(minutes=60)
### Need instance of sc
sc = pyana.pep.pep_attitude.NadirLookingYawSteeringSc()
rg = 2634.1
for fovt in fovtlist:
scpos_gany = js.get_position(fovt, origin='GANYMEDE', frame='IAU_GANYMEDE')
scvel_gany = js.get_velocity(fovt, origin='GANYMEDE', frame='IAU_GANYMEDE')
### For steering angle GUESS,
lon, lat = vec2lonlat(scpos_gany)
print(lon, lat)
yaw = get_yaw(lat)
print('Lat', lat)
print('Yaw', yaw)
sc.set_posvel(scpos_gany, scvel_gany, yaw)
for ch in range(0, 7):
thetalist = np.array([-3., 0, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, -3, -3, -3, -3, -3, -3, -3, -3, -3], dtype=np.float)
philist = -70 + 20 * ch + np.array([0., 0, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 20, 20, 17.5, 15, 12.5, 10, 7.5, 5, 2.5, 0], dtype=np.float)
jnacoords = np.array([frame.angles2jna(t, p) for t, p in zip(thetalist, philist)])
pepcoords = frame.jna2nsc(jnacoords.T).T # This is PEP (virtual NSC) frame.
nsc2ysc = sc.get_matrix_nsc2ysc() # Now, YSC frame, the actially NSC frame
sccoords = nsc2ysc.dot(pepcoords.T).T
nsc2ref = sc.get_matrix_nsc2ref()
print('nsc2ref', nsc2ref)
ganycoords = nsc2ref.dot(sccoords.T).T
xsect = np.ma.ones(ganycoords.shape) * np.nan
for idx, vec in enumerate(ganycoords):
v2 = pyana.util.intersection.los_limb_sphere(scpos_gany, vec, [0, 0, 0], rg)
if v2 != None:
xsect[idx, :] = np.array([v2.x, v2.y, v2.z])
xsect = np.ma.masked_invalid(xsect)
xsectn = xsect / rg
lat = np.arcsin(xsectn[:, 2].compressed()) * 180. / np.pi
lon = np.unwrap(np.arctan2(xsectn[:, 1].compressed(), xsectn[:, 0].compressed()))
print(lon)
lon = lon * 180. / np.pi
x, y = m(lon, lat)
m.plot(x, y, 'r')
if output == None:
fig.savefig('flyover_yaw_1.png')
else:
fig.savefig(output)
### Only for mapping now.
return
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
jnatime = t0 + datetime.timedelta(minutes=2) * np.arange(90)
jnacnt = np.zeros((len(jnatime), 7)) + np.nan
import pyana.jdc.fluxsimulator
jdcsim = pyana.jdc.fluxsimulator.SimpleSimulator()
jdcsim.set_species(16.0, 1.0)
pp = pyana.pep.mhddata.PlasmaParameter1205()
jdcflx = np.zeros((len(jnatime), 128)) + np.nan
for ilat, t in enumerate(jnatime):
print(ilat, t)
scpos_gany = js.get_position(t, origin='GANYMEDE', frame='IAU_GANYMEDE')
scvel_gany = js.get_velocity(t, origin='GANYMEDE', frame='IAU_GANYMEDE')
# Convert to plasma frame
scpos_gany = geo2pla(scpos_gany) / rg
scvel_gany = geo2pla(scvel_gany)
sc.set_posvel(scpos_gany, scvel_gany)
for ch in (1, 2, 3, 4, 5): # Only center 5 channels.
# Center direction
az = jnafov.azim_pix_center(ch)
el = jnafov.elev_pix_center(ch)
vecjna = jnaframe.angles2jna(el, az)
vecnsc = jnaframe.jna2nsc(vecjna)
vecgan = sc.convert_to_ref(vecnsc) # Looking vector in Ganymede coordinates.
v2 = pyana.util.intersection.los_sphere(scpos_gany, vecgan, [0, 0, 0], 1)
if v2 == None:
continue
longit = np.arctan2(v2[1], v2[0]) * 180. / np.pi
if longit < 0: longit += 360
latit = np.arcsin(v2[2]) * 180. / np.pi
print('LONLAT', longit, latit)
try:
idxlon, idxlat = fmap.nearest_neighbor(longit, latit)
except IndexError as e:
from pudb import set_trace; set_trace()
print('Strange, index error. Anyway, skip it.', longit, latit)
continue
jnacnt[ilat, ch] = flxlist[idxlon, idxlat] * 0.2 / (2 * np.pi) * 2e-6 # in /s. Indeed , count rate
# 0.2 comes from 20% from CENA result.
# Sensor g-factor is taken from proposal v0.3, 2e-6 cm2 sr/pixel.
# print jnacnt[ilat, ch]
### JDC
n, vx, vy, vz, temp, pth = pp.interpolate3d(scpos_gany[0], scpos_gany[1], scpos_gany[2])
vth = pyana.pep.mhddata.t2vth(temp, mass=16.) # in m/s
n *= 1e6 # in /m3
vx *= 1e3; vy *= 1e3; vz *= 1e3 # in m/s
jdcsim.set_spacecraft_pvat(scpos_gany, scvel_gany)
jdcsim.set_plasmaparameter(n, vx, vy, vz, vth)
for ie in range(128):
jdcflx[ilat, ie] = jdcsim.get_countrate(ie, 31, 0)
print(jdcflx[ilat, :])
# xax = [(j - jnatime[0]).seconds for j in jnatime]
import matplotlib.dates as mpldates
xax = [mpldates.date2num(j) for j in jnatime]
yax = np.arange(8)
x, y = np.meshgrid(xax, yax)
# img = ax2.pcolor(x, y, np.ma.log10(jnacnt).T)
img = ax2.pcolor(x, y, np.ma.masked_invalid(jnacnt).T, vmin=0, vmax=10)
ax2.set_ylabel('JNA: Direction')
cb = plt.colorbar(img, ax=ax2)
cb.set_label('JNA count rate (log)')
cb.set_label('JNA count rate')
ax2.xaxis.set_major_formatter(mpldates.DateFormatter('%H:%M'))
# locator = mpldates.AutoDateLocator()
# ax2.xaxis.set_major_locator(locator)
ax2.xaxis.set_ticks([mpldates.date2num(j) for j in tick_time])
ax2.set_xlim(xax[0], xax[-1])
import pyana.jdc.energy0 as energy
yax = energy.getEnergy()
x, y = np.meshgrid(xax, yax)
img = ax1.pcolor(x, y, np.ma.log10(jdcflx).T, vmin=0.1, vmax=4)
ax1.set_ylabel('JDC: Energy [eV]')
ax1.set_yscale('log')
ax1.set_ylim(yax.min(), yax.max())
cb = plt.colorbar(img, ax=ax1)
cb.set_label('JDC count rate (log)')
ax1.xaxis.set_ticks([mpldates.date2num(j) for j in tick_time])
ax1.xaxis.set_ticklabels('')
# ax1.xaxis.set_major_locator(locator)
ax1.set_xlim(xax[0], xax[-1])
fig.autofmt_xdate()
fig.savefig('flyover_yaw_2.png')
if __name__ == "__main__":
main(get_yaw=lambda x: 0, output="flyover_noyaw.png") # Constantly 70 deg.
main(get_yaw=lambda x: 12, output="flyover_yaw_12const.png") # Constantly 70 deg.
main(get_yaw=lambda x: 12 / 90. * x, output="flyover_yaw_12steer.png") # Constantly 70 deg.
main(get_yaw=lambda x: 28, output="flyover_yaw_28const.png") # Constantly 70 deg.
main(get_yaw=lambda x: 28 / 90. * x, output="flyover_yaw_28steer.png") # Constantly 70 deg.
main(get_yaw=lambda x: 70, output="flyover_yaw_70const.png") # Constantly 70 deg.
main(get_yaw=lambda x: 70 / 90. * x, output="flyover_yaw_70steer.png") # Constantly 70 deg.
# main(type='footprint')