Python 基本語法與科學計算套件的使用： Python科學計算套件(二)

In [1]:

    

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline


    

In [2]:

    

time = np.linspace(0, 10, 200)
counts = 50 * np.sin(2 * np.pi * 1. / 2.5 * time) + 100  +  np.random.normal(0, 5., len(time))
plt.plot(time, counts, 'k.')
counts_err = 4 * np.random.rand(len(time)) + 1
plt.errorbar(time, counts, yerr=counts_err, fmt="r")


    

In [3]:

    

def lc_model(time, amplitude, period, dc):
    return amplitude * np.sin(2 * np.pi * 1. / period * time) + dc


    

In [4]:

    

popt, pcov = curve_fit(lc_model, time, counts,  p0=[50, 2.5, 100], sigma=counts_err)
print(popt)
print(pcov)


    

    




[ 49.85353949   2.50166026  99.89739415]
[[  2.22598077e-01   5.26041656e-06   1.51224999e-02]
 [  5.26041656e-06   2.45107993e-06   2.79458011e-05]
 [  1.51224999e-02   2.79458011e-05   1.12455146e-01]]

In [5]:

    

perr = np.sqrt(np.diag(pcov))
print("Amplitude =", popt[0], "+/-", perr[0])
print("Period =", popt[1], "+/-", perr[1])
print("DC =", popt[2], "+/-", perr[2])


    

    




Amplitude = 49.8535394854 +/- 0.471803006579
Period = 2.50166026043 +/- 0.0015655925175
DC = 99.8973941468 +/- 0.335343325886

In [6]:

    

plt.errorbar(time, counts, yerr=counts_err, fmt="none", label='Data')
plt.plot(time, lc_model(time, popt[0], popt[1], popt[2]),'r-', label='Model')
plt.xlabel('Time (hour)')
plt.ylabel('Counts')
plt.legend()


    

In [7]:

    

from astropy.io import fits
hdulist = fits.open("XTE_J1550_564_30191011500A_2_13kev_0.01s_0_2505s.fits")
hdulist.info()


    

    




Filename: XTE_J1550_564_30191011500A_2_13kev_0.01s_0_2505s.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      30   ()              
1    RATE        BinTableHDU     82   250500R x 5C   [D, E, E, E, 1D]   
2    STDGTI      BinTableHDU     41   1R x 2C      [1D, 1D]   

In [8]:

    

primary = hdulist[0]
primary.header


    

    
Out[8]:





SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                    8 / number of bits per data pixel                  
NAXIS   =                    0 / number of data axes                            
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
CONTENT = 'LIGHT CURVE'        / light curve file                               
ORIGIN  = 'NASA/GSFC/XTE/GOF'  / origin of fits file                            
DATE    = '2013-07-22T06:57:10' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
TELESCOP= 'XTE     '           / Telescope (mission) name                       
INSTRUME= 'PCA     '           / Instrument used for observation                
MJDREFI =                49353 / Integer part of MJDREF                         
MJDREFF = 6.965740740000000E-04 / Fractional part of MJDREF                     
DATE-OBS= '29/09/98'           / EARLIEST observation date of files             
TIME-OBS= '05:15:13'           / EARLIEST time of all input files               
DATE-END= '29/09/98'           / LATEST observation date of files               
TIME-END= '05:57:07'           / LATEST time of all input files                 
TIMESYS = 'TT      '           / The time system used was                       
TIMEUNIT= 's       '           / unit for TSTARTI/F and TSTOPI/F, TIMEZERO      
TSTART  = 1.49663661378431E+08 / Observation start time                         
TSTOP   = 1.49666166378431E+08 / Observation stop time                          
OBJECT  = 'XTE_J1550-564'      / OBJECT from the FIRST input file               
RA_PNT  =       2.37669998E+02 / RA of First input object                       
DEC_PNT =      -5.64599991E+01 / DEC of First input object                      
EQUINOX =              2000.00 / Equinox of the FIRST object                    
RADECSYS= 'FK5     '           / Co-ordinate frame used for equinox             
TIMVERSN= 'OGIP/93-003'        / OGIP memo number for file format               
CREATOR = 'SAEXTRCT version 4.2g' / Program name that produced this file        
CHECKSUM= 'WhNSWfKRWfKRWfKR'   / HDU checksum updated 2013-07-22T06:57:10       
DATASUM = '         0'         / data unit checksum updated 2013-07-22T06:57:10 

In [9]:

    

print(primary.header['OBJECT '])


    

    




XTE_J1550-564

In [10]:

    

rate = hdulist[1]
rate.data


    

    
Out[10]:





FITS_rec([(149663661.37843132, 17682.732, 1329.1155, 1.0009766, 149663476.18321466),
       (149663661.38843131, 20479.963, 1430.3768, 1.0009766, 149663476.193214),
       (149663661.39843133, 19281.191, 1387.8818, 1.0009766, 149663476.20321333),
       ...,
       (149666166.34843132, 1498.5366, 386.91968, 1.0009766, 149665980.97182214),
       (149666166.35843131, 1398.6373, 373.80069, 1.0009766, 149665980.98182136),
       (149666166.36843133, 1618.1729, 404.54321, 0.98876953, 149665980.9918206)], 
      dtype=[('TIME', '>f8'), ('RATE', '>f4'), ('ERROR', '>f4'), ('FRACEXP', '>f4'), ('BARYTIME', '>f8')])

In [11]:

    

print(rate.data.shape)
print(rate.data.field(1))


    

    




(250500,)
[ 17682.73242188  20479.96289062  19281.19140625 ...,   1498.53662109
   1398.6373291    1618.17285156]

In [12]:

    

time = rate.data.field(0)[:500]
counts = rate.data.field(1)[:500]
plt.plot(time, counts)
plt.xlabel('Time (s)')
plt.ylabel('Rate (count/s)')


    

In [13]:

    

hdulist.close()