We have presented the procedure which can obtain dynamical properties such as radius and mass of cluster from the kinematics and CCD photometry. The membership stars of NGC 7142 were selected from the proper motions of the open cluster that we measured in 4 epochs of observations. The proper motion of stars was corrected from the effect of parallax. The membership probabilities belonging to NGC 7142 were derived from maximum-likelihood method with proper motion error and there are 774 stars have been selected as membership stars. The interstellar reddening and extinction correction were found by fitting the TCD red giant branch line and the result showed the average reddening E(B−V) = 0.42, E(U−B) = 0.45 and mean extinction Av = 1.30. while Z = 0.010 or [M/H] = -0.28. By fitting the standard evolution curve we got the distance modulus (m−M) = 11.60 ± 0.1 and logarithm age = 9.5 ~ 9.8, or distance ~ 2100 pc while age = 3.2 ~ 6.3 Gyr. The tidal radius from King model fitting is 36.37 arcmin, or about 22 pc. The mass of NGC 7142 is about 4100 Mʘ.
There still have several problems going to be solved in this research.
Although the observations passed a decade is enough to derived proper motion, the 4 epochs are still a little insufficient to fit parallax properly.
Secondly, the NGC 7142 area have uneven distributing interstellar medium that cause the reddening correction more complex and difficult, but here we still using general correction method and that might have more uncertainty on reddening and extinction correction. From the V images of 2015O and 2015N, there are circular ring-shape nebula seems locate across the NGC 7142. The reddening correction should consider
52
these nebula but we are not able to correct it for now. The concept of correcting these reddening and extinction causing by the nebula is that to fit these ring and correct the reddening and extinction ring by ring separately. The King model profile fitting shows the data point are concentrate at core part and density at outer part of the cluster have large error. It is probably due to small field of view in 2003 and 2012 images.
The further observations with better FOV are required in the future work.
The procedure we used is currently aimed to one open cluster, NGC7142. In the future, we hope there will be more open cluster observation in different age, location, distance, metallicity or even spiral arm to know how the open cluster evolve, how they relate to the evolution of Milky Way Galaxy and mapping 3 dimensional distribution of open cluster.
53
References
Abell, G. O., 1959, ASPL, 8, 121
Balaguer-Núňez, L., Tian, K. P. & Zhao J. L., 1998, A&ASS, 133, 387 Blanco, V, M., 1956, ApJ, 123, 64
Cardelli, J. A., Clayton, G. C. & Mathis, J. S., 1989, 345, 245 Crinklaw, G & Talbert, F. D., 1991, PASP103, 536
Feast, M. & Whitelock, P., 1997, MNRAS, 291,683 Hiltner, W. A. & Johnson, H. L., 1956, ApJ, 124, 367 Hoag, A. A. et al., 1961, PUSNO, 13, 343
Hsieh, H.-Y., 2013, NTNU master thesis Huang, S.-T., 2007, NTNU master thesis
Jennens, P. A. & Helfer, H. L., 1975, MNRAS, 172, 681 Johnson, H. L. & Morgan, W. W., 1953, ApJ, 117, 313 Johnson, H. L. & Morgan, W. W., 1955, ApJ, 122, 142 Johnson, H. L. et al., 1961, LowOB, 5, 133
King, I. R., 1962, AJ, 67, 471 King, I. R., 1966, AJ, 71, 64
Lasker, B. M. et al, 1990, AJ, 99, 2019 Lindholm, E. H., 1957, ApJ, 126, 588
Morgan, W. W., Harris, D. L. & Johnson, H. L., 1953, ApJ, 118, 92 Oort, F. H., 1927, BAN, 3, 275
Palmer, J. & Davenhall, 2011, StarC, 6 Punanova, A. F. et al. 2012, AASP, 2, 11 Sanders, W. L., 1971, A&A, 14, 226
Savage, B. D. & Mathis, J. S., 1979, ARA&A, 17, 73 Sharov, A. S., 1965, SvA, 8, 780
Sharov, A. S., 1968, SvA, 12, 550
54
Slovak, M. H., 1977, AJ, 82, 818
Smart, W. A., 1949, Text-book on Spherical Astronomy, Cambrige universty press, UK
Straižys, V., Sūdžius, J & Kurilinenė, G. 1976, A&A, 50, 431 Straižys, V. et al, 2014, MNRAS, 437, 1628
van den Bergh, S. 1962, JRASC, 56, 41
van den Bergh, S. & Heeringa, R., 1970, A&A, 9, 209
Vasilevskis, S., Klemola, A. & Perston, G., 1958, AJ, 63, 387
Vasilevskis, S., Sanders, W. L. & van Altena, W. F., 1965, AJ, 70, 806 Wang, J.-H., 2003, NTNU master thesis
Wu, C.-K., 2004, NTNU master thesis
Zhao, J. L. & He, Y. P., 1990, A&A, 237, 54 WEBDA, 2006, http://www.univie.ac.at/webda/
CMD 2.8, 2016, http://stev.oapd.inaf.it/cgi-bin/cmd
55
Appendixes
Appendix A: IRAF information
We use the DAOPHOT package of IRAF program to progress the position and photometry data. There are the parameters in DAOPHOT of each image frame.
File 1. 20030807V.fits IRAF file name: 20030807V
IMEXAM: STDDEV: 6.65 6.78 7.69 6.98 7.56 7.21 6.88 7.04 7.58 6.42 avg 7.04 ≈ 7.0 IMEXAM: FWHM: 2.54 2.56 2.31 2.75 2.67 2.58 2.9 2.47 2.43 2.79 avg 2.60
DAOFIND:
FWHMPSF = 2.6 SIGMA = 7.0 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 10 DANNULUS = 10 APERTURE = 2.6 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 11 FITRAD = 6
20 psf stars. par1: 1.250444. par2: 1.311973 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 10 DANNULUS = 10 MAXGROUP = 60 txdump > 20030807V.als id,xcen,ycen,mag,merr
56 file 2. 20120708V.fits
IRAF file name: 20120708V
IMEXAM: STDDEV: 35.4 49.6 42.5 35.8 31.2 47.5 39.2 34.5 42.1 46.7 avg 40.45 ≈ 40.5 IMEXAM: FWHM: 2.84 2.56 3.11 2.89 2.97 2.58 2.94 2.47 3.01 2.68 avg 2.8
DAOFIND:
FWHMPSF = 2.8 SIGMA = 40.5 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 12 DANNULUS = 12 APERTURE = 2.8 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 12 FITRAD = 6
18 psf stars. par1: 1.441058 par2: 1.540893 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 12 DANNULUS = 6 MAXGROUP = 60 txdump > 20120708V.als id,xcen,ycen,mag,merr
57 file 3. 20151014V.fits
IRAF file name: 20151014V
IMEXAM: STDDEV: 43.5 39.6 42.5 45.8 41.2 47.5 39.2 54.2 42.1 46.7 avg 44.5 ≈ 45 IMEXAM: FWHM: 2.54 2.33 2.22 2.47 2.41 2.36 2.42 2.51 2.36 2.34 avg 2.39
DAOFIND:
FWHMPSF = 2.39 SIGMA = 45 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 10 DANNULUS = 10 APERTURE = 2.4 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 11 FITRAD = 7
14 psf stars. par1: 1.278506 par2: 1.17063 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 10 DANNULUS = 10 MAXGROUP = 60 txdump > 20151014V.als id,xcen,ycen,mag,merr
58 file 4. 20151014B.fits
IRAF file name: 20151014B
IMEXAM: STDDEV: 30.2 36.5 25.0 26.7 31.2 27.6 32.6 29.8 24.1 35.1 avg 29.9 ≈ 30 IMEXAM: FWHM: 2.54 2.77 2.36 2.76 2.84 2.68 2.73 2.61 2.46 2.34 avg 2.61
DAOFIND:
FWHMPSF = 2.61 SIGMA = 30 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 12 DANNULUS = 7 APERTURE = 2.6 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 12 FITRAD = 7
17 psf stars. par1: 1.393874 par2: 1.326062 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 12 DANNULUS = 7 MAXGROUP = 60 txdump > 20151014B.als id,xcen,ycen,mag,merr
59 file 5. 20151014U.fits
IRAF file name: 20151014U
IMEXAM: STDDEV: 20.3 24.7 25.0 19.3 21.3 24.2 18.9 19.6 24.1 21.5 avg 21.9 ≈ 22 IMEXAM: FWHM: 2.54 2.77 2.36 2.76 2.84 2.68 2.73 2.61 2.46 2.34 avg 3.2
DAOFIND:
FWHMPSF = 3.2 SIGMA = 22 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 18 DANNULUS = 6 APERTURE = 3.2 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 15 FITRAD = 9
15 psf stars. par1: 1.940688 par2: 1.454519 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 18 DANNULUS = 6 MAXGROUP = 60 txdump > 20151014U.als id,xcen,ycen,mag,merr
60 file 6. 20151014U.fits
IRAF file name: 20151014U
IMEXAM: STDDEV: 20.3 24.7 25.0 19.3 21.3 24.2 18.9 19.6 24.1 21.5 avg 21.9 ≈ 22 IMEXAM: FWHM: 2.98 3.31 3.57 3.34 2.84 3.44 2.73 2.89 3.46 3.61 avg 3.2
DAOFIND:
FWHMPSF = 3.2 SIGMA = 22 THRESHOLD = 4 PHOT:
CALGORITHM = none SALGORITHM = mode ANNULUS = 18 DANNULUS = 6 APERTURE = 3.2 PSF:
FUNCTION = gauss VARORDER = 1 PSFRAD = 15 FITRAD = 9
15 psf stars. par1: 1.940688 par2: 1.454519 ALLSTAR:
RECENTER = yes GRPSKY = yes FITSKY = yes ANNULUS = 18 DANNULUS = 6 MAXGROUP = 60 txdump > 20151014U.als id,xcen,ycen,mag,merr
61
Appendix B: Python source code
Here we present the Python code made by our own. There are 6 programs and they are sorted by the order of the sequence.
1. photometry.py
fin.append(open(path+"20151014"+fltr[i]+".als", "r").readlines()) star.append([])
for j in range(len(fin[i])):
star[i].append(np.array(fin[i][j].split(),dtype=np.float))
fin=[]
ref=[]
for i in range(len(fltr)):
fin.append(open(path+"20151014BV_"+fltr[i]+".photref", "r").readlines()[2:]) ref.append([])
for j in range(len(fin[i])):
ref[i].append(np.array([fin[i][j].split()[0],fin[i][j].split()[1]],dtype=np.float))
#match ref in star and get info rstar=[]
62
*bvcalv[0]+bvcalv[1])+Vfit[1]), phot[0],"k.",label=r"$V_{cal}$") plt.xlabel(r"$v\ &\ V_{cal}$",fontsize=14)
plt.text(1, 5.3, "V-v = %.3f (B-V) + %.3f"%(Vfit[0], Vfit[1]), color="r") plt.subplot(212)
plt.plot(rstar[1][2]-rstar[0][2],phot[0]-rstar[0][2], "ko") plt.plot(range(3),vfit[0]*np.array(range(3))+vfit[1],"r-") plt.xlabel(r"$(b-v)$",fontsize=14)
plt.ylabel(r"$V - v$",fontsize=14)
plt.text(1, 5.3, "V-v = %.3f (b-v) + %.3f"%(vfit[0], vfit[1]), color="r") plt.tight_layout()
plt.text(0.25, 5.3, "B-b = %.3f (B-V) + %.3f"%(Bfit[0], Bfit[1]), color="r") plt.subplot(212)
plt.plot(rstar[1][2]-rstar[0][2],phot[1]-rstar[1][2], "ko") plt.plot(range(3),bfit[0]*np.array(range(3))+bfit[1],"r-") plt.xlabel(r"$(b-v)$",fontsize=14)
plt.ylabel(r"$B - b$",fontsize=14)
plt.text(0.25, 5.3, "B-b = %.3f (b-v) + %.3f"%(bfit[0], bfit[1]), color="r")
63
fout.writelines("%4d %5.2f %5.2f\n"%(bvmatch[0][5][i], bcal[i], vcal[i])) fout.close()
fout=open(path2+"Thesis_bvxy.txt", "w")
fout.writelines("%4s %7s %7s %7s\n"%("ID", "V", "X", "Y")) for i in range(len(vcal)):
fout.writelines("%4d %7.3f %7.3f %7.3f\n"%(bvmatch[0][5][i], vcal[i], bvmatch[0][1][i], bvmatch[0][2][i]))
plt.title("Stars field of NGC7142 in 2015O image",fontsize=20) for i in range(len(vcal)):
plt.plot(bvmatch[0][1][i], bvmatch[0][2][i], "ko", ms=10**((-0.1)*vcal[i]+2.5))
64
fin.append(open(path+"20151014"+fltr[i]+".als", "r").readlines()) star.append([])
for j in range(len(fin[i])):
65
star[i].append(np.array(fin[i][j].split(),dtype=np.float))
fin=[]
ref=[]
for i in range(len(fltr)):
fin.append(open(path+"20151014"+fltr[i]+".photref", "r").readlines()[2:]) ref.append([])
#match ref in star and get info rstar=[]
66
plt.text(1.5, 1.7, "U-u = %.3f (U-B) + %.3f"%(Ufit[0], Ufit[1]), color="r") plt.subplot(212)
plt.plot(rstar[1][2]-rstar[0][2],phot[2]-rstar[1][2], "ko") plt.plot(range(3,7),ufit[0]*np.array(range(3,7))+ufit[1],"r-") plt.xlabel(r"$(u-b)$",fontsize=14)
plt.ylabel(r"$U - u$",fontsize=14)
plt.text(4.5, 1.7, "U-u = %.3f (u-b) + %.3f"%(ufit[0], ufit[1]), color="r") plt.tight_layout()
fout.writelines("%4s %5s %5s %5s\n"%("ID", "U", "B", "V")) for i in range(len(ucal)):
67
fout.writelines("%4d %5.2f %5.2f %5.2f\n"%(ubvmatch[0][5][i], ucal[i], bcal[i], vcal[i]))
fout.close()
fout=open(path2+"Thesis_ubvxy.txt", "w")
fout.writelines("%4s %7s %7s %7s\n"%("ID", "V", "X", "Y")) for i in range(len(vcal)):
fout.writelines("%4d %7.3f %7.3f %7.3f\n"%(ubvmatch[0][5][i], vcal[i], ubvmatch[0][1][i], ubvmatch[0][2][i]))
68 for i in range(len(vcal)):
if ucal[i]-bcal[i] < -0.2:
plt.plot((ubvmatch[2][1][i]-coeffx[1])/coeffx[0], (ubvmatch[2][2][i]-coeffy[1])/coeffy[0],"ro",fillstyle="none",mew=1,ms=6)
plt.xlim(0,2048) plt.ylim(0,2048)
plt.savefig(path2+"position_ub_lt-0.2.png")
69
2. match.py
# -*- coding: utf-8 -*-
"""
This program read the IRAF all star file and reference stars list, match stars and calculate the coordinates coversion coefficient
@author: Yi-Hsiang Hsu
#allstar files, reference star file name
iraflist=["20030807V","20120708V","20151014V","20151111V"]
#obs data in JD
obstime=np.array([2452858.95903,2456117.95139,2457310.0609028,2457338.17847222]) imgsize=np.array([1340,1300]) #XYsize
pixscl=0.51 #pixel scale in arcsec, in the ref frame
imgctr=np.array([SkyCoord("21h45m13s", "+65d46m25s").ra.degree,\
fin.append(open(path+"\\"+dirlist[i]+"\\"+iraflist[i]+".ref", "r"))
#read the files
GMLS=(280.46646+Jc*(36000.76983 + Jc*0.0003032))%360. #Geom Mean Long Sun (deg) GMAS=357.52911+Jc*(35999.05029 - 0.0001537*Jc) #Geom Mean Anom Sun (deg)
SEqC=np.sin(np.deg2rad(GMAS))*(1.914602-Jc*(0.004817+0.000014*Jc))+\
np.sin(np.deg2rad(2*GMAS))*(0.019993-0.000101*Jc)+\
np.sin(np.deg2rad(3*GMAS))*0.000289 #Sun equation of center STL=GMLS+SEqC #sun true long (deg)
STA=GMAS+SEqC #sun true anom
Ecc=0.016708634-Jc*(0.000042037 + 0.0000001267*Jc) #Eccent Earth Orbit AU=(1.000001018*(1-Ecc*Ecc))/(1+Ecc*np.cos(np.deg2rad(STA)))
70
#read in star files fin=[]
allstar=[]
for i in range(len(iraflist)):
fin.append(open(path+"\\"+dirlist[i]+"\\"+iraflist[i]+".als", "r")) allstar.append(fin[i].readlines())
xref[i].append(float(datalist[i][j].split()[0]))
yref[i].append(float(datalist[i][j].split()[1])) xref[i]=np.array(xref[i])
yref[i]=np.array(yref[i])
#define coefficients calculater
#calculate the convert coefficient with least-square def coeff(x,y,cnt):
coeffx=[]
coeffy=[]
for i in range(len(x)):
coeffx.append(np.linalg.lstsq(np.array([x[i],y[i],\
np.ones(len(x[i]))]).T,x[1])[0])
coeffy.append(np.linalg.lstsq(np.array([y[i],x[i],\
np.ones(len(x[i]))]).T,y[1])[0])
#write the result in to a file
fout=open(path+"coorcal\\Thesis_coeff"+str(cnt)+".txt","w") fout.writelines("%8s %10s %10s %10s %10s %10s %10s\n"%\
("Obs","x=ax'","+by'","+c","y=dy'","+ex'","+f"))
plt.title("Reference stars in 2012 image",fontsize=18) plt.plot(x[1],y[1],"ro",fillstyle="none",mew=1.5,ms=8)
plt.title(r"$X_{%s}\ --\ X_{%s}$"%(iraflist[1],iraflist[i])) plt.xlim(0,1340)
plt.subplot(222)
plt.plot(y[1],y[i], "ko")
plt.title(r"$Y_{%s}\ --\ Y_{%s}$"%(iraflist[1],iraflist[i]))
71 plt.xlim(0,1300)
plt.subplot(223)
plt.title(r"$Residual\ of\ X_{%s}\ to\ X_{%s}$"%\
(iraflist[1],iraflist[i]))
plt.plot(x[1],coeffx[i][0]*x[i]+coeffx[i][1]*y[i]+coeffx[i][2]-x[1],\
"ko")
plt.plot(y[1],coeffy[i][0]*y[i]+coeffy[i][1]*x[i]+coeffy[i][2]-y[1],\
"ko")
plt.plot([0,1300],[0,0], "b-") plt.xlim(0,1300)
plt.savefig(path+"coorcal\\"+iraflist[i]+"coor_cal"+str(cnt)+".png") return coeffx,coeffy
#define match star subrutine def matchstar(cx,cy,cnt):
#convert coordinate of each star list
xycoor=[] #original coordinates xycotf=[] #transfer coordinates for i in range(len(iraflist)):
xycotf.append([cx[i][0]*x+cx[i][1]*y+cx[i][2],cy[i][0]*y+cy[i][1]*x+\
cy[i][2]])
72
star.append([float(fin[i].split()[3]),float(fin[i].split()[4]),\
#calculate dispacemaent of each stars dx=[]
dy=[]
for i in range(len(star)):
dx.append([star[i][0]-star[i][2],star[i][4]-star[i][2],star[i][6]-star[i][2]])
dy.append([star[i][1]-star[i][3],star[i][5]-star[i][3],star[i][7]-star[i][3]])
dx,dy=-np.array(dx)*pixscl*1000, np.array(dy)*pixscl*1000
73
pmpxx.append(np.linalg.lstsq(X.T,dx[i])[0]) pmpxy.append(np.linalg.lstsq(Y.T,dy[i])[0]) varx.append(np.linalg.lstsq(X.T,dx[i])[1][0]) vary.append(np.linalg.lstsq(Y.T,dy[i])[1][0]) pmpxx=np.array(pmpxx).T
plt.suptitle("Histogram of parallax", fontsize=16) plt.subplot(211)
plt.subplots_adjust(top=0.9)
plt.xlabel(r"$\pi_{\alpha}\ (mas)$", fontsize=14)
74 plt.ylabel("N star", fontsize=12)
plt.hist(pmpxx[1],bins=(np.linspace(-2000,2000,80))) plt.subplot(212)
plt.xlabel(r"$\pi_{\delta}\ (mas)$", fontsize=14) plt.ylabel("N star", fontsize=12)
plt.hist(pmpxy[1],bins=(np.linspace(-2000,2000,80))) plt.xlabel(r"$Error (mas)$", fontsize=14)
plt.ylabel("N star", fontsize=12)
plt.savefig(path+"coorcal\\errhist"+str(cnt)+".png")
star.append([float(fin[i].split()[3]),float(fin[i].split()[4]),\
#calculate dispacemaent of each stars dx=[]
dy=[]
for i in range(len(star)):
dx.append([star[i][0]-star[i][4],star[i][2]-star[i][4],star[i][6]-star[i][4]])
dy.append([star[i][1]-star[i][5],star[i][3]-star[i][5],star[i][7]-star[i][5]])
dx,dy=-np.array(dx)*pixscl*1000, np.array(dy)*pixscl*1000
75 plt.xlabel(r"$Error (mas)$", fontsize=14)
plt.ylabel("N star", fontsize=12)
plt.savefig(path+"coorcal\\pmerrhist"+str(cnt)+".png") return star, varx, vary
#find new coordinates reference star def findref(m,cnt): float(fin[i].split()[2]),float(fin[i].split()[3])]) #pick new ref
76
77
if len(np.where(bvid == id15[i])[0]) > 0 and np.abs(pmaall[i]) < 100 and np.abs(pmdall[i]) <100:
78
plt.title("VPD of NGC 7142 region (selected)") plt.xlim(-50,50)
79
if sum_sigma_m > 0.01: sigma_m=sigma_m+0.01 if sum_sigma_m < -0.01: sigma_m=sigma_m-0.01
if sum_ux_0 > 0.01: ux_0=ux_0+0.01 if sum_ux_0 < -0.01: ux_0=ux_0-0.01 if sum_uy_0 > 0.01: uy_0=uy_0+0.01 if sum_uy_0 < -0.01: uy_0=uy_0-0.01
if sum_sigma_bx > 0.01: sigma_bx=sigma_bx+0.01 if sum_sigma_bx < -0.01: sigma_bx=sigma_bx-0.01 if sum_sigma_by > 0.01: sigma_by=sigma_by+0.01 if sum_sigma_by < -0.01: sigma_by=sigma_by-0.01
iterpro.append([i+1,n_m,ua_0,ud_0,sigma_m,ux_0,uy_0,sigma_bx,sigma_by]) itersum.append([i+1,sum_n_m,sum_ua_0,sum_ud_0,sum_sigma_m\
,sum_ux_0,sum_uy_0,sum_sigma_bx,sum_sigma_by]) iterpro=np.array(iterpro).T
itersum=np.array(itersum).T
80
81
plt.plot([-1.57*50, 1.57*50], [-0.76*50, 0.76*50], "g-") #galatic plane plt.arrow(-20*np.cos(-thita)+pmx0, -20*np.sin(-thita)+pmy0\
, -20*np.cos(-thita), -20*np.sin(-thita), width=0.1,color="brown")
#plt.text(25, -2, "Rotated pricipal axis", color="brown", rotation=-10) plt.arrow(-55.345000*0.3, -35.768944*0.3, -55.345000*0.35, -35.768944*0.35\
+((pmdr-ud_0)**2)/(sigma_m**2+err**2)))
alpha=np.exp(-0.5*(((pmar-ux_0)**2)/(sigma_bx**2+err**2)\
+((pmdr-uy_0)**2)/(sigma_by**2+err**2)))
psi=(((len(pma)-n_m)/(2*np.pi*np.sqrt(sigma_bx**2+err**2)\ plt.title("Histogram of membership probability") plt.xlabel("Percentage")
for i in np.where(prob >= 85)[0]:
plt.plot(pma[i],pmd[i],"r.", ms=3) plt.xlim(-50,50)
for i in np.where(prob < 85)[0]:
82 plt.plot(pma[i],pmd[i],"k.", ms=3) plt.savefig(path+"VPD_f.png")
plt.clf()
plt.gca().set_aspect('equal', adjustable='box') plt.xlim(-50,50)
plt.ylim(-50,50)
plt.xlabel(r"$\mu_{\alpha}^{*}\ (mas/yr)$",fontsize=16) plt.ylabel(r"$\mu_{\delta}\ (mas/yr)$",fontsize=16) for i in np.where(prob >= 85)[0]:
plt.plot(pma[i],pmd[i],"r.", ms=3) plt.savefig(path+"VPD_m.png")
fout=open(path+"Thesis_prob.txt","w") for i in range(len(ids)):
fout.writelines("%4d %5.2f %5.2f %5.2f %4d\n"%\
(ids[i],prob[i], vs[i], bvs[i], id15s[i])) fout.close()
plt.clf() plt.cla()
plt.figure(figsize=[8,8]) plt.plot(bvs, vs, "k.", ms=2)
plt.title("CMD in NGC 7142 region",size=14) plt.xlabel("B-V",fontsize=12)
plt.ylabel("V",fontsize=12) plt.tick_params(labelsize=12) plt.ylim(22,10)
plt.xlim(0,2.0)
plt.savefig(path+"CMD_sel.png")
83
float(uphot[i].split()[1])-float(uphot[i].split()[2]) v,bv,ub=np.array(v),np.array(bv),np.array(ub)
84 +4*ks*(-kr*bv+fit[0]*bv+fit[1]-isofit[1])))/(2*ks) ub0=ub-(bv-bv0)*(kr+ks*bv)
85
plt.title("TCD of NGC7142, Z=0.010, age=5 Gyr") plt.plot(fitbv,fitub)
fout.writelines("%s %5.2f %5.2f %5.2f\n"%\
(phot[i].strip("\n"), v0[i],bv0[i],ub0[i])) fout.close()
86 age.append(float(data[i].split()[1])) isou.append(float(data[i].split()[8])) isob.append(float(data[i].split()[9])) isov.append(float(data[i].split()[10])) stage.append(int(data[i].split()[23])) age=np.array(age)
isou=np.array(isou) isob=np.array(isob) isov=np.array(isov)
stage=np.array(stage,dtype=np.int) isobv,isoub=isob-isov,isou-isob plt.clf()
dm=11.6
for i in range(len(bv0)):
if prob[i] > 85:
plt.plot(bv0[i],v0[i],"k.",ms=4) for i in [950,965,980]:
fitv=[]
fitbv=[]
for j in range(len(isov)):
if age[j] == round(i*.01,4) : fitv.append(isov[j]) fitbv.append(isobv[j]) fitv=np.array(fitv)
fitbv=np.array(fitbv)
plt.title("CMD of NGC7142, Z=0.010,(m-M)=%5.2f"%(dm)) plt.plot(fitbv,fitv+dm,label="age = %.2f"%(i*0.01),lw=1.5) plt.ylim(22,10)
plt.xlim(-0.5,2)
plt.xlabel(r"$(B-V)_0$",fontsize=14) plt.ylabel(r"$V_0$",fontsize=14) plt.tick_params(labelsize=12) plt.legend(loc=2,fontsize=10) plt.savefig(path+"CMD_fit.png")
print 10**((dm+5.1)/5),10**((dm+5)/5),10**((dm+5-.1)/5)
87
5. king.py
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 20 13:33:47 2016
@author: Yi-Hsiang Hsu
plt.plot(np.log10(rrc[:r]), np.log10(f), "k-") plt.text(max(np.log10(rrc[:r]))+0.025 \
,min(np.log10(f)), str(np.log10(rtrc[i]))) plt.xlabel(r"$log\ r/r_c$",fontsize=12)
88
coeff=open(path+"\\coorcal\\Thesis_coeff10.txt","r").readlines()[3]
fit=[float(coeff.split()[1]),float(coeff.split()[2])\
,float(coeff.split()[3]),float(coeff.split()[4])\
,float(coeff.split()[5]),float(coeff.split()[6])]
xcal=fit[0]*x+fit[1]*y+fit[2]
89 from scipy.stats import norm
from scipy.optimize import curve_fit
#plot ring plt.clf()
#do gauss fitting to find center of cluster cx=norm.fit(xm)[0]
np.sin(cir)*(r[i])*(60/0.512)+cy, "b-") plt.xlim(0,1340)
90
#plot datapoint on king standard curve plt.clf()
*(((1/np.sqrt(1+rrc[rad]**2))-(1/np.sqrt(1+rtrc[i]**2)))**2)) rad=rad+1
plt.plot(np.log10(rrc[:rad]), np.log10(f),"-",\
label="log rt/rc = %.1f"%(np.log10(rtrc[i]))) plt.xlabel(r"$log\ r/r_c$",fontsize=12)
91
#do gauss fitting to find center of cluster cx=norm.fit(xs)[0]
92 plt.gca().set_aspect('equal', adjustable='box') for i in range(len(xs)):
plt.plot(xs[i],ys[i], "ko", ms=10**((-0.1)*vs[i]+2.0)) plt.plot(cx,cy,"b+")
np.sin(cir)*(r[i])*(60/0.512)+cy, "b-") plt.xlabel("X pix",fontsize=12)
fout.writelines("%d %.2f %.2f %d %.3f %.3f\n"%(rad+1,r[rad]-rw,r[rad],sc[rad],area[rad],d[rad]))
93
*(((1/np.sqrt(1+rrc[rad]**2))-(1/np.sqrt(1+rtrc[i]**2)))**2)) rad=rad+1
plt.plot(np.log10(rrc[:rad]), np.log10(f),"-",\
label="log rt/rc = %.2f"%(np.log10(rtrc[i]))) plt.xlabel(r"$log\ r/r_c$",fontsize=12)
94
for i in np.where(prob >= 85)[0]:
vm.append(v[i])
for i in np.where(prob < 85)[0]:
vf.append(v[i])
plt.title("Membership star CMD in NGC 7142 region",size=14) plt.xlabel("B-V",fontsize=12)
95 plt.plot(bvm, vm, "r.", ms=3)
plt.title("Membership star CMD in NGC 7142 region",size=14) plt.xlabel("B-V",fontsize=12)
plt.title("Field star CMD in NGC 7142 region",size=14) plt.xlabel("B-V",fontsize=12) pmxid.append(int(fin[i].split()[4])) pmxid=np.array(pmxid)
96 for i in range(len(fin)):
x03.append(float(fin[i].split()[6])) y03.append(float(fin[i].split()[7])) xyid.append(int(fin[i].split()[8])) xyid=np.array(xyid)
fout=open(path+"Thesis_totallist.txt","w")
fout.writelines("%4s %8s %8s %8s %8s %8s %8s %8s %8s %8s\n" \ %("No.","x_12", "y_12", "pma", "pmd", "err"\
,"U0", "B0", "V0", "prob")) count=1
for i in range(len(xyid)):
j,k,l=[],[],[]
if xyid[i] < 9999:
j=np.where(pmxid == xyid[i])[0]
l=np.where(id15== xyid[i])[0]
k=np.where(id0 == xyid[i])[0]
if len(j) > 0 and len(l) > 0 and len(k) > 0:
j,k,l=j[0],k[0],l[0]
if u0[k] < 90 and b0[k] != np.nan:
fout.writelines("%4d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f
%8.2f %8.2f\n" \ %
(count,x03[i],y03[i],pma[j],pmd[j],err[j],u0[k],b0[k],v0[k],prob[l])) count=count+1
elif b0[k] != np.nan:
fout.writelines("%4d %8.2f %8.2f %8.2f %8.2f %8.2f %8s %8.2f %8.2f
%8.2f\n" \
%
(count,x03[i],y03[i],pma[j],pmd[j],err[j],"---",b0[k],v0[k],prob[l])) count=count+1 fout.close()
97
Appendix C: Stars List of NGC 7142
No. x_12 y_12 μα μδ ε U0 B0 V0 prob
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
No. x_12 y_12 μα μδ ε U0 B0 V0 prob
1091 649.88 907.28 13.45 3.82 3.66 12.81 12.96 12.58 3.69 1092 36.59 1113.93 11.37 6.86 7.13 16.54 16.37 15.63 75.85 1093 1244.45 1271.22 -11.32 11.17 6.81 --- 18.50 17.90 82.14 1094 156.00 53.19 22.59 -23.07 8.90 16.49 16.49 15.88 10.21 1095 736.43 224.78 -2.57 -0.61 2.76 13.54 12.36 10.81 97.78 1096 61.83 516.46 11.92 -3.93 5.68 16.37 16.31 15.72 61.61 1097 779.28 555.48 1.52 -21.96 7.82 --- 19.06 17.49 45.31 1098 784.21 601.80 -0.54 2.08 2.68 13.60 13.04 12.00 98.04 1099 108.16 903.84 8.29 2.88 4.46 16.86 17.01 16.38 79.20 1100 786.04 1007.27 -4.38 -8.45 4.97 14.07 14.04 13.22 84.93 1101 596.31 1179.04 -0.63 -1.62 1.77 13.98 12.87 11.65 98.88 1102 414.36 112.97 8.21 -7.71 5.52 17.31 17.56 17.05 69.96 1103 435.15 565.85 -4.61 -7.89 3.52 --- 18.88 18.20 70.11 1104 560.21 1212.06 4.35 4.01 5.84 13.30 12.88 12.91 91.20 1105 437.61 154.67 -0.17 -7.80 4.51 13.10 13.36 12.97 86.94 1106 350.35 602.27 33.87 -4.56 7.26 --- 13.48 12.50 0.34 1107 401.73 794.54 -2.68 -4.68 3.51 13.79 13.36 12.36 93.62 1108 1185.01 815.71 1.95 0.66 2.44 --- 17.96 17.13 98.02 1109 1107.06 1047.68 -3.44 -0.60 1.54 15.17 15.34 14.76 94.49 1110 462.93 109.34 -2.93 0.28 4.79 --- 18.46 17.94 95.78 1111 103.50 487.47 10.18 -3.41 5.11 14.95 15.00 14.36 68.60 1112 761.88 570.25 -0.94 3.65 2.19 13.17 11.88 10.45 96.18 1113 824.20 890.52 5.27 -4.42 3.50 12.64 12.91 12.56 83.51 1114 980.56 984.15 1.89 3.23 3.10 13.59 12.51 11.27 96.07 1115 123.00 1136.60 1.24 4.02 5.34 16.43 16.24 15.50 93.97 1116 994.54 1137.17 -6.71 0.95 4.65 15.99 16.39 15.74 93.59 1117 51.11 927.24 14.78 7.37 6.56 12.31 12.16 11.81 52.64 1118 320.35 735.57 -13.63 -11.85 6.57 11.71 11.93 11.60 64.66 1119 781.66 558.75 1.52 -21.96 7.82 --- 19.06 17.49 45.31 1120 208.58 247.90 15.29 -5.12 6.71 --- 18.69 17.96 50.12 1121 320.44 738.91 -13.63 -11.85 6.57 11.71 11.93 11.60 64.66 1122 552.69 797.14 -6.44 -5.05 3.63 16.82 17.19 16.46 84.00 1123 545.67 809.31 -2.54 4.33 4.58 12.98 12.33 11.13 94.96 1124 840.82 208.73 3.36 -5.59 9.39 11.51 11.78 11.20 87.03