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))