Amasouverts¶

M67 = NGC¶

Data : vizier search M67 r = 30'
Isophotes :http://stev.oapd.inaf.it/cgi-bin/cmd Photometric system Gaia EDR3

In [27]:
import os
import pandas as pd
from matplotlib import pyplot as plt
import numpy as np

def ellipse (x0,y0,a,b):
    theta = np.linspace( 0 , 2 * np.pi , 360 )
    x = x0 + a*np.cos(theta)
    y = y0 + b*np.cos(theta)
    print(x)
    return x, y


def pc2al (d_pc):
    d_al = d_pc * 3.26
    return d_al

def plx2d_pc(plx):
    d_pc = 1/plx
    return d_pc


    



path = r'D:\DUAO\Amas'
os.chdir(path)
df0 = pd.read_csv("M67.csv",sep = ";")


_describe = 0
_display = 0

df1=df0.dropna(subset=['pmRA','pmDE',"Plx","Gmag","BPmag","RPmag"])
if _describe :display(df1.describe())
if _display : display(df1.sample(100))
In [28]:
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2)
plt.rcParams['font.size'] = 6

ax1.plot(abs(df1["pmRA"]),df1["e_pmRA"],'o',markersize = 1)
ax1.set_title("e = f(pmRA)")
ax2.plot(abs(df1["pmDE"]),df1["e_pmDE"],'o',markersize = 1)
ax2.set_title("e )= f(pmDE)")
ax3.plot(df1["Plx"],df1["e_Plx"],'o',markersize = 1)
ax3.set_title("e = f(Plx)")
ax4.plot(abs(df1["Gmag"]),df1["e_Gmag"],'o',markersize = 1)
ax4.set_title("e = f(Gmag")
Out[28]:
Text(0.5, 1.0, 'e = f(Gmag')
In [29]:
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2)
plt.rcParams['font.size'] = 6

ax1.hist(df1["e_pmRA"],bins = 30)
ax1.set_title("e = f(pmRA)")
ax2.hist(df1["e_pmDE"],bins = 30)
ax2.set_title("e = f(pmDE)")
ax3.hist(df1["e_Plx"],bins = 30)
ax3.set_title("e = f(Plx)")
ax4.hist(df1["e_Gmag"],bins = 30)
ax4.set_title("e = f(Gmag)")
Out[29]:
Text(0.5, 1.0, 'e = f(Gmag)')
In [30]:
display(df1["pmRA"].describe(),df1["pmDE"].describe(),df1["Plx"].describe(),df1["Gmag"].describe())

count    5286.000000
mean       -5.401819
std        10.173552
min      -229.438000
25%       -10.893750
50%        -4.727500
75%        -0.546750
max       163.074000
Name: pmRA, dtype: float64
count    5286.000000
mean       -4.804093
std        11.224085
min      -387.849000
25%        -5.264750
50%        -2.937500
75%        -1.490500
max        49.365000
Name: pmDE, dtype: float64
count    5286.000000
mean        0.979845
std         1.186904
min        -9.358600
25%         0.417275
50%         1.000250
75%         1.246400
max        22.044800
Name: Plx, dtype: float64
count    5286.000000
mean       17.861055
std         2.451274
min         7.543855
25%        16.211755
50%        18.553976
75%        19.882385
max        21.007792
Name: Gmag, dtype: float64
In [31]:
percentile = [.0,0.5,0.75,0.9,0.95,0.99]
display(df1["e_pmRA"].describe(percentiles = percentile),df1["e_pmDE"].describe(percentiles = percentile),df1["e_Plx"].describe(percentiles = percentile),df1["e_Gmag"].describe(percentiles = percentile))

count    5286.000000
mean        0.401080
std         0.508151
min         0.013000
0%          0.013000
50%         0.210000
75%         0.524000
90%         1.047000
95%         1.525750
99%         2.423500
max         3.210000
Name: e_pmRA, dtype: float64
count    5286.000000
mean        0.289461
std         0.359850
min         0.010000
0%          0.010000
50%         0.157000
75%         0.386750
90%         0.743000
95%         1.074000
99%         1.716300
max         2.416000
Name: e_pmDE, dtype: float64
count    5286.000000
mean        0.367427
std         0.442563
min         0.013300
0%          0.013300
50%         0.204200
75%         0.503575
90%         0.946600
95%         1.331075
99%         1.972975
max         3.576100
Name: e_Plx, dtype: float64
count    5286.000000
mean        0.004217
std         0.002141
min         0.002758
0%          0.002758
50%         0.003262
75%         0.004814
90%         0.007180
95%         0.008852
99%         0.011810
max         0.028395
Name: e_Gmag, dtype: float64
In [32]:
fig, ax = plt.subplots()
ax.hist(df1["Plx"],bins = 30)
Out[32]:
(array([1.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 4.000e+00, 2.000e+00,
        2.800e+01, 7.900e+01, 4.970e+02, 2.412e+03, 1.818e+03, 2.760e+02,
        9.100e+01, 3.300e+01, 2.200e+01, 7.000e+00, 4.000e+00, 3.000e+00,
        3.000e+00, 0.000e+00, 0.000e+00, 2.000e+00, 0.000e+00, 1.000e+00,
        0.000e+00, 2.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 1.000e+00]),
 array([-9.3586 , -8.31182, -7.26504, -6.21826, -5.17148, -4.1247 ,
        -3.07792, -2.03114, -0.98436,  0.06242,  1.1092 ,  2.15598,
         3.20276,  4.24954,  5.29632,  6.3431 ,  7.38988,  8.43666,
         9.48344, 10.53022, 11.577  , 12.62378, 13.67056, 14.71734,
        15.76412, 16.8109 , 17.85768, 18.90446, 19.95124, 20.99802,
        22.0448 ]),
 <BarContainer object of 30 artists>)
In [33]:
#Critères
Plx_max = 0.13
e_pmRA_max = 1
e_pmDE_max = 1
e_Gmag_max = 0.1
Gmag_max = 18
In [34]:
print(df1["Plx"].count())

Plx_min = 0
df2 = df1.loc[df1["Plx"] > 0.2]
print("Nb Stars Plx >" + str(Plx_min) +": ",df2["Plx"].count())

df2 = df2.loc[df2["Gmag"] < Gmag_max]
print("Nb Stars Gmag <" + str(Gmag_max) +": ",df2["Plx"].count())
df2 = df2.loc[df1["e_Plx"] < 0.2]
print(df2["Plx"].count())
df2 = df2.loc[df1["e_pmRA"] < 1]
print(df2["Plx"].count())
df2 = df2.loc[df1["e_pmDE"] < 1]
print(df2["Plx"].count())

df2["Grandeur"] = Gmag_max - df2.Gmag
stats = df2.Grandeur.describe()
print(stats)
delta_mag = stats[3] - stats[7]
df2.Grandeur = df2.Grandeur.apply(lambda x: 1.1**x)
stats = df2.Grandeur.describe()
print(stats)

df2.Plx.describe()

5286
Nb Stars Plx >0:  4449
Nb Stars Gmag <18:  2186
2162
2162
2162
count    2162.000000
mean        2.567459
std         1.837016
min         0.002792
25%         0.995230
50%         2.295671
75%         3.821713
max        10.456145
Name: Grandeur, dtype: float64
count    2162.000000
mean        1.297827
std         0.243906
min         1.000266
25%         1.099500
50%         1.244583
75%         1.439431
max         2.708993
Name: Grandeur, dtype: float64
Out[34]:
count    2162.000000
mean        1.230843
std         1.203307
min         0.202500
25%         0.784975
50%         1.129950
75%         1.201750
max        22.044800
Name: Plx, dtype: float64
In [35]:
fig, ax = plt.subplots()
#sizes = df2.Grandeur
ax.plot(df2["RA_ICRS"],df2["DE_ICRS"], 'o',markersize = 1)
ax.set_title("Champ M67 d = 30'",fontsize = 10)
ax.set_xlabel("RA [°]",fontsize = 10)
ax.set_ylabel("RA [°]",fontsize = 10)
#figsize, axis equal, Grandeur...
Out[35]:
Text(0, 0.5, 'RA [°]')
In [36]:
fig, (ax1,ax2) = plt.subplots(1,2)
fig.set_figheight(4)
fig.set_figwidth(10)


ax1.plot(df2["pmRA"],df2["pmDE"], 'o',markersize = 1)
#ax1.set_title("Champ M67 d = 30'",fontsize = 10)
ax1.set_xlabel("pmRA [mas]",fontsize = 10)
ax1.set_ylabel("pmDE [mas]",fontsize = 10)
ax1.set_xlim(-40, 20)
ax1.set_ylim(-40,20)
#ax.set_aspect( 1 )

x0 = -11
y0 = -2.9
a =2
b=2
c=1# ajsutement de la dimension du crop


theta = np.linspace( 0 , 2 * np.pi , 360 )
x = x0 + a*np.cos(theta)
y = y0 + b*np.sin(theta)

#ellipse (x0,y0,a,b)
ax1.plot(x,y)

ax2.plot(df2["pmRA"],df2["pmDE"], 'o',markersize = 1)
ax2.set_title("crop",fontsize = 10)
ax2.set_xlabel("pmRA [mas]",fontsize = 10)
ax2.set_ylabel("pmDE [mas]",fontsize = 10)
ax2.set_xlim(x0-c, x0+c)
ax2.set_ylim(y0-c,y0+c)

plt.plot([x0,x0],[y0-c,y0+c],color = "g",linewidth = 0.5)
plt.plot([x0-c,x0+c],[y0,y0],color = "g",linewidth = 0.5)
plt.text(x0+c-0.5, y0+c-0.1,"pmRA ~ "+ str(x0) + " mas")
plt.text(x0+c-0.5, y0+c-0.2,"pmDE ~ "+ str(y0) + " mas")
Out[36]:
Text(-10.5, -2.1, 'pmDE ~ -2.9 mas')
In [37]:
### add densité
pm_r = 1
df2["pm_d"] = np.sqrt((df2["pmRA"] - x0)**2 + (df2["pmDE"] - y0)**2)
df3= df2.loc[df2["pm_d"] < pm_r]
In [38]:
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2)
plt.rcParams['font.size'] = 6

ax1.hist(df3["pmRA"],bins = 30)
ax1.set_title("pmRA)")
ax2.hist(df3["pmDE"],bins = 30)
ax2.set_title("pmDE")
ax3.hist(df3["Plx"],bins = 30)
ax3.set_title("Plx)")
ax4.hist(df3["Gmag"],bins = 30)
ax4.set_title("Gmag")
Out[38]:
Text(0.5, 1.0, 'Gmag')
In [39]:
df3.describe()
Out[39]:
RA_ICRS e_RA_ICRS DE_ICRS e_DE_ICRS Plx e_Plx PM pmRA e_pmRA pmDE e_pmDE Gmag e_Gmag BPmag e_BPmag RPmag e_RPmag BP-RP _RAJ2000 _DEJ2000 Grandeur pm_d
count 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000 963.000000
mean 132.849697 0.038487 11.818487 0.020768 1.151889 0.046642 11.351873 -10.969757 0.048394 -2.911672 0.035827 14.910676 0.002813 15.425139 0.006548 14.264629 0.004770 1.160510 132.849747 11.818500 1.363698 0.269416
std 0.175392 0.028489 0.172438 0.015429 0.086118 0.034785 0.224924 0.224209 0.036219 0.228344 0.026905 1.847972 0.000169 2.043984 0.006709 1.698319 0.002121 0.425569 0.175393 0.172438 0.247898 0.175506
min 132.349112 0.011300 11.317111 0.006100 0.696600 0.013300 10.457000 -11.799000 0.013000 -3.814000 0.010000 8.222472 0.002758 9.105073 0.002802 7.300713 0.003779 -0.102004 132.349162 11.317124 1.000266 0.009055
25% 132.742186 0.016850 11.719098 0.009000 1.120450 0.020250 11.210500 -11.091000 0.021000 -3.051000 0.016000 13.492197 0.002763 13.796419 0.002880 13.034657 0.003825 0.788301 132.742237 11.719111 1.167364 0.148600
50% 132.845627 0.027400 11.815933 0.014700 1.152800 0.033400 11.345000 -10.974000 0.034000 -2.918000 0.025000 14.985489 0.002770 15.392708 0.003308 14.402118 0.003979 1.034775 132.845675 11.815946 1.332842 0.235926
75% 132.957440 0.050350 11.911461 0.027400 1.185400 0.060950 11.483000 -10.827000 0.064000 -2.783000 0.048000 16.376373 0.002815 17.046806 0.006722 15.609744 0.004882 1.462490 132.957490 11.911474 1.536704 0.335306
max 133.333803 0.167500 12.302884 0.090100 1.862500 0.193800 12.231000 -10.058000 0.209000 -1.942000 0.167000 17.997208 0.006879 19.217754 0.043493 17.028120 0.038080 2.323196 133.333852 12.302897 2.539324 0.982153

Tri complémentaire sur la base de la parallaxe

In [40]:
#Vérification Parallaxe
Plx_med = 1.153

d_mean = "{:.3}".format(1/Plx_med)
print(d_mean)

0.867
In [41]:
df4 = df3.loc[(df3["Plx"]>1.153-0.08) & (df3["Plx"]<1.153+0.08)]
In [42]:
df4.Plx.describe()
Out[42]:
count    763.000000
mean       1.154127
std        0.034779
min        1.074900
25%        1.130450
50%        1.153500
75%        1.177700
max        1.231900
Name: Plx, dtype: float64
In [43]:
fig, (ax1,ax2) = plt.subplots(1,2)
fig.set_figheight(4)
fig.set_figwidth(10)
#sizes = df2.Grandeur
ax1.plot(df4["RA_ICRS"],df4["DE_ICRS"], 'o',markersize = 1)
ax1.set_title("Champ M67 d = 30'",fontsize = 10)
ax1.set_xlabel("RA [°]",fontsize = 10)
ax1.set_ylabel("RA [°]",fontsize = 10)
ax1.set_xlim(132, 134)
ax1.set_ylim(11,12.5)
ax2.quiver(df4["RA_ICRS"],df4["DE_ICRS"],df4["pmRA"],df4["pmDE"],scale = 200)



#figsize, axis equal, Grandeur...
Out[43]:
<matplotlib.quiver.Quiver at 0x2dda71a16a0>
In [44]:
fig, (ax1,ax2,ax3) = plt.subplots(1,3)
fig.set_figheight(4)
fig.set_figwidth(10)
ax1.hist(df4["PM"],bins = 30)
ax2.hist(df4["pmRA"],bins = 30)
ax3.hist(df4["pmDE"],bins = 30)
Out[44]:
(array([  3.,   1.,   2.,   1.,   4.,   3.,   4.,  18.,  15.,  25.,  60.,
         68.,  84., 103.,  90.,  93.,  54.,  49.,  30.,  21.,   8.,   5.,
          8.,   6.,   1.,   3.,   2.,   0.,   1.,   1.]),
 array([-3.746     , -3.68586667, -3.62573333, -3.5656    , -3.50546667,
        -3.44533333, -3.3852    , -3.32506667, -3.26493333, -3.2048    ,
        -3.14466667, -3.08453333, -3.0244    , -2.96426667, -2.90413333,
        -2.844     , -2.78386667, -2.72373333, -2.6636    , -2.60346667,
        -2.54333333, -2.4832    , -2.42306667, -2.36293333, -2.3028    ,
        -2.24266667, -2.18253333, -2.1224    , -2.06226667, -2.00213333,
        -1.942     ]),
 <BarContainer object of 30 artists>)
In [45]:
pd.options.mode.chained_assignment = None  # default='warn'
df4["B-R"] = df4.BPmag - df4.RPmag
fig, ax = plt.subplots()
#sizes = df4.Grandeur
ax.plot(df4["B-R"],df4["Gmag"], 'o',markersize = 1)
ax.set_title("Diagramme HR M67",fontsize = 10)
ax.set_xlabel("B - R",fontsize = 10)
ax.set_ylabel("G",fontsize = 10)
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(18,8)



#figsize, axis equal, Grandeur...
Out[45]:
(18.0, 8.0)

Magnitude absolue Choix Extinction interstellaire A https://irsa.ipac.caltech.edu/applications/DUST/

In [46]:
A = 0.1
def abs_mag(m,d,A):
    # m = magnitude relative ; d= distance en parsec ; A = absorption interstellaire en magnitude
    M = m - 5*np.log10(d) + 5 - A
    return M

absMag = abs_mag(8,883,0)
print(absMag)

-1.7298035178878433
In [47]:
pd.set_option('display.max_columns', None)    
d = 883    
    
df4["absB"]  = df4["BPmag"].apply(lambda x: x - 5*np.log10(883) + 5 - .1)
df4["absG"]  = df4["Gmag"].apply(lambda x: x - 5*np.log10(883) + 5 - .1)
df4["absR"]  = df4["RPmag"].apply(lambda x: x - 5*np.log10(883) + 5 - .1)
In [48]:
pd.options.mode.chained_assignment = None  # default='warn'
df4["B-R"] = df4.absB - df4.absR
df4["absG"].describe()
Out[48]:
count    763.000000
mean       4.759315
std        1.738028
min       -1.607332
25%        3.512710
50%        4.711729
75%        6.023683
max        8.157878
Name: absG, dtype: float64
In [49]:
fig, ax = plt.subplots()
#sizes = df4.Grandeur
ax.plot(df4["B-R"],df4["absG"], 'o',markersize = 1)
ax.set_title("Diagramme HR M67 Magnitude absolue",fontsize = 10)
ax.set_xlabel("B - R",fontsize = 10)
ax.set_ylabel("G",fontsize = 10)
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(8,-2)


x = 0.74
ax.plot([x,x],[10,-5],color = "g",linewidth = 0.5)


#figsize, axis equal, Grandeur...
Out[49]:
[<matplotlib.lines.Line2D at 0x2dda76156a0>]

Ajustement isochrone

In [50]:
fig, ax = plt.subplots()
#sizes = df4.Grandeur
ax.plot(df4["B-R"],df4["absG"], 'o',markersize = 1)
ax.set_title("Diagramme HR M67 Magnitude absolue",fontsize = 10)
ax.set_xlabel("B - R",fontsize = 10)
ax.set_ylabel("G",fontsize = 10)
ax.set_xlim(-0.5, 2.5)
ax.set_ylim(8,-2)

df_iso0 = pd.read_csv(r"D:\DUAO\Amas\isochrones1.csv")


df_iso1 = df_iso0[df_iso0["logAge"] == 9.6]

df_iso1 = df_iso1[df_iso1["B-R"] < 2.5]
ax.plot(df_iso1["B-R"],df_iso1["Gmag"], 'o',markersize = 1)
Out[50]:
[<matplotlib.lines.Line2D at 0x2dda54ca730>]