notebook.community

Edit and run

如何自訂函式、模組及類別?

參考資料

準備工作：範例所需檔案下載及說明

  • MAXI_lc_1day-Cyg_X-1.csv, MAXI_lc_1orbit-Cyg_X-1.csv, MAXI_lc_1day-GX_339-4.csv, MAXI_lc_1orbit-GX_339-4.csv, MAXI_lc_1day-LMC_X-4.csv, MAXI_lc_1orbit-LMC_X-4.csv, MAXI_lc_1day-4U_1323-619.csv, MAXI_lc_1orbit-4U_1323-619.csv, MAXI_lc_1day-Cyg_X-3.csv, MAXI_lc_1orbit-Cyg_X-3.csv (都已在files4examples資料夾中無須下載)：

這幾個X-ray雙星系統的光變曲線資料下載自X射線望遠鏡MAXI的官網，檔名中1day和1orbit分別表示以1天和90分鐘取樣。 檔案內容為不同能量區段的光變曲線資料，資料的詳細格式請見MAXI官網的說明。

  • 其他範例檔案(待補)

範例1-1：自訂專門畫MAXI某段能量區段光變觀測資料的函數





In [1]:



    


def plot_maxi_lc(filename, band='all', ploterr=True):
    """    
    輸入檔名，畫出MAXI某段能量區段觀測資料的光變曲線。
    
    Parameters
    ----------------   
    filename: 字串
    band: 字串 
          能量區段， 2-20 keV:  'all'(預設)
                   2-4 keV:   'low'
                   4-10 keV:  'mid'
                   10-20 keV: 'high'
    ploterr: 布林值
             是否也把error bar畫出來(預設為True)   
    """
    from astropy.io import ascii
    import matplotlib.pyplot as plt
    lc = ascii.read(filename)
    time = lc['col1']
    rate_2_20keV = lc['col2']
    err_2_20keV = lc['col3']
    rate_2_4keV = lc['col4']
    err_2_4keV = lc['col5']
    rate_4_10keV = lc['col6']
    err_4_10keV = lc['col7']
    rate_10_20keV = lc['col8']
    err_10_20keV = lc['col9']
    rate = {
        'all': [rate_2_20keV, err_2_20keV, '(2-20 keV)'],
        'low': [rate_2_4keV, err_2_4keV, '(2-4 keV)'],
        'mid': [rate_4_10keV, err_4_10keV, '(4-10 keV)'],
        'high': [rate_10_20keV, err_10_20keV, '(10-20 keV)'],        
    }
    if ploterr == True:
        plt.errorbar(time, rate.get(band)[0], yerr=rate.get(band)[1], fmt='.b')
    else:
        plt.plot(time, rate.get(band)[0], '.k')        
    plt.xlabel('Time (MJD)', fontsize=12)
    plt.ylabel('Count rate', fontsize=12)
    plt.title('Cyg X-1 light curve observed by MAXI ' + rate.get(band)[2], fontsize=12)
    plt.show()



    











In [2]:



    


# 查看函數用途
plot_maxi_lc?



    











In [3]:



    


# 此程式區塊為畫出個別能量區段的光變曲線 (每天1個資料點)
%matplotlib inline
lc_1day = "../files4examples/MAXI_lc_1day-Cyg_X-1.csv"
plot_maxi_lc(lc_1day)
plot_maxi_lc(lc_1day, band='low')
plot_maxi_lc(lc_1day, band='mid', ploterr=False)
plot_maxi_lc(lc_1day, band='high')



    


















    

















    

















    

















    
























In [4]:



    


# 此程式區塊為畫出個別能量區段的光變曲線 (90分鐘1個資料點)
lc_1orbit = "../files4examples/MAXI_lc_1orbit-Cyg_X-1.csv"
plot_maxi_lc(lc_1orbit)
plot_maxi_lc(lc_1orbit, band='low')
plot_maxi_lc(lc_1orbit, band='mid', ploterr=False)
plot_maxi_lc(lc_1orbit, band='high')

範例1-2：自訂專門畫MAXI所有能量區段光變觀測資料的函數





In [5]:



    


def plot_maxi_lcs(filename, subplot=False):
    """    
    輸入檔名，將MAXI所有能量區段的光變曲線畫在同一張圖。
    
    Parameters
    ----------------   
    filename: 字串
    subplot: 布林值
             True: 以subplot的型式畫圖
             Flase(預設): 直接把所有光變曲線疊起來               
    """
    from astropy.io import ascii
    import matplotlib.pyplot as plt
    lc = ascii.read(filename)
    time = lc['col1']
    rate_2_20keV = lc['col2']
    err_2_20keV = lc['col3']
    rate_2_4keV = lc['col4']
    err_2_4keV = lc['col5']
    rate_4_10keV = lc['col6']
    err_4_10keV = lc['col7']
    rate_10_20keV = lc['col8']
    err_10_20keV = lc['col9']
    if subplot == False:
        plt.errorbar(time, rate_2_20keV, yerr=err_2_20keV, fmt='.r')
        plt.hold(True)
        plt.errorbar(time, rate_2_4keV, yerr=err_2_4keV, fmt='.b')
        plt.errorbar(time, rate_4_10keV, yerr=err_4_10keV, fmt='.k')
        plt.errorbar(time, rate_10_20keV, yerr=err_10_20keV, fmt='.m')
        plt.hold(False)
        plt.legend(('2-20 keV', '2-4 keV', '4-10 keV', '10-20 keV'))
        plt.ylabel('Count rate', fontsize=12)
        plt.title('Cyg X-1 light curves observed by MAXI', fontsize=12)
    else:
        plt.subplot(411)
        plt.errorbar(time, rate_2_20keV, yerr=err_2_20keV, fmt='.r')
        plt.title('Cyg X-1 light curves observed by MAXI')
        plt.ylabel('2-20 keV')
        plt.xticks([])
        plt.subplot(412)
        plt.errorbar(time, rate_2_4keV, yerr=err_2_4keV, fmt='.b')
        plt.ylabel('2-4 keV')
        plt.xticks([])
        plt.subplot(413)
        plt.errorbar(time, rate_4_10keV, yerr=err_4_10keV, fmt='.k')
        plt.ylabel('4-10 keV')
        plt.xticks([])
        plt.subplot(414)
        plt.errorbar(time, rate_10_20keV, yerr=err_10_20keV, fmt='.m')
        plt.ylabel('10-20 keV')        
    plt.xlabel('Time (MJD)', fontsize=12)    
    plt.show()



    











In [6]:



    


# 查看函數用途
plot_maxi_lcs?



    











In [7]:



    


# 此程式區塊為在一張圖上畫出所有能量區段的光變曲線 (每天1個資料點)
%matplotlib inline
lc_1day = "../files4examples/MAXI_lc_1day-Cyg_X-1.csv"
plot_maxi_lcs(lc_1day)
plot_maxi_lcs(lc_1day, True)



    


















    

















    
























In [8]:



    


# 此程式區塊為在一張圖上畫出所有能量區段的光變曲線 (90分鐘1個資料點)
lc_1orbit = "../files4examples/MAXI_lc_1orbit-Cyg_X-1.csv"
plot_maxi_lcs(lc_1orbit)
plot_maxi_lcs(lc_1orbit, True)

範例1-3：將上例兩個函數寫進一個檔案中，成為可引入的模組





In [9]:



    


output_string = '''
"""用來畫出MAXI光變觀測資料的模組

可用函數：
- plot_maxi_lc: 畫出MAXI單一能量區段的光變曲線
- plot_maxi_lcs: 畫出MAXI所有能量區段的光變曲線
"""

from astropy.io import ascii
import matplotlib.pyplot as plt

def plot_maxi_lc(filename, band='all', ploterr=True):
    """    
    輸入檔名，畫出MAXI某段能量區段觀測資料的光變曲線。
    
    Parameters
    ----------------   
    filename: 字串
    band: 字串 
          能量區段， 2-20 keV:  'all'(預設)
                   2-4 keV:   'low'
                   4-10 keV:  'mid'
                   10-20 keV: 'high'
    ploterr: 布林值
             是否也把err bar畫出來(預設為True)   
    """
    lc = ascii.read(filename)
    time = lc['col1']
    rate_2_20keV = lc['col2']
    err_2_20keV = lc['col3']
    rate_2_4keV = lc['col4']
    err_2_4keV = lc['col5']
    rate_4_10keV = lc['col6']
    err_4_10keV = lc['col7']
    rate_10_20keV = lc['col8']
    err_10_20keV = lc['col9']
    rate = {
        'all': [rate_2_20keV, err_2_20keV, '(2-20 keV)'],
        'low': [rate_2_4keV, err_2_4keV, '(2-4 keV)'],
        'mid': [rate_4_10keV, err_4_10keV, '(4-10 keV)'],
        'high': [rate_10_20keV, err_10_20keV, '(10-20 keV)'],        
    }
    if ploterr == True:
        plt.errorbar(time, rate.get(band)[0], yerr=rate.get(band)[1], fmt='.b')
    else:
        plt.plot(time, rate.get(band)[0], '.k')        
    plt.xlabel('Time (MJD)', fontsize=12)
    plt.ylabel('Count rate', fontsize=12)
    plt.title('Cyg X-1 light curve observed by MAXI ' + rate.get(band)[2], fontsize=12)
    plt.show()
    
    
def plot_maxi_lcs(filename, subplot=False):
    """    
    輸入檔名，將MAXI所有能量區段的光變曲線畫在同一張圖。
    
    Parameters
    ----------------   
    filename: 字串
    subplot: 布林值
             True: 以subplot的型式畫圖
             Flase(預設): 直接把所有光變曲線疊起來               
    """
    lc = ascii.read(filename)
    time = lc['col1']
    rate_2_20keV = lc['col2']
    err_2_20keV = lc['col3']
    rate_2_4keV = lc['col4']
    err_2_4keV = lc['col5']
    rate_4_10keV = lc['col6']
    err_4_10keV = lc['col7']
    rate_10_20keV = lc['col8']
    err_10_20keV = lc['col9']
    if subplot == False:
        plt.errorbar(time, rate_2_20keV, yerr=err_2_20keV, fmt='.r')
        plt.hold(True)
        plt.errorbar(time, rate_2_4keV, yerr=err_2_4keV, fmt='.b')
        plt.errorbar(time, rate_4_10keV, yerr=err_4_10keV, fmt='.k')
        plt.errorbar(time, rate_10_20keV, yerr=err_10_20keV, fmt='.m')
        plt.hold(False)
        plt.legend(('2-20 keV', '2-4 keV', '4-10 keV', '10-20 keV'))
        plt.ylabel('Count rate', fontsize=12)
        plt.title('Cyg X-1 light curves observed by MAXI', fontsize=12)
    else:
        plt.subplot(411)
        plt.errorbar(time, rate_2_20keV, yerr=err_2_20keV, fmt='.r')
        plt.title('Cyg X-1 light curves observed by MAXI')
        plt.ylabel('2-20 keV')
        plt.xticks([])
        plt.subplot(412)
        plt.errorbar(time, rate_2_4keV, yerr=err_2_4keV, fmt='.b')
        plt.ylabel('2-4 keV')
        plt.xticks([])
        plt.subplot(413)
        plt.errorbar(time, rate_4_10keV, yerr=err_4_10keV, fmt='.k')
        plt.ylabel('4-10 keV')
        plt.xticks([])
        plt.subplot(414)
        plt.errorbar(time, rate_10_20keV, yerr=err_10_20keV, fmt='.m')
        plt.ylabel('10-20 keV')        
    plt.xlabel('Time (MJD)', fontsize=12)    
    plt.show()
'''

output_file = open('maxi_lc.py','w')
output_file.write(output_string)
output_file.close()



    











In [10]:



    


# 引入模組並查看該模組的相關說明
import maxi_lc as mlc

mlc?
#mlc.plot_maxi_lc?
#mlc.plot_maxi_lcs?



    











In [11]:



    


# 使用模組中的函式
lc_1day = "../files4examples/MAXI_lc_1day-Cyg_X-1.csv"
lc_1orbit = "../files4examples/MAXI_lc_1orbit-Cyg_X-1.csv"
%matplotlib inline
mlc.plot_maxi_lc(lc_1day)
mlc.plot_maxi_lc(lc_1day, band='low')
mlc.plot_maxi_lc(lc_1day, band='mid', ploterr=False)
mlc.plot_maxi_lc(lc_1day, band='high')
mlc.plot_maxi_lc(lc_1orbit)
mlc.plot_maxi_lc(lc_1orbit, band='low')
mlc.plot_maxi_lc(lc_1orbit, band='mid', ploterr=False)
mlc.plot_maxi_lc(lc_1orbit, band='high')
mlc.plot_maxi_lcs(lc_1day)
mlc.plot_maxi_lcs(lc_1day, True)
mlc.plot_maxi_lcs(lc_1orbit)
mlc.plot_maxi_lcs(lc_1orbit, True)



    


















    

















    

















    

















    

















    

















    

















    

















    

















    

















    

















    

















    
























In [12]:



    


# 刪除模組及其相關檔案
import os
import shutil
os.remove('maxi_lc.py')
shutil.rmtree('__pycache__')

範例1-4：自訂MAXI光變曲線類別





In [13]:



    


class MaxiLc():
    """ MAXI光變曲線類別 """

    
    def __init__(self, source_name, lc_1day_filename, lc_1orbit_filename):
        from astropy.io import ascii
        lc_1day = ascii.read(lc_1day_filename, names=["TIME", "RATE (2-20keV)", "err (2-20keV)", 
                                                      "RATE (2-4keV)", "err (2-4keV)", 
                                                      "RATE (4-10keV)", "err (4-10keV)", 
                                                      "RATE (10-20keV)", "err (10-20keV)"] )
        lc_1orbit = ascii.read(lc_1orbit_filename, names=["TIME", "RATE (2-20keV)", "err (2-20keV)", 
                                                      "RATE (2-4keV)", "err (2-4keV)", 
                                                      "RATE (4-10keV)", "err (4-10keV)", 
                                                      "RATE (10-20keV)", "err (10-20keV)"] )
        self.source_name = source_name
        self.lc_1day = lc_1day
        self.lc_1orbit = lc_1orbit

        
    def plot_maxi_lc(self, dt, band='all', ploterr=True):
        """    
        輸入觀測時間間隔，畫出MAXI某段能量區段觀測資料的光變曲線。
    
        Parameters
        ----------------   
        dt: 字串
            觀測時間間隔("day"或"orbit")
        band: 字串 
              能量區段， 2-20 keV:  'all'(預設)
                       2-4 keV:   'low'
                       4-10 keV:  'mid'
                       10-20 keV: 'high'
        ploterr: 布林值
                 是否也把error bar畫出來(預設為True)
        """
        import matplotlib.pyplot as plt
        if dt == 'day':
            lc = self.lc_1day    
        elif dt == 'orbit':
            lc = self.lc_1orbit        
        else:
            print('請輸入正確的觀測時間間隔("day"或"orbit")')    
        time = lc['TIME']
        rate_2_20keV = lc['RATE (2-20keV)']
        err_2_20keV = lc['err (2-20keV)']
        rate_2_4keV = lc['RATE (2-4keV)']
        err_2_4keV = lc['err (2-4keV)']
        rate_4_10keV = lc['RATE (4-10keV)']
        err_4_10keV = lc['err (4-10keV)']
        rate_10_20keV = lc['RATE (10-20keV)']
        err_10_20keV = lc['err (10-20keV)']
        rate = {
            'all': [rate_2_20keV, err_2_20keV, '(2-20 keV)'],
            'low': [rate_2_4keV, err_2_4keV, '(2-4 keV)'],
            'mid': [rate_4_10keV, err_4_10keV, '(4-10 keV)'],
            'high': [rate_10_20keV, err_10_20keV, '(10-20 keV)'],        
        }
        if ploterr == True:
            plt.errorbar(time, rate.get(band)[0], yerr=rate.get(band)[1], fmt='.b')
        else:
            plt.plot(time, rate.get(band)[0], '.k')        
        plt.xlabel('Time (MJD)', fontsize=12)
        plt.ylabel('Count rate', fontsize=12)
        plt.title(self.source_name + ' light curve observed by MAXI ' + rate.get(band)[2], fontsize=12)
        plt.show()



    











In [14]:



    


# 利用MAXI光變曲線類別建立不同天體的光變曲線(物件)
lc_cyg_x1 = MaxiLc("Cyg X-1", "../files4examples/MAXI_lc_1day-Cyg_X-1.csv", "../files4examples/MAXI_lc_1orbit-Cyg_X-1.csv")
lc_gx339_4 =  MaxiLc("GX 339-4", "../files4examples/MAXI_lc_1day-GX_339-4.csv", "../files4examples/MAXI_lc_1orbit-GX_339-4.csv") 
lc_lmc_x4 =  MaxiLc("LMC X-4", "../files4examples/MAXI_lc_1day-LMC_X-4.csv", "../files4examples/MAXI_lc_1orbit-LMC_X-4.csv") 
lc_4u_1323_619 =  MaxiLc("4U 1323-619", "../files4examples/MAXI_lc_1day-4U_1323-619.csv", "../files4examples/MAXI_lc_1orbit-4U_1323-619.csv")
lc_cyg_x3 = MaxiLc("Cyg X-3", "../files4examples/MAXI_lc_1day-Cyg_X-3.csv", "../files4examples/MAXI_lc_1orbit-Cyg_X-3.csv")



    











In [15]:



    


# 此程式區塊為示範使用光變曲線物件的"屬性"
print(lc_cyg_x1.source_name)
print(lc_cyg_x1.lc_1day)



    


















    






Cyg X-1
    TIME     RATE (2-20keV) err (2-20keV) ... RATE (10-20keV) err (10-20keV)
------------ -------------- ------------- ... --------------- --------------
     55058.5       0.917442      0.034054 ...        0.172521       0.020003
     55059.5       1.043909      0.036054 ...        0.184611       0.021066
     55060.5       1.164664      0.033378 ...         0.17648       0.018508
     55061.5       1.334263       0.04054 ...        0.216942       0.022606
     55062.5       1.333345      0.041085 ...        0.224568       0.022853
     55063.5       1.306505      0.052945 ...        0.246716       0.029038
     55064.5       1.138778      0.039837 ...        0.229176       0.022868
     55065.5       1.172685      0.088554 ...         0.20063       0.043301
     55066.5       1.162694      0.023896 ...        0.252453       0.013336
     55067.5       0.970559      0.020609 ...        0.159804       0.011458
         ...            ...           ... ...             ...            ...
57579.499977       1.119939      0.033945 ...        0.198253       0.016014
57580.499977       1.023292      0.033558 ...        0.190985       0.016068
57581.499977       1.218075      0.037125 ...        0.203354       0.016809
57582.499977       1.120713      0.035425 ...        0.192787       0.016299
57583.499977        1.21908      0.035061 ...        0.215413       0.016227
57584.499977       1.145633      0.036962 ...        0.195677       0.017052
57585.499977       1.042149      0.035383 ...        0.174579       0.016525
57586.499977       1.189613      0.060465 ...        0.248028       0.029518
57587.499977       1.153689      0.034472 ...        0.185533       0.015477
57588.499977        1.18946      0.036611 ...        0.237238       0.018198
57589.499988       1.282976      0.043708 ...        0.182316        0.01899
Length = 1625 rows

















In [16]:



    


# 此程式區塊為示範使用光變曲線物件的"方法"
%matplotlib inline
lc_cyg_x1.plot_maxi_lc('day')
lc_gx339_4.plot_maxi_lc('day', band='mid')
lc_lmc_x4.plot_maxi_lc('day', band='low', ploterr=False)
lc_4u_1323_619.plot_maxi_lc('orbit', band='mid')
lc_cyg_x3.plot_maxi_lc('orbit', band='high', ploterr=False)

其他範例：尚待補充有函數回傳值的例子

Content source: Astrohackers-TW/IANCUPythonAdventure

Similar notebooks:

notebook.community | gallery | about