''' Calculate precipitation flux from JIAMHD data
Calculate precipitation following the scheme below.
Iterate over the surface
1. Specify the point at surface.
2. Trace back to the point at 1.2 Rg along B field.
3. Obtain the plasma parameters (n, vxyz, vth, Bxyz) at the point.
4. Calculate the velocity vector along magnetic field.
5. Calculate the precipitation flux along the magnetic field.
The output is seven figures.
The most important one is the precipitating flux map
.. image:: ../../../src/scripts/apps120811_flux_recalc/precipitation_1.png
'''
import os
import sys
import logging
logging.basicConfig()
logger = logging.getLogger('project_at_500km.py')
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
from pyana.util import streamline
from pyana.util import maxwell
from pyana.pep import mhddata
def derivepoint(p0, p1, r=1.2):
r''' Derive the point where the distance is ``r`` between p0 and p1
p0 and p1 should be a np.array with a shape of (3,).
The equation is simple enough.
.. math::
(p_1 - p_0)^2 k^2 + 2 p_0 (p_1 - p_0) k + p_0^2 - r^2 = 0
The returned k (two roots) is selected to be between 0 and 1.
Then, you can get the point :math:`p = p_0 + k (p_1 - p_0)`.
'''
a = ((p1 - p0) ** 2).sum()
b = ((p1 - p0) * p0).sum()
c = (p0 ** 2).sum() - r ** 2
k1 = (-b + np.sqrt(b * b - a * c)) / a
k2 = (-b - np.sqrt(b * b - a * c)) / a
k = np.array([k1, k2])
k = np.ma.masked_outside(k, 0, 1).compressed()
if k.size == 1:
return p0 + k[0] * (p1 - p0)
elif k.size == 2:
return p0 + k[1] * (p1 - p0)
else:
return np.zeros([3]) + np.nan
def traceto(lond, latd, mf, r=1.2):
if r <= 1.0:
raise IndexError('Given r (=%f) should be >1')
r2 = r * r
### Bfield function
def bfield(pos):
''' Return the bfield direction in x-z plane
'''
if (pos[0] ** 2 + pos[1] ** 2 + pos[2] ** 2).sum() < .999:
raise IndexError('Inside sphere')
bx, by, bz = mf.interpolate3d(pos[0], pos[1], pos[2])
babs = np.sqrt(bx * bx + by * by + bz * bz)
return np.array((bx / babs, by / babs, bz / babs))
### Position at surface
lonr = np.deg2rad(lond)
latr = np.deg2rad(latd)
pos = np.array([np.cos(latr) * np.cos(lonr), np.cos(latr) * np.sin(lonr), np.sin(latr)])
tracer = streamline.SimpleTracing(pos, bfield, 0.01)
for t in range(1000):
try:
tracer.step_forward()
x, y, z = tracer.x, tracer.y, tracer.z
l2 = x * x + y * y + z * z
if l2 > r2:
return derivepoint(np.array((x0, y0, z0)), np.array((x, y, z)), r=r)
else:
x0 = x; y0 = y; z0 = z
except IndexError:
logger.warn('Index error. Stop calculation')
break
tracer2 = streamline.SimpleTracing(pos, bfield, -0.01)
for t in range(1000):
try:
tracer2.step_forward()
x, y, z = tracer2.x, tracer2.y, tracer2.z
l2 = x * x + y * y + z * z
if l2 > r2:
return derivepoint(np.array((x0, y0, z0)), np.array((x, y, z)), r=r)
else:
x0 = x; y0 = y; z0 = z
except IndexError:
logger.warn('Index error. Stop calculation')
break
# If no data available (e.g. L-value inside give r), nan is returend.
return np.zeros([3]) + np.nan
def main(mf=None, pp=None, dlat=10, dlon=45):
if mf == None:
mf = mhddata.MagneticField1201()
mf.set_interpolate_strategy('project') # Using scaling projection.
if pp == None:
pp = mhddata.PlasmaParameter1205()
lats = np.arange(-90, 91, dlat)
lons = np.arange(0, 361, dlon)
projlonlist = np.zeros([len(lons), len(lats)]) + np.nan
projlatlist = np.zeros([len(lons), len(lats)]) + np.nan
projnlist = np.zeros([len(lons), len(lats)]) + np.nan
projvlist = np.zeros([3, len(lons), len(lats)]) + np.nan
projvthlist = np.zeros([len(lons), len(lats)]) + np.nan
projblist = np.zeros([3, len(lons), len(lats)]) + np.nan
projvblist = np.zeros([len(lons), len(lats)]) + np.nan
projfluxlist = np.zeros([len(lons), len(lats)]) + np.nan
# Suppress warnings.
logging.getLogger('MagneticField1201').setLevel(logging.WARN)
for ilon in range(len(lons)):
lon = lons[ilon]
for ilat in range(len(lats)):
lat = lats[ilat]
# Trace back to r=1.2.
pos1 = traceto(lon, lat, mf, r=1.2)
if np.isnan(pos1).any():
continue
r1 = np.sqrt((pos1 ** 2).sum())
# Longitude and latitude
projlon = np.rad2deg(np.arctan2(pos1[1], pos1[0]))
projlat = 90 - np.rad2deg(np.arccos(pos1[2] / r1))
projlonlist[ilon, ilat] = projlon
projlatlist[ilon, ilat] = projlat
# Parameters
print(pos1)
n, vx, vy, vz, temp, pth = pp.interpolate3d(pos1[0], pos1[1], pos1[2])
projnlist[ilon, ilat] = n * 1e6
projvlist[:, ilon, ilat] = vx * 1e3, vy * 1e3, vz * 1e3 # In m/s
vth = mhddata.t2vth(temp, mass=16.) # In m/s
projvthlist[ilon, ilat] = vth
# Bfield
bx, by, bz = mf.interpolate3d(pos1[0], pos1[1], pos1[2])
b = np.array([bx, by, bz])
# Change oxrientation toward surface!!
if (b * pos1).sum() > 0: # if b vector pointing upward
b = -b
projblist[:, ilon, ilat] = b
# Velocity component along b.
bn = b / np.sqrt((b * b).sum())
vb = (np.array([vx, vy, vz]) * bn).sum()
projvblist[ilon, ilat] = vb
# Precipitating flux
precip = maxwell.flux_precip(n * 1e6, vb, vth) * 1e-4 # In /cm2 s
projfluxlist[ilon, ilat] = precip
print(lon, lat, projlon, projlat, n, vx, vy, vz, temp, projvthlist[ilon, ilat], bx, by, bz, precip)
return mf, pp, lons, lats, projfluxlist, projlonlist, projlatlist, projnlist, projvlist, projvthlist, projblist, projvblist,
if __name__ == "__main__":
# mf, pp, lons, lats, projfluxlist, projlonlist, projlatlist, projnlist, projvlist, projvthlist, projblist, projvblist, = main()
retval = main()
lvalue = mhddata.LValue1201()
ocb = lvalue.opencloseboundaries()
plt.figure()
plt.pcolor(retval[2], retval[3], np.ma.log10(retval[4].T), vmin=7, vmax=8.5)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Downgoing Flux (log) /cm2 s')
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.colorbar()
plt.savefig('precipitation_1.png')
plt.figure()
plt.pcolor(retval[2], retval[3], retval[7].T / 1e6, vmin=0)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Density /cm3')
plt.colorbar()
plt.savefig('precipitation_2.png')
plt.figure()
v = np.sqrt((retval[8] ** 2).sum(0))
plt.pcolor(retval[2], retval[3], v.T / 1e3, vmin=0)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Velocity km/s')
plt.colorbar()
plt.savefig('precipitation_3.png')
plt.figure()
plt.pcolor(retval[2], retval[3], retval[9].T / 1e3, vmin=0)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Vth km/s')
plt.colorbar()
plt.savefig('precipitation_4.png')
plt.figure()
plt.pcolor(retval[2], retval[3], retval[11].T / 1e3, vmin=0)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('V// km/s')
plt.colorbar()
plt.savefig('precipitation_5.png')
plt.figure()
plt.pcolor(retval[2], retval[3], retval[5].T)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Longitude at project')
plt.colorbar()
plt.savefig('precipitation_6.png')
plt.pcolor(retval[2], retval[3], retval[6].T)
plt.plot(lvalue.lonlist(), ocb[0], 'w:', lw=2)
plt.plot(lvalue.lonlist(), ocb[1], 'w:', lw=2)
plt.title('Latitude at project')
plt.colorbar()
plt.savefig('precipitation_7.png')