Skip to content

Home

Jump to bottom
Oscar Barragán edited this page Nov 26, 2021 · 8 revisions

First steps

Dependencies

Before starting to use the code, be sure that you have the following compilers and python libraries in your machine

  • gfortran (or your favorite fortran compiler)
  • OpenMP
  • numpy
  • matplotlib
  • seaborn (optional)

We tested this code using Linux, including Fedora and Ubuntu. This code has also been being built on MAC computers.

Test run for pyaneti

Just clone or download pyaneti.

git clone https://github.com/oscaribv/pyaneti

The advantage about cloning the repository is the possibility to follow the changes to this package easily with git pull (learn more about git at https://git-scm.com/).

The next step is to enter the pyaneti folder and see what we can find inside it

cd pyaneti
ls
  LICENSE  src inpy  makefile  pyaneti.py  README.md

Now compile the code typing make (you can also compile the code in parallel following these steps)

make

If you have all the dependencies installed, the making process should end without error.

You know pyaneti has compiled if there is a *.so file created in your main pyaneti directory (example pyaneti.cpython-36m-x86_64-linux-gnu.so or pyaneti.so)

If you have problems compiling the code, check this.

Now you are ready to run the code for the first time, just type

./pyaneti.py test

or

python pyaneti.py test

The program will start. You will see something like:

==============================
------------------------------
    INITIAL CONFIGURATION
------------------------------
Star           = test
No. planets    = 1
------------------------------
iter max       = 100000000
thin factor    = 1
nconv          = 500
nwalkers       = 100
------------------------------
fit RV         = [True]
fit Transit    = [True]
------------------------------
LC data        = free
cadence time   =  30.000 min
n rebinning    = 1
Stellar priors = [True]
------------------------------
  PLANET testb
------------------------------
        PRIOR RANGES
------------------------------
T0 = u[ 2448285.0500 , 2448285.1500 ]
P  = u[ 365.2060 , 365.3060 ]
ew1= f[ 0.0000 , 1.0000 ]
ew2= f[ 0.0000 , 1.0000 ]
b  = f[ 0.0000 , 1.0000 ]
a  = g[ 200.0000 , 250.0000 ]
rp = u[ 0.0000 , 0.1000 ]
K  = u[ 0.0000 , 0.0010 ]
------------------------------
 Other parameter priors
------------------------------
q1 = g[ 0.3464 , 0.0500 ]
q2 = g[ 0.2839 , 0.0500 ]
Super_telescope = u[ 21.0719 , 23.0721 ]
==============================
 CREATING RANDOM SEED
 CREATING CHAINS

 STARTING MCMC CALCULATION
 RV datapoints  =           36
 TR datapoints  =          112
 No. parameters =            8
 dof            =          140
 ==================================
      Chain statistics
 ==================================
 chain |  reduced chi^2
 best  :   64667.93
 worst : **********
 mean  : **********
 ==================================
 STARTING INFINITE LOOP!
 ==================================
      Chain statistics
 ==================================
 chain |  reduced chi^2
 best  :       0.95
 worst :       1.05
 mean  :       0.99
 ==================================
 ==================================
      PERFOMING GELMAN-RUBIN
       TEST FOR CONVERGENCE
 ==================================
 ==================================
 CHAINS HAVE NOT CONVERGED YET!
         500  ITERATIONS MORE!
 ==================================
 ==================================
      Chain statistics
 ==================================
 chain |  reduced chi^2
 best  :       0.95
 worst :       1.09
 mean  :       0.99
 ==================================
 ==================================
      PERFOMING GELMAN-RUBIN
       TEST FOR CONVERGENCE
 ==================================
 ==================================
       CHAINS HAVE CONVERGED
 ==================================
    CREATING OUTPUT DATA FILES
 ==================================
STARTING CHAIN CLUSTERING
Initial number of chains: 100
Final number of chains: 100

--------------------------------------------------------------
Summary:
N_chains         =      100
N_iter           =      500
thin_factor      =        1
N_rv_data        =       36
N_tr_data        =      112
N_data           =      148
N_pars           =        8
chi2_rv          = 7.2765
chi2_tr          = 125.9964
chi2             = 133.2729
dof              =      140
chi2/dof         = 0.9519
ln likelihood_rv = 377.7453
ln likelihood_tr = 1131.1755
ln likelihood    = 1508.9208
BIC              = -2977.8640
AIC              = -3001.8417
--------------------------------------------------------------
             INPUT STELLAR PARAMETERS
--------------------------------------------------------------
M_*     = 1.0000000 - 0.1000000 + 0.1000000 solar masses
R_*     = 1.0000000 - 0.1000000 + 0.1000000 solar radii
T_*     = 5772.0000000 - 100.0000000 + 100.0000000 K
--------------------------------------------------------------
                   Parameters test b
-------------------------Fitted-------------------------------
        T0 = 2448285.0955436 - 0.0031192 + 0.0031409     days
         P = 365.2508520 - 0.0049167 + 0.0048356     days
      ew 1 = 0.0000000 - 0.0000000 + 0.0000000
      ew 2 = 0.0000000 - 0.0000000 + 0.0000000
         b = 0.0000000 - 0.0000000 + 0.0000000
      a/R* = 215.4640517 - 2.3708930 + 2.5122541
     rp/R* = 0.0093235 - 0.0000712 + 0.0000775
         K = 0.0881821 - 0.0024223 + 0.0023856      m/s
-------------------------Derived------------------------------
        Mp = 0.9841032 - 0.0704355 + 0.0718640  M_earth
        Rp = 1.0173252 - 0.1028264 + 0.1013333  R_earth
         e = 0.0000000 - 0.0000000 + 0.0000000
         w = 0.0000000 - 0.0000000 + 0.0000000      deg
         i = 90.0000000 - 0.0000000 + 0.0000000      deg
         a = 1.0019439 - 0.1014392 + 0.1005043       AU
Insolation = 0.9967687 - 0.0707513 + 0.0741733  F_Earth
      rho* = 1.4183711 - 0.0463035 + 0.0502306 g/cm^3 (transit)
      rho* = 1.4065628 - 0.3676615 + 0.5586840 g/cm^3 (stellar paramters)
     rho_p = 5.1298925 - 1.3074050 + 1.9887038   g/cm^3
       g_p = 936.8970215 - 30.5861395 + 32.3670771 cm/s^2 (K and Rp/R*)
       g_p = 933.8572371 - 172.6787114 + 234.1233896 cm/s^2 (planet parameters)
     Tperi = 2448193.7827742 - 0.0039813 + 0.0041039     days
       Teq = 278.1026529 - 5.0720559 + 5.0352570 K (albedo=0)
     T_tot = 13.0713846 - 0.1507021 + 0.1448215    hours
    T_full = 12.8294707 - 0.1483685 + 0.1432300    hours
--------------------------------------------------------------
--------------------  Other parameters -----------------------
        q1 = 0.3760168 - 0.0418577 + 0.0442539
        q2 = 0.3157777 - 0.0460379 + 0.0468480
        u1 = 0.3855446 - 0.0574350 + 0.0598140
        u2 = 0.2246861 - 0.0583548 + 0.0612362
Sys. vel. Super_telescope = 22.0719865 - 0.0000015 + 0.0000016      m/s
--------------------------------------------------------------

Creating  outpy/test_out/test_posterior.pdf
Creating  outpy/test_out/test_correlations.pdf
Creating  outpy/test_out/testb_tr.pdf
Creating  outpy/test_out/test_rv_all.pdf
Creating  outpy/test_out/testb_rv.pdf

If you see this output it means that pyaneti ended succesfully.

Now let us check the plots (I use evince but you can use your favourite pdf reader)

evince outpy/test_out/testb_tr.pdf outpy/test_out/testb_rv.pdf

You will see some nice plots that look like this

testb_tr testb_rv

Let me explain you briefly what this test fit was about:

If you were an advanced alien civilization with really high technology, and "lucky" enough to see an Earth-like planet crossing in front of a Sun-like star, this is how the Earth would look like to you.

Look at those well-known parameters:

  • 1 Earth Mass
  • 1 Earth radii
  • Period of 365 days
  • 1 AU semi-major axis
  • Density of ~5.5 g/cm^2,
  • Gravity of ~10 m/s^2.

Of course, you would need a spectrograph with a precision of a few cm/s and also a very nice photometer, and an excellent way to remove stellar intrinsic signals in your data, but that is another story.

If you are at this point, you learned two things. First, with good data you can obtain really nice planet parameters and second, you learned how to run pyaneti.

Play with test - Joint radial velocity and transit light curve fitting.

  • There is a directory called inpy, inside this folder, you will find a second directory called test.

  • This directory contains the input and data files to perform the test fit.

  • You can create an input director inside inpy for each of your systems!

  • We encourage you to start to play with the code.

Create your own test directory and copy all the files from inpy/test folder.

mkdir inpy/my_test
cp inpy/test/* inpy/my_test

Now you are ready to run my_test

./pyaneti.py my_test

You will see an output similar to that of the test case. Now the output files are inside outpy/my_test_out. You will notice that inside this directory you will find an extra file with the posterior distribution and correlation plots of the fitted parameters.

Now open the file inpy/my_test/input_fit.py and start to play with it. The file is commented. Let us change the priors for some parameters. Uncomment lines 56, 91 and 92 to fit for the scaled semi-major axis and impact factor. Save the changes and re-run the code.

./pyaneti.py my_test

Now you can see that the fitted parameters are different in comparison with the values are given by test.

If you have some RV and/or transit data you only have to put the name of your data files, change the prior ranges, and start to fit your data!

Clone this wiki locally