Plotting Color Magnitude Diagrams

using plotly.graph_objects

The example uses HST data and MS isochrone from Bertelli

  1. Import packages
  2. Import data
    • Trim data if needed
    • Define B_V (color) and V (magnitude)
    • Define name of cluster/object, and x/yrange for the plot
    • Define ZAMS (zero age main sequence) data and mass (if needed)
  3. Plot and save
In [183]:
import pandas as pd
import numpy as np
from scipy.interpolate import interp1d,interp2d
import plotly.graph_objects as go

Simple CMD

Globular Clusters
In [184]:
xrange = [0.1,1.2]
yrange = [22,12]
name = 'M12'
data = pd.read_csv('./Labs/CMD/Globulars/M12/M12.dat',sep=' ',comment='#',index_col=False,header=None)
data.drop(columns=0,inplace=True)
columns = ['id','x','y','v','b','F555W','F439W','errf555w','errf439w','V_nr','B_nr','F555W_nr','F439W_nr','chi','sha','file_num']
data.columns = columns
V_0 = data['v']
B_0 = data['b']
B_V_0 = B_0 - V_0

# if name is 'M12':
dist_mod = 12.92
dist_mod_orig = 14.02
trim = 1.75
EB_V = 0.0
EB_V_orig = 0.0
# bert_ref = 10.06
# if name is 'M71':
# dist_mod = 12.0
# dist_mod_orig = 13.75 
# trim = 1.75
# EB_V = 0.06
# EB_V_orig = 0.0   
# bert_ref = 10.06

bert_ref = 6.6

AB_EB_V = 4.315
AV_EB_V = 3.315

data['B_V'] = B_V_0 - EB_V
data['V'] = V_0 - AV_EB_V * EB_V
data['M_V'] = V_0 - dist_mod

# Trimming data with Bertelli main sequence 
table = pd.read_csv('./Labs/CMD/Bertelli/table5.dat',comment='#',sep=' ')
table.drop(table.columns[0], axis=1, inplace=True)

masslim = 20
tab = table[(table['lgage']==bert_ref) & (table['mass'] < masslim)]
M_V_bert = tab['V_abs']
B_V_bert = tab['B-V']
interped = interp1d(x=B_V_bert,y=(M_V_bert + dist_mod + trim) , fill_value='extrapolate')
Vlim = interped(data['B_V'])
trimmed = data[(data['V'] < Vlim)]

# Ultimately, 
B_V = trimmed['B_V']
V = trimmed['V']

/Users/anikaslizewski/anaconda3/lib/python3.7/site-packages/scipy/interpolate/interpolate.py:609: RuntimeWarning:

divide by zero encountered in true_divide
In [185]:
fig = go.Figure()

fig.add_trace(go.Scatter(
        x=B_V,
        y=V,
        mode='markers',
        marker=dict(
            color='red',
            size=5,
            opacity=0.5,
            symbol='circle-open',line_width=1.5),name='')) 

# Layout:
fig.update_layout(title={'text': name,
            'y':0.86,
            'x':0.5},
        xaxis={'title':"Color (B - V)",
            'range': xrange,
            'linecolor': 'rgba(0, 0, 0, 1)',
            'mirror': True, 
            'zeroline':False},
        yaxis={'title':'Apparent Magnitude (V)',
            'range':yrange,
            'linecolor': 'rgba(0, 0, 0, 1)',
            'mirror': True, 
            'zeroline':False},
        plot_bgcolor = 'rgba(0, 0, 0, 0)',
        showlegend=False,
        width=600,
        height=600)

# fig.write_html('./'+name+'cmd.html')
fig.show();

CMD with ZAMS

Open Clusters
In [186]:
xrange = [-0.5,1.7]
yrange = [21,7]
name = 'NGC2168'
masslim = 15

data_read = pd.read_csv('./Labs/CMD/NGC2168/webda_2168.dat',sep='\s+',comment='#',index_col=False,header=None,usecols=(1,2,3))
columns = ['ref','V_0','B_V_0']
data_read.columns = columns

# if "2264"
# refnum = 168
# dist_mod = 9.28
# dist_mod_orig = 9.28
# EB_V = 0.051
# EB_V_orig = 0.051
# lgage = 6.954
# trim = 1.5
# bert_ref = 7.0

# # if "2547"
# refnum = 235
# dist_mod = 8.0
# dist_mod_orig = 8.42
# EB_V = 0.041
# EB_V_orig = 0.041
# lgage = 7.557
# trim = 1.5
# bert_ref = 7.5

# # if "2168"
refnum = 209
dist_mod = 9.6
dist_mod_orig = 10.37
EB_V = 0.262
EB_V_orig = 0.262
lgage = 7.979
trim = 0.7
# bert_ref = 8.0
bert_ref = 6.6

# # if "2099"
# refnum = 285
# dist_mod = 10.8
# dist_mod_orig = 11.64
# EB_V = 0.375
# EB_V_orig = 0.302
# lgage = 8.540
# trim = 1.
# bert_ref = 8.5

# # if "5822"
# refnum = 236
# dist_mod =   9.6
# dist_mod_orig = 10.28
# EB_V = 0.25
# EB_V_orig = 0.150
# lgage = 8.821
# trim = 1.5
# bert_ref = 8.8

# # if "2682"
# refnum = 54
# dist_mod = 9.7
# dist_mod_orig = 9.97
# EB_V = 0.059
# EB_V_orig = 0.059
# lgage = 9.409
# trim = 0.75
# bert_ref = 9.4

# # if "0188"
# refnum = 240
# dist_mod = 11.15
# dist_mod_orig = 11.81
# EB_V = 0.12
# EB_V_orig = 0.082
# lgage = 9.632
# trim = 1
# bert_ref = 9.6

# # if "6791"
# refnum = 46
# dist_mod = 12.6
# dist_mod_orig = 14.20
# EB_V = 0.217
# EB_V_orig = 0.117
# lgage = 9.643
# lgage_orig = 9.643
# trim = 0.75
# bert_ref = 9.9
# bert_ref_orig = 9.6

# correct for extinction (using Schlegel table, landolt B,V, R_V=3.1)
AB_EB_V = 4.315
AV_EB_V = 3.315

data = data_read[data_read['ref'] == refnum]
B_V_0 = data['B_V_0']
V_0 = data['V_0']

data['B_V'] = B_V_0 - EB_V
data['V'] = V_0 - AV_EB_V * EB_V
data['M_V'] = V_0 - dist_mod

table = pd.read_csv('./Labs/CMD/Bertelli/table5.dat',comment='#',sep=' ')
tab = table[(table['lgage']==6.6) & (table['mass'] < masslim)]


M_V_bert = tab['V_abs']
B_V_bert = tab['B-V']

interped = interp1d(x=B_V_bert,y=(M_V_bert + dist_mod + trim),fill_value='extrapolate')
Vlim = interped(data['B_V'])
trimmed = data[(data['V'] < Vlim)]
B_V = trimmed['B_V']
V = trimmed['V']
# ZAMS stuff:
zams_B_V = B_V_bert
zams_M_V = M_V_bert
zams_mass = tab['mass']
zams_pts_mass = [0.75,1,1.5,2,3,4,5,7,10]
zams_pts_mass.reverse()
foo = np.sort(zams_B_V )
goo = np.sort(zams_M_V) 
joo = np.sort(zams_mass)
interped1 = interp1d(zams_mass,zams_B_V,fill_value='extrapolate')
interped2 = interp1d(zams_mass,zams_M_V,fill_value='extrapolate')
zams_pts_B_V = interped1(zams_pts_mass)
zams_pts_M_V = interped2(zams_pts_mass)

/Users/anikaslizewski/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:101: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy

/Users/anikaslizewski/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:102: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy

/Users/anikaslizewski/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:103: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
In [187]:
fig = go.Figure()

fig.add_trace(go.Scatter(
        x = B_V,
        y = V,
        mode = 'markers',
        marker = dict(
            color = 'red',
            size = 5,
            opacity = 0.5,
            symbol = 'circle-open',
            line_width = 1.5),name='')) 

# ZAMS line and numbers:
fig.add_trace(go.Scatter(
        x = zams_pts_B_V,
        y = zams_pts_M_V + dist_mod,
        mode = 'lines',
        line_color = 'black',
        line_width = 1,
        name = '',hoverinfo='skip'))
for i in range(len(zams_pts_B_V)):
    fig.add_annotation(
        x = zams_pts_B_V[i],
        y = zams_pts_M_V[i] + dist_mod,
        text = zams_pts_mass[i],
        ax = -10, ay = 30)
# Layout:
fig.update_layout(title = {'text': name,
            'y' : 0.86,
            'x' : 0.5},
        xaxis = {'title':"Color (B - V)",
            'range': xrange,
            'linecolor': 'rgba(0, 0, 0, 1)',
            'mirror': True, 
            'zeroline': False},
        yaxis = {'title':'Apparent Magnitude (V)',
            'range': yrange,
            'linecolor': 'rgba(0, 0, 0, 1)',
            'mirror': True, 
            'zeroline': False},
        plot_bgcolor = 'rgba(0, 0, 0, 0)',
        showlegend = False,
        width = 600,
        height = 600)

# fig.write_html('./'+name+'cmd.html')
fig.show();