%matplotlib inline
%config InlineBackend.figure_format = 'retina'

from __future__ import print_function
import sys

import numpy as N
import libstempo as T
import libstempo.plot as LP, libstempo.toasim as LT

T.data = T.__path__[0] + '/data/' # example files

print("Python version   :",sys.version.split()[0])
print("libstempo version:",T.__version__)
print("Tempo2 version   :",T.libstempo.tempo2version())

"""
Explanation: libstempo tutorial: simulating residuals with toasim
Michele Vallisneri, vallis@vallis.org, 2014/10/31
This notebook demonstrates the libstempo module toasim, which allows the simple simulation of various kinds of noise.
End of explanation
"""


psr = T.tempopulsar(parfile = T.data + 'B1953+29_NANOGrav_dfg+12.par',
                    timfile = T.data + 'B1953+29_NANOGrav_dfg+12.tim')
LP.plotres(psr)

"""
Explanation: We open up a NANOGrav par/tim file combination with libstempo, and plot the residuals.
End of explanation
"""


LT.make_ideal(psr)
LP.plotres(psr)

"""
Explanation: We now remove the computed residuals from the TOAs, obtaining (in effect) a perfect realization of the deterministic timing model. The pulsar parameters will have changed somewhat, so make_ideal calls fit() on the pulsar object.
End of explanation
"""


#LT.add_line(psr,f=10**6.5,A=1e-5)
LT.add_efac(psr,efac=1.0,seed=1234)
LP.plotres(psr)

"""
Explanation: We now add a single line of noise at $10^{6.5}$ Hz, with an amplitude of 10 us. We also put back radiometer noise, with rms amplitude equal to 1x the nominal TOA errors.
All the noise-generating commands take an optional argument seed that will reseed the numpy pseudorandom-number generator, so you are able to reproduce the same instance of noise. However, if you issue several noise-generating commands in sequence, you should use different seeds.
End of explanation
"""


LT.add_rednoise(psr,1e-12,3)
LP.plotres(psr)

"""
Explanation: We could also add EQUAD quadrature noise (with add_equad) or its coarse-grained version (with add_jitter), but instead we prefer some red noise of "GW-like" amplitude $10^{-12}$ and spectral slope $\gamma = -3$.
End of explanation
"""


LT.add_gwb(psr,flow=1e-8,gwAmp=5e-12)
LP.plotres(psr)

help(LT.add_gwb)

LT.createGWB([psr],Amp=5e-15,gam=13./3.)
LP.plotres(psr)

"""
Explanation: Or, we may add a GW background as simulated by the tempo2 GWbkgrd plugin (see the docstring below).
End of explanation
"""


psr.fit()
LP.plotres(psr)

"""
Explanation: Refitting will remove some of the power.
End of explanation
"""


psr.savepar('B1953+29-simulate.par')
psr.savetim('B1953+29-simulate.tim')

"""
Explanation: All done! We can save the resulting par and tim file, and analyze them with a favorite pipeline.
End of explanation
"""


T.purgetim('B1953+29-simulate.tim')

"""
Explanation: Note that currently the tim file that is output by tempo2 has a spurious "MODE 1" line that tempo2 does not like upon reloading. To erase it, you can do
End of explanation
"""


psr2 = T.tempopulsar(parfile = 'B1953+29-simulate.par',
                     timfile = 'B1953+29-simulate.tim')
LP.plotres(psr2)

"""
Explanation: And if we reload the files we get pack the same thing...
End of explanation
"""


psr = LT.fakepulsar(parfile=T.data+'B1953+29_NANOGrav_dfg+12.par',
                    obstimes=N.arange(53000,54800,30)+N.random.randn(60),  # observe every 30+-1 days
                    toaerr=0.1)

LT.add_efac(psr,efac=1.0,seed=1234)
LP.plotres(psr)

help(LT.fakepulsar)

"""
Explanation: It's also possible to obtain a perfect realization of the timing model described in a par file without a tim file, by specifying a new set of observation times (in MJD) and errors (in us). The observation frequency, observatory, and flags can also be specified (see the docstring below).
End of explanation
"""


