1 files changed, 195 insertions, 0 deletions
diff --git a/mass_profile/analyze_mass_profile.py b/mass_profile/analyze_mass_profile.py
new file mode 100755
index 0000000..3ee128c
--- /dev/null
+++ b/mass_profile/analyze_mass_profile.py
@@ -0,0 +1,195 @@
+#!/usr/bin/env python
+
+import sys
+import numpy
+import scipy.interpolate
+
+confidence_level=.68
+def read_file(param):
+    delta=float(param[0])
+
+    file_mass_center=open("mass_int_center.qdp").readlines();
+    file_delta_center=open("overdensity_center.qdp").readlines();
+    
+    center_r=0
+    center_m=0
+    center_gm=0
+    center_gf=0
+
+    
+    for i in range(0,len(file_mass_center)):
+        lm=file_mass_center[i].strip();
+        ld=file_delta_center[i].strip();
+        r,m=lm.split()
+        r,d=ld.split()
+        r=float(r)
+        d=float(d)
+        m=float(m)
+        if m<1e11:
+            continue
+        if d<delta:
+            center_r=r
+            center_m=m
+            for j in open("gas_mass_int_center.qdp"):
+                rgm,gm=j.strip().split()
+                rgm=float(rgm)
+                gm=float(gm)
+                if rgm>r:
+
+                    center_gm=gm
+                    center_gf=gm/m
+                    break
+            break
+    if len(param)>1 and param[1]=='c':
+        #print("%s(<r%d)=%E solar mass"%("mass",delta,center_m))
+        #print("%s%d=%E kpc"%("r",delta,center_r))
+        #print("%s(<r%d)=%E solar mass"%("gas mass",delta,center_gm))
+        #print("%s(<r%d)=%E"%("gas fraction",delta,center_gf))
+        return center_m,center_r,center_gm,center_gf,None,None,None,None
+    
+
+#print(center_gm,center_gf)
+    file_mass=open('summary_mass_profile.qdp').readlines()
+    file_delta=open('summary_overdensity.qdp').readlines()
+    file_gm=open('summary_gas_mass_profile.qdp')
+
+
+    flag=True
+    rlist=[]
+    mlist=[]
+    gmlist=[]
+    gflist=[]
+    old_m=0
+    invalid_count=0
+    for i in range(0,len(file_mass)):
+        lm=file_mass[i].strip()
+        ld=file_delta[i].strip()
+        if lm[0]=='n':
+            flag=True
+            old_m=0
+            continue
+        if not flag:
+            continue
+        r,m=lm.split()
+        m=float(m)
+        if m<1e12:
+            continue
+        if m<old_m:
+            invalid_count+=1
+            flag=False
+            continue
+        r,d=ld.split()
+        d=float(d)
+        r=float(r)
+
+        if d<delta:
+            #print("%s %e"%(d,m))
+            mlist.append(m)
+            rlist.append(r)
+            flag1=True
+            while True:
+                lgm=file_gm.readline().strip()
+                if lgm[0]=='n':
+                    break
+                rgm,gm=lgm.split()
+                rgm=float(rgm)
+                gm=float(gm)
+                if rgm>r and flag1:
+                    gmlist.append(gm)
+
+                    flag1=False
+                    gflist.append(gm/mlist[-1])
+                #print(gm,gflist[-1])
+            flag=False
+        old_m=m
+    print("%d abnormal data dropped"%(invalid_count))
+
+
+    return center_m,center_r,center_gm,center_gf,mlist,rlist,gmlist,gflist
+#center_m=numpy.mean(mlist)
+#center_r=numpy.mean(rlist)
+
+center_m,center_r,center_gm,center_gf,mlist,rlist,gmlist,gflist=read_file(sys.argv[1:])
+delta=float(sys.argv[1])
+
+if len(sys.argv)>2 and sys.argv[2]=='c':
+    print("%s(<r%d)=%E solar mass"%("mass",delta,center_m))
+    print("%s%d=%E kpc"%("r",delta,center_r))
+    print("%s(<r%d)=%E solar mass"%("gas mass",delta,center_gm))
+    print("%s(<r%d)=%E"%("gas fraction",delta,center_gf))
+    sys.exit(0)
+
+
+mlist.sort()
+rlist.sort()
+gflist.sort()
+gmlist.sort()
+
+m_idx=-1
+r_idx=-1
+gm_idx=-1
+gf_idx=-1
+delta=float(sys.argv[1])
+for i in range(len(mlist)-1):
+    if (center_m-mlist[i])*(center_m-mlist[i+1])<=0:
+        m_idx=i
+        break
+
+for i in range(len(rlist)-1):
+    if (center_r-rlist[i])*(center_r-rlist[i+1])<=0:
+        r_idx=i
+        break
+
+for i in range(len(gmlist)-1):
+    if (center_gm-gmlist[i])*(center_gm-gmlist[i+1])<=0:
+        gm_idx=i
+        break
+
+for i in range(len(gflist)-1):
+    if (center_gf-gflist[i])*(center_gf-gflist[i+1])<=0:
+        gf_idx=i
+        break
+
+
+if m_idx==-1 or r_idx==-1 or gf_idx==-1 or gm_idx==-1:
+    print("Error, the center value is not enclosed by the Monte-Carlo realizations, please check the result!")
+    print("m:%E %E %E"%(center_m,mlist[0],mlist[-1]))
+    print("gm:%E %E %E"%(center_gm,gmlist[0],gmlist[-1]))
+    print("gf:%E %E %E"%(center_gf,gflist[0],gflist[-1]))
+    print("r:%E %E %E"%(center_r,rlist[0],rlist[-1]))
+    sys.exit(1)
+
+
+mlidx=int(m_idx*(1-confidence_level))
+muidx=m_idx-1+int((len(mlist)-m_idx)*confidence_level)
+
+
+rlidx=int(r_idx*(1-confidence_level))
+ruidx=r_idx-1+int((len(rlist)-r_idx)*confidence_level)
+
+gmlidx=int(gm_idx*(1-confidence_level))
+gmuidx=gm_idx-1+int((len(gmlist)-gm_idx)*confidence_level)
+
+gflidx=int(gf_idx*(1-confidence_level))
+gfuidx=gf_idx-1+int((len(gflist)-gf_idx)*confidence_level)
+
+
+merr1=mlist[mlidx]-center_m
+merr2=mlist[muidx]-center_m
+
+rerr1=rlist[rlidx]-center_r
+rerr2=rlist[ruidx]-center_r
+
+gmerr1=gmlist[gmlidx]-center_gm
+gmerr2=gmlist[gmuidx]-center_gm
+
+gferr1=gflist[gflidx]-center_gf
+gferr2=gflist[gfuidx]-center_gf
+
+#print("%d %d %d"%(mlidx,m_idx,muidx))
+#print("%d %d %d"%(rlidx,r_idx,ruidx))
+
+print("m%d=\t%e\t %e/+%e solar mass (1 sigma)"%(delta,center_m,merr1,merr2))
+print("gas_m%d=\t%e\t %e/+%e solar mass (1 sigma)"%(delta,center_gm,gmerr1,gmerr2))
+print("gas_fraction%d=\t%e\t %e/+%e (1 sigma)"%(delta,center_gf,gferr1,gferr2))
+print("r%d=\t%d\t %d/+%d kpc (1 sigma)"%(delta,center_r,rerr1,rerr2))