====== Modeling Brown Dwarf Evolution Influenced by Close M-dwarf Companion ======
=== RudModel.ipynb ===
This Jupyter notebook outlines a model to simulate the evolution of a brown dwarf when influenced by a close M-dwarf companion. The model includes importing necessary libraries such as numpy, scipy, and matplotlib for numerical computations and plotting, along with a specific library for astronomical simulations. It details the initialization of constants like solar luminosity, radius, and mass, alongside brown dwarf parameters. It further describes the loading of an MDwarf model, the creation and normalization of various functions for the model, and the establishment of classes for modeling the interaction between a brown dwarf and its M-dwarf companion, specifically focusing on the aspects of flaring events, luminosity calculations, and temperature variations over time due to these stellar interactions. The model attempts to quantify how such interactions, particularly flaring events from the M-dwarf, could affect the temperature evolution of a brown dwarf, including detailed methods for calculating incident luminosity, orbital parameters, and resulting surface parameters and temperatures over significant time scales.


```python
import numpy as np
import scipy
import matplotlib.pyplot as plt
from pysep.dsep import load_model
from scipy.interpolate import interp1d
from tqdm.notebook import tqdm
from scipy.optimize import curve_fit

try:
    plt.style.use(pubStyle)
except NameError:
    pass
```


```python
SLUM = 3.846e33 # ergs / s
SRAD = 6.957e+10 # cm
SMASS = 1.989e+33 # g
BDLUM = 1/100_000 # L/Lsolar
BDTEMP = 1000 # K
BDRAD = 1.001e+10 # cm
sb = 5.670374419e-5 # cgs
kb =  1.3807e-16#  cm2 g s-2 K-1 (cgs
BDATMMFRAC = 0.00001  #0.000001
G = 6.67430e-8 # cgs
```


```python
basedMDwarfModel = load_model("model.dsep")
```


```python
def PosNormal(mean, sigma):
    x = np.random.normal(loc = mean, scale = sigma, size=1)[0]
    return(x if x>=0 else PosNormal(mean,sigma))
```


```python
def normalize(a):
    if a.max() == a.min():
        return a
    return a - a.min()/(a.max()-a.min())
```


```python
def interpTRK(targetAge, df):
    ages = df.AGE
    closestGreater = ages[ages > targetAge].min()
    closestLess = ages[ages < targetAge].max()
    row_closest_less = df[df['AGE'] == closestLess]
    row_closest_greater = df[df['AGE'] == closestGreater]
    fraction = (targetAge - closestLess) / (closestGreater - closestLess)
    interpolated_row = row_closest_less.drop(columns=['AGE']).reset_index(drop=True) + \
                   fraction * (row_closest_greater.drop(columns=['AGE']).reset_index(drop=True) - row_closest_less.drop(columns=['AGE']).reset_index(drop=True))

    return interpolated_row

```


```python
def generate_flare_events_array(timeArraySize, flareProp):
    randomNumbers = np.random.rand(timeArraySize)
    flareEventsArray = randomNumbers < flareProp
    
    return flareEventsArray
```


```python
def build_flare(t12):
    trise = t12[t12 <= 0]
    tdecay = t12[t12 > 0]
    frise = 1 + 1.941*trise - 0.175*trise**2 - 2.246*trise**3 - 1.1225 * trise**4
    fdecay = 0.6890 * np.exp(-1.6 * tdecay) + 0.3030*np.exp(-0.2783 * tdecay)
    return np.hstack([frise, fdecay])
```


```python
class TriangulatedSphere:
    _SLUM = 3.846e33
    _SRAD = 6.957e+10

    def __init__(self, nTheta, nPhi, times, model, modelAge=5, noiseScale=0.1):
        trk = interpTRK(modelAge, model.trk)
        self.radius = 10**trk.log_R.values[0] * self._SRAD
        self.lum = 10**trk.log_L.values[0] * self._SLUM
        
        self.nTheta = nTheta
        self.nPhi = nPhi
        self.times = times
        self.fluxMap = np.zeros((times.shape[0], nTheta, nPhi))

        
        self.thetaValues = np.linspace(0, np.pi, nTheta )
        self.phiValues = np.linspace(0, 2 * np.pi, nPhi)
        self.surfaceArea = 4 * np.pi * self.radius**2
        self.triangleArea = self.surfaceArea / (nTheta * nPhi)
        self.meanFlux = (self.lum)/(4*np.pi * (self.radius)**2)
        self.set_overall_flux(self.meanFlux)
        self.noiseLevel = noiseScale*self.meanFlux

    def set_overall_flux(self, flux):
        """Set the same flux value for the entire surface."""
        self.fluxMap[:, :, :] = flux

    def set_overall_flux_at_time(self, flux, time):
        tID = self.time_to_index(time)
        self.fluxMap[tID, :, :] = flux

    def incriment_overall_flux_at_time(self, df, time):
        tID = self.time_to_index(time)
        self.fluxMap[tID, :, :] += df
    
    def set_flux_at(self, theta, phi, flux, time):
        iTheta, iPhi = self.theta_phi_to_index(theta, phi)
        tID = self.time_to_index(time)
        self.fluxMap[tID, iTheta, iPhi] = flux
        
    def get_flux_at(self, theta, phi, time):
        iTheta, iPhi = self.theta_phi_to_index(theta, phi)
        tID = self.time_to_index(time)
        return self.fluxMap[tID, iTheta, iPhi]

    def get_luminosity_over_solid_angle(self, thetaCenter, phiCenter, omega, time, solar=False):
        # Convert time to indices
        if isinstance(time, float):
            tIDs = np.array([self.time_to_index(time)])
        elif isinstance(time, slice):
            tIDs = np.array([self.time_to_index(t) for t in self.times[time]])
        elif isinstance(time, str) and time == '*':
            tIDs = np.arange(len(self.times))

        lums = np.zeros(len(tIDs))
        
        thetaGrid, phiGrid = np.meshgrid(self.thetaValues, self.phiValues, indexing='ij')
        distance = self.angular_distance_vectorized(thetaGrid, phiGrid, thetaCenter, phiCenter)
        withinSolidAngle = distance < np.sqrt(omega / np.pi)

        for i, tID in enumerate(tIDs):
            if solar:
                lums[i] = np.sum(self.fluxMap[tID, withinSolidAngle] * self.triangleArea / self._SLUM)
            else:
                lums[i] = np.sum(self.fluxMap[tID, withinSolidAngle] * self.triangleArea)

        return lums

    @staticmethod
    def angular_distance_vectorized(theta1, phi1, theta2, phi2):
        """Vectorized calculation of angular distance."""
        return np.arccos(np.sin(theta1) * np.sin(theta2) + np.cos(theta1) * np.cos(theta2) * np.cos(phi1 - phi2))

    def index_to_theta_phi(self, iTheta, iPhi):
        """Convert grid indices to theta and phi values."""
        return self.thetaValues[iTheta], self.phiValues[iPhi]
    
    def theta_phi_to_index(self, theta, phi):
        iTheta = np.argmin(np.abs(self.thetaValues - theta))
        iPhi = np.argmin(np.abs(self.phiValues - phi))
        return iTheta, iPhi

    def index_to_time(self, tID):
        return self.times[tID]

    def time_to_index(self, time):
        return np.argmin(np.abs(self.times - time))

    def evolve(self, flareProb, flareMeanDuration, flareMeanAmplitude, flareMeanLat, flareDurationScale, flareAmplitudeScale, flareLatScale):
        flareDecisions = np.random.uniform(0, 1, size=len(self.times)) < flareProb
        flareDurations = np.abs(np.random.normal(loc=flareMeanDuration, scale=flareDurationScale, size=flareDecisions.sum()))
        flareAmplitudes = np.abs(np.random.normal(loc=flareMeanAmplitude, scale=flareAmplitudeScale, size=flareDecisions.sum()))
        
        flarePhis = np.radians(np.random.choice([-1, 1], size=flareDecisions.sum())) * np.abs(np.random.normal(loc=flareMeanLat, scale=flareLatScale, size=flareDecisions.sum()))
        flareThetas = np.random.uniform(0, 2*np.pi, size=flareDecisions.sum())
        
        noise = np.random.normal(loc=0, scale=self.noiseLevel, size=len(self.times))
        
        fID = 0
        for tid, (time, doFlare, s) in tqdm(enumerate(zip(self.times, flareDecisions, noise)), total=len(self.times), desc="Evolving"):
            if doFlare:
                timeCond = (self.times >= time) & (self.times <= time + flareDurations[fID])
                flareTime = np.linspace(-1, 7, timeCond[timeCond==True].shape[0])
                flareData = build_flare(flareTime)              
                timeCond = (self.times >= time) & (self.times <= time + flareDurations[fID])
                flare = self.meanFlux + (self.meanFlux * flareAmplitudes[fID] * normalize(flareData))
                for flareTimeId, (flareFlux, dft) in enumerate(zip(flare, self.times[timeCond])):
                    self.set_flux_at(flareThetas[fID], flarePhis[fID], flareFlux, time + dft)
                fID += 1

            self.incriment_overall_flux_at_time(s, time)
```


```python
estimateParticlesBBN = lambda mass: mass * (0.75 / 1.67e-27 + 0.25 / 6.64e-27 + 1e-9 / 1.15e-26)
```


```python
def incident_luminosity(luminosityA, aR, bR, d):
    crossSectionalAreaB = np.pi * bR**2
    surfaceAreaD = 4 * np.pi * d**2
    fractionIncident = crossSectionalAreaB / surfaceAreaD
    incidentLuminosityB = luminosityA * fractionIncident
    
    return incidentLuminosityB

```


```python
def orbital_params(m1, m2, p):
    Msolar = 1.989e+33 # g
    mt = m1 + m2
    mt = Msolar * mt
    G = 6.67430e-8 # cgs
    a = ((p**2 * G * mt)/(4*np.pi**2))**(1/3) # cm
    return a
```

$$ L = 4\pi R^{2}\sigma T^{4} $$
$$ U = \frac{3}{2} N k_{B} T $$


```python
def calculate_new_surface_params(iL, temp, mass, dt):
    N = estimateParticlesBBN(mass*BDATMMFRAC) # assume all energy is distributed through the atmosphere
    Q = iL * dt
    dT = Q / ((3/2) * N * kb)
    return temp + dT
```


```python
def cool_blackbody(radius, initialTemperature, mass, delta_time, internal_energy):
    N = estimateParticlesBBN(mass*BDATMMFRAC) # assume all energy is distributed through the atmosphere
    Ui = 3/2 * N * kb * initialTemperature
    
    initial_luminosity = 4 * np.pi * radius**2 * sb * (initialTemperature)**4
    energy_lost = initial_luminosity * delta_time

    Uf = (Ui - energy_lost) + internal_energy
    finalTemperature = (2/3 * Uf)/(N*kb)
    return finalTemperature
```

All units are kept in cgs (roughly, using pi 10 ^ 7 for sec - > year conversion)


```python
# Define how many time steps are used. The cadence right now has a probelmatically large effect on the outcome. This needs
#  to be smoothed out.
period = 6 * 3600
orbits = 50
tmax = period * orbits
cadence = 30
times = np.arange(0, tmax, cadence)
```

One major issue right now is that the flaring model predicts flares based on the cadence of the data (i.e. a 5% flare rate means
that every time step there is a 5 percent change of a flare event. This means that if you increase the observational
frequency the number of flares will go up. This needs to be chaged to a flare rate / time)


```python
sphere = TriangulatedSphere(model=basedMDwarfModel, times=times, nTheta=5, nPhi=5, noiseScale=0.001)
sphere.evolve(0.05, 1000*60, 0.25, 30, 5*60, 1, 10)
#sphere.evolve(0.05, 100*60, 0.01, 30, 5*60, 0.1, 10)
#flareProb, flareMeanDuration, flareMeanAmplitude, flareMeanLat, flareDurationScale, flareAmplitudeScale, flareLatScale):

# Assuming that the components are tidally locked we only need to subtend one hemisphere to get the incident luminosity at any given time
hemiLuminosity = sphere.get_luminosity_over_solid_angle(0, 0, omega=2*np.pi, time='*')
```


    Evolving:   0%|          | 0/36000 [00:00<?, ?it/s]



```python
fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(times, hemiLuminosity)
ax.set_xlabel("Time [s]")
ax.set_ylabel("Luminsoity [erg s$^{-1}$]")
```




    Text(0, 0.5, 'Luminsoity [erg s$^{-1}$]')




    
![png](media:7a650a09f81d2cb895880d1e4f53652cc84c7d4139a21af15f69598e43dfc9d0:output_19_1.png)
    



```python
# The conversion from solar mass to g is done in the orbital params function since dsep tracks
# everyging in either solar units or log solar uints
a = orbital_params(basedMDwarfModel.mass, 0.075, period)

dt = (times[1] - times[0]) * 60

BrownDwarfTemps = list()
bdtemp = BDTEMP
bdmass = 0.075 * SMASS

# Here is a handwavey thing. I am basically assuming that all the energy
#  dumped onto the star stays in the atmosphere (and that the atmosphre is 
#  and idea monatomic gas (f=3))
N = estimateParticlesBBN(bdmass*BDATMMFRAC)

# this is used to figure out how much energy the star radiates away so
#  that some interneral energy fudge factor can be added to keep the
#  atmosphere from rapidly cooling to 0
initialCooling = bdtemp - cool_blackbody(BDRAD, bdtemp, bdmass, dt, 0)

# Due to the high temporal resolutions required to resolve flares and the long
# baseline needed to heat the star I run this section to just fine the mean 
# temperature change per orbit.
porbitID = 0
temps = list()
pOrbitTemp = bdtemp
for tid, (t, l) in enumerate(zip(times, hemiLuminosity)):
    orbitID = t // period
    if orbitID != porbitID:
        temps.append(bdtemp - pOrbitTemp)
        pOrbitTemp = bdtemp
        porbitID = orbitID
    incidentUponBD = incident_luminosity(l, sphere.radius, 1.001e+10, a)
    bdtemp = calculate_new_surface_params(incidentUponBD, bdtemp, bdmass, dt)
    BrownDwarfTemps.append(bdtemp)
```

### Mean Temperature Change Per Orbit


```python
line = lambda x, m, b: m*x + b
fit, ax = curve_fit(line, range(len(temps)), temps)
plt.plot(temps)
plt.plot(range(len(temps)), line(range(len(temps)), *fit))
```




    [<matplotlib.lines.Line2D at 0x7f75cf909370>]




    
![png](media:7a650a09f81d2cb895880d1e4f53652cc84c7d4139a21af15f69598e43dfc9d0:output_22_1.png)
    


Now we will assume that this mean temperature change per orbit (on the hot side) holdes true for ever and evolve this over many
many more orbits. We make the assumption that the hot side temperature is independent of the cold side temperature so we first calculate
just the hot side temperature|. This allows us to more easily use a lookback time for the cold side temperature latter


```python
meanTempChangePerOrbit = np.mean(temps)

morbits = 1e9
orbitJump = 5e2

orbits = np.arange(0, morbits, orbitJump)

hotSideTemp = np.zeros_like(orbits)
coldSideTemp = np.zeros_like(orbits)

hottemp = BDTEMP
coldtemp = BDTEMP

period = 6*3600 #seconds
longTimes = np.linspace(0, morbits*period, orbits.shape[0])

bdmass = 0.075 * SMASS
intrinsicCoolingCoef = 0.9
initialCooling = hottemp - cool_blackbody(BDRAD, hottemp, bdmass, period*orbitJump, 0)
internalEnergy = (3/2) * N * kb * initialCooling * intrinsicCoolingCoef


for oID, orbit in tqdm(enumerate(orbits), total=orbits.shape[0]):
    hottemp += (orbitID + 1) * meanTempChangePerOrbit
    hottemp = cool_blackbody(BDRAD, hottemp, bdmass, orbitJump * period, internalEnergy) 
    hotSideTemp[oID] = hottemp
```


      0%|          | 0/2000000 [00:00<?, ?it/s]



```python
fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(longTimes/(np.pi*1e7 * 1e6), hotSideTemp)
ax.set_xlabel("Time [Myr]", fontsize=27)
ax.set_ylabel("Hot Side Teff [K]", fontsize=27)
```




    Text(0, 0.5, 'Hot Side Teff [K]')




    
![png](media:7a650a09f81d2cb895880d1e4f53652cc84c7d4139a21af15f69598e43dfc9d0:output_25_1.png)
    



```python
# we then build an interpolator so that we can look back on the hot side temp from before
hst = interp1d(longTimes, hotSideTemp)
```


```python
tauCO = 70*24*3600 # s (wright et al. 2018 -- calibrated for mid M dwarfs, likeley not correct for Brown Dwarfs)
offset = 300 # to allow for thermal transfer time # need to connect this to the convective overturn timescale or better yet dont be stupid about it
coldSideTemp = np.zeros_like(orbits) + BDTEMP
for oID, orbit in tqdm(enumerate(orbits[offset:]), total=orbits[offset:].shape[0]):
    delayedHotSideTemp = hst(oID*period)
    coldSideTemp[oID + offset] = delayedHotSideTemp
    coldSideTemp[oID + offset] = cool_blackbody(BDRAD, coldSideTemp[oID + offset], bdmass, orbitJump * period, internalEnergy) 
```


      0%|          | 0/1999700 [00:00<?, ?it/s]



```python
fig, ax = plt.subplots(1,1,figsize=(10, 7))
ax.plot(longTimes/(np.pi*1e7*1e6), hotSideTemp, label="Hot Side")
ax.plot(longTimes/(np.pi*1e7*1e6), coldSideTemp, label="Cold Side")
ax.set_xlabel("Time [Myr]", fontsize=27)
ax.set_ylabel("Teff [K]", fontsize=27)
ax.legend()
```




    <matplotlib.legend.Legend at 0x7f75cfd272b0>




    
![png](media:7a650a09f81d2cb895880d1e4f53652cc84c7d4139a21af15f69598e43dfc9d0:output_28_1.png)
    



```python

```


```python

```