''' Solar wind related module.
'''
import numpy as np
import datetime
import spiceypy as spice
import irfpy.mexpvat.mexspice as ms
ms.init('full')
import irfpy.mima.fov as ifov
[docs]def aberrate(t, v_mso, frame='MSO'):
''' Calculate the aberrated solar wind velocity vector.
:param t: Time
:param v_mso: Non-aberrated solar wind velocity vector
in the MSO frame.
:keyword frame: The frame.
Originally developed in :mod:`appl_130717_swparam_statistics.solarwind_direction`.
Now ported with a slihgt modification of interface..
The MSO/non-aberrated solar wind is converted to the aberrated solar wind velocity
in any frame. The aberration of the spacecraft velocity is calculated inside this
function by subtracting the spacecraft velocity in an inertia frame (J2000 in this case).
Therefore, the algorithm is
1. The non-aberrated SW velocity (``v_mso``) is converted to the J2000 inertial frame,
in which frame the solar wind flow is "fixed".
2. The non-aberrated SW velocity in J2000 frame is aberrated by considering the spacecraft
velocity in J2000 frame.
3. The aberrated SW velocity (J2000) is converted to the frame given by the parameter ``frame``.
>>> t = datetime.datetime(2007, 12, 10)
>>> vsw_mso_noab = np.array([-400, 0, 0])
>>> if not ms.isdb():
... print("SPICE not loaded. Check MEX/SPICE module (mexspice)") #doctest: +SKIP
... else:
... print('{0[0]:.2f} {0[1]:.2f} {0[2]:.2f}'.format(aberrate(t, vsw_mso_noab)))
-398.11 22.19 -0.31
'''
et = ms.spice.str2et(t.strftime('%FT%T.%f'))
### First, covert to J2000 frame.
M = np.array(ms.spice.pxform('MSO', 'J2000', et))
Vj2000_noab = M.dot(v_mso)
### Then, calculate the MEX velocity wrt SUN in J2000
vsc = ms.get_posvel(t, target='MEX', origin='SUN', correction='None', frame='J2000')[3:]
### This is used to aberrate the solar wind velocity
Vj2000_ab = Vj2000_noab - vsc
### Convert from J2000 to given frame
M=np.array(ms.spice.pxform('J2000', frame, et))
return M.dot(Vj2000_ab)
[docs]def vsw(swvel, theta_deg, phi_deg):
''' Calculate the solar wind velocity vector.
The definition of th_d and ph_d is slightly different from the usual practice.
The theta is defined from the ecliptic plane of Mars.
The phi is defined from the ``-x``, namely the solar wind flowing direction.
``+y`` direction is positive for phi.
Thus, the conversion is
``[-swvel * cos t * cos p, swvel * cos t * sin p, swvel * sin t]``
:param swvel: Solar wind velocity (positive)
:param theta_deg: Angle of theta.
:param phi_deg: Angle of phi.
>>> print(vsw(400, 0, 0))
[-400. 0. 0.]
>>> print(vsw(400, 10, 0)) # 10 degree "north"
[-393.9231012 0. 69.45927107]
>>> print(vsw(400, 0, 10)) # 10 degree "+y"
[-393.9231012 69.45927107 0. ]
'''
t, p = np.deg2rad([theta_deg, phi_deg])
st, sp = np.sin([t, p])
ct, cp = np.cos([t, p])
vmso_noab = np.array([-swvel * ct * cp, swvel * ct * sp, swvel * st])
return vmso_noab
[docs]def direction_imachannel(t, swvel=400., theta_deg=0., phi_deg=0.):
''' Return the direction of SW in IMA channel.
>>> t0 = datetime.datetime(2010, 9, 12, 5, 20)
>>> if not ms.isdb():
... print("SPICE not loaded. Check MEX/SPICE module (mexspice)") #doctest: +SKIP
... else:
... print('{0[0]:.2f} {0[1]:.2f} {0[2]:.2f}'.format(direction_imachannel(t0, swvel=420)))
422.13 1.20 7.46
'''
vsw_mso_noab = vsw(swvel, theta_deg, phi_deg)
vsw_imas_ab = aberrate(t, vsw_mso_noab, frame='MEX_ASPERA_IMAS')
view_imas_ab = -vsw_imas_ab # Viewing direction from IMAS
fov = ifov.SimpleFOV()
az, el = fov.numbering_imasframe(view_imas_ab)
return np.linalg.norm(vsw_imas_ab), az, el
[docs]def direction_imachannel_noab(t, swvel=400., theta_deg=0., phi_deg=0.):
''' Return the direction of SW in IMA channel.
>>> t0 = datetime.datetime(2010, 9, 12, 5, 20)
>>> if not ms.isdb():
... print("SPICE not loaded. Check MEX/SPICE module (mexspice)") #doctest: +SKIP
... else:
... print('{0[0]:.2f} {0[1]:.2f} {0[2]:.2f}'.format(direction_imachannel_noab(t0, swvel=420)))
420.00 1.35 7.50
'''
vsw_mso_noab = vsw(swvel, theta_deg, phi_deg)
et = spice.str2et(t.strftime('%FT%T.%f'))
matx_mso2imas = spice.pxform('MSO', 'MEX_ASPERA_IMAS', et)
vsw_imas_noab = np.dot(matx_mso2imas, vsw_mso_noab)
fov = ifov.SimpleFOV()
az, el = fov.numbering_imasframe(-vsw_imas_noab)
return np.linalg.norm(vsw_imas_noab), az, el
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')