apps120811_flux_recalc.project_at_500km

''' Start from a point at surface, calculate the 500 km point
'''


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.pep import mhddata


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 np.array((x0, y0, z0)), np.array((x, y, z))
            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 np.array((x0, y0, z0)), np.array((x, y, z))
            else:
                x0 = x; y0 = y; z0 = z

        except IndexError:
            logger.warn('Index error. Stop calculation')
            break

    return np.zeros([3]), np.zeros([3])



def main(mf=None):

    if mf == None:
        mf = mhddata.MagneticField1201()
    mf.set_interpolate_strategy('project')

#    lats = np.arange(-90, 91, 1)
#    lons = np.arange(0, 361, 2)
    lats = np.arange(-90, 91, 10)
    lons = np.arange(0, 361, 45)

    projlonlist = np.zeros([len(lons), len(lats)])
    projlatlist = np.zeros([len(lons), len(lats)])

    logging.getLogger('MagneticField1201').setLevel(logging.WARN)

    for ilon in range(len(lons)):
        lon = lons[ilon]
        for ilat in range(len(lats)):
            lat = lats[ilat]

            pos0, pos1 = traceto(lon, lat, mf, r=1.2)
            r1 = np.sqrt((pos1 ** 2).sum())

            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

            print(lon, lat, projlon, projlat)

    return mf, lons, lats, projlonlist, projlatlist

if __name__ == "__main__":
    
    mf = main()