Problem 2

CMD

Extract data from ascii files and overplot star evolutionary tracks


In [10]:
#Extract data form ascii files
#SUN-MS
# open file
f = open('./ms_1.00.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL100=[]
logTef100=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL100.append(sou1)
    logTef100.append(sou2)
f.close()

In [11]:
#SUN-AGB
# open file
f = open('./agb_1.00.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL_agb100=[]
logTef_agb100=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL_agb100.append(sou1)
    logTef_agb100.append(sou2)
f.close()

In [12]:
#SUN-HB
# open file
f = open('./hb_1.00.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL_hb100=[]
logTef_hb100=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL_hb100.append(sou1)
    logTef_hb100.append(sou2)
f.close()

In [13]:
#Extract data form ascii files
#vega-MS
# open file
f = open('./ms_2.10.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL210=[]
logTef210=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL210.append(sou1)
    logTef210.append(sou2)
f.close()

In [14]:
#vega-AGB
# open file
f = open('./agb_2.10.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL_agb210=[]
logTef_agb210=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL_agb210.append(sou1)
    logTef_agb210.append(sou2)
f.close()

In [16]:
#vega-HB
# open file
f = open('./hb_2.00.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL_hb200=[]
logTef_hb200=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL_hb200.append(sou1)
    logTef_hb200.append(sou2)
f.close()

In [26]:
#Extract data form ascii files
#regulus-MS
# open file
f = open('./ms_3.50.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL350=[]
logTef350=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL350.append(sou1)
    logTef350.append(sou2)
f.close()

In [28]:
#Extract data form ascii files
#regulues-agb
# open file
f = open('./ms_3.50.dat', 'r')
#remove file's header
header1 = f.readline()
header2 = f.readline()
#create new lists
logL_agb350=[]
logTef_agb350=[]

#extract columns
for line in f:
    line = line.strip()
    columns = line.split()
    sou1=float(columns[2])
    sou2=float(columns[3])
    logL_agb350.append(sou1)
    logTef_agb350.append(sou2)
f.close()

H-R & Evolutionary tracks


In [33]:
#plot H-R map
import matplotlib.pyplot as plt
from astropy.io import fits
import numpy as np
import math

#figure
fig, ax = plt.subplots()
ax.invert_xaxis()

ax.set_title('HR Diagram',fontweight='bold',fontsize='20')
ax.set_xlim(4.2,3.4)
ax.set_ylim(-4,7)
ax.set_xlabel('log(Tef/K)',fontsize=20)
ax.set_ylabel('log(L/Lsun)',fontsize=20)

plt.xticks(fontsize=20)
plt.yticks(fontsize=20)

#major ticks every 2, minor ticks every0.5
major_ticks = np.arange(-4, 7, 2)                                              
minor_ticks = np.arange(-4, 7, 0.5)                                               
ax.set_yticks(major_ticks)                                                       
ax.set_yticks(minor_ticks, minor=True)                                           



#plot H-R map use hippacos star data
hdu3=fits.open('./hippacos.fit')
datas3=hdu3[1].data
mask = (datas3['B-V']!=0)
datas3 = datas3[mask]
color=datas3['B-V']
vmag=datas3['Vmag']

#absolute magnitude
pc=list(map(abs,1000/datas3['Plx']))# abs of Plx(mas)--distance(pc)
Mav=vmag+5*(1-np.log10(pc))
# covert color index B-V into Tef
Tef=4600*(1/((0.92*color)+1.7)+1/((0.92*color)+0.62))
L=-0.4*(Mav-4.83)
Tef=np.log10(Tef)
ax.scatter(Tef,L,2,c='k',alpha=0.5,label='stars')




#plot 20 brightes stars 
a=np.arange(0,20)
strnam=['Achernar.fit','Acrux.fit','Aldebaran.fit','Altair.fit','Antares.fit',
        'Betelgeuse.fit','Canopus.fit','Capella.fit','Deneb.fit','Fomalhaut.fit',
        'Pollux.fit','Procyon.fit','Regulus.fit','Rigel.fit','Shaula.fit',
        'Vega.fit','Sirius.fit','Spica.fit','Mimosa.fit','Alnilam.fit']

for each in a:
    hdu4=fits.open(strnam[each])
    datas4=hdu4[1].data
    color=datas4['B-V']#star colors
    vmagstar=datas4['Vmag']
    mas=datas4['Plx']#million second arc mas
    dpc=list(map(abs,1000/mas))# distance in pc unit
    Mv=vmagstar+5*(1-np.log10(dpc))#absolute magnitude
    Tef=4600*(1/((0.92*color)+1.7)+1/((0.92*color)+0.62))
    #L=np.log10(10**(-0.4*Mv))
    L=-0.4*(Mv-4.83)
    Tef=np.log10(Tef)
    ax.scatter(Tef,L,60,c='orange',alpha=1,label=strnam[each][:-4])#substract from strnam by ".fit"
    plt.text(Tef,L,strnam[each][:-4],fontweight='bold',fontsize=13,color='red',rotation=45)

   

#sun
ax.scatter(np.log10(5772),0,60,c='orange',alpha=1,label='Sun')
plt.text(np.log10(5772),0.3,'Sun',fontweight='bold',fontsize=18,color='red',rotation=45)


#set figure size
fig.set_size_inches(20, 15)
#legend = ax.legend(loc='lower left', shadow=True, fontsize='x-large')
#legend.get_frame().set_facecolor('#00FFCC')


#plot star evolutionary tracks
#sun-ms
plt.plot(logTef100,logL100,'#ADFF2F',linewidth=3.0)

#sun-agb
plt.plot(logTef_agb100,logL_agb100,'#ADFF2F',linestyle='--',linewidth=4.0)
#sun-hb
plt.plot(logTef_hb100,logL_hb100,'#ADFF2F',linestyle='-.',linewidth=4.0)



#vega-ms
plt.plot(logTef210,logL210,'#ADFF2F',linewidth=3.0)
#vega-hb
plt.plot(logTef_hb200,logL_hb200,'#ADFF2F',linewidth=3.0)
#vega-agb
plt.plot(logTef_agb210,logL_agb210,'#ADFF2F',linewidth=3.0)


#reguls-ms
plt.plot(logTef350,logL350,'#ADFF2F',linewidth=3.0)
#regulus-agb
plt.plot(logTef_agb350,logL_agb350,'#ADFF2F',linewidth=3.0)





plt.text(logTef100[-1],logL100[-1]+0.5,'z017y26-1M',fontweight='bold',fontsize=15,color='g',rotation=45)
plt.text(logTef210[-1]-0.1,logL_agb210[-1]+1,'z017y26-2.10M',fontweight='bold',fontsize=15,color='g',rotation=45)
plt.text(logTef350[-1],logL_agb350[-1]+1,'z017y26-3.50M',fontweight='bold',fontsize=15,color='g',rotation=45)






plt.show()


/usr/local/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:37: RuntimeWarning: divide by zero encountered in true_divide

About CMD and Evolutionary

The star masses are critical to star evolutionary,actually, lifetimes of stars are as  function of their masses.The more massive stars are, more energy they lose.

The massive stars which are very luminous and thus have very rapid stellar winds, lose mass so rapidly due to radiation pressure that they tend to strip off their own envelopes before they can expand to become red supergiants, and thus retain extremely high surface temperatures (and blue-white color) from their main-sequence time onwards.

A 0.1 solar mass star's lifetime is about 1000 billion years,but a 60 solar mass can only exist 3 million years!The most massve stars' evolutionary stage are mainly out of mainsequence,for example,the superblue star Alnilam,still on main sequence,with 30 solar mass,its age is 5.7 Myrs,very young star.

a)

Large mass stars evolve rapidly and die young,with brilliant explosive ends,on the contrary,low mass star evolve slowly,die quitely,with gloomy ends. 

b)

The fraction of the brightest 20 stars are on the main sequence versus in a post-main sequence phase is 6:4,which is smaller than the real ratio 9:1. The resonable explanation is the we obtained small samples.

Problem 3

a)


In [6]:
import matplotlib.pyplot as plt
from astropy.io import fits
import numpy as np

plt.figure(1)
plt.figure(figsize=(30,20))
ax1=plt.subplot(111)
hdu=fits.open('./lumfun.fit')
hdudat=hdu[1].data
mask1=(hdudat['Plx']>1)#the distance out to stars whoes distance to the Sun is smaller than 1000pc.
mask2=(hdudat['Plx']<1000)# greater than 1pc

hdudat=hdudat[mask1*mask2]

pc=list(map(abs,1000/hdudat['Plx']))# abs of Plx(mas)--distance(pc)
Mav=hdudat['Vmag']+5*(1-np.log10(pc))


#plot histgrams and its envelope curve.
returns=plt.hist(Mav,bins='auto',align='left')
lum=returns[1][:-1]
num=returns[0][:]
plt.plot(lum,num,'k--')
#plt.xlim(-10,1)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.xlabel('absolute magitude',fontsize=20)
plt.ylabel('number of stars',fontsize=20)

plt.show()


<matplotlib.figure.Figure at 0x1131c3d30>
The Hippacos data selection effect!

There is very strong selection effect in hippacos data.So,the luminosity function is almost a normal distribution... Because of limit of time,I didn't find complete star catalogue..