<jupyter_start><jupyter_text>Demonstration of OrbitObjectsByHP - 20190312 - MJP - How to use the OrbitObjectsByHP routine to get the list of objects which touch a given healpix on a specified (integer) night - Healpix must be specified from a center: if nothing is specified, the default will be the geocenter, but otherwise this can be supplied as an obscode<jupyter_code>%matplotlib inline

# Standard imports 
import numpy as np 
import matplotlib.pyplot as plt 
import healpy as hp 

import MPC_library as MPCL
import phys_const as PHYS
import mpchecker 
import angles as ANG
import params<jupyter_output>/Users/matthewjohnpayne/anaconda/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)
/Users/matthewjohnpayne/anaconda/lib/python3.5/importlib/_bootstrap.py:222: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility. Expected 96, got 88
  return f(*args, **kwds)<jupyter_text>The results returned below will depend on the assumed HP scale for the preprocessing - Finer-scale HP will return smaller numbers of objects-per-HP - At the time of evaluation the assumed parameters are ...<jupyter_code>print("Healpix: Nside", params.hp_nside)
print("Healpix: Npix", params.hp_npix)
print("Healpix: Area [deg^2]", params.hp_area)
print("Healpix: Side [deg]", params.hp_side)<jupyter_output>Healpix: Nside 128
Healpix: Npix 196608
Healpix: Area [deg^2] 0.20982341130279172
Healpix: Side [deg] 0.45806485490898746<jupyter_text>Get a list of "OrbitObjects" which are in HP 348 & 375 - This currently takes 4 seconds ... much better than the ~2mins it used to take when directly reading "Objects" using pickle<jupyter_code>%%time
HPIndexList = [348,375]
JDint = 2458360
JDhalfint = 2458360.5
JDintStr = '2458360_5'
obsCode = '500'
OrbitObjectList = mpchecker.OrbitObjectsByHP(HPIndexList, JDintStr, obsCode)<jupyter_output>Opening '/Users/matthewjohnpayne/Dropbox/mpchecker/preprocessed_data/OrbitNumberByHP_500_2458360_5.npy'
Opening '/Users/matthewjohnpayne/Dropbox/mpchecker/preprocessed_data/orbitObject_states_2458360_5.npy'
states[0] [ 2.95912208e-04  2.45836050e+06  2.04408462e+00 -2.32335426e+00
  2.65748670e-01  6.56987317e-03  5.54146168e-03  3.55313810e-04]
Opening '/Users/matthewjohnpayne/Dropbox/mpchecker/preprocessed_data/orbitObject_desigs_2458360_5.npy'
desigs[0] J95O12J
CPU times: user 3.91 s, sys: 242 ms, total: 4.16 s
Wall time: 4.16 s<jupyter_text>Examine the returned data ...<jupyter_code>lenList = len(OrbitObjectList)
print("# Returned Objects = %d"%lenList )

OO = OrbitObjectList[0]
print("Type of first item in list = ", type(OO) )
print("Attributes of first item in list ...")
print("\t OO.cart_state = ", OO.cart_state)
print("\t OO.GMcenter = ", OO.GMcenter)
print("\t OO.epoch = ", OO.epoch)
print("\t OO.desig = ", OO.desig)<jupyter_output># Returned Objects = 7
Type of first item in list =  <class 'classes.OrbitObject'>
Attributes of first item in list ...
	 OO.cart_state =  <universal_kepler.universal_kepler.CartState object at 0xb1d302620>
	 OO.GMcenter =  0.00029591220828559115
	 OO.epoch =  2458360.5
	 OO.desig =  A0405<jupyter_text>Use RADEC_from_state_C to calculate the RA,DEC, (and hence the Unit Vector  & the healpix) from the Geocenter- I.e. we are double-checking that the returned objects make sense as possibilities for the input healpix<jupyter_code># When & Where are the simulated observations from?
obsCode = '500'
array_of_observation_times_UTC = np.array( [JDhalfint - 0.00080074] ) # HORRIBLE HACK ...
array_of_observation_times_TT = np.array( [MPCL.EOP.jdTT(array_of_observation_times_UTC[0])])
posEQ_OBS_wrt_HELIO = MPCL.Observatory().getObservatoryPosition(obsCode, array_of_observation_times_UTC[0])
array_of_observatory_positions_wrt_center_ECLIPTIC = np.dot(np.transpose(PHYS.rot_mat), posEQ_OBS_wrt_HELIO )

print("UTC",   array_of_observation_times_UTC)
print("TT",    array_of_observation_times_TT)
print("posEQ", posEQ_OBS_wrt_HELIO)
print("posEC", array_of_observatory_positions_wrt_center_ECLIPTIC)

# Advance the OrbitObjects from the OrbitObjectList above ...
RA_, DEC_, HP_ = [],[],[]
for OO in OrbitObjectList:
    r, d = OO.get_RADEC(array_of_observation_times_TT, array_of_observatory_positions_wrt_center_ECLIPTIC)
    # Save the RA & Dec
    RA_.append(r[0])
    DEC_.append(d[0])
    # Calculate the UV & HP
    UV = ANG.radec2unit(r,d)[0]
    HP_.append(hp.vec2pix(params.hp_nside, *UV, nest=params.hp_nested))
    
    
# Display the RA,DEC and calculate the UV & healpix
for i in range(5):
    print(RA_[i],DEC_[i], HP_[-1])
# Check that all of the healpix are as expected
uniqueHP_ = list(set(HP_))
print("uniqueHP_, len(uniqueHP_) = " , uniqueHP_, len(uniqueHP_) )
print("All unique output healpix are in input list: ", np.all([ uHP in HPIndexList for uHP in uniqueHP_ ]))
print("All input healpix are found in unique output list: ", np.all([ HP in uniqueHP_ for HP in HPIndexList ]))<jupyter_output>UTC [2458360.49919926]
TT [2458360.5]
posEQ [ 0.92484883 -0.37188159 -0.16121244]
posEC [ 9.24848826e-01 -4.05321309e-01  1.64062931e-05]
54.76705919525536 10.129123719735848 375
54.56340692149976 9.934946828054201 375
55.018831421836566 9.985703545848724 375
54.06997815604782 11.132633262624594 375
54.218404780871026 11.308403785940689 375
uniqueHP_, len(uniqueHP_) =  [348, 375] 2
All unique output healpix are in input list:  True
All input healpix are found in unique output list:  True<jupyter_text>Use RADEC_from_state_C to advance the returned OrbitObjects to a different time in the night - I.e. to the time of a specific pointing- Note that RADEC_from_state_C has been implemented as a method on the OrbitObjects<jupyter_code># When & Where are the simulated observatons from?
obsCode = 'G96'
array_of_observation_times_UTC = np.array( [JDint + 0.1] )
array_of_observation_times_TT = np.array( [MPCL.EOP.jdTT(array_of_observation_times_UTC[0])])
posEQ_OBS_wrt_HELIO = MPCL.Observatory().getObservatoryPosition(obsCode, array_of_observation_times_UTC[0])
array_of_observatory_positions_wrt_center_ECLIPTIC = np.dot(np.transpose(PHYS.rot_mat), posEQ_OBS_wrt_HELIO )

print("UTC",   array_of_observation_times_UTC)
print("TT",    array_of_observation_times_TT)
print("posEQ", posEQ_OBS_wrt_HELIO)
print("posEC", array_of_observatory_positions_wrt_center_ECLIPTIC)

# Advance the first OrbitObject from above ...
RA_, DEC_ = OO.get_RADEC(array_of_observation_times_TT, array_of_observatory_positions_wrt_center_ECLIPTIC)
for RA,DEC in zip(RA_, DEC_): print(RA,DEC)<jupyter_output>UTC [2458360.1]
TT [2458360.10080074]
posEQ [ 0.92218657 -0.37758264 -0.16367669]
posEC [ 9.22186573e-01 -4.11532150e-01  2.32497329e-05]
54.092385594043606 11.158747845008978
{'fraction_non_alphanumeric': 0.0988072448829333, 'fraction_numerical': 0.10351936386393756, 'mean_word_length': 5.778443113772455, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '152266', 'license': 'MIT', 'n_tokens_mistral': 2988, 'n_tokens_neox': 2523, 'n_words': 612, 'path': 'notebooks/zzz_old/Demonstrate_OrbitObjectsByHP.ipynb', 'repo_name': 'matthewjohnpayne/MPChecker2', 'size': '6791'}