diff options
author | Aaron LI <aaronly.me@outlook.com> | 2017-02-20 12:11:25 +0800 |
---|---|---|
committer | Aaron LI <aaronly.me@outlook.com> | 2017-02-20 12:11:25 +0800 |
commit | cba36462e9e70f45341e432074c726dda5e31092 (patch) | |
tree | 0bc9885d0ff86b563e0fb5321c56d51d5687ba70 /mass_profile/fg_2500_500.py | |
parent | 2a069ed00d6f1c83153be9174c296e5540f37d30 (diff) | |
download | chandra-acis-analysis-cba36462e9e70f45341e432074c726dda5e31092.tar.bz2 |
Move shell/python scripts from 'mass_profile/' to 'bin/'
Diffstat (limited to 'mass_profile/fg_2500_500.py')
-rwxr-xr-x | mass_profile/fg_2500_500.py | 153 |
1 files changed, 0 insertions, 153 deletions
diff --git a/mass_profile/fg_2500_500.py b/mass_profile/fg_2500_500.py deleted file mode 100755 index 67a1a11..0000000 --- a/mass_profile/fg_2500_500.py +++ /dev/null @@ -1,153 +0,0 @@ -#!/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) - -if len(sys.argv)>1: - center_m2500,center_r2500,center_gm2500,center_gf2500,mlist2500,rlist2500,gmlist2500,gflist2500=read_file([2500,sys.argv[1]]) - center_m500,center_r500,center_gm500,center_gf500,mlist500,rlist500,gmlist500,gflist500=read_file([500,sys.argv[1]]) -else: - center_m2500,center_r2500,center_gm2500,center_gf2500,mlist2500,rlist2500,gmlist2500,gflist2500=read_file([2500]) - center_m500,center_r500,center_gm500,center_gf500,mlist500,rlist500,gmlist500,gflist500=read_file([500]) - -if mlist2500!=None and len(mlist2500)!=len(mlist500): - raise Exception("Something wrong, the number of 2500 and 500 data are different") - - -if mlist2500==None: - print("gas fraction between r2500 and r500 is %E"%((center_gm500-center_gm2500)/(center_m500-center_m2500))) - sys.exit(0) - -gf_2500_500=[] - -for i in range(0,len(mlist500)): - if mlist500[i]-mlist2500[i]<=0: - continue - gf_2500_500.append((gmlist500[i]-gmlist2500[i])/(mlist500[i]-mlist2500[i])) - -gf_2500_500.sort(); - - -center_gf_2500_500=(center_gm500-center_gm2500)/(center_m500-center_m2500) -gf_idx=-1 - -for i in range(len(gf_2500_500)-1): - if (center_gf_2500_500-gf_2500_500[i])*(center_gf_2500_500-gf_2500_500[i+1])<=0: - gf_idx=i - break -if gf_idx==-1: - raise Exception("Something wrong!") - -gflidx=int(gf_idx*(1-confidence_level)) -gfuidx=gf_idx-1+int((len(gf_2500_500)-gf_idx)*confidence_level) - -gferr1=gf_2500_500[gflidx]-center_gf_2500_500 -gferr2=gf_2500_500[gfuidx]-center_gf_2500_500 - -print("gas_fraction between r2500 and r500=\t%e\t %e/+%e (1 sigma)"%(center_gf_2500_500,gferr1,gferr2)) |