''' WIND module
WIND is earth L1 point spacecraft measuring the solar wind parameters.
Supporse WIND observed the solar wind velocity ``(vx0, vy0, vz0)``
in the GSE frame at t=t0. You may convert the vector in any frame
via spice, to be ``(vx, vy, vz)``.
.. todo::
Detailed explanation of conversion to Mars should be done.
If the WIND is at (x0, y0, z0) at t0 in IAU_SUN and
Mars, for example, at t0 (X0, Y0, Z0).
The reason of using IAU_SUN is that the best angular separation
is addressed to the Sun rotation.
Assume the polar coordinates (r:distance, ph:longitude, th:latitude)
and then the WIND position (r0, ph0, th0) at t0 and
the Mars position (R0, Ph0, Th0) at t0.
The solar wind velocity vector to be (vx0, vy0, vz0) in IAU_SUN at t0.
.. note::
For all the calculation, SPICE correction should be ``none``.
No light time correction is suitable.
**(To be written)**
*Data file generation*
Download all the needed dataset.
.. code-block:: sh
for yr in {2004..2013}; do wget -nc -np -r ftp://space.mit.edu/pub/plasma/wind/kp_files/$yr; done
ln -s space.mit.edu/pub/plasma/wind/kp_files/ raw
Then, start calculating the data.
.. code-block:: sh
for yr in {2004..2013}; do cd raw/$yr; python -c 'from irfpy.asperacommon import wind; wind.run()'; cd ../..; done
You may find ``.mars`` file in the directories.
The directory should be pointed by [asperacommon] windatmarsbase entry in .irfpyrc
.. todo::
MAG data propagation in a coherent way.
.. todo::
Extend this module to Venus in coherent way. It is found not very simple...
Also, extension to other planets may be useful for future.
In this sense, you may outdate this module, and put the module to another project or sub-package.
'''
import sys
import logging
import datetime
from io import StringIO
import numpy as np
import irfpy.mexpvat.mexspice as ms
ms.init_noatt() # Attitude not needed.
[docs]def vec2rtp(vec):
''' Return r, theta(latitude) and phi in radians
'''
r = np.sqrt((np.array(vec)[:3] ** 2).sum())
v = np.array(vec[:3]) / r
# print v
t = np.arcsin(v[2])
p = np.arctan2(v[1], v[0])
return np.array([r, t, p])
[docs]def rtp2vec(r, t, p):
x = r * np.cos(p) * np.cos(t)
y = r * np.sin(p) * np.cos(t)
z = r * np.sin(t)
[docs]def showvec(vec, fp=sys.stderr):
print('VEC= %g %g %g' % tuple(vec), file=fp)
r, t, p = vec2rtp(vec)
print(' %g (%f %f)' % (r, np.rad2deg(t), np.rad2deg(p)), file=fp)
[docs]def approach(x_ref, x_tar, v_sw, w_lon=2.666253e-6, logfile=None):
''' To calculate the WIND observation toward Mars.
Suppose the solar wind velocity vector ``vsw=(vx, vy, vz)``
measured at x_ref, the radial and the longitudal propagations
will make the solar wind velocity vector near to the
target postion, ``x_tar``. Not exactly matches, becuse
``x_tar`` will be changed.
This function will return the tuple, stating
(t_diff, x_ref_new), where t_diff is the propagation time
and x_ref_new is the new position of the reference point
(closest to x_tar at t=0).
The coordinate system should be IAU_SUN.
This is because the longitudal propagation we assume the
Sun rotation.
>>> if not ms.isdb(): # Check if SPICE is initialized
... print("SPICE is not loaded. Check if MEX/SPICE module is enabled") # doctest: +SKIP
... else:
... xr = [4, 1, 2]
... xt = [-1, 5, 1]
... vsw = [1, 1, 0.1]
... t, xr2, vsw = approach(xr, xt, vsw)
... print('%.3f' % t)
541313.542
... print '%.3f %.3f %.3f' % tuple(xr2)
-0.936 4.681 2.052
If the y location is behind, minus will be returned.
>>> if not ms.isdb(): # Check if SPICE is initialized
... print("SPICE is not loaded. Check if MEX/SPICE module is enabled") # doctest: +SKIP
... else:
... xr = [4, 1, 2]
... xt = [-1, -5, 1]
... vsw = [1, -1, 0.1]
... t, xr2, vsw = approach(xr, xt, vsw)
... print('%.3f' % t)
-682286.883
... print '%.3f %.3f %.3f' % tuple(xr2)
-0.934 -4.671 2.076
'''
logstr = StringIO()
xref = np.array(x_ref)[:3]
xtar = np.array(x_tar)[:3]
vsw = np.array(v_sw)[:3]
logstr.write('### START\n')
logstr.write('#Reference_point ')
logstr.write('%.2f %.2f %.2f\n' % tuple(xref))
logstr.write('#Target_point ')
logstr.write('%.2f %.2f %.2f\n' % tuple(xtar))
logstr.write('#SWvel ')
logstr.write('%.2f %.2f %.2f\n' % tuple(vsw))
rref, tref, pref = vec2rtp(xref)
rtar, ttar, ptar = vec2rtp(xtar)
vrsw, vtsw, vpsw = vec2rtp(vsw)
xref_vsw = np.dot(xref, vsw)
rref2 = rref ** 2
rtar2 = rtar ** 2
vrsw2 = vrsw ** 2
### Radial direction propagation.
t_radial = (- xref_vsw + np.sqrt(xref_vsw **2 - vrsw2 * (rref2 - rtar2))) / vrsw2
logging.debug('Radial time: {}'.format(t_radial))
logstr.write('#RadialTime ')
logstr.write('%.8f\n' % t_radial)
xref = xref + t_radial * vsw
rref, tref, pref = vec2rtp(xref) # Update the coordinate
logging.debug('xref={} vsw={} xtar={} xref.vsw={}'.format(xref, vsw, xtar, np.dot(xref, vsw)))
logstr.write('#RadialPoint ')
logstr.write('%.2f %.2f %.2f\n' % tuple(xref))
### Longitudal propgataion
dlon = ptar - pref
while dlon > np.pi: dlon -= (2 * np.pi)
while dlon < -np.pi: dlon += (2 * np.pi)
t_longitude = dlon / w_lon
logging.debug('dlon={} dlon(deg)={} t={} t(day)={}'.format(dlon, np.rad2deg(dlon), t_longitude, t_longitude / 86400.))
### Longitudal rotation
matx = np.array([[np.cos(dlon), -np.sin(dlon), 0],
[np.sin(dlon), np.cos(dlon), 0],
[0, 0, 1]])
xref = matx.dot(xref)
vsw = matx.dot(vsw)
logstr.write('#LongitudeTime ')
logstr.write('%.8f\n' % t_longitude)
logstr.write('#LongitudePoint ')
logstr.write('%.2f %.2f %.2f\n' % tuple(xref))
logstr.write('#SW_newframe ')
logstr.write('%.2f %.2f %.2f\n' % tuple(vsw))
logging.debug('xref={} vsw={} xtar={} xref.vsw={}'.format(xref, vsw, xtar, np.dot(xref, vsw)))
logstr.write('### END\n')
if logfile is not None:
lf = open(logfile, 'a')
lf.write(logstr.getvalue())
lf.close()
return t_radial + t_longitude, xref, vsw
[docs]def wind2mars(t0, wind_gse, v_sw_gse, tlimsec=1, logfile=None):
''' WIND observation converted to the Mars position.
Use of WIND, not ACE, as abberation is included for ACE.
>>> if not ms.isdb(): # Check if SPICE is initialized
... print("SPICE is not loaded. Check if MEX/SPICE module is enabled") # doctest: +SKIP
... else:
... t0 = datetime.datetime(2009, 3, 1, 12, 0, 0)
... wind = [1257103.8, -102048., -127687.56]
... vsw = [-437.6, -0.5192, -6.457]
... tm, xm, distrate = wind2mars(t0, wind, vsw)
... print(tm.strftime('%FT%T'))
2009-03-14T13:40:32
... print '{:.5e} {:.5e} {:.5e}'.format(xm[0], xm[1], xm[2])
4.33638e+06 2.78700e+06 -4.21764e+07
'''
logstr = StringIO()
logstr.write('### START\n')
# t0 = datetime.datetime(2009, 3, 1, 12, 0, 0) # DOY=60
# ### Data of WIND from rdP2009060.01
# winddat = '2009 60 0.49996424 437.7 -0.06798 -0.8452 2.737 28.84 -1.000e+31 -437.6 -0.5192 -6.457 197.1 -16.00 -20.02'.split()
#
# v_sw_gse = np.array([float(i) for i in winddat[-6:-3]])
# wind_gse = np.array([float(i) for i in winddat[-3:]]) * 6378.
logging.debug(str([v_sw_gse, wind_gse]))
et0 = ms.spice.str2et(t0.strftime('%FT%T.%f'))
#### GSE -> Sun-centered GSE.
## Sun in GSE.
suningse = np.array(ms.spice.spkezr('Sun', et0, 'GSE', 'NONE', 'Earth')[0])[:3]
## Earth in SCGSE
earth_insgse = np.array(ms.spice.spkezr('Earth', et0, 'GSE', 'NONE', 'Sun')[0])[:3]
## WIND in SCGSE
wind_insgse = - suningse + wind_gse
dist0sq = (wind_insgse ** 2).sum()
rotmatx = np.array(ms.spice.pxform('GSE', 'IAU_SUN', et0))
## WIND/Earth in IAU_SUN
wind_iausun = rotmatx.dot(wind_insgse)
earth_iausun = rotmatx.dot(earth_insgse)
earth_iausun_spice = np.array(ms.spice.spkezr('Earth', et0, 'IAU_SUN', 'NONE', 'Sun')[0])[:3]
vsw_iausun = rotmatx.dot(v_sw_gse)
mars_iausun = np.array(ms.spice.spkezr('Mars', et0, 'IAU_SUN', 'NONE', 'Sun')[0])[:3]
## For debugging, J2000 positions are needed
rotmatxj2000 = np.array(ms.spice.pxform('IAU_SUN', 'ECLIPJ2000', et0))
wind_j2000 = rotmatxj2000.dot(wind_iausun)
mars_j2000 = rotmatxj2000.dot(mars_iausun)
vsw_j2000 = rotmatxj2000.dot(vsw_iausun)
logstr.write('START 0\n')
logstr.write('REF0 %g %g %g\n' % tuple(wind_j2000))
logstr.write('TAR0 %g %g %g\n' % tuple(mars_j2000))
logstr.write('VEL0 %g %g %g\n' % tuple(vsw_j2000))
logstr.write('END 0\n')
logging.debug(str(['*** Original ***']))
logging.debug(str(['t0=', t0]))
logging.debug(str(['Position of WIND [1.5e8 km ~au]']))
# showvec(wind_iausun / 1.5e8)
logging.debug(str(['Position of Mars [1.5e8 km ~au]']))
# showvec(mars_iausun / 1.5e8)
logging.debug(str(['Vsw [km/s]']))
# showvec(vsw_iausun)
t_t0 = t0
et_t0 = et0
x_t0 = wind_iausun
vsw_t0 = vsw_iausun
mars_t0 = mars_iausun
for i in range(1, 11):
logging.debug(str(['*** %d-th iteration***' % i]))
t_prop, x_prop, vsw_prop = approach(x_t0, mars_t0, vsw_t0, logfile=logfile)
### x_prop is the reference point at t1, but the coordinate IAUSUN is defined t0.
### Thus, it is first converted to any inertia frame at t0.
t_t1 = t_t0 + datetime.timedelta(seconds=t_prop)
et1 = ms.spice.str2et(t_t1.strftime('%FT%T.%f'))
logging.debug(str(['*** after %d-th iteration ***' % i]))
logging.debug(str(['Propagation time', t_prop, 's = ', (t_prop / 3600), 'h = ', (t_prop / 86400), 'd']))
logging.debug(str(['t1 = ', t_t1]))
logging.debug(str(['Propagation point [1.5e8 km ~au] in t0-IAUSUN']))
# showvec(x_prop / 1.5e8)
logging.debug(str(['Propagation point [1.5e8 km ~au] in ECLIPJ2000']))
rotmatx = np.array(ms.spice.pxform('IAU_SUN', 'ECLIPJ2000', et_t0))
x_prop_j2000 = rotmatx.dot(x_prop)
# showvec(x_prop_j2000 / 1.5e8)
vsw_prop_j2000 = rotmatx.dot(vsw_prop)
logging.debug(str(['Mars@t0 [1.5e8 km ~au] in ECLIPJ2000']))
mars_j2000 = np.array(ms.spice.spkezr('Mars', et_t0, 'ECLIPJ2000', 'NONE', 'Sun')[0])[:3]
# showvec(mars_j2000 / 1.5e8)
logging.debug(str(['Mars@t1 [1.5e8 km ~au] in ECLIPJ2000']))
mars_j2000 = np.array(ms.spice.spkezr('Mars', et1, 'ECLIPJ2000', 'NONE', 'Sun')[0])[:3]
# showvec(mars_j2000 / 1.5e8)
mars_iausun = np.array(ms.spice.spkezr('Mars', et1, 'IAU_SUN', 'NONE', 'Sun')[0])[:3]
rotmatx = np.array(ms.spice.pxform('ECLIPJ2000', 'IAU_sun', et1))
x_prop = rotmatx.dot(x_prop_j2000)
vsw_prop = rotmatx.dot(vsw_prop_j2000)
logging.debug(str(['Position of Mars @t1 [1.5e8 km ~au] t1-IAUSUN']))
# showvec(mars_iausun / 1.5e8)
logging.debug(str(['Position of refpoint [1.5e8 km ~au]']))
# showvec(x_prop / 1.5e8)
wind_j2000 = x_prop_j2000
mars_j2000 = mars_j2000
vsw_j2000 = vsw_prop_j2000
logstr.write('START %d\n' % i)
logstr.write('REF%d %g %g %g\n' % (i, wind_j2000[0], wind_j2000[1], wind_j2000[2]))
logstr.write('TAR%d %g %g %g\n' % (i, mars_j2000[0], mars_j2000[1], mars_j2000[2]))
logstr.write('VEL%d %g %g %g\n' % (i, vsw_j2000[0], vsw_j2000[1], vsw_j2000[2]))
logstr.write('END %d\n' % i)
if t_prop <= tlimsec:
rotmatx = np.array(ms.spice.pxform('IAU_SUN', 'MSO', et1))
marspos = np.array(ms.spice.spkezr('Mars', et1, 'MSO', 'NONE', 'Sun')[0])[:3]
logstr.write('### END\n')
if logfile is not None:
with open(logfile, 'w') as fp:
fp.write(logstr.getvalue())
dist1sq = (x_prop ** 2).sum()
return t_t1, rotmatx.dot(x_prop) - marspos, dist1sq / dist0sq
t_t0 = t_t1
et_t0 = et1
x_t0 = x_prop
vsw_t0 = vsw_prop
mars_t0 = mars_iausun
raise RuntimeError('Too many iteration. Stop calculating.')
[docs]def wind2mars_fromfile(filename, outfile):
fp2 = open(outfile, 'w')
with open(filename) as fp:
for l in fp.readlines():
winddat = l.split()
if len(winddat) !=15:
continue
yr = int(winddat[0])
doi = int(winddat[1])
tod = float(winddat[2]) * 86400 # seconds
obsdat = datetime.datetime(yr, 1, 1)
obsdat = obsdat + datetime.timedelta(days=doi - 1)
obsdat = obsdat + datetime.timedelta(seconds = tod)
v_sw_gse = np.array([float(i) for i in winddat[-6:-3]])
wind_gse = np.array([float(i) for i in winddat[-3:]]) * 6378.
try:
tmars, pos, r1r0_sq = wind2mars(obsdat, wind_gse, v_sw_gse)
pos = np.sqrt((pos * pos).sum()) / 3397.
except RuntimeError:
tmars, pos = datetime.datetime.max, -1e20
print(l.strip(), obsdat.strftime('%FT%T.%f'), tmars.strftime('%FT%T.%f'), r1r0_sq, pos, file=fp2)
[docs]def run():
import glob
fnams = glob.glob('./rdP???????.??')
for f in fnams:
print(f)
wind2mars_fromfile(f, f + '.mars')
[docs]class Wind:
''' WIND data base.
Now, the WIND data is mass-prepared by :func:`run` function,
and the ``*.mars`` files should be placed under the directory
specified by [asperacommon]-windatmarsbase.
From the directory below, the data is stored into a sqlite database.
The database name is {``winddatamarsbase``}/``wind.db``
'''
filetabledefinition = """FileTable
(Filename TEXT,
DataDate TEXT,
DataVersion INTEGER,
DataParseDate TEXT,
FileCreationDate TEXT)"""
datatabledefinition = """SwParamTable
(Filename TEXT,
Time TEXT,
TimeMars TEXT,
ErrorMars REAL,
Velocity REAL,
Density REAL,
Temperature REAL)"""