<jupyter_start><jupyter_text>Calcualting the stellar wind<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
import os
import glob
from astropy.table import Table
from astropy.io import ascii
import astropy.units as u
import astropy.constants as const
from scipy.interpolate import interpolate
from craftroom import resample
from astropy.convolution import convolve, Box1DKernel
from scipy.io.idl import readsav
from astropy import constants as const

#matplotlib set up
%matplotlib inline
from matplotlib import rcParams
rcParams["figure.figsize"] = (14, 5)
rcParams["font.size"] = 20

rates = Table.read('tables/pceb_fit_results_rates.csv')
rates
els = ['si', 'c', 'o', 's']
el_mass = [28.085, 12.011, 15.999, 32.06]
sol_frac = [6.649e-4, 2.365e-3,5.733e-3, 3.093e-4]
def find_v(a, Rrd, Mrd, Mwd, p):
    """
    velocity of the stellar wind at distance a
    """
    a, Rrd, Mrd, Mwd, p = a.to(u.m), Rrd.to(u.m), Mrd.to(u.kg), Mwd.to(u.kg), p.to(u.s)
    vorb = (2*np.pi/p) * a * (Mrd/(Mrd+Mwd)) #orbital velocity of wd
#     print('vorb=', vorb)
    vwind = ((2*const.G*Mrd)/Rrd)**0.5 #approximate as escape velocity
#     print('vwind=' , vwind)
    vrel = (vwind**2+vorb**2)**0.5
    return vrel

def mass_loss(mdot, Mwd, a, v): 
    """mass loss rate of mdwarf to produce an accretion rate mdot onto a white dwarf of mass Mwd at separation a"""
#     m_loss = (((2*a*Mrd)/(Rrd*Mwd))**2) * mdot
#     print(Mwd, Mrd)
    mdot, Mwd, a, v =  mdot.to(u.g/u.s), Mwd.to(u.kg), a.to(u.m), v.to(u.m/u.s)
#     print(v)
  #  v= 600000*u.m/u.s
#     print(a.to(u.Rsun))
   # print(v)
#     print (mdot, Mwd, a, v)
    scale = (v**4*a**2)/(const.G**2*Mwd**2)
#     print('scale=',scale)
    m_loss = mdot*scale
    return m_loss

def total_infer(mdot, el):
#     print(el)
    els = np.genfromtxt('solar_abd.tsv', names=True, delimiter='\t', dtype=None, encoding=None)
    mf = els['mass_fraction'][np.where(els['Atom_X']==el)][0]
    total_mdot = mdot/mf
    return total_mdot

def p_to_a(m1, m2, p):
    m1, m2, p = m1.to(u.kg), m2.to(u.kg), p.to(u.s)
    a = ((const.G*(m1+m2) * p**2)/(4*np.pi**2))**(1/3)
#     print(a.to(u.Rsun))
    return a
                                   
def find_wind(mdot, Mwd,p, Rrd, Mrd, el):
#     print(el)
#     print(mdot)
    a = p_to_a(Mwd, Mrd, p)
    v = find_v(a, Rrd, Mrd, Mwd, p)
    rate = mass_loss(mdot, Mwd, a, v)
#     print(rate)
#     print(rate)
    total = total_infer(rate, el.capitalize())
#     print(total)
    return(total)
                                   
wind_rates = []
wind_errors = []
masses = []
sptypes = [] 
for row in rates:
#     if row['Target'] == 'WD0710+741':
    if row['Teff'] < 30000:
        masses.append(row['mass_est'])
        sptypes.append(row['spt'])
        mdots = []
        for i, el in enumerate(els):
            if row['{}_e'.format(el)] > 0.0:
                acc = row['F{}'.format(el)]*(u.g/u.s)
                Mwd, Mrd = row['mass']*u.M_sun, row['mass_est']*u.M_sun
                Rrd = row['rad_est']*u.R_sun 
                p = row['porb_min']*u.min
#                 print(p)
                mdot = find_wind(acc, Mwd,p, Rrd, Mrd, el)
                mdots.append(mdot.value)
      #  print(mdots)
#         mdots = np.array(mdots)*u.g/u.s
     #   print(np.mean(mdots))
        wind_rates.append(np.mean(mdots))
        wind_errors.append(np.std(mdots))

# print(wind_rates)
wind_rates = np.array(wind_rates)*(u.g/u.s).to(u.M_sun/u.yr)
wind_errors = np.array(wind_errors)*(u.g/u.s).to(u.M_sun/u.yr)

log_wind_rates = np.log10(wind_rates)
log_wind_errors = 0.434 * (wind_errors/wind_rates)
# print(np.log10(wind_rates))
masses= np.array(masses)

plt.errorbar(masses, log_wind_rates, yerr=log_wind_errors, marker='o', ls='none')
plt.xlim(0.49, 0.051)
plt.axhline(np.log10(2e-14), c='C1', ls='--')
plt.xlabel(r'Mass (M$_{\odot}$)')
plt.ylabel(r'\.M (M$_{\odot}$yr$^{-1}$)')
plt.tight_layout()
#  print(len(masses))
#print(masses)

rates_tab= np.genfromtxt('wind_measurements.csv', delimiter=',', dtype=None, encoding=None, names=True)
sol_n = 2e-14
sol = np.log10(2e-14)

# new_types = [4, 5, 3, 4.0, 4.5, 3.5, 2.5, 4.5, 3.0, 3.5, 5.0, 6, 4.0, 4.5]
# row_1 = [4, 5, 3, 4.5, 3.5, 2.5,  3.5, 6, 8]
# row_2 = [4.0, 4.5,3.0, 5.0,]
# row_3 = [4.0, 4.5]

data_m = rates_tab[rates_tab['e_mdot'] != 0]
data_u = rates_tab[rates_tab['e_mdot'] == 0]

fig, ax =plt.subplots(figsize=(10,6))
plt.errorbar(data_m['spn'][data_m['method'] =='pceb'], data_m['log_mdot'][data_m['method'] =='pceb'], 
             yerr= data_m['e_mdot'][data_m['method'] =='pceb'], marker='o', ls='none', label='Binaries')
plt.errorbar(data_m['spn'][data_m['method'] !='pceb'], data_m['log_mdot'][data_m['method'] !='pceb'], 
             yerr= data_m['e_mdot'][data_m['method'] !='pceb'], marker='o', ls='none', label='Single Stars')
plt.errorbar(data_u['spn'][data_u['method'] !='pceb'], data_u['log_mdot'][data_u['method'] !='pceb'], yerr=0.5, marker='o', ls='none', uplims=True, c='C1')

newtypes = [float(spt[1:]) for spt in sptypes]
plt.errorbar(newtypes+np.random.normal(0, 0.1, len(newtypes)), log_wind_rates, yerr=log_wind_errors, marker='o', ls='none', c='C3', label ='This work')


plt.xlabel('Spectral Type (M$n$)')
plt.ylabel('log\,[\.M] ($\mathrm{M}_{\odot}\mathrm{yr}^{-1}$)')
ax.tick_params(
    axis='x',          # changes apply to the x-axis
    which='minor',      # both major and minor ticks are affected
    top=False,
    bottom=False)      # ticks along the bottom edge are off 
plt.xlim(1.5, 8.5)
plt.axhline(sol, ls='--', c='C2', alpha =0.7, label='Solar')
plt.ylim(-16.9, -13.1)
# plt.scatter(row_1, np.full(len(row_1), -13), marker='x', c='C3', label = 'New Systems')
# plt.scatter(row_2, np.full(len(row_2), -13.2), marker='x', c='C3')
# plt.scatter(row_3, np.full(len(row_3), -12.8), marker='x', c='C3')



plt.legend(frameon=True, loc=4)
plt.tight_layout()
# plt.savefig('../plots/wind_sota_150.pdf', dpi=150)
[float(spt[1:]) for spt in sptypes]
a[0].upper()
const.G

a= 0.7*u.Rsun
Mwd = 0.389*u.Msun
vs = np.arange(50, 601, 50)*u.km/u.s
scales = np.array([(v**4*a**2)/(const.G**2*Mwd**2) for v in vs])
plt.plot(vs, scales)
plt.axhline(1, c='C1')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20798454603992272, 'fraction_numerical': 0.03976175144880876, 'mean_word_length': 3.482683982683983, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '504926', 'license': 'MIT', 'n_tokens_mistral': 2892, 'n_tokens_neox': 2716, 'n_words': 517, 'path': '.ipynb_checkpoints/windv1-checkpoint.ipynb', 'repo_name': 'davidjwilson/pceb', 'size': '6212'}