aboutsummaryrefslogtreecommitdiffstats
path: root/mass_profile/analyze_entropy_profile.py
diff options
context:
space:
mode:
authorAaron LI <aaronly.me@gmail.com>2016-05-27 22:47:24 +0800
committerAaron LI <aaronly.me@gmail.com>2016-05-27 22:47:24 +0800
commitffd178e0bd72562a3c2cff9747b6e656edc881dc (patch)
tree8800b7b5b2e8bc3df1a6760df5cd54eaaa686702 /mass_profile/analyze_entropy_profile.py
parent5c35fad9240fb42c1371c721e0b2af7379bd9ea0 (diff)
downloadchandra-acis-analysis-ffd178e0bd72562a3c2cff9747b6e656edc881dc.tar.bz2
Add mass_profile tools
* These tools are mainly use to calculate the total gravitational mass profile, as well as the intermediate products (e.g., surface brightness profile fitting, gas density profile, NFW fitting, etc.) * There are additional tools for calculating the luminosity and flux. * These tools mainly developed by Junhua GU, and contributed by Weitian (Aaron) LI, and Zhenghao ZHU.
Diffstat (limited to 'mass_profile/analyze_entropy_profile.py')
-rwxr-xr-xmass_profile/analyze_entropy_profile.py54
1 files changed, 54 insertions, 0 deletions
diff --git a/mass_profile/analyze_entropy_profile.py b/mass_profile/analyze_entropy_profile.py
new file mode 100755
index 0000000..c123681
--- /dev/null
+++ b/mass_profile/analyze_entropy_profile.py
@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+
+import sys
+import numpy
+import scipy.interpolate
+
+center_entropy_file=open('entropy_center.qdp')
+entropy_file=open('summary_entropy.qdp')
+confidence_level=.68
+rout=float(sys.argv[1])
+
+center_s=0
+for i in center_entropy_file:
+ r,s=i.split()
+ r=float(r)
+ s=float(s)
+ if r>rout:
+ center_s=s
+ break
+
+new_data=True
+
+
+s_list=[]
+for i in entropy_file:
+ if i[0]=='n':
+ new_data=True
+ continue
+ if new_data==False:
+ continue
+ r,s=i.split()
+ r=float(r)
+ s=float(s)
+ if r>rout:
+ new_data=False
+ s_list.append(s)
+
+s_idx=-1
+
+s_list.sort()
+for i in range(len(s_list)-1):
+ if (center_s-s_list[i])*(center_s-s_list[i+1])<=0:
+ m_idx=i
+ break
+
+
+slidx=int(s_idx*(1-confidence_level))
+suidx=s_idx-1+int((len(s_list)-s_idx)*confidence_level)
+
+
+serr1=s_list[slidx]-center_s
+serr2=s_list[suidx]-center_s
+
+print("S=\t%e\t %e/+%e keV cm^2 (1 sigma)"%(center_s,serr1,serr2))