Skip to content

Instantly share code, notes, and snippets.

@cjtu

cjtu/ast_501_pset1.ipynb

Last active October 16, 2019 01:47
Show Gist options
  • Save cjtu/8ade72b627d6e32aa7eb017d96cc5f66 to your computer and use it in GitHub Desktop.
Save cjtu/8ade72b627d6e32aa7eb017d96cc5f66 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
ast_501_pset1.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# AST501: Problem Set 1\n",
"\n",
"## Q1: Minimum Mass Trappist-1 Nebula\n",
"\n",
"Compute the minimum mass of the Trappist-1 Nebula based on the known planets in [Grimm et al. (2018)](https://arxiv.org/abs/1802.01377).\n",
"\n",
"### a) Planet surface densities\n",
"\n",
"Take midpoint between each planet and use the formula:\n",
"\n",
"$$\\frac{m_{planet}}{Area}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R\n",
"name \n",
"b 0.011548 1.020 1.010\n",
"c 0.015815 1.180 1.060\n",
"d 0.022280 0.281 0.697\n",
"e 0.029283 0.766 0.913\n",
"f 0.038534 0.926 0.986\n",
"g 0.046877 1.140 1.050\n",
"h 0.061935 0.313 0.719"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"# Planet properties taken from Grimm et al. (2018)\n",
"d = {\n",
" 'name': ['b','c','d','e','f','g','h'],\n",
" 'a': [0.01154775, 0.01581512, 0.02228038, 0.02928282, 0.03853361, 0.04687692, 0.06193488],\n",
" 'M': [1.02, 1.18, 0.281, 0.766, 0.926, 1.14, 0.313],\n",
" 'R': [1.01, 1.06, 0.697, 0.913, 0.986, 1.05, 0.719]\n",
"}\n",
"df = pd.DataFrame(d)\n",
"df.set_index('name', inplace=True)\n",
"df.head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the midpoints of the semi-major axes"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" <th>ann_min</th>\n",
" <th>ann_max</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" <td>0.009414</td>\n",
" <td>0.013681</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" <td>0.013681</td>\n",
" <td>0.019048</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" <td>0.019048</td>\n",
" <td>0.025782</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" <td>0.025782</td>\n",
" <td>0.033908</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" <td>0.033908</td>\n",
" <td>0.042705</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" <td>0.042705</td>\n",
" <td>0.054406</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" <td>0.054406</td>\n",
" <td>0.069464</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R ann_min ann_max\n",
"name \n",
"b 0.011548 1.020 1.010 0.009414 0.013681\n",
"c 0.015815 1.180 1.060 0.013681 0.019048\n",
"d 0.022280 0.281 0.697 0.019048 0.025782\n",
"e 0.029283 0.766 0.913 0.025782 0.033908\n",
"f 0.038534 0.926 0.986 0.033908 0.042705\n",
"g 0.046877 1.140 1.050 0.042705 0.054406\n",
"h 0.061935 0.313 0.719 0.054406 0.069464"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"midpoints = df['a'][:-1].values + (df['a'][1:].values - df['a'][:-1].values)/2\n",
"inner_pt = 2*df['a'][0] - midpoints[0]\n",
"outer_pt = 2*df['a'][-1] - midpoints[-1]\n",
"df['ann_min'] = np.insert(midpoints, 0, inner_pt)\n",
"df['ann_max'] = np.concatenate([midpoints, [outer_pt]])\n",
"df.head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute the area of annuli based on the midpoints."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" <th>ann_min</th>\n",
" <th>ann_max</th>\n",
" <th>ann_area</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" <td>0.009414</td>\n",
" <td>0.013681</td>\n",
" <td>0.000310</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" <td>0.013681</td>\n",
" <td>0.019048</td>\n",
" <td>0.000552</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" <td>0.019048</td>\n",
" <td>0.025782</td>\n",
" <td>0.000948</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" <td>0.025782</td>\n",
" <td>0.033908</td>\n",
" <td>0.001524</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" <td>0.033908</td>\n",
" <td>0.042705</td>\n",
" <td>0.002117</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" <td>0.042705</td>\n",
" <td>0.054406</td>\n",
" <td>0.003570</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" <td>0.054406</td>\n",
" <td>0.069464</td>\n",
" <td>0.005860</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R ann_min ann_max ann_area\n",
"name \n",
"b 0.011548 1.020 1.010 0.009414 0.013681 0.000310\n",
"c 0.015815 1.180 1.060 0.013681 0.019048 0.000552\n",
"d 0.022280 0.281 0.697 0.019048 0.025782 0.000948\n",
"e 0.029283 0.766 0.913 0.025782 0.033908 0.001524\n",
"f 0.038534 0.926 0.986 0.033908 0.042705 0.002117\n",
"g 0.046877 1.140 1.050 0.042705 0.054406 0.003570\n",
"h 0.061935 0.313 0.719 0.054406 0.069464 0.005860"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['ann_area'] = np.pi*(df['ann_max']**2 - df['ann_min']**2)\n",
"df.head(7)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" <th>ann_min</th>\n",
" <th>ann_max</th>\n",
" <th>ann_area</th>\n",
" <th>rho_surf</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" <td>0.009414</td>\n",
" <td>0.013681</td>\n",
" <td>0.000310</td>\n",
" <td>87906.146493</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" <td>0.013681</td>\n",
" <td>0.019048</td>\n",
" <td>0.000552</td>\n",
" <td>57065.994962</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" <td>0.019048</td>\n",
" <td>0.025782</td>\n",
" <td>0.000948</td>\n",
" <td>7906.552869</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" <td>0.025782</td>\n",
" <td>0.033908</td>\n",
" <td>0.001524</td>\n",
" <td>13412.989686</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" <td>0.033908</td>\n",
" <td>0.042705</td>\n",
" <td>0.002117</td>\n",
" <td>11670.124142</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" <td>0.042705</td>\n",
" <td>0.054406</td>\n",
" <td>0.003570</td>\n",
" <td>8521.832518</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" <td>0.054406</td>\n",
" <td>0.069464</td>\n",
" <td>0.005860</td>\n",
" <td>1425.343835</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R ann_min ann_max ann_area rho_surf\n",
"name \n",
"b 0.011548 1.020 1.010 0.009414 0.013681 0.000310 87906.146493\n",
"c 0.015815 1.180 1.060 0.013681 0.019048 0.000552 57065.994962\n",
"d 0.022280 0.281 0.697 0.019048 0.025782 0.000948 7906.552869\n",
"e 0.029283 0.766 0.913 0.025782 0.033908 0.001524 13412.989686\n",
"f 0.038534 0.926 0.986 0.033908 0.042705 0.002117 11670.124142\n",
"g 0.046877 1.140 1.050 0.042705 0.054406 0.003570 8521.832518\n",
"h 0.061935 0.313 0.719 0.054406 0.069464 0.005860 1425.343835"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"au2cm = 1.496e+13\n",
"M_earth = 5.972e27 # grams\n",
"df['rho_surf'] = df['M']*M_earth / (df['ann_area']*au2cm**2)\n",
"df.head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### b) Correction factor"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" <th>ann_min</th>\n",
" <th>ann_max</th>\n",
" <th>ann_area</th>\n",
" <th>rho_surf</th>\n",
" <th>rho</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" <td>0.009414</td>\n",
" <td>0.013681</td>\n",
" <td>0.000310</td>\n",
" <td>87906.146493</td>\n",
" <td>5.458137</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" <td>0.013681</td>\n",
" <td>0.019048</td>\n",
" <td>0.000552</td>\n",
" <td>57065.994962</td>\n",
" <td>5.462265</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" <td>0.019048</td>\n",
" <td>0.025782</td>\n",
" <td>0.000948</td>\n",
" <td>7906.552869</td>\n",
" <td>4.575266</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" <td>0.025782</td>\n",
" <td>0.033908</td>\n",
" <td>0.001524</td>\n",
" <td>13412.989686</td>\n",
" <td>5.549129</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" <td>0.033908</td>\n",
" <td>0.042705</td>\n",
" <td>0.002117</td>\n",
" <td>11670.124142</td>\n",
" <td>5.325846</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" <td>0.042705</td>\n",
" <td>0.054406</td>\n",
" <td>0.003570</td>\n",
" <td>8521.832518</td>\n",
" <td>5.429319</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" <td>0.054406</td>\n",
" <td>0.069464</td>\n",
" <td>0.005860</td>\n",
" <td>1425.343835</td>\n",
" <td>4.642651</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R ann_min ann_max ann_area rho_surf \\\n",
"name \n",
"b 0.011548 1.020 1.010 0.009414 0.013681 0.000310 87906.146493 \n",
"c 0.015815 1.180 1.060 0.013681 0.019048 0.000552 57065.994962 \n",
"d 0.022280 0.281 0.697 0.019048 0.025782 0.000948 7906.552869 \n",
"e 0.029283 0.766 0.913 0.025782 0.033908 0.001524 13412.989686 \n",
"f 0.038534 0.926 0.986 0.033908 0.042705 0.002117 11670.124142 \n",
"g 0.046877 1.140 1.050 0.042705 0.054406 0.003570 8521.832518 \n",
"h 0.061935 0.313 0.719 0.054406 0.069464 0.005860 1425.343835 \n",
"\n",
" rho \n",
"name \n",
"b 5.458137 \n",
"c 5.462265 \n",
"d 4.575266 \n",
"e 5.549129 \n",
"f 5.325846 \n",
"g 5.429319 \n",
"h 4.642651 "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"R_earth = 6.371e8 # cm\n",
"V = (4/3)*np.pi*(df['R']*R_earth)**3\n",
"df['rho'] = df['M']*M_earth / V\n",
"df.head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fit a line to the correction factors of Earth and Neptune, compute correction factor of each planet from `corr(rho)`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>a</th>\n",
" <th>M</th>\n",
" <th>R</th>\n",
" <th>ann_min</th>\n",
" <th>ann_max</th>\n",
" <th>ann_area</th>\n",
" <th>rho_surf</th>\n",
" <th>rho</th>\n",
" <th>corr</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>0.011548</td>\n",
" <td>1.020</td>\n",
" <td>1.010</td>\n",
" <td>0.009414</td>\n",
" <td>0.013681</td>\n",
" <td>0.000310</td>\n",
" <td>87906.146493</td>\n",
" <td>5.458137</td>\n",
" <td>594.632937</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>0.015815</td>\n",
" <td>1.180</td>\n",
" <td>1.060</td>\n",
" <td>0.013681</td>\n",
" <td>0.019048</td>\n",
" <td>0.000552</td>\n",
" <td>57065.994962</td>\n",
" <td>5.462265</td>\n",
" <td>595.162212</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>0.022280</td>\n",
" <td>0.281</td>\n",
" <td>0.697</td>\n",
" <td>0.019048</td>\n",
" <td>0.025782</td>\n",
" <td>0.000948</td>\n",
" <td>7906.552869</td>\n",
" <td>4.575266</td>\n",
" <td>481.444375</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>0.029283</td>\n",
" <td>0.766</td>\n",
" <td>0.913</td>\n",
" <td>0.025782</td>\n",
" <td>0.033908</td>\n",
" <td>0.001524</td>\n",
" <td>13412.989686</td>\n",
" <td>5.549129</td>\n",
" <td>606.298574</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>0.038534</td>\n",
" <td>0.926</td>\n",
" <td>0.986</td>\n",
" <td>0.033908</td>\n",
" <td>0.042705</td>\n",
" <td>0.002117</td>\n",
" <td>11670.124142</td>\n",
" <td>5.325846</td>\n",
" <td>577.672583</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>0.046877</td>\n",
" <td>1.140</td>\n",
" <td>1.050</td>\n",
" <td>0.042705</td>\n",
" <td>0.054406</td>\n",
" <td>0.003570</td>\n",
" <td>8521.832518</td>\n",
" <td>5.429319</td>\n",
" <td>590.938284</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>0.061935</td>\n",
" <td>0.313</td>\n",
" <td>0.719</td>\n",
" <td>0.054406</td>\n",
" <td>0.069464</td>\n",
" <td>0.005860</td>\n",
" <td>1425.343835</td>\n",
" <td>4.642651</td>\n",
" <td>490.083466</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" a M R ann_min ann_max ann_area rho_surf \\\n",
"name \n",
"b 0.011548 1.020 1.010 0.009414 0.013681 0.000310 87906.146493 \n",
"c 0.015815 1.180 1.060 0.013681 0.019048 0.000552 57065.994962 \n",
"d 0.022280 0.281 0.697 0.019048 0.025782 0.000948 7906.552869 \n",
"e 0.029283 0.766 0.913 0.025782 0.033908 0.001524 13412.989686 \n",
"f 0.038534 0.926 0.986 0.033908 0.042705 0.002117 11670.124142 \n",
"g 0.046877 1.140 1.050 0.042705 0.054406 0.003570 8521.832518 \n",
"h 0.061935 0.313 0.719 0.054406 0.069464 0.005860 1425.343835 \n",
"\n",
" rho corr \n",
"name \n",
"b 5.458137 594.632937 \n",
"c 5.462265 595.162212 \n",
"d 4.575266 481.444375 \n",
"e 5.549129 606.298574 \n",
"f 5.325846 577.672583 \n",
"g 5.429319 590.938284 \n",
"h 4.642651 490.083466 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m = (600-100)/(5.5-1.6) # slope\n",
"b = 100 - m*1.6 # y-intercept\n",
"def corr(rho):\n",
" return m*rho + b\n",
"\n",
"df['corr'] = corr(df['rho'])\n",
"df.head(7)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>rho_surf</th>\n",
" <th>corr</th>\n",
" <th>rho_surf_corr</th>\n",
" </tr>\n",
" <tr>\n",
" <th>name</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>b</th>\n",
" <td>87906.146493</td>\n",
" <td>594.632937</td>\n",
" <td>5.227189e+07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>c</th>\n",
" <td>57065.994962</td>\n",
" <td>595.162212</td>\n",
" <td>3.396352e+07</td>\n",
" </tr>\n",
" <tr>\n",
" <th>d</th>\n",
" <td>7906.552869</td>\n",
" <td>481.444375</td>\n",
" <td>3.806565e+06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>e</th>\n",
" <td>13412.989686</td>\n",
" <td>606.298574</td>\n",
" <td>8.132277e+06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>f</th>\n",
" <td>11670.124142</td>\n",
" <td>577.672583</td>\n",
" <td>6.741511e+06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>g</th>\n",
" <td>8521.832518</td>\n",
" <td>590.938284</td>\n",
" <td>5.035877e+06</td>\n",
" </tr>\n",
" <tr>\n",
" <th>h</th>\n",
" <td>1425.343835</td>\n",
" <td>490.083466</td>\n",
" <td>6.985374e+05</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" rho_surf corr rho_surf_corr\n",
"name \n",
"b 87906.146493 594.632937 5.227189e+07\n",
"c 57065.994962 595.162212 3.396352e+07\n",
"d 7906.552869 481.444375 3.806565e+06\n",
"e 13412.989686 606.298574 8.132277e+06\n",
"f 11670.124142 577.672583 6.741511e+06\n",
"g 8521.832518 590.938284 5.035877e+06\n",
"h 1425.343835 490.083466 6.985374e+05"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df['rho_surf_corr'] = df['rho_surf'] * df['corr']\n",
"df[['rho_surf', 'corr', 'rho_surf_corr']].head(7)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### c) Surf density vs Semimajor Axis"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"df.plot(kind='scatter', x='a', y='rho_surf_corr', logx=True, logy=True,\n",
" title='Surface density vs semimajor axis')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### d) Power law fit\n",
"\n",
"Do log linear fit using relationship:\n",
"\n",
"$$y = ax^b$$\n",
"$$ln(y) = b \\cdot ln(x) + ln(a)$$\n",
"$$ln(y) = b \\cdot ln(x) + c$$\n",
"\n",
"Where $c = ln(a)$, so $a = e^{c}$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.optimize import leastsq\n",
"\n",
"powerlaw = lambda x, a, b: a * x**b\n",
"\n",
"# Linear fit funcs\n",
"fitfunc = lambda p, x: p[0] * x + p[1]\n",
"errfunc = lambda p, x, y: (y - fitfunc(p, x))\n",
"\n",
"# Initial guess\n",
"p0 = [1, -3]\n",
"\n",
"logx = np.log(df['a'].values)\n",
"logy = np.log(df['rho_surf_corr'].values)\n",
"\n",
"pfinal, covar, _, _, _ = leastsq(errfunc, p0, args=(logx, logy), full_output=True)\n",
"\n",
"b, c = pfinal\n",
"a = np.exp(c)\n",
"\n",
"print(\"Power Law ax^b, where a={:.3f}, b={:.3f}\".format(a, b))\n",
"\n",
"x = np.linspace(min(df['a']), max(df['a']), 1000)\n",
"\n",
"ax = df.plot(kind='scatter', x='a', y='rho_surf_corr', logx=True, logy=True,\n",
" title='Surface density vs semimajor axis')\n",
"ax.plot(x, powerlaw(x, a, b), label='a:{:.3f}\\nb{:.3f}'.format(a, b))\n",
"ax.set_xlim((0.01,0.1))\n",
"ax.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### e) Compare Trappist-1 to the Solar System\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import matplotlib.pyplot as plt\n",
"f = os.path.join(os.getcwd(), \"ss_densities.png\")\n",
"plt.imshow(plt.imread(f))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The surface densities are $10^7 gcm^{-1}$, much higher than those of our solar system. This could be because the host star's much smaller solar mass, potentially leaving more material in the disk to form large planets relatively near to the star.\n",
"\n",
"The overall shape of the Trappist-1 plot is similar to our solar system, with an overal trend of decreasing surface density from the star and an anomalously low point in the middle (corresponding to the asteroid belt in our SS, and planet d in Trappist-1).\n",
"\n",
"Orbital evolution also controls the computed surface densities, since they are based upon the observed objects remaining in the time since formation. Material ejected from the system or planet migration could explain the deviations from a smooth curve of decreasing surface density with distance from the star."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Q2: Star Masses\n",
"\n",
"### Part 1: Minimum Mass Derivation\n",
"\n",
"#### a) Derive M\n",
"\n",
"Set the radiated energy (luminosity) of the collapsing cloud to potential energy of collapsing cloud.\n",
"\n",
"$$L = PE_{sphere}$$\n",
"\n",
"Where $L = 4 \\pi R^2 \\sigma T^4$, and $PE = \\frac{GM^2}{Rt_{ff}}$\n",
"\n",
"$$4 \\pi R^2 \\sigma T^4 = \\frac{GM^2}{Rt_{ff}}$$\n",
"\n",
"Solving for M:\n",
"\n",
"$$M = \\frac{3 \\sigma T^4}{(G \\rho)^{3/2}}$$\n",
"\n",
"#### b) Derive density\n",
"\n",
"We know the Jeans mass is given by:\n",
"\n",
"$$M_{Jeans} = \\left(\\frac{5k_bT}{Gm}\\right)^{3/2} \\left(\\frac{3}{4 \\pi \\rho}\\right)^{1/2}$$\n",
"\n",
"Setting it equal to our mass expression and solving for density:\n",
"\n",
"$$\\left(\\frac{5k_bT}{Gm}\\right)^{3/2} \\left(\\frac{3}{4 \\pi \\rho}\\right)^{1/2} = \\frac{3 \\sigma T^4}{(G \\rho)^{3/2}}$$\n",
"\n",
"Solving for $\\rho$:\n",
"\n",
"$$\\rho = \\sigma T^{5/2} (\\frac{4 \\pi}{3})^{1/2} (\\frac{m}{5k_B})^{3/2}$$\n",
"\n",
"#### c) Plug density into mimimum mass equation\n",
"\n",
"$$M = \\frac{3 \\sigma T^4}{(G \\rho)^{3/2}}$$\n",
"\n",
"$$M = \\frac{3 \\sigma T^4}{G^{3/2}} \\frac{1}{(\\sigma T^{5/2} (\\frac{4 \\pi}{3})^{1/2} (\\frac{m}{5k_B})^{3/2})^{3/2}}$$\n",
"\n",
"$$M = \\frac{3 \\sigma T^4}{G^{3/2}} \\frac{1}{\\sigma^{3/2} T^{15/4} (\\frac{4 \\pi}{3})^{3/4} (\\frac{m}{5k_B})^{9/4}}$$\n",
"\n",
"$$M = \\frac{(5k_B)^{9/4} (3T)^{1/4}}{(4 \\pi)^{3/4} \\sigma^{1/2} G^{3/2} m^{9/4}}$$\n",
"\n",
"\n",
"### Part 2: Min and Max mass stars\n",
"\n",
"#### a) Compute the minimum and maximum mass of star after formation\n",
"\n",
"Compute min mass given $T = 10 K, n = 10^3 cm^{-3}$\n",
"Compute max mass given $T = 100 K$\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Given\n",
"T_min = 100 # K\n",
"n_max = 1e3*1e6 # m^-3\n",
"T_max = 10 # K\n",
"\n",
"# Constants\n",
"sig = 5.670373e-8 # W /(m^2 K^4)\n",
"G = 6.674e-11 # m^3 / (kg s^2)\n",
"kb = 1.3806e-23 # J/K\n",
"M_sun = 1.989e30 # kg\n",
"m_H = 2*1.6735575e-27 # kg\n",
"\n",
"def minmass(T):\n",
" num = (5 * kb)**2.25 * (3 * T)**0.25\n",
" denom = (4 * np.pi)**0.75 * sig**0.5 * G**1.5 * m_H**2.25\n",
" \n",
" return (1/M_sun) * num / denom\n",
"\n",
"def jeansmass(T, rho):\n",
" num = ((5 * kb * T)**1.5) * (3**0.5)\n",
" denom = ((G * m_H)**1.5) * ((4 * np.pi * rho)**0.5)\n",
" return (1/M_sun) * num / denom\n",
"\n",
"\n",
"rho_max = n_max * m_H # kg / m^3\n",
"\n",
"# Compute the min and max mass solar nebulae\n",
"M_min = minmass(T_min) # Msun\n",
"M_max = jeansmass(T_max, rho_max) # Msun\n",
"\n",
"print(\"Min mass star: {:.3f} Msun\".format(M_min))\n",
"print(\"Max mass star: {:.3f} Msun\".format(M_max))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### b) Plot min and max mass on m22 mass distribution"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m22 = pd.read_csv(\"./M22data.txt\", delimiter=\" \")\n",
"p = m22.plot(kind='scatter', x='Color Index', y='Mass', \n",
" title='Mass vs Color Index (logy)')\n",
"\n",
"x = np.linspace(np.min(m22['Color Index']), np.max(m22['Color Index']))\n",
"p.plot(x, [M_max]*len(x), label='Max {:.3f}'.format(M_max))\n",
"p.plot(x, [M_min]*len(x), label='Min {:.3f}'.format(M_min))\n",
"plt.yscale(\"log\")\n",
"p.legend()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Part 3: Distribution of mass\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The masses we computed were the initial minimum and maximum star masses, whereas M22 is an old globular cluster. In the 12 Ga since M22 was formed, most of the larger, faster burning stars would have burned out, leaving stars ~ 1 solar mass or smaller which can remain on the main sequence for ~$10^{10}$ years. Stars remaining which are ~5 solar masses usually burn out in ~70 million years are likely next-generation stars formed in nebulae of stars which previously went supernova. Stars smaller than 0.08 solar masses are brown dwarfs and do not have enough mass for nuclear fusion, therefore would also not be resolved by our observations of M22."
]
}
],
"metadata": {
"file_extension": ".py",
"gist": {
"data": {
"description": "ast501-pset1.ipynb",
"public": false
},
"id": ""
},
"kernelspec": {
"display_name": "Python [conda env:astro]",
"language": "python",
"name": "conda-env-astro-py"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.8"
},
"mimetype": "text/x-python",
"name": "python",
"npconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": 3
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment