aboutsummaryrefslogtreecommitdiffstats
path: root/IC_Helpers/genic_astorb.py
blob: 9667e9d23dfa7b1523934e281617505e9ae88fa3 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
"""
Generate Genga Initial Conditions File from Astorb Entries.

Usage:
$ python ./genic_astorb.py --fname_in astorb.hdf5 --crossers > initial.dat

Changes:

* 15 May 2016 / Volker Hoffmann <volker@cheleb.net>
  Limit Number of Objects (--limit)

* 15 May 2016 / Volker Hoffmann <volker@cheleb.net>
  Add Perturbed Clone Generator (--clones)

* 15 May 2016 / Volker Hoffmann <volker@cheleb.net>
  Initial Release
"""

import sys
import pandas as pd
import kepler_helpers as kh
import constants as C
import argparse
import numpy as np


# Parse Arguments
parser = argparse.ArgumentParser()
parser.add_argument('--fname_in', required=True, \
                    help='Name of Astorb Source File.')
parser.add_argument('--clones', action='store_true', \
                   help="Generate Perturbed Clones for Each Crosser.")
parser.add_argument('--limit', type=int, default=-1, \
                   help="Limit Number of Objects.")
group = parser.add_mutually_exclusive_group(required=True)
group.add_argument('--crossers', action='store_true', \
                   help="Use Earth Crossers (~1'000).")
group.add_argument('--all', action='store_true', \
                   help="Use All Asteroid (~600'000).")
args = parser.parse_args()

# Load Data
sys.stderr.write('// Loading Data\n')
df = pd.read_hdf("%s" % args.fname_in, 'df')

# Select Earth Crossers
# Cf. ftp://cdsarc.u-strasbg.fr/pub/cats/B/astorb/ReadMe
if args.crossers:
    sys.stderr.write('// Select Earth Crossers\n')
    df = df[df.Xflg==1]

# Limit Number of Objects
if args.limit > 0:
    sys.stderr.write("// Limiting to %i Objects\n" % args.limit)
    df = df.head(args.limit)

# Generate IC Lines
if args.clones:
    sys.stderr.write('// Generating Genga IC Lines, With Clones\n')
    pid_base = 10000; lines = []

    # Semi-Major Axis Offsets
    # a = a * [ 1,
    #           1-10^{-16}, 1-10^{-15}, ..., 1-10^{-1}
    #           1+10^{-1}, ..., 1+10^{-15}, 1-10^{-16} ]

    da_lo = 1.0 - np.logspace(-16, -1, 16)
    da_hi = 1.0 + np.logspace(-16, -1, 16)
    da_xx = np.concatenate([np.array([1.0]), da_lo, da_hi])

    for ii, [ irow, row ] in enumerate(df.iterrows()):
        for jj, da in enumerate(da_xx):
            # Convert to Cartesian State Vector
            x, v = \
                kh.kep2cart(row.a * da, row.e, row.i * C.d2r, \
                            row.Omega * C.d2r, row.omega * C.d2r, \
                            row.M * C.d2r, \
                            0.0, 1.0)
            v *= C.kms_to_genga

            # Format IC Line
            line = "0.0 %06d %.16e %.16e " % ( pid_base + jj, 0.0, \
                                               row.Diameter / C.au2km )
            line += "%+.16e %+.16e %+.16e " % ( x[0], x[1], x[2] )
            line += "%+.16e %+.16e %+.16e " % ( v[0], v[1], v[2] )
            line += "0.0 0.0 0.0"
            lines.append(line)

        # Base Particle ID Counter
        pid_base += 100

else:
    sys.stderr.write('// Generating Genga IC Lines\n')
    pid = 10000; lines = []
    for ii, [ irow, row ] in enumerate(df.iterrows()):
        # Convert to Cartesian State Vector
        x, v = \
            kh.kep2cart(row.a, row.e, row.i * C.d2r, \
                        row.Omega * C.d2r, row.omega * C.d2r, row.M * C.d2r, \
                        0.0, 1.0)
        v *= C.kms_to_genga

        # Format IC Line
        line = "0.0 %06d %.16e %.16e " % ( pid, 0.0, row.Diameter / C.au2km )
        line += "%+.16e %+.16e %+.16e " % ( x[0], x[1], x[2] )
        line += "%+.16e %+.16e %+.16e " % ( v[0], v[1], v[2] )
        line += "0.0 0.0 0.0"
        lines.append(line)

        # Particle IC Counter
        pid += 1

# Print Lines
sys.stderr.write('// Printing Genga IC Lines\n')
for line in lines:
    print line