Source code for irfpy.earth.magnetosphere

import os,sys
import math
import logging

class MsModelShue:
	''' Magnetosphere model by Shue et al., 1997.  Assuming that Bz=0/Dp=3 '''

	def __init__(self, Bz=0, Dp=3):
		self._re=6378.
		if Bz >= 0:
			self.r0 = (11.4+0.013 * Bz) * Dp ** (-1./6.6)
		else:
			self.r0 = (11.4+0.14 * Bz) * Dp ** (-1./6.6)

		logging.debug('R0 = %f'%self.r0)
		self.alpha = (0.58 - 0.01*Bz) * (1+0.010*Dp)
		logging.debug('Alpha = %f'%self.alpha)

	def isInside(self, x, y, z):
		''' Return true if the position (x,y,z)[Rm] is in or on the bow shock. '''
		r = math.sqrt(x*x+y*y+z*z)
		t = math.atan2(math.sqrt(y*y+z*z), x)
		rb = self.r0 * ((2/(1+math.cos(t)))  ** self.alpha)
		logging.debug('MS Distance for angle %f = %f'%((t*180./math.pi), rb))
		return r<rb

	def isInsideKm(self, x, y, z):
		''' Return true if the position (x,y,z)[Km] is in or on the bow shock. '''
		rx = x/self._re
		ry = y/self._re
		rz = z/self._re
		return self.isInside(rx,ry,rz)

	def getBounds2D(self):
		bnd=[]
		for td in range(180):
			t=td*math.pi/180.
			rb = self.r0 * ((2/(1+math.cos(t)))  ** self.alpha)
			x=rb*math.cos(t)
			y=rb*math.sin(t)
			bnd.append([x,y])
		return bnd
		
def getboundary():
    ''' Return the boundary in Rm.
    '''
    model = MsModelShue()
    return model.getBounds2D()
			
import unittest
import doctest


def doctests():
    return unittest.TestSuite((
        doctest.DocTestSuite(),
        ))
if __name__ == '__main__':
    unittest.main(defaultTest='doctests')