Skip to content

Commit 2a5f9ce

Browse files
Jolien128Jolien128
committed
updates AccrDisk notebooks and script
Add functions and notebook to make plots of keplerian velocity, tangential velocity, and mean radial velocity ifo radius
1 parent 813b847 commit 2a5f9ce

File tree

5 files changed

+350
-452
lines changed

5 files changed

+350
-452
lines changed

‎docs/src/1_examples/2_accretion_disk/accrDisks_calcHRM_Plons.ipynb

Lines changed: 56 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,9 @@
1010
{
1111
"cell_type": "code",
1212
"execution_count": 1,
13-
"metadata": {},
13+
"metadata": {
14+
"tags": []
15+
},
1416
"outputs": [],
1517
"source": [
1618
"import plons\n",
@@ -26,7 +28,9 @@
2628
{
2729
"cell_type": "code",
2830
"execution_count": 2,
29-
"metadata": {},
31+
"metadata": {
32+
"tags": []
33+
},
3034
"outputs": [],
3135
"source": [
3236
"'''\n",
@@ -51,8 +55,11 @@
5155
" thetaArray = np.linspace(thetaMin,(numberOfThetas-1)/numberOfThetas*thetaMax,numberOfThetas)\n",
5256
"\n",
5357
" crit = 0.3\n",
54-
" lowerR = 0.1\n",
55-
" r_step = 0.1\n",
58+
" # lowerR = 0.1 #to test\n",
59+
" # r_step = 0.1 #to test\n",
60+
" \n",
61+
" lowerR = 0.02\n",
62+
" r_step = 0.01\n",
5663
"\n",
5764
" # Calculate arrays with radii, disk mass, relative addes mass, and scale heights\n",
5865
" r_array,totalMassDisk,rel_mass_added,SH_array = ad.get_SH_diskMass_radius(lowerR,r_step,thetaArray,dumpData,maxH,n_grid,run,crit,xH,phiQuadr = False, printOut = printOut)\n",
@@ -106,7 +113,9 @@
106113
{
107114
"cell_type": "code",
108115
"execution_count": 3,
109-
"metadata": {},
116+
"metadata": {
117+
"tags": []
118+
},
110119
"outputs": [],
111120
"source": [
112121
"'''\n",
@@ -128,7 +137,8 @@
128137
"\n",
129138
" # Choose which and how many r and theta values you want to use to calculate scale heights and disk mass and for plot of rho vs h\n",
130139
" n_grid = 10000\n",
131-
" maxH = 1.5 * cgs.au \n",
140+
" # maxH = 1.5 * cgs.au \n",
141+
" maxH = 3 * cgs.au \n",
132142
" numberOfThetas = 10 #Number of thetas withing one array \n",
133143
" \n",
134144
" '''\n",
@@ -161,8 +171,11 @@
161171
" thetaArray4 = np.linspace(thetaMin4,thetaMax4,numberOfThetas)\n",
162172
"\n",
163173
" crit = 0.3\n",
164-
" lowerR = 0.1\n",
165-
" r_step1 = 0.1\n",
174+
" # lowerR = 0.1 #to test\n",
175+
" # r_step1 = 0.1 #to test\n",
176+
" \n",
177+
" lowerR = 0.02\n",
178+
" r_step1 = 0.01\n",
166179
" \n",
167180
" if run == '/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/v10e00_T3000_res8_racc01/':\n",
168181
" r_step3 = 0.015 # Was nescessary for my model v10e00\n",
@@ -206,7 +219,9 @@
206219
{
207220
"cell_type": "code",
208221
"execution_count": 4,
209-
"metadata": {},
222+
"metadata": {
223+
"tags": []
224+
},
210225
"outputs": [],
211226
"source": [
212227
"'''\n",
@@ -358,6 +373,37 @@
358373
"main(run,dump,['1H'],full=True,quadrants=True,printOut=True)"
359374
]
360375
},
376+
{
377+
"cell_type": "code",
378+
"execution_count": 1,
379+
"metadata": {
380+
"tags": []
381+
},
382+
"outputs": [
383+
{
384+
"data": {
385+
"text/plain": [
386+
"'\\nExample\\n'"
387+
]
388+
},
389+
"execution_count": 1,
390+
"metadata": {},
391+
"output_type": "execute_result"
392+
}
393+
],
394+
"source": [
395+
"'''\n",
396+
"Example\n",
397+
"'''\n",
398+
"\n",
399+
"dir = '/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/'\n",
400+
"\n",
401+
"#Non-eccentric runs --> use final full dump; do calculations for both scale height options and for both the full thetaregion and the four quadrants seperately to study asymmetries\n",
402+
"run = os.path.join(dir, 'v20e00_T3000_res8_racc01/')\n",
403+
"dump = 292\n",
404+
"main(run,dump,['2H'],full=True,quadrants=True)"
405+
]
406+
},
361407
{
362408
"cell_type": "code",
363409
"execution_count": null,
@@ -443,7 +489,7 @@
443489
"name": "python",
444490
"nbconvert_exporter": "python",
445491
"pygments_lexer": "ipython3",
446-
"version": "3.11.5"
492+
"version": "3.10.13"
447493
}
448494
},
449495
"nbformat": 4,

‎docs/src/1_examples/2_accretion_disk/accrDisks_plotsHRM_Plons.ipynb

Lines changed: 22 additions & 61 deletions
Large diffs are not rendered by default.

‎docs/src/1_examples/2_accretion_disk/accrDisks_vrvrPlots_Plons.ipynb

Lines changed: 36 additions & 379 deletions
Large diffs are not rendered by default.

‎docs/src/1_examples/2_accretion_disk/accrDisks_vtKepl_Plons_final.ipynb

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "code",
5+
"execution_count": 1,
6+
"id": "9310d776-8884-41a6-a8e8-f71f94a75ad1",
7+
"metadata": {},
8+
"outputs": [],
9+
"source": [
10+
"import math as math\n",
11+
"import matplotlib.pyplot as plt\n",
12+
"import numpy as np\n",
13+
"# import necessary plons scripts\n",
14+
"import plons.AccrDisk as ad"
15+
]
16+
},
17+
{
18+
"cell_type": "code",
19+
"execution_count": 2,
20+
"id": "b3e09372-3c74-479a-8f07-58f34c9cb739",
21+
"metadata": {},
22+
"outputs": [
23+
{
24+
"name": "stderr",
25+
"output_type": "stream",
26+
"text": [
27+
"/lhome/jolienm/anaconda3/envs/plons/lib/python3.10/site-packages/numpy/core/fromnumeric.py:3464: RuntimeWarning: Mean of empty slice.\n",
28+
" return _methods._mean(a, axis=axis, dtype=dtype,\n",
29+
"/lhome/jolienm/anaconda3/envs/plons/lib/python3.10/site-packages/numpy/core/_methods.py:192: RuntimeWarning: invalid value encountered in scalar divide\n",
30+
" ret = ret.dtype.type(ret / rcount)\n"
31+
]
32+
}
33+
],
34+
"source": [
35+
"%matplotlib inline \n",
36+
"\n",
37+
"dump = 292\n",
38+
"models = ['v05e00','v10e00','v20e00']\n",
39+
"rs = [0.94, 0.75, 0.36] #for crit03\n",
40+
"colors = ['firebrick','goldenrod','navy']\n",
41+
"\n",
42+
"#rvalue to start from + also determines step for rr array\n",
43+
"rstep = 0.02\n",
44+
"rmin = rstep/2\n",
45+
"# rmin = 0.005\n",
46+
"fig, ax = plt.subplots()\n",
47+
"\n",
48+
"i = 0\n",
49+
"for model in models:\n",
50+
" run = '/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/'+str(model)+'_T3000_res8_racc01/'\n",
51+
" dumpData,setup = ad.loadDataForSmoothing(run,dump)\n",
52+
" rmax = rs[i]\n",
53+
" ad.plot_vt_divided_vtKepl(ax,rmax,rmin,rstep,dumpData,model,colors,i)\n",
54+
" i = i+1\n",
55+
"ax.legend(fontsize = 12)\n",
56+
"ax.set_xlabel(r'$r$ [au]',fontsize = 12)\n",
57+
"ax.set_ylabel(r'$v_{\\rm t}/v_{\\rm kepl}$ [km/s]',fontsize = 12,rotation = 90)\n",
58+
"ax.tick_params(axis='x', labelsize=12)\n",
59+
"ax.tick_params(axis='y', labelsize=12)\n",
60+
"# fig.show()\n",
61+
"plt.savefig('/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/plotsAnalysis/vt:vtKepl.png')\n",
62+
"plt.close()"
63+
]
64+
},
65+
{
66+
"cell_type": "code",
67+
"execution_count": 3,
68+
"id": "c90e4805-f642-4a55-9211-6239e81a569f",
69+
"metadata": {
70+
"tags": []
71+
},
72+
"outputs": [],
73+
"source": [
74+
"%matplotlib inline \n",
75+
"\n",
76+
"dump = 292\n",
77+
"fig, ax = plt.subplots()\n",
78+
"rstep = 0.02\n",
79+
"rmin = rstep/2\n",
80+
"models = ['v05e00','v10e00','v20e00']\n",
81+
"colors = ['firebrick','goldenrod','navy']\n",
82+
"lineStyles=['-','-','-','-']\n",
83+
"rs = [0.94, 0.75, 0.36] #for crit03\n",
84+
"k_vt = np.array([])\n",
85+
"r_vt_max = np.array([])\n",
86+
"i = 0\n",
87+
"for model in models:\n",
88+
" run = '/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/'+str(model)+'_T3000_res8_racc01/'\n",
89+
" dumpData,setup = ad.loadDataForSmoothing(run,dump)\n",
90+
" rmax = rs[i]#-rstep/2\n",
91+
" r_vt_max = ad.plot_vt_vKepl(ax,rmax,rmin,rstep,dumpData,model,k_vt,r_vt_max,colors,lineStyles,i)\n",
92+
" i = i+1\n",
93+
"ax.legend(fontsize = 12)\n",
94+
"ax.set_xlabel(r'$r$ [au]',fontsize = 12)\n",
95+
"ax.set_ylabel(r'$v_t$ [km/s]',fontsize = 12,rotation = 90)\n",
96+
"ax.set_ylim(1,np.max(r_vt_max)+40)\n",
97+
"ax.tick_params(axis='x', labelsize=12)\n",
98+
"ax.tick_params(axis='y', labelsize=12)\n",
99+
"# fig.show()\n",
100+
"plt.savefig('/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/plotsAnalysis/vt_vtKepl.png')\n",
101+
"plt.close()"
102+
]
103+
},
104+
{
105+
"cell_type": "code",
106+
"execution_count": 4,
107+
"id": "ca9f0a06-1dd2-4f0e-ae2a-7d8282b91f36",
108+
"metadata": {},
109+
"outputs": [],
110+
"source": [
111+
"%matplotlib inline \n",
112+
"\n",
113+
"dump = 292\n",
114+
"fig, ax = plt.subplots()\n",
115+
"rstep = 0.01\n",
116+
"rmin = rstep/2\n",
117+
"models = ['v05e00','v10e00','v20e00']\n",
118+
"colors = ['firebrick','goldenrod','navy']\n",
119+
"rs = [0.94, 0.75, 0.36] \n",
120+
"\n",
121+
"i = 0\n",
122+
"for model in models:\n",
123+
" run = '/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/'+str(model)+'_T3000_res8_racc01/'\n",
124+
" dumpData,setup = ad.loadDataForSmoothing(run,dump)\n",
125+
" rmax = rs[i]#-rstep/2\n",
126+
" ad.plot_vrmean(ax,rmax,rmin,rstep,dumpData,model,colors,i)\n",
127+
" i = i+1\n",
128+
"\n",
129+
"ax.hlines(0,0,np.max(rs),linewidth = 0.7,linestyle = '--',color = 'k') \n",
130+
"ax.legend(fontsize = 12)\n",
131+
"ax.set_xlabel(r'$r$ [au]',fontsize = 12)\n",
132+
"ax.set_ylabel(r'mean $v_r$ [km/s]',fontsize = 12,rotation = 90)\n",
133+
"ax.set_ylim(-10,2)\n",
134+
"ax.tick_params(axis='x', labelsize=12)\n",
135+
"ax.tick_params(axis='y', labelsize=12)\n",
136+
"# fig.show()\n",
137+
"plt.savefig('/lhome/jolienm/Documents/TierModels/R_Aql/cooling/binariesInPaper/finalAccrDisks/plotsAnalysis/mean_vr.png')\n",
138+
"plt.close()"
139+
]
140+
},
141+
{
142+
"cell_type": "code",
143+
"execution_count": null,
144+
"id": "ded43fa6-1de7-4a3c-af37-1fa0a5b44814",
145+
"metadata": {},
146+
"outputs": [],
147+
"source": []
148+
}
149+
],
150+
"metadata": {
151+
"kernelspec": {
152+
"display_name": "Python 3 (ipykernel)",
153+
"language": "python",
154+
"name": "python3"
155+
},
156+
"language_info": {
157+
"codemirror_mode": {
158+
"name": "ipython",
159+
"version": 3
160+
},
161+
"file_extension": ".py",
162+
"mimetype": "text/x-python",
163+
"name": "python",
164+
"nbconvert_exporter": "python",
165+
"pygments_lexer": "ipython3",
166+
"version": "3.10.13"
167+
}
168+
},
169+
"nbformat": 4,
170+
"nbformat_minor": 5
171+
}

‎src/plons/AccrDisk.py

Lines changed: 65 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,13 @@ def getRadTanVelocity(x,y,v_x,v_y):
6565

6666
return v_rad,v_tan
6767

68+
'''
69+
Calculate keplerian velocity from companion mass and radius
70+
'''
71+
def kepler_vt(r,dumpData):
72+
kep_vt = np.sqrt(cgs.G * dumpData['massComp'] / (r*cgs.au))/cgs.kms
73+
return kep_vt
74+
6875
def loadDataForSmoothing(run,dump):
6976
'''
7077
Load in data
@@ -181,6 +188,62 @@ def plot_vrvtRho(axs, smooth,zoom,r):
181188
ax3.add_patch(circle3)
182189

183190

191+
192+
# Plot of tangential velocity divided by keplerian velocity; in function of radius
193+
def plot_vt_divided_vtKepl(ax,rmax,rmin,rstep,dumpData,model,colors,i):
194+
# set number of rvalues for which we will plot velocities (rmin = rstep)
195+
rnum = int(rmax /rstep)
196+
# make rr array
197+
rr = np.linspace(rmin,rmax,rnum)
198+
# make empty tangential velocity array
199+
r_vt = np.array([])
200+
# make empty keplerian velocity array
201+
k_vt = np.array([])
202+
# for each r in the array, calculate Keplerian velocities
203+
for r in rr:
204+
k_vt = np.append(k_vt,kepler_vt(r,dumpData))
205+
# Calculate the mean tangential velocity of all data in r-step
206+
filter = (dumpData['new_r']/cgs.au<(r+(rmin/2))) & (dumpData['new_r']/cgs.au>(r-(rmin/2)))
207+
r_vt = np.append(r_vt, np.mean(dumpData['new_vt'] [filter]))
208+
ax.plot(rr,r_vt/k_vt,c=colors[i],linewidth = 1.2,label = str(model))
209+
210+
# Plot of tangential velocity and keplerian velocity, in function of radius
211+
def plot_vt_vKepl(ax,rmax,rmin,rstep,dumpData,model,k_vt,r_vt_max,colors,lineStyles,i):
212+
# set number of rvalues for which we will plot velocities
213+
rnum = int(rmax /rstep)
214+
# make rr array
215+
rr = np.linspace(rmin,rmax,rnum)
216+
# make empty tangential velocity array
217+
r_vt = np.array([])
218+
# for each r in the array, calculate tangential and keplerian velocities (Kepl only for longest model --> for i==0)
219+
for r in rr:
220+
if i == 0:
221+
k_vt = np.append(k_vt,kepler_vt(r,dumpData))
222+
filter = (dumpData['new_r']/cgs.au<(r+(rmin/2))) & (dumpData['new_r']/cgs.au>(r-(rmin/2)))
223+
r_vt = np.append(r_vt, np.mean(dumpData['new_vt'] [filter]))
224+
if i ==0:
225+
ax.plot(rr,k_vt,c='k',linestyle = '--',linewidth = 1,label = 'Keplerian')
226+
ax.plot(rr,k_vt*0.9,c='k',linestyle = ':',linewidth = 1,label = '90% (sub)Keplerian')
227+
ax.plot(rr,r_vt,c=colors[i],linestyle = lineStyles[i],linewidth = 1.2,label = str(model))
228+
r_vt_max = np.append(r_vt_max,np.nanmax(r_vt))
229+
return r_vt_max
230+
231+
'''
232+
Plot of mean of radial velocity (not absolute value!); in function of radius
233+
'''
234+
def plot_vrmean(ax,rmax,rmin,rstep,dumpData,model,colors,i):
235+
# set number of rvalues for which we will plot velocities
236+
rnum = int(rmax /rstep)
237+
# make rr array
238+
rr = np.linspace(rmin,rmax,rnum)
239+
r_vr = np.array([])
240+
for r in rr:
241+
filter = (dumpData['new_r']/cgs.au<(r+(rmin/2))) & (dumpData['new_r']/cgs.au>(r-(rmin/2)))
242+
r_vr = np.append(r_vr, np.mean(dumpData['new_vr'] [filter]))
243+
ax.plot(rr,r_vr,c=colors[i],linestyle = '-',linewidth = 1.2,label = str(model))
244+
245+
246+
184247
'''
185248
Adapted function from SmoothingKernelScript, to get pixCoord needed for function getSmoothingKernelledPix
186249
'''
@@ -346,14 +409,14 @@ def getSH_andPlot_one_r(r,thetaArray,dumpData,maxH,n_grid,run,xH):
346409

347410
def get_SH_diskMass_radius(lowerR,r_step,thetaArray,dumpData,maxH,n_grid,run,crit,xH,phiQuadr = False, printOut = False):
348411
# Make array to store the SH estimate for each r (0 at r=0.01 au (racc))
349-
SH_array = np.array([0.1])
412+
SH_array = np.array([lowerR-r_step])
350413
# Make array to store the total mass of the disk at each r (0 at r=0)
351414
totalMassDisk = np.array([0])
352415
# Make array to store the relative added mass of the disk at each r (1 at r=0)
353416
rel_mass_added = np.array([1])
354417

355418
# For each radius ri in the given array, but starting from r=racc
356-
r_array = np.append([0.1*cgs.au],lowerR*cgs.au)
419+
r_array = np.append([(lowerR-r_step)*cgs.au],lowerR*cgs.au)
357420
i=0
358421
# Calculate for each r the scale height and disk mass, until the relative mass added exceeds 0.01
359422
while rel_mass_added[-1]/r_step>crit:

0 commit comments

Comments
 (0)