''' Pitch angle distribution analysis (preliminary)
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import math
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pyana.pep.ganymede_magfield
def main():
# Specify the phase angle (longitude).
phase_angle = 70
phase_angle = phase_angle * np.pi / 180.
radius = 2631.
# height = 500.
height = 200.
# height = radius * 2
### Define the field model.
# fm = pyana.pep.ganymede_magfield.SimpleDipole()
fm = pyana.pep.ganymede_magfield.TiltedDipole()
### The FoV model.
# fov = fov_circle
# title_trail = '360 fov'
# fov = fov_nadir_circle
# title_trail = '180 nadir'
fov = fov_zenith_circle
title_trail = '180 zenith'
# Polar angle every 0.5 degrees
npol = 72
ndiv = 180
polar_angles = np.linspace(0, 2 * np.pi, npol)
pitchangle_distribution = []
losscones = []
bx = []
by = []
bz = []
bt = []
px = []
py = []
pz = []
for polar_angle in polar_angles:
pos = np.array([np.sin(polar_angle) * np.cos(phase_angle),
np.sin(polar_angle) * np.sin(phase_angle),
np.cos(polar_angle)])
pos_norm = pos
print('POS=', pos)
t0 = np.array([-np.sin(phase_angle),
np.cos(phase_angle),
0])
vel = np.cross(t0, pos)
print('VEL=', vel)
pos = pos * (height + radius)
px.append(pos[0])
py.append(pos[1])
pz.append(pos[2])
# A magnetic field vector
bvec = fm.get_field(pos)
print('PS,B', pos, bvec)
bx.append(bvec[0])
by.append(bvec[1])
bz.append(bvec[2])
btot = np.linalg.norm(bvec)
bvec_norm = bvec / btot
bt.append(btot)
# Surface footpoint B field
### The local latitude lat_loc in the magnetic moment coordinate
mvec_norm = fm.dipole_moment / np.linalg.norm(fm.dipole_moment)
colat_loc = np.arccos(mvec_norm.dot(pos_norm))
print('coLAT (local)=', colat_loc * 180. / np.pi)
print('polar angle', polar_angle * 180. / np.pi)
lat_loc = np.pi / 2. - colat_loc
print('LAT (local)=', lat_loc * 180. / np.pi)
### Absorpotion point
# r / cos^2 lat is constant.
# cos lat_abs = cos lat_local * sqrt(r_absorption / r_local)
lat_abs = np.arccos(np.cos(lat_loc) * np.sqrt(radius / (radius + height)))
print('LAT (absorption)=', lat_abs * 180. / np.pi)
dip = np.linalg.norm(fm.dipole_moment)
print('DIP', dip)
mag_abs = 1e-7 * dip / ((radius * 1e3) ** 3) * np.sqrt(1 + 3 * np.sin(lat_abs) ** 2) * 1e9
print('MAG (absorption)=', mag_abs)
### Loss cone angle, alpha = arcsin(sqrt(mag_local / mag_absorp))
losscone = np.arcsin(np.sqrt(btot / mag_abs)) * 180. / np.pi
print('LS', losscone)
losscones.append(losscone)
### The latitude: B = M / r^3 sqrt(1 + 3 sin^2 lat)
# FoV vectors
fov_vectors = fov(-pos, vel, ndiv=ndiv)
# Pitchangles
pitchangles = []
for fov_vec in fov_vectors:
pa = fov_vec.dot(bvec_norm)
pitchangles.append(pa) # So far, save just a inner dot.
pitchangles = -np.array(pitchangles) # ABove is FoV, now velocity
pitchangles = np.arccos(pitchangles) * 180. / np.pi # Make to the pitchangle.
print(pitchangles)
pitchangle_distribution.append(pitchangles)
pitchangle_distribution = np.array(pitchangle_distribution)
# Shape is (npol, ndiv)
fig = plt.figure()
ax = fig.add_subplot(111)
for ipol in range(npol):
polar_angle = polar_angles[ipol] * 180. / np.pi
p = polar_angle * np.ones_like(pitchangle_distribution[ipol, :])
ax.plot(p, pitchangle_distribution[ipol, :], 'b.')
ax.plot(polar_angles * 180. / np.pi, losscones, 'r-')
ax.plot(polar_angles * 180. / np.pi, 180 - np.array(losscones), 'r-')
ax.set_xlim(0, 360)
ax.set_xlabel('Polar angle (deg) = Time')
ax.set_ylabel('Pitch angle')
ax.set_title('Pitch angle coverage (longitude=%g) %s' % (phase_angle * 180. / np.pi, title_trail))
fig.savefig('analysis_pitchangle_%.1f.eps' % height)
fig = plt.figure()
ax = fig.add_subplot(211)
ax.plot(polar_angles * 180. / np.pi, bx)
ax.plot(polar_angles * 180. / np.pi, by)
ax.plot(polar_angles * 180. / np.pi, bz)
ax.set_xlim(0, 360)
ax.set_xlabel('Polar angle (deg) = Time')
ax.set_ylabel('B-field [nT]')
ax = fig.add_subplot(212)
ax.plot(polar_angles * 180. / np.pi, px)
ax.plot(polar_angles * 180. / np.pi, py)
ax.plot(polar_angles * 180. / np.pi, pz)
ax.set_xlabel('Polar angle (deg) = Time')
ax.set_ylabel('Position [km]')
return pitchangle_distribution
def fov_zenith_circle(nadir, vvec, ndiv=360):
fovs = fov_circle(nadir, vvec, ndiv=ndiv)
nadir_vec = np.array(nadir)
valid_fov = []
for i in range(ndiv):
fov = fovs[i, :]
# Only nadir pointing is returned.
if np.dot(nadir_vec, fov) <= 0:
valid_fov.append(fov)
else:
valid_fov.append([np.nan, np.nan, np.nan])
return np.array(valid_fov)
def fov_nadir_circle(nadir, vvec, ndiv=360):
fovs = fov_circle(nadir, vvec, ndiv=ndiv)
nadir_vec = np.array(nadir)
valid_fov = []
for i in range(ndiv):
fov = fovs[i, :]
# Only nadir pointing is returned.
if np.dot(nadir_vec, fov) >= 0:
valid_fov.append(fov)
else:
valid_fov.append([np.nan, np.nan, np.nan])
return np.array(valid_fov)
def fov_circle(nadir, vvec, ndiv=360):
''' FoV model with nadir-pointing circular plane with 360 degrees.
Returned is np.array with shape of (ndiv, 3).
All the vector is normalized.
'''
logger = logging.getLogger('fov_circle')
logger.setLevel(logging.WARN)
# First calculate "a" vector that is perpendicualr to vvec
vvec = np.array(vvec)
vmin = np.argmin(np.abs(vvec))
vx, vy, vz = vvec
if vmin == 0:
nvec = np.array([0, vz, -vy])
elif vmin == 1:
nvec = np.array([vz, 0, -vx])
elif vmin == 2:
nvec = np.array([vy, -vx, 0])
else:
raise ValueError('')
nvec_len = np.sqrt((nvec ** 2).sum())
nvec = nvec / nvec_len
tvec = np.cross(vvec, nvec)
tvec_len = np.linalg.norm(tvec)
tvec = tvec / tvec_len
logger.info('Velocity vector = %s' % str(vvec))
logger.info('1st orth vector = %s (%e)' % (str(nvec), vvec.dot(nvec)))
logger.info('2nd orth vector = %s (%e) (%e)' % (str(tvec), vvec.dot(tvec), nvec.dot(tvec)))
angles = np.linspace(0, 2 * np.pi, ndiv)
vecs = []
for angle in angles:
vec = nvec * np.cos(angle) + tvec * np.sin(angle)
logger.debug('(%f) %s (dot=%e)' % (angle * 180. / np.pi,
str(vec), vec.dot(vvec)))
vec_len = np.linalg.norm(vec)
vec = vec / vec_len
vecs.append(vec)
return np.array(vecs)
if __name__ == '__main__':
pa = main()