''' ENA flux calculation, batch script version
'''
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 scipy.interpolate
import pyana.pep.iotorus
from .app03_enaflux import getenaflux
def main():
iotorus = pyana.pep.iotorus.TwoPeakModel()
scpos = [0, 15.9 * np.cos(np.deg2rad(30)), 15.9 * np.sin(np.deg2rad(30))]
scpos = np.array(scpos)
print('# POS', scpos[0], scpos[1], scpos[2])
# ax0 = [0, -1, 0] # Boresight direction, theta=0
ax0 = np.array([5.88, 0, 0]) - scpos # Look at (close to io)
# ax0 = -scpos # Look at Jupiter
ax0 = ax0 / np.sqrt(np.dot(ax0, ax0))
ax0 = np.array(ax0)
ax0 = ax0 / np.sqrt((ax0 * ax0).sum())
print('# AX0', ax0[0], ax0[1], ax0[2])
ax1 = [1, 0, 0] # Phi origin. theta=90 (or whatever), phi=0
ax1 = np.array(ax1)
ax1 = ax1 / np.sqrt((ax1 * ax1).sum())
ax2 = np.cross(ax0, ax1)
ax1 = np.cross(ax2, ax0)
print('# AX1', ax1[0], ax1[1], ax1[2])
print('# AX2', ax2[0], ax2[1], ax2[2])
thetalist = list(range(0, 41, 3))
thetalist[0] = 0.1 # Avoid zero.
philist = list(range(0, 360, 3))
nall = len(thetalist) * len(philist)
ndone = 0
logger = logging.getLogger('app03_enaflux.py')
logger.setLevel(logging.INFO)
import time
t0 = time.time()
for phi in philist:
for theta in thetalist:
t = np.deg2rad(theta)
p = np.deg2rad(phi)
lookdir = ax0 * np.cos(t) + ax1 * np.sin(t) * np.cos(p) + ax2 * np.sin(t) * np.sin(p)
flx = getenaflux(sc=scpos, los=lookdir, iotorus_model=iotorus)
print(theta, phi, lookdir[0], lookdir[1], lookdir[2], flx)
logger.info( 'th=%f ph=%f look=%f, %f, %f, flx=%e' % (theta, phi, lookdir[0], lookdir[1], lookdir[2], flx))
ndone += 1.0
t1 = time.time()
logger.info('%d/%d, %ds. Left %ds. Total %ds.' % (
ndone, nall, t1 - t0, (t1 - t0) / ndone * (nall - ndone),
(t1 - t0) / ndone * nall))
if __name__ == "__main__":
main()