6.4 Fitting a periodic lightcurve ########################################## Reference script ------------------------------------------- The least-square fit can be made with the following script. A similar script is provided in the file "06_cdr\\Lightcurve_fitperiod.py". :: from aspylib import cdr, astro import numpy as np import matplotlib.pyplot as plt #------ input files -------------------------------- Folder = u"C:\\Mesures\\" Filelist = ["JC4_01Sept2011.txt", "JC4_14Nov2011.txt", "JC4_15Nov2011.txt"] Filelist = [Folder + x for x in Filelist] #------ loading data ------------------------------- alldata = cdr.load_magdata(Filelist) #------ time correction ---------------------------- alpha = astro.hours_to_decdegrees('00:08:29.67') delta = astro.degrees_to_decdegrees('+64:19:46.9') alldata = cdr.time_correction(alldata, 2, [alpha,delta]) #------ performing the fit ------------------------- T0 = 2455800 period = 0.04982 #initial guess for the period, in days popt, pcov, period, dperiod, amplitude, maxdate, mindate = cdr.fit_CDR(alldata, period, T0, nb_harmo = 2, show_results = True, show_plot = True) plt.ylim([-0.061,0.061]) plt.title('Lightcurve with its Fourier fit') plt.show() Results for star JC4 ------------------------------------------- For the variable star JC4 we find the following results. With the parameter "period" set to 0.04982 we deliberately forced the fit to find the same period as CourbRot. .. image:: fit.png .. image:: fit2.png On the CdR-CdL website we find: .. image:: JC4.png As we see we recover very similar results in both cases: * period = 0.0498230 days (AsPyLib) / 0.0498228 days (CdR-CdL) * period uncertainty = +/- 0.0000002 (AsPyLib) / +/- 0.0000004 days (CdR-CdL) which probably indicates the use of different noise estimates * amplitude = 0.0289 (AsPyLib) / 0.0291 (CdR-CdL) * heliocentric date of min magnitude = 01.9463 September 2011 in both cases. We can note the very small value of period uncertainty, which comes directly from the very narrow width of the oscillation feature in which the fitting algorithm finally ends (see the analysis perfored `here `_). Since it's not clear if the minimum found by the algorithm is the right one, this uncertainty value is not yet representative. Calculating a phase ephemeris ------------------------------------------- It is easy to calculate a phase ephemeris. This can be done by adding the following code to the script above: :: #------ phase ephemeris ---------------------------- JJ = np.linspace(2455822.5, 2455842.5, 201) print print "------ phase ephemerids ------" phase = (JJ-T0)/popt[0]-np.floor((JJ-T0)/popt[0]) for i,dd in enumerate(JJ): print dd, phase[i] raw_input("Press Enter")