apps120911_iotorus.app04_eneflux_batchΒΆ

''' 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()