Running vplanet

Overview

Let’s run vplanet on the CircumbinaryOrbit example.

import vplanet

output = vplanet.run("examples/CircumbinaryOrbit/vpl.in")
output

<vplanet.Output: kepler16>

The variable output is an instance of vplanet.Output. We can peek at its contents by accessing the members attribute:

output.members

['sysname', 'bodies', 'log', 'path', 'primary', 'secondary', 'cbp', 'earth']

The properties sysname, bodies, and path are just strings, but log is a vplanet.Log instance and the remaining four are vplanet.Body instances.

The log file

Let’s peek at log:

log = output.log
log

<vplanet.Log: kepler16>

Again, let’s look at the members attribute:

log.members

['sysname', 'path', 'header', 'initial', 'final']

The header, initial, and final properties are instances of vplanet.LogStage, containing information about specific stages of the simulation. Let’s check out everything in the header of the log file:

for member in log.header.members:
    print("{}: {}".format(member, getattr(log.header, member)))

Executable: /usr/share/miniconda/envs/vplot/bin/vplanet
Version: Unknown
SystemName: kepler16
PrimaryInputFile: vpl.in
BodyFile1: primary.in
BodyFile2: secondary.in
BodyFile3: cbp.in
BodyFile4: earth.in
AllowFilesToBeOverwitten: True
MassUnits: Grams
LengthUnits: Meters
TimeUnits: Seconds
AngleUnits: Radians
VerbosityLevel: 5
CrossoverDecadeForScientificNotation: 4
NumberOfDigitsAfterDecimal: 6
IntegrationMethod: Runge-Kutta4
Direction: Forward
TimeStep: 31557600.0
StopTime: 3155760000.0
OutputInterval: 315576.0
UseVariableTimestep: True
dEta: 0.01
MinimumValueOfEccAndObl: 1e-10

Note that all of these can, of course, be accessed as regular properties:

log.header.Executable

'/usr/share/miniconda/envs/vplot/bin/vplanet'

Next, let’s inspect the initial stage of the log file:

log.initial.members

['system', 'primary', 'secondary', 'cbp', 'earth']

There is one attribute per body, plus one for the system. Here’s what’s in earth:

for member in log.initial.earth.members:
    print("{}: {}".format(member, getattr(log.initial.earth, member)))

ActiveModules: BINARY
ModuleBitSum: 1025
Color: 0
Mass: 5.972186e+24 kg
Radius: 6378100.0 m
RadGyra: 0.5
RotAngMom: 4.416946e+33 kg m2 / sec
BodyType: 0.0
Density: 5495.038549 kg / m3
HZLimitDryRunaway: -1.0 m
HZLimRecVenus: -1.0
HZLimRunaway: -1.0
HZLimMoistGreenhouse: -1.0
HZLimMaxGreenhouse: -1.0
HZLimEarlyMars: -1.0
Instellation: -1.0 kg / sec3
Eccentricity: 0.009302
OrbEnergy: 0.0 kg m2 / sec2
MeanMotion: 1.877507e-07 1 / sec
OrbPeriod: 33905420.0 sec
SemiMajorAxis: 150905800000.0 m
CriticalSemiMajorAxis: -1.0 m
COPP: 0.0
OrbAngMom: 0.0 kg m2 / sec
ArgP: 3.970693 rad
Inc: 0.005425 rad
LongA: 3.112206 rad
LongP: 0.799715 rad
TotOrbEnergy: 0.0 kg m2 / sec2
OrbPotEnergy: 0.0 kg m2 / sec2
FreeEcc: 0.03
FreeInc: 0.005381 rad
LL13N0: 1.88385e-07 sec
LL13K0: 1.870849e-07 sec
LL13V0: 1.896763e-07 sec
CBPR: 149950500000.0 m
CBPZ: 0.03614 m
CBPPhi: 6.253799 rad
CBPRDot: -191.811805 m / sec
CBPZDot: -152.682529 m / sec
CBPPhiDot: 1.87673e-07 1 / sec
R0: 149597900000.0 m
CBPInsol: -1.0 F/F_Earth
BinPriR: 8366404000.0 m
BinPriPhi: 0.658853 rad
BinSecR: 28406370000.0 m
BinSecPhi: 3.800446 rad
OutputOrder: Time[year] CBPR[au] ArgP[deg] LongA[deg] Eccentricity[] Inc[deg] LongP[deg]
GridOutputOrder: None

The other attributes similarly list all of the initial properties of each body (or the system).

Before we move on, note the type of each of these properties:

type(log.initial.earth.Inc)

vplanet.quantity.VPLANETQuantity

This isn’t a regular number, but an instance of vplanet.Quantity, which is a subclasss of astropy.units.Quantity. These objects have numerical values (either scalars or arrays) and units; plus, they automatically handle unit conversions when operated on.

inc = log.initial.earth.Inc
print(inc)
print(inc.to("deg"))

0.005425 rad
0.31082960385847164 deg

mass = log.initial.earth.Mass
print(mass)
print(mass.to("Mearth"))

5.972186e+24 kg
1.000003036118378 earthMass

Check out the astropy documentation for more information on supported units.

The bodies

Going back to the output object, let’s inspect its earth member:

output.earth

<vplanet.Body: earth>

This object contains the information of the arrays computed during the simulation, which are specified in the saOutputOrder line of the input files.

output.earth.members

['name',
 'infile',
 'fwfile',
 'bwfile',
 'climfile',
 'Time',
 'CBPR',
 'ArgP',
 'LongA',
 'Eccentricity',
 'Inc',
 'LongP']

Let’s check out one of them:

output.earth.Inc

$[0.310855,~0.311045,~0.311036,~\dots,~0.310957,~0.310998,~0.311019] \; \mathrm{{}^{\circ}}$

As before, these are vplanet.Quantity instances, which behave like numpy arrays except they have units attached.

print(output.earth.Inc.to("rad"))

[0.00542544 0.00542876 0.0054286  ... 0.00542722 0.00542794 0.00542831] rad

That’s it for this tutorial! Check out the plotting tutorial for information on how to easily visualize all of these quantities.