import urllib
cat = 'J/A+A/415/241/table1'
uri = 'http://vizier.u-strasbg.fr/viz-bin/asu-binfits/VizieR?&-source={}&-out.all'.format(cat)
url = urllib.URLopener()
filen,msg = url.retrieve(uri,filename='cat.fits')
url.close()


import pyfits
hdus = pyfits.open('cat.fits')
hdus
hdus.info()


table = hdus[1].data






hdus[0].header


hdus[1].header


hdus[1].columns.names


hjd = hdus[1].data.field('HJD')
vmag = hdus[1].data.field('Vmag')
hdus.close()


hjd.shape
vmag.dtype


hjd.ptp()
vmag.mean()


import matplotlib.pyplot as plt
plt.plot(hjd,vmag,'ko')


import numpy as np
a = np.array([10,20,30,40])
b = np.arange(4)
a + b**2  # elementwise operation


f = np.ones([3, 4])                 # 3 x 4 float array of ones
f
g = np.zeros([2, 3, 4], dtype=int)  # 3 x 4 x 5 int array of zeros
i = np.zeros_like(f)                # array of zeros with same shape/type as f


w = np.arange(12)
w.shape = [3, 4]       # does not modify the total number of elements
w
x = np.arange(5)
y = x.reshape(5, 1)
y = x.reshape(-1, 1)  # Same thing but NumPy figures out correct length


x = np.linspace(0,1,3)
y = np.linspace(0,10,3)
np.hstack([x,y])         # shape (6,)
np.vstack([x,y])         # shape (2,3)
np.column_stack([x,y])   # shape (3,2)
xgrid, ygrid = np.meshgrid(x,y)         # 2x shape (3,3)
print(xgrid)
print(ygrid)


x = np.linspace(-20,20,160)
y = np.linspace(-20,20,160)
x,y = np.meshgrid(x,y)        # converts x and y into 2D arrays representing coordinates
r = np.sqrt(x**2 + y**2)
z = np.cos(r) / (r + 5)
import matplotlib.pyplot as plt
plt.imshow(z)


import scipy.signal


df = 0.2 / hjd.ptp()
freq_interval = np.arange(6.2, 7.2, df)*2*np.pi


vmag = vmag - vmag.mean()


hjd = np.array(hjd,float)
vmag = np.array(vmag,float)
pgram = scipy.signal.lombscargle(hjd,vmag,freq_interval)


Ndata = len(hjd)         # equivalent to Ndata = hjd.shape[0]
norm = np.sqrt(4.0 * pgram / Ndata)
freq_interval_cy = freq_interval/(2*np.pi)
plt.plot(freq_interval_cy,norm,'k-')


freq = freq_interval[np.argmax(pgram)]
print(freq/(2*np.pi))


C = np.column_stack([np.sin(freq*hjd), np.cos(freq*hjd)])
X, sumsq, rank, s = np.linalg.lstsq(C,vmag)


amplitude = np.sqrt(X[0]**2 + X[1]**2)
phase = np.arctan2(X[1],X[0])
fit = amplitude * np.sin(freq*hjd + phase)


phased = 2*np.pi * (hjd % (2*np.pi/freq))      # frequency in rad/s, but time in s!
plt.plot(phased,vmag,'ko');


hjd
phased


x = np.linspace(0,1,11)
print(x)
index = np.array([0,1,2,3,4])
print(x[index])
index = np.array([0,3,1,4,2])
print(x[index])


x = np.random.uniform(size=5)
x
index = np.argsort(x)
index
x[index]


x = np.linspace(0,1,11)
x
index = np.array([True,False,True,True])
x[index]


x = np.random.normal(size=5)
x
is_positive = x > 0
is_positive
x[is_positive]


sortarray = np.argsort(phased)
select = (fit[sortarray] < 0) & (phased[sortarray] < 0.5)
plt.plot(phased[sortarray], fit[sortarray], 'r-',lw=2);
plt.plot(phased[sortarray][select], fit[sortarray][select], 'b-',lw=4);


J = np.column_stack([np.sin(freq*hjd + phase), amplitude * np.cos(freq*hjd + phase)])
cov = sumsq / (Ndata - 2) * np.linalg.inv(np.dot(J.T, J))
e_amplitude = np.sqrt(cov.diagonal()[0])
e_phase = np.sqrt(cov.diagonal()[1])
print("Amplitude = {:.4f} +/- {:.4f} mag".format(amplitude,e_amplitude))
print("Phase = {:.2f} +/- {:.2f} rad".format(phase,e_phase))