import sys
#sys.path.insert(0,'D:/Python1/pydocument/seniorproject_quenching2/practice/_pycache_/session_1')
sys.path.insert(0,'D:/Python1/pydocument/seniorproject_quenching2/practice')
##导入自己编写的脚本，需要加入这两句，一句声明符号应用，然后声明需要引入文-
#件的路径具体到文件夹
#from read_m16 import read_m16_ds
#from read_m16 import read_m16_mass
#from read_m16 import test_read_m16_mass
from read_m16 import test_read_m16_ds
#导入的质量区间为11~15dex-Msolar
import os.path
#该句表示路径包调用
import numpy as np
import matplotlib.pyplot as plt
#import pandas as pd
from scipy.stats import t as st
#画误差棒函数须引入scipy库
def dolad_data(m,hz):
    _mh_path ='D:/Python1/pydocument/seniorproject_quenching2/practice/'
    fname = os.path.join(_mh_path,'data_mh1.txt')
    #上述两句表示文件调用并打开数据
    m = np.loadtxt(fname,delimiter=',',dtype=np.float,unpack=True)    
    print('m=',m)
    fname = os.path.join(_mh_path,'z_eff.txt')
    hz = np.loadtxt(fname,delimiter=',',dtype=np.float,unpack=True)
    #print('hz=',hz)
    #该处读入的红移为实际的红移数据
    #plt.plot(m[0,:],m[1,:],label='Doload-successfully')
    #plt.legend()
    #plt.show()
    lmw = len(m[:,0])
    lml = len(m[0,:])
    m_dex = np.linspace(11,15,lmw)
    # print('lmw=',lmw)
    # print('lml=',lml)
    # print('dex=',m_dex)
    return(m,m_dex,lmw,lml)
#dolad_data(m=True,hz=True)
#sction2:数据处理
##下面把deltasigma的计算转为Rp的函数关系，以方便调用
#part1:#下面这一部分为共动坐标情况（z=0）
def calcu_sigma(Rpc,m_,x):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
#a flat CDM model,omegak=0;omegalamd=0
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67*10**(-11)
    G = G_
  #  global ms_
  #  ms_ = 1.989*10**(30)
  #  ms = ms_  
    global c_
    c_ = 6
    ##调试结果c=3.25~3.75之间比较合适
    c = c_
    #下面开始计算
    #对导入的数据做单位转化转为:太阳质量/h
    m1 = h*m_*10**x
    L = len(Rpc)
    Rps = np.zeros(L,dtype=np.float)
    g_x = np.zeros(L,dtype=np.float)
    deltaSigma = np.zeros(L,dtype=np.float)
    Sigma = np.zeros(L,dtype=np.float)
    # rs =  np.zeros(LL,dtype=np.float)
    # r_200 =  np.zeros(LL,dtype=np.float)
    # rho_0 = np.zeros(LL,dtype=np.float)
    #单位修正因子
    #Qc = 0.1412#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    #Q200 = 0.7
    ##重新对单位修正因子订正如下：
    Qc = 3.084*10**(-2)/(1.989*h**2)#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    Q200 = 1#对r200的单位修正（此时对rs也做了修正）
    H = h*100
    #修正
    #E = np.sqrt((1+z)**3*omegam)
    #H = h*100*E
    ##修正可以选择两处，对于计算过程涉及的物理量修正，或者对于常数修正
    rhoc = Qc*(3*H**2)/(8*np.pi*G)
    rhom = 200*rhoc*omegam
    r_200 = Q200*(3*m1/(4*np.pi*200*rhoc))**(1/3)
    rs = r_200/c
    rho_0 = m1/((np.log(1+c)-c/(1+c))*4*np.pi*rs**3)
    for t in range(0,L):
        Rps[t] = Rpc[t]
        #引入中间函数
        f0 = Rps[t]/rs#这是考虑把R参数化的情况
        if Rps[t]<rs:
           f1 = np.arctanh(np.sqrt((1-f0)/(1+f0)))
           f2 = np.log(f0/2)
           f4 = f0**2*np.sqrt(1-f0**2)
           f5 = (f0**2-1)*np.sqrt(1-f0**2)
           g_x[t] = 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
           delta_segma = rs*rho_0*g_x[t]
          # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
          # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
           deltaSigma[t] = delta_segma*10**(-12)
           Sigma[t] = (2*rs*rho_0/(f0**2-1))*(1-2*f1/(np.sqrt(1-f0**2)))
        elif Rps[t]==rs:
             g_x[t] = 10/3+4*np.log(1/2)
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0/3
        else:
             f1 = np.arctan(np.sqrt((f0-1)/(f0+1)))
             f2 = np.log(f0/2)
             f4 = f0**2*np.sqrt(f0**2-1)
             f5 = (f0**2-1)**(3/2)
             g_x[t]= 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0*(1-2*f1/np.sqrt(f0**2-1))/(f0**2-1)
    return(Rpc,Sigma,deltaSigma)
#calcu_sigma(Rpc=True,m_=True,x=True)   
#part2:#下面这部分为考虑红移的情况
def calcu_sigmaz(Rpc,m_,x,z):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
#a flat CDM model,omegak=0;omegalamd=0
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67*10**(-11)
    G = G_
    #global ms_
    #ms_ = 1.989*10**(30)
    #ms = ms_  
    global c_
    c_ = 6
    ##调试结果c=3.25~3.75之间比较合适
    c = c_
    #下面开始计算
    #对导入的数据做单位转化转为:太阳质量/h
    m1 = h*m_*10**x
    L = len(Rpc)
    Rps = np.zeros(L,dtype=np.float)
    g_x = np.zeros(L,dtype=np.float)
    deltaSigma = np.zeros(L,dtype=np.float)
    Sigma = np.zeros(L,dtype=np.float)
    # rs =  np.zeros(LL,dtype=np.float)
    # r_200 =  np.zeros(LL,dtype=np.float)
    # rho_0 = np.zeros(LL,dtype=np.float)
    #单位修正因子
    #Qc = 0.1412#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    #Q200 = 0.7
    ##重新对单位修正因子订正如下：
    Qc = 3.084*10**(-2)/(1.989*h**2)#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    Q200 = 1#对r200的单位修正（此时对rs也做了修正）
    H = h*100
    #修正
    #E = np.sqrt((1+z)**3*omegam)
    #H = h*100*E
    ##修正可以选择两处，对于计算过程涉及的物理量修正，或者对于常数修正
    rhoc = (Qc*(3*H**2)/(8*np.pi*G))/(1/(1+z))**3
    rhom = 200*rhoc*omegam
    r_200 = Q200*(3*m1/(4*np.pi*200*rhoc))**(1/3)
    rs = r_200/c
    rho_0 = m1/((np.log(1+c)-c/(1+c))*4*np.pi*rs**3)
    for t in range(0,L):
        Rps[t] = Rpc[t]
        #引入中间函数
        f0 = Rps[t]/rs#这是考虑把R参数化的情况
        if Rps[t]<rs:
           f1 = np.arctanh(np.sqrt((1-f0)/(1+f0)))
           f2 = np.log(f0/2)
           f4 = f0**2*np.sqrt(1-f0**2)
           f5 = (f0**2-1)*np.sqrt(1-f0**2)
           g_x[t] = 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
           delta_segma = rs*rho_0*g_x[t]
          # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
          # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
           deltaSigma[t] = delta_segma*10**(-12)
           Sigma[t] = (2*rs*rho_0/(f0**2-1))*(1-2*f1/(np.sqrt(1-f0**2)))
        elif Rps[t]==rs:
             g_x[t] = 10/3+4*np.log(1/2)
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0/3
        else:
             f1 = np.arctan(np.sqrt((f0-1)/(f0+1)))
             f2 = np.log(f0/2)
             f4 = f0**2*np.sqrt(f0**2-1)
             f5 = (f0**2-1)**(3/2)
             g_x[t]= 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0*(1-2*f1/np.sqrt(f0**2-1))/(f0**2-1)
    return(Rpc,Sigma,deltaSigma,rs)
#calcu_sigmaz(Rpc=True,m_=True,x=True,z=True)
#单独调用模块时需要对此句修改，引入直接的Rp和质量
def fig_f(inm,Rp,r,rs,r_200,rho_R,sigma,deltasigma):
#def fig_f(ff):
#inm,Rs,r,rs,r_200,rho_R,sigma,deltasigma=semula_tion(m_,x)
    lm = len(inm)
    plt.figure()
    plt.subplot(121)
    # delta1 = np.zeros(lm,dtype=np.float)
    # delta2 = np.zeros(lm,dtype=np.float)
    for k in range(0,lm):
        plt.loglog(r[k,:],rho_R[k])
        '''
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
        '''
    plt.xlabel(r'$r-Mpc/h$')
    plt.ylabel(r'$\rho(R)-M_\odot/h$')
    plt.grid()
    print('mh_inte=',inm)
    plt.subplot(122)
    #下面引入具体的deltasigma数据的观测值,和所给质量比较
    test_read_m16_ds()
    delta1 = np.zeros(lm,dtype=np.float)
    delta2 = np.zeros(lm,dtype=np.float)
    for k in range(0,lm):
        plt.loglog(Rp[k,:],deltasigma[k,:])
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')
    plt.tight_layout()
    plt.show()
    out=[inm,rho_R,deltasigma,sigma]
    return(out)
#fig_f(inm=True,Rs=True,r=True,rs=True,r_200=True,rho_R=True,sigma=True,deltasigma=True)
#fig_f(ff=True)
##定义一个五个质量最佳预测和观测的对比图象
def fig_ff(Rpc,ds_sim,lmw,r):
    test_read_m16_ds()
    for k in range(0,lmw):
        pa = r[k,1]
        plt.plot(Rpc[:],ds_sim[k,pa,:],'-*')
        nn = ds_sim[k,pa,:].size#求样本大小
        x_mean = np.mean(ds_sim[k,pa,:])#求算术平均值
        x_std = np.std(ds_sim[k,pa,:])#求标准偏差
        x_se = x_std/np.sqrt(nn)#求标准误差
        dof = nn-1#自由度计算
        alpha = 1.0-0.95
        conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
        np.sqrt(1.+1./nn)#设置置信区间
        plt.errorbar(Rpc[:],ds_sim[k,pa,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_ff(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_ff(f=True)
##定义一个质量最佳预测和观测的对比图象
def fig_fff1(Rpc,ds_sim,k):
    test_read_m16_ds()
    plt.plot(Rpc[:],ds_sim[k,:],'-*')
    nn = ds_sim[k,:].size#求样本大小
    x_mean = np.mean(ds_sim[k,:])#求算术平均值
    x_std = np.std(ds_sim[k,:])#求标准偏差
    x_se = x_std/np.sqrt(nn)#求标准误差
    dof = nn-1#自由度计算
    alpha = 1.0-0.95
    conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
    np.sqrt(1.+1./nn)#设置置信区间
    plt.errorbar(Rpc[:],ds_sim[k,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_fff1(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_fff1(f=True)
def fig_fff2(Rpc,ds_sim,k,t):
    test_read_m16_ds()
    plt.plot(Rpc[:],ds_sim[k,t,:],'-*')
    nn = ds_sim[k,t,:].size#求样本大小
    x_mean = np.mean(ds_sim[k,t,:])#求算术平均值
    x_std = np.std(ds_sim[k,t,:])#求标准偏差
    x_se = x_std/np.sqrt(nn)#求标准误差
    dof = nn-1#自由度计算
    alpha = 1.0-0.95
    conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
    np.sqrt(1.+1./nn)#设置置信区间
    plt.errorbar(Rpc[:],ds_sim[k,t,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_fff2(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_fff2(f=True)
##
def fit_datar(y):
    h = 0.673
    #先找出观测值对应的rp
    rsa,dssa,ds_errsa = test_read_m16_ds()
    # print('ds=',dssa)
    # print('ds_err=',ds_errsa)
    #print('rp=',rsa)  
    a = np.shape(dssa)
    print('size a=',a)
    ##注意到观测数值选取的是物理坐标，需要靠里尺度银子修正，修正如下
    z_r = 0.230
    #z_r = 0.105
    a_r = 1/(1+z_r)
    ##此时对于预测的R也要做修正
    rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2])
    #rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=5]) 
    #rp = rsa[0,:]
    #该句直接对数组筛选，计算2Mpc以内（保证NFW模型适用的情况下）信号
    print(rp)
    b = len(rp)
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    #下面做两组数据的方差比较,注意不能对观测数据进行插值
    #比较方法，找出观测的sigma数值对应的Rp,再根据模型计算此时的模型数值Sigma（该步以完成）    
    ds_simr = np.zeros((lmw,lml,b),dtype=np.float)
    rr = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            m_ = m[k,t]
            x = m_dex[k]
            #Rpc = rp
            #Rpc,Sigma,deltaSigma = calcu_sigma(Rpc,m_,x)    
            #ds_simr[k,t,:] = deltaSigma
            #对比文献，加入尺度因子修正如下,把物理的距离转为共动的距离
            #计算模型在对应的投射距离上的预测信号取值
            Rpc = rp
            Rpc,Sigma,deltaSigma,rs = calcu_sigmaz(Rpc,m_,x,z_r)
            #对预测信号做相应的变化，把共动的转为物理的
            ds_simr[k,t,:] = deltaSigma
            rr[k,t] = rs
    yy = np.shape(ds_simr)
    #print('rr=',rr)#输出查看ds_sim的维度，即模型预测下的透镜信号    
    #比较观测的sigma和预测的Sigma,比较结果用fitr和fitb记录比较结果
    fitr = np.zeros((lmw,2),dtype=np.int)
    delta_r = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            d_r = 0
            for n in range(0,yy[2]):
                d_r = d_r+((ds_simr[k,t,n]-dssa[0,n])/ds_errsa[0,n])**2
            delta_r[k,t] = d_r
    #print(delta_r)
    #下面这段求每个质量级对应的方差最小值
    for k in range(0,lmw):
        #先把np的定义数组转为列表，再寻找索引号，tolist()表示内置函数调用
        aa = delta_r[k,:].tolist()
        xa = aa.index(min(delta_r[k,:]))
        fitr[k,0] = k
        fitr[k,1] = xa
    print(fitr)
    #下面做图比较几个最佳预测值与观测的对比情况
    Rpc = rp
    fig_ff(Rpc,ds_simr,lmw,fitr)
    plt.title('R-galaxy')
    plt.show() 
    #比较提取出的五个最小值的最小值为最后结果
    deltar = np.zeros(lmw,dtype=np.float)
    for k in range(0,lmw):
        deltar[k] = delta_r[fitr[k,0],fitr[k,1]]
    #下面用bestr和bestb记录最后比较结果
    bestfr = np.zeros(2,dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa,:]  
    bestr = m[bestfr[0],bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    #下面两句是为与文献对比做的修正
    dexr = m_dex[bestfr[0]]+np.log10(bestr*h) 
    #作图对比最佳情况
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp,ds_simr,k,t)
    plt.title('Red')
    #plt.savefig('fit_r.png',dpi=600)
    plt.show()
    print('x^2=',min(deltar))
    print('mr=',best_mr) #此句质量计算正常
    print('positionr=',bestfr)
    print('dexr=',dexr1)
    print('co_dexr=',dexr)  
    print('corr_mr=',10**dexr)
    return rp,best_mr,dexr1,a,b,yy,dexr,a_r,delta_r 
#fit_datar(y=True)
def fit_datab(y):
    h = 0.673
    #先找出观测值对应的rp
    rsa,dssa,ds_errsa = test_read_m16_ds()
    # print('ds=',dssa)
    # print('ds_err=',ds_errsa)
    #print('rp=',rsa)  
    a = np.shape(dssa)
    print('size a=',a)
    ##注意到观测数值选取的是物理坐标，需要靠里尺度银子修正，修正如下
    z_b = 0.246
    #z_b = 0.124
    a_b = 1/(1+z_b)
    ##此时对于预测的R也要做修正
    rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2]) 
   # rp = [rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2]
   # rp = rsa[0,:]
    #该句直接对数组筛选，计算2Mpc以内（保证NFW模型适用的情况下）信号
    print(rp)
    b = len(rp)
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    #下面做两组数据的方差比较,注意不能对观测数据进行插值
    #比较方法，找出观测的sigma数值对应的Rp,再根据模型计算此时的模型数值Sigma（该步以完成）    
    ds_simb = np.zeros((lmw,lml,b),dtype=np.float)
    rb = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            m_ = m[k,t]
            x = m_dex[k]
            #Rpc = rp
            #Rpc,Sigma,deltaSigma = calcu_sigma(Rpc,m_,x,z_b)
            #ds_simb[k,t,:] = deltaSigma
            #对比文献，把物理距离转为共动距离
            Rpc = rp
            #计算模型在对应的投射距离上的预测信号取值
            Rpc,Sigma,deltaSigma,rs = calcu_sigmaz(Rpc,m_,x,z_b)
            #预测信号也做修正
            ds_simb[k,t,:] = deltaSigma
            rb[k,t] = rs
    yy = np.shape(ds_simb)
    #print('rb=',rb)#输出查看ds_sim的维度，即模型预测下的透镜信号    
    #比较观测的sigma和预测的Sigma,比较结果用fitr和fitb记录比较结果
    fitb = np.zeros((lmw,2),dtype=np.int)
    delta_b = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            d_b = 0
            for n in range(0,yy[2]):
                d_b = d_b+((ds_simb[k,t,n]-dssa[1,n])/ds_errsa[1,n])**2
            delta_b[k,t] = d_b
    #print(delta_b)
    #下面这段求每个质量级对应的方差最小值
    for k in range(0,lmw):
        #先把np的定义数组转为列表，再寻找索引号，tolist()表示内置函数调用
        bb = delta_b[k,:].tolist()
        xb = bb.index(min(delta_b[k,:]))
        fitb[k,0] = k
        fitb[k,1] = xb
    print(fitb)
    #下面做图比较几个最佳预测值与观测的对比情况
    Rpc = rp
    fig_ff(Rpc,ds_simb,lmw,fitb)
    plt.title('B-galaxy')
    plt.show() 
    #比较提取出的五个最小值的最小值为最后结果
    deltab = np.zeros(lmw,dtype=np.float)
    for k in range(0,lmw):
        deltab[k] = delta_b[fitb[k,0],fitb[k,1]]
    #下面用bestr和bestb记录最后比较结果
    bestfb = np.zeros(2,dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb,:]    
    bestb = m[bestfb[0],bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    #下面是与文献对比做的修正
    dexb = m_dex[bestfb[0]]+np.log10(bestb*h)    
    #作图对比最佳情况
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp,ds_simb,k,t)
    plt.title('Blue')
    #plt.savefig('fit_b.png',dpi=600)
    plt.show()
    print('x^2=',min(deltab))    
    print('mb=',best_mb) #此句质量计算正常
    print('positionb=',bestfb)
    print('dexb=',dexb1)
    print('co_dexb=',dexb) 
    print('corr_mb=',10**dexb)
    return rp,best_mb,dexb1,a,b,yy,dexb,a_b,delta_b 
#fit_datab(y=True)
#下面做图比较最佳预测值与观测的对比情况
def fig_(T):
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    rp,best_mr,dexr1,a,b,yy,dexr,a_r,delta_r = fit_datar(y=True)
    rp,best_mb,dexb1,a,b,yy,dexb,a_b,delta_b = fit_datab(y=True)
    plt.subplot(121)
    plt.plot(m_dex,np.log10(delta_r))
    plt.title('R-galaxy')
    plt.xlabel(r'$log(\frac{M_h}{M_\odot})$')
    plt.ylabel(r'$log(\chi^2)$')
    plt.grid()
    plt.subplot(122)
    plt.plot(m_dex,np.log10(delta_b))
    plt.xlabel(r'$log(\frac{M_h}{M_\odot})$')
    plt.ylabel(r'$log(\chi^2)$')
    plt.title('B-galaxy')
    plt.grid()
    plt.tight_layout()
    plt.show()  
    print('mr=',dexr)
    print('mb=',dexb)
    return dexr,dexb
fig_(T=True)
'''
if __name__ == "__main__":
    dolad_data(m=True,hz=True)
    pass
'''
4106f11b53126636f268da277d2d5e2ad945c385
<mask token>


def fig_ff(Rpc, ds_sim, lmw, r):
    test_read_m16_ds()
    for k in range(0, lmw):
        pa = r[k, 1]
        plt.plot(Rpc[:], ds_sim[k, pa, :], '-*')
        nn = ds_sim[k, pa, :].size
        x_mean = np.mean(ds_sim[k, pa, :])
        x_std = np.std(ds_sim[k, pa, :])
        x_se = x_std / np.sqrt(nn)
        dof = nn - 1
        alpha = 1.0 - 0.95
        conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 
            1.0 / nn)
        plt.errorbar(Rpc[:], ds_sim[k, pa, :], yerr=x_std, fmt='-',
            linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff1(Rpc, ds_sim, k):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, :], '-*')
    nn = ds_sim[k, :].size
    x_mean = np.mean(ds_sim[k, :])
    x_std = np.std(ds_sim[k, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff2(Rpc, ds_sim, k, t):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, t, :], '-*')
    nn = ds_sim[k, t, :].size
    x_mean = np.mean(ds_sim[k, t, :])
    x_std = np.std(ds_sim[k, t, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, t, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fit_datar(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_r = 0.23
    a_r = 1 / (1 + z_r)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simr = np.zeros((lmw, lml, b), dtype=np.float)
    rr = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_r)
            ds_simr[k, t, :] = deltaSigma
            rr[k, t] = rs
    yy = np.shape(ds_simr)
    fitr = np.zeros((lmw, 2), dtype=np.int)
    delta_r = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_r = 0
            for n in range(0, yy[2]):
                d_r = d_r + ((ds_simr[k, t, n] - dssa[0, n]) / ds_errsa[0, n]
                    ) ** 2
            delta_r[k, t] = d_r
    for k in range(0, lmw):
        aa = delta_r[k, :].tolist()
        xa = aa.index(min(delta_r[k, :]))
        fitr[k, 0] = k
        fitr[k, 1] = xa
    print(fitr)
    Rpc = rp
    fig_ff(Rpc, ds_simr, lmw, fitr)
    plt.title('R-galaxy')
    plt.show()
    deltar = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltar[k] = delta_r[fitr[k, 0], fitr[k, 1]]
    bestfr = np.zeros(2, dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa, :]
    bestr = m[bestfr[0], bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    dexr = m_dex[bestfr[0]] + np.log10(bestr * h)
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp, ds_simr, k, t)
    plt.title('Red')
    plt.show()
    print('x^2=', min(deltar))
    print('mr=', best_mr)
    print('positionr=', bestfr)
    print('dexr=', dexr1)
    print('co_dexr=', dexr)
    print('corr_mr=', 10 ** dexr)
    return rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r


def fit_datab(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_b = 0.246
    a_b = 1 / (1 + z_b)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simb = np.zeros((lmw, lml, b), dtype=np.float)
    rb = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_b)
            ds_simb[k, t, :] = deltaSigma
            rb[k, t] = rs
    yy = np.shape(ds_simb)
    fitb = np.zeros((lmw, 2), dtype=np.int)
    delta_b = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_b = 0
            for n in range(0, yy[2]):
                d_b = d_b + ((ds_simb[k, t, n] - dssa[1, n]) / ds_errsa[1, n]
                    ) ** 2
            delta_b[k, t] = d_b
    for k in range(0, lmw):
        bb = delta_b[k, :].tolist()
        xb = bb.index(min(delta_b[k, :]))
        fitb[k, 0] = k
        fitb[k, 1] = xb
    print(fitb)
    Rpc = rp
    fig_ff(Rpc, ds_simb, lmw, fitb)
    plt.title('B-galaxy')
    plt.show()
    deltab = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltab[k] = delta_b[fitb[k, 0], fitb[k, 1]]
    bestfb = np.zeros(2, dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb, :]
    bestb = m[bestfb[0], bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    dexb = m_dex[bestfb[0]] + np.log10(bestb * h)
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp, ds_simb, k, t)
    plt.title('Blue')
    plt.show()
    print('x^2=', min(deltab))
    print('mb=', best_mb)
    print('positionb=', bestfb)
    print('dexb=', dexb1)
    print('co_dexb=', dexb)
    print('corr_mb=', 10 ** dexb)
    return rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b


<mask token>

<mask token>


def dolad_data(m, hz):
    _mh_path = 'D:/Python1/pydocument/seniorproject_quenching2/practice/'
    fname = os.path.join(_mh_path, 'data_mh1.txt')
    m = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    print('m=', m)
    fname = os.path.join(_mh_path, 'z_eff.txt')
    hz = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    lmw = len(m[:, 0])
    lml = len(m[0, :])
    m_dex = np.linspace(11, 15, lmw)
    return m, m_dex, lmw, lml


def calcu_sigma(Rpc, m_, x):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G)
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma


def calcu_sigmaz(Rpc, m_, x, z):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G) / (1 / (1 + z)) ** 3
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma, rs


def fig_f(inm, Rp, r, rs, r_200, rho_R, sigma, deltasigma):
    lm = len(inm)
    plt.figure()
    plt.subplot(121)
    for k in range(0, lm):
        plt.loglog(r[k, :], rho_R[k])
        """
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
        """
    plt.xlabel('$r-Mpc/h$')
    plt.ylabel('$\\rho(R)-M_\\odot/h$')
    plt.grid()
    print('mh_inte=', inm)
    plt.subplot(122)
    test_read_m16_ds()
    delta1 = np.zeros(lm, dtype=np.float)
    delta2 = np.zeros(lm, dtype=np.float)
    for k in range(0, lm):
        plt.loglog(Rp[k, :], deltasigma[k, :])
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k], ls='--', linewidth=0.5, color='red')
        plt.axvline(delta2[k], ls='--', linewidth=0.5, color='blue')
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    plt.tight_layout()
    plt.show()
    out = [inm, rho_R, deltasigma, sigma]
    return out


def fig_ff(Rpc, ds_sim, lmw, r):
    test_read_m16_ds()
    for k in range(0, lmw):
        pa = r[k, 1]
        plt.plot(Rpc[:], ds_sim[k, pa, :], '-*')
        nn = ds_sim[k, pa, :].size
        x_mean = np.mean(ds_sim[k, pa, :])
        x_std = np.std(ds_sim[k, pa, :])
        x_se = x_std / np.sqrt(nn)
        dof = nn - 1
        alpha = 1.0 - 0.95
        conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 
            1.0 / nn)
        plt.errorbar(Rpc[:], ds_sim[k, pa, :], yerr=x_std, fmt='-',
            linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff1(Rpc, ds_sim, k):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, :], '-*')
    nn = ds_sim[k, :].size
    x_mean = np.mean(ds_sim[k, :])
    x_std = np.std(ds_sim[k, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff2(Rpc, ds_sim, k, t):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, t, :], '-*')
    nn = ds_sim[k, t, :].size
    x_mean = np.mean(ds_sim[k, t, :])
    x_std = np.std(ds_sim[k, t, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, t, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fit_datar(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_r = 0.23
    a_r = 1 / (1 + z_r)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simr = np.zeros((lmw, lml, b), dtype=np.float)
    rr = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_r)
            ds_simr[k, t, :] = deltaSigma
            rr[k, t] = rs
    yy = np.shape(ds_simr)
    fitr = np.zeros((lmw, 2), dtype=np.int)
    delta_r = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_r = 0
            for n in range(0, yy[2]):
                d_r = d_r + ((ds_simr[k, t, n] - dssa[0, n]) / ds_errsa[0, n]
                    ) ** 2
            delta_r[k, t] = d_r
    for k in range(0, lmw):
        aa = delta_r[k, :].tolist()
        xa = aa.index(min(delta_r[k, :]))
        fitr[k, 0] = k
        fitr[k, 1] = xa
    print(fitr)
    Rpc = rp
    fig_ff(Rpc, ds_simr, lmw, fitr)
    plt.title('R-galaxy')
    plt.show()
    deltar = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltar[k] = delta_r[fitr[k, 0], fitr[k, 1]]
    bestfr = np.zeros(2, dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa, :]
    bestr = m[bestfr[0], bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    dexr = m_dex[bestfr[0]] + np.log10(bestr * h)
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp, ds_simr, k, t)
    plt.title('Red')
    plt.show()
    print('x^2=', min(deltar))
    print('mr=', best_mr)
    print('positionr=', bestfr)
    print('dexr=', dexr1)
    print('co_dexr=', dexr)
    print('corr_mr=', 10 ** dexr)
    return rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r


def fit_datab(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_b = 0.246
    a_b = 1 / (1 + z_b)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simb = np.zeros((lmw, lml, b), dtype=np.float)
    rb = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_b)
            ds_simb[k, t, :] = deltaSigma
            rb[k, t] = rs
    yy = np.shape(ds_simb)
    fitb = np.zeros((lmw, 2), dtype=np.int)
    delta_b = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_b = 0
            for n in range(0, yy[2]):
                d_b = d_b + ((ds_simb[k, t, n] - dssa[1, n]) / ds_errsa[1, n]
                    ) ** 2
            delta_b[k, t] = d_b
    for k in range(0, lmw):
        bb = delta_b[k, :].tolist()
        xb = bb.index(min(delta_b[k, :]))
        fitb[k, 0] = k
        fitb[k, 1] = xb
    print(fitb)
    Rpc = rp
    fig_ff(Rpc, ds_simb, lmw, fitb)
    plt.title('B-galaxy')
    plt.show()
    deltab = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltab[k] = delta_b[fitb[k, 0], fitb[k, 1]]
    bestfb = np.zeros(2, dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb, :]
    bestb = m[bestfb[0], bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    dexb = m_dex[bestfb[0]] + np.log10(bestb * h)
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp, ds_simb, k, t)
    plt.title('Blue')
    plt.show()
    print('x^2=', min(deltab))
    print('mb=', best_mb)
    print('positionb=', bestfb)
    print('dexb=', dexb1)
    print('co_dexb=', dexb)
    print('corr_mb=', 10 ** dexb)
    return rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b


def fig_(T):
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r = fit_datar(y=True)
    rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b = fit_datab(y=True)
    plt.subplot(121)
    plt.plot(m_dex, np.log10(delta_r))
    plt.title('R-galaxy')
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.grid()
    plt.subplot(122)
    plt.plot(m_dex, np.log10(delta_b))
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.title('B-galaxy')
    plt.grid()
    plt.tight_layout()
    plt.show()
    print('mr=', dexr)
    print('mb=', dexb)
    return dexr, dexb


<mask token>

<mask token>
sys.path.insert(0, 'D:/Python1/pydocument/seniorproject_quenching2/practice')
<mask token>


def dolad_data(m, hz):
    _mh_path = 'D:/Python1/pydocument/seniorproject_quenching2/practice/'
    fname = os.path.join(_mh_path, 'data_mh1.txt')
    m = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    print('m=', m)
    fname = os.path.join(_mh_path, 'z_eff.txt')
    hz = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    lmw = len(m[:, 0])
    lml = len(m[0, :])
    m_dex = np.linspace(11, 15, lmw)
    return m, m_dex, lmw, lml


def calcu_sigma(Rpc, m_, x):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G)
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma


def calcu_sigmaz(Rpc, m_, x, z):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G) / (1 / (1 + z)) ** 3
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma, rs


def fig_f(inm, Rp, r, rs, r_200, rho_R, sigma, deltasigma):
    lm = len(inm)
    plt.figure()
    plt.subplot(121)
    for k in range(0, lm):
        plt.loglog(r[k, :], rho_R[k])
        """
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
        """
    plt.xlabel('$r-Mpc/h$')
    plt.ylabel('$\\rho(R)-M_\\odot/h$')
    plt.grid()
    print('mh_inte=', inm)
    plt.subplot(122)
    test_read_m16_ds()
    delta1 = np.zeros(lm, dtype=np.float)
    delta2 = np.zeros(lm, dtype=np.float)
    for k in range(0, lm):
        plt.loglog(Rp[k, :], deltasigma[k, :])
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k], ls='--', linewidth=0.5, color='red')
        plt.axvline(delta2[k], ls='--', linewidth=0.5, color='blue')
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    plt.tight_layout()
    plt.show()
    out = [inm, rho_R, deltasigma, sigma]
    return out


def fig_ff(Rpc, ds_sim, lmw, r):
    test_read_m16_ds()
    for k in range(0, lmw):
        pa = r[k, 1]
        plt.plot(Rpc[:], ds_sim[k, pa, :], '-*')
        nn = ds_sim[k, pa, :].size
        x_mean = np.mean(ds_sim[k, pa, :])
        x_std = np.std(ds_sim[k, pa, :])
        x_se = x_std / np.sqrt(nn)
        dof = nn - 1
        alpha = 1.0 - 0.95
        conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 
            1.0 / nn)
        plt.errorbar(Rpc[:], ds_sim[k, pa, :], yerr=x_std, fmt='-',
            linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff1(Rpc, ds_sim, k):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, :], '-*')
    nn = ds_sim[k, :].size
    x_mean = np.mean(ds_sim[k, :])
    x_std = np.std(ds_sim[k, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff2(Rpc, ds_sim, k, t):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, t, :], '-*')
    nn = ds_sim[k, t, :].size
    x_mean = np.mean(ds_sim[k, t, :])
    x_std = np.std(ds_sim[k, t, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, t, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fit_datar(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_r = 0.23
    a_r = 1 / (1 + z_r)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simr = np.zeros((lmw, lml, b), dtype=np.float)
    rr = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_r)
            ds_simr[k, t, :] = deltaSigma
            rr[k, t] = rs
    yy = np.shape(ds_simr)
    fitr = np.zeros((lmw, 2), dtype=np.int)
    delta_r = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_r = 0
            for n in range(0, yy[2]):
                d_r = d_r + ((ds_simr[k, t, n] - dssa[0, n]) / ds_errsa[0, n]
                    ) ** 2
            delta_r[k, t] = d_r
    for k in range(0, lmw):
        aa = delta_r[k, :].tolist()
        xa = aa.index(min(delta_r[k, :]))
        fitr[k, 0] = k
        fitr[k, 1] = xa
    print(fitr)
    Rpc = rp
    fig_ff(Rpc, ds_simr, lmw, fitr)
    plt.title('R-galaxy')
    plt.show()
    deltar = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltar[k] = delta_r[fitr[k, 0], fitr[k, 1]]
    bestfr = np.zeros(2, dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa, :]
    bestr = m[bestfr[0], bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    dexr = m_dex[bestfr[0]] + np.log10(bestr * h)
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp, ds_simr, k, t)
    plt.title('Red')
    plt.show()
    print('x^2=', min(deltar))
    print('mr=', best_mr)
    print('positionr=', bestfr)
    print('dexr=', dexr1)
    print('co_dexr=', dexr)
    print('corr_mr=', 10 ** dexr)
    return rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r


def fit_datab(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_b = 0.246
    a_b = 1 / (1 + z_b)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simb = np.zeros((lmw, lml, b), dtype=np.float)
    rb = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_b)
            ds_simb[k, t, :] = deltaSigma
            rb[k, t] = rs
    yy = np.shape(ds_simb)
    fitb = np.zeros((lmw, 2), dtype=np.int)
    delta_b = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_b = 0
            for n in range(0, yy[2]):
                d_b = d_b + ((ds_simb[k, t, n] - dssa[1, n]) / ds_errsa[1, n]
                    ) ** 2
            delta_b[k, t] = d_b
    for k in range(0, lmw):
        bb = delta_b[k, :].tolist()
        xb = bb.index(min(delta_b[k, :]))
        fitb[k, 0] = k
        fitb[k, 1] = xb
    print(fitb)
    Rpc = rp
    fig_ff(Rpc, ds_simb, lmw, fitb)
    plt.title('B-galaxy')
    plt.show()
    deltab = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltab[k] = delta_b[fitb[k, 0], fitb[k, 1]]
    bestfb = np.zeros(2, dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb, :]
    bestb = m[bestfb[0], bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    dexb = m_dex[bestfb[0]] + np.log10(bestb * h)
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp, ds_simb, k, t)
    plt.title('Blue')
    plt.show()
    print('x^2=', min(deltab))
    print('mb=', best_mb)
    print('positionb=', bestfb)
    print('dexb=', dexb1)
    print('co_dexb=', dexb)
    print('corr_mb=', 10 ** dexb)
    return rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b


def fig_(T):
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r = fit_datar(y=True)
    rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b = fit_datab(y=True)
    plt.subplot(121)
    plt.plot(m_dex, np.log10(delta_r))
    plt.title('R-galaxy')
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.grid()
    plt.subplot(122)
    plt.plot(m_dex, np.log10(delta_b))
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.title('B-galaxy')
    plt.grid()
    plt.tight_layout()
    plt.show()
    print('mr=', dexr)
    print('mb=', dexb)
    return dexr, dexb


fig_(T=True)
<mask token>

import sys
sys.path.insert(0, 'D:/Python1/pydocument/seniorproject_quenching2/practice')
from read_m16 import test_read_m16_ds
import os.path
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t as st


def dolad_data(m, hz):
    _mh_path = 'D:/Python1/pydocument/seniorproject_quenching2/practice/'
    fname = os.path.join(_mh_path, 'data_mh1.txt')
    m = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    print('m=', m)
    fname = os.path.join(_mh_path, 'z_eff.txt')
    hz = np.loadtxt(fname, delimiter=',', dtype=np.float, unpack=True)
    lmw = len(m[:, 0])
    lml = len(m[0, :])
    m_dex = np.linspace(11, 15, lmw)
    return m, m_dex, lmw, lml


def calcu_sigma(Rpc, m_, x):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G)
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma


def calcu_sigmaz(Rpc, m_, x, z):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67 * 10 ** -11
    G = G_
    global c_
    c_ = 6
    c = c_
    m1 = h * m_ * 10 ** x
    L = len(Rpc)
    Rps = np.zeros(L, dtype=np.float)
    g_x = np.zeros(L, dtype=np.float)
    deltaSigma = np.zeros(L, dtype=np.float)
    Sigma = np.zeros(L, dtype=np.float)
    Qc = 3.084 * 10 ** -2 / (1.989 * h ** 2)
    Q200 = 1
    H = h * 100
    rhoc = Qc * (3 * H ** 2) / (8 * np.pi * G) / (1 / (1 + z)) ** 3
    rhom = 200 * rhoc * omegam
    r_200 = Q200 * (3 * m1 / (4 * np.pi * 200 * rhoc)) ** (1 / 3)
    rs = r_200 / c
    rho_0 = m1 / ((np.log(1 + c) - c / (1 + c)) * 4 * np.pi * rs ** 3)
    for t in range(0, L):
        Rps[t] = Rpc[t]
        f0 = Rps[t] / rs
        if Rps[t] < rs:
            f1 = np.arctanh(np.sqrt((1 - f0) / (1 + f0)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(1 - f0 ** 2)
            f5 = (f0 ** 2 - 1) * np.sqrt(1 - f0 ** 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / (f0 ** 2 - 1) * (1 - 2 * f1 / np.
                sqrt(1 - f0 ** 2))
        elif Rps[t] == rs:
            g_x[t] = 10 / 3 + 4 * np.log(1 / 2)
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 / 3
        else:
            f1 = np.arctan(np.sqrt((f0 - 1) / (f0 + 1)))
            f2 = np.log(f0 / 2)
            f4 = f0 ** 2 * np.sqrt(f0 ** 2 - 1)
            f5 = (f0 ** 2 - 1) ** (3 / 2)
            g_x[t] = 8 * f1 / f4 + 4 * f2 * f0 ** -2 - 2 / (f0 ** 2 - 1
                ) + 4 * f1 / f5
            delta_segma = rs * rho_0 * g_x[t]
            deltaSigma[t] = delta_segma * 10 ** -12
            Sigma[t] = 2 * rs * rho_0 * (1 - 2 * f1 / np.sqrt(f0 ** 2 - 1)) / (
                f0 ** 2 - 1)
    return Rpc, Sigma, deltaSigma, rs


def fig_f(inm, Rp, r, rs, r_200, rho_R, sigma, deltasigma):
    lm = len(inm)
    plt.figure()
    plt.subplot(121)
    for k in range(0, lm):
        plt.loglog(r[k, :], rho_R[k])
        """
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
        """
    plt.xlabel('$r-Mpc/h$')
    plt.ylabel('$\\rho(R)-M_\\odot/h$')
    plt.grid()
    print('mh_inte=', inm)
    plt.subplot(122)
    test_read_m16_ds()
    delta1 = np.zeros(lm, dtype=np.float)
    delta2 = np.zeros(lm, dtype=np.float)
    for k in range(0, lm):
        plt.loglog(Rp[k, :], deltasigma[k, :])
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k], ls='--', linewidth=0.5, color='red')
        plt.axvline(delta2[k], ls='--', linewidth=0.5, color='blue')
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    plt.tight_layout()
    plt.show()
    out = [inm, rho_R, deltasigma, sigma]
    return out


def fig_ff(Rpc, ds_sim, lmw, r):
    test_read_m16_ds()
    for k in range(0, lmw):
        pa = r[k, 1]
        plt.plot(Rpc[:], ds_sim[k, pa, :], '-*')
        nn = ds_sim[k, pa, :].size
        x_mean = np.mean(ds_sim[k, pa, :])
        x_std = np.std(ds_sim[k, pa, :])
        x_se = x_std / np.sqrt(nn)
        dof = nn - 1
        alpha = 1.0 - 0.95
        conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 
            1.0 / nn)
        plt.errorbar(Rpc[:], ds_sim[k, pa, :], yerr=x_std, fmt='-',
            linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff1(Rpc, ds_sim, k):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, :], '-*')
    nn = ds_sim[k, :].size
    x_mean = np.mean(ds_sim[k, :])
    x_std = np.std(ds_sim[k, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fig_fff2(Rpc, ds_sim, k, t):
    test_read_m16_ds()
    plt.plot(Rpc[:], ds_sim[k, t, :], '-*')
    nn = ds_sim[k, t, :].size
    x_mean = np.mean(ds_sim[k, t, :])
    x_std = np.std(ds_sim[k, t, :])
    x_se = x_std / np.sqrt(nn)
    dof = nn - 1
    alpha = 1.0 - 0.95
    conf_region = st.ppf(1 - alpha / 2.0, dof) * x_std * np.sqrt(1.0 + 1.0 / nn
        )
    plt.errorbar(Rpc[:], ds_sim[k, t, :], yerr=x_std, fmt='-', linewidth=0.5)
    plt.xlabel('$R-Mpc/h$')
    plt.ylabel('$\\Delta\\Sigma-hM_\\odot/{pc^2}$')
    return ()


def fit_datar(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_r = 0.23
    a_r = 1 / (1 + z_r)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simr = np.zeros((lmw, lml, b), dtype=np.float)
    rr = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_r)
            ds_simr[k, t, :] = deltaSigma
            rr[k, t] = rs
    yy = np.shape(ds_simr)
    fitr = np.zeros((lmw, 2), dtype=np.int)
    delta_r = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_r = 0
            for n in range(0, yy[2]):
                d_r = d_r + ((ds_simr[k, t, n] - dssa[0, n]) / ds_errsa[0, n]
                    ) ** 2
            delta_r[k, t] = d_r
    for k in range(0, lmw):
        aa = delta_r[k, :].tolist()
        xa = aa.index(min(delta_r[k, :]))
        fitr[k, 0] = k
        fitr[k, 1] = xa
    print(fitr)
    Rpc = rp
    fig_ff(Rpc, ds_simr, lmw, fitr)
    plt.title('R-galaxy')
    plt.show()
    deltar = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltar[k] = delta_r[fitr[k, 0], fitr[k, 1]]
    bestfr = np.zeros(2, dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa, :]
    bestr = m[bestfr[0], bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    dexr = m_dex[bestfr[0]] + np.log10(bestr * h)
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp, ds_simr, k, t)
    plt.title('Red')
    plt.show()
    print('x^2=', min(deltar))
    print('mr=', best_mr)
    print('positionr=', bestfr)
    print('dexr=', dexr1)
    print('co_dexr=', dexr)
    print('corr_mr=', 10 ** dexr)
    return rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r


def fit_datab(y):
    h = 0.673
    rsa, dssa, ds_errsa = test_read_m16_ds()
    a = np.shape(dssa)
    print('size a=', a)
    z_b = 0.246
    a_b = 1 / (1 + z_b)
    rp = np.array([rsa[0, k] for k in range(0, len(rsa[0, :])) if rsa[0, k] *
        h <= 2])
    print(rp)
    b = len(rp)
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    ds_simb = np.zeros((lmw, lml, b), dtype=np.float)
    rb = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            m_ = m[k, t]
            x = m_dex[k]
            Rpc = rp
            Rpc, Sigma, deltaSigma, rs = calcu_sigmaz(Rpc, m_, x, z_b)
            ds_simb[k, t, :] = deltaSigma
            rb[k, t] = rs
    yy = np.shape(ds_simb)
    fitb = np.zeros((lmw, 2), dtype=np.int)
    delta_b = np.zeros((lmw, lml), dtype=np.float)
    for k in range(0, lmw):
        for t in range(0, lml):
            d_b = 0
            for n in range(0, yy[2]):
                d_b = d_b + ((ds_simb[k, t, n] - dssa[1, n]) / ds_errsa[1, n]
                    ) ** 2
            delta_b[k, t] = d_b
    for k in range(0, lmw):
        bb = delta_b[k, :].tolist()
        xb = bb.index(min(delta_b[k, :]))
        fitb[k, 0] = k
        fitb[k, 1] = xb
    print(fitb)
    Rpc = rp
    fig_ff(Rpc, ds_simb, lmw, fitb)
    plt.title('B-galaxy')
    plt.show()
    deltab = np.zeros(lmw, dtype=np.float)
    for k in range(0, lmw):
        deltab[k] = delta_b[fitb[k, 0], fitb[k, 1]]
    bestfb = np.zeros(2, dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb, :]
    bestb = m[bestfb[0], bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    dexb = m_dex[bestfb[0]] + np.log10(bestb * h)
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp, ds_simb, k, t)
    plt.title('Blue')
    plt.show()
    print('x^2=', min(deltab))
    print('mb=', best_mb)
    print('positionb=', bestfb)
    print('dexb=', dexb1)
    print('co_dexb=', dexb)
    print('corr_mb=', 10 ** dexb)
    return rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b


def fig_(T):
    m, m_dex, lmw, lml = dolad_data(m=True, hz=True)
    rp, best_mr, dexr1, a, b, yy, dexr, a_r, delta_r = fit_datar(y=True)
    rp, best_mb, dexb1, a, b, yy, dexb, a_b, delta_b = fit_datab(y=True)
    plt.subplot(121)
    plt.plot(m_dex, np.log10(delta_r))
    plt.title('R-galaxy')
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.grid()
    plt.subplot(122)
    plt.plot(m_dex, np.log10(delta_b))
    plt.xlabel('$log(\\frac{M_h}{M_\\odot})$')
    plt.ylabel('$log(\\chi^2)$')
    plt.title('B-galaxy')
    plt.grid()
    plt.tight_layout()
    plt.show()
    print('mr=', dexr)
    print('mb=', dexb)
    return dexr, dexb


fig_(T=True)
<mask token>

import sys
#sys.path.insert(0,'D:/Python1/pydocument/seniorproject_quenching2/practice/_pycache_/session_1')
sys.path.insert(0,'D:/Python1/pydocument/seniorproject_quenching2/practice')
##导入自己编写的脚本，需要加入这两句，一句声明符号应用，然后声明需要引入文-
#件的路径具体到文件夹
#from read_m16 import read_m16_ds
#from read_m16 import read_m16_mass
#from read_m16 import test_read_m16_mass
from read_m16 import test_read_m16_ds
#导入的质量区间为11~15dex-Msolar
import os.path
#该句表示路径包调用
import numpy as np
import matplotlib.pyplot as plt
#import pandas as pd
from scipy.stats import t as st
#画误差棒函数须引入scipy库
def dolad_data(m,hz):
    _mh_path ='D:/Python1/pydocument/seniorproject_quenching2/practice/'
    fname = os.path.join(_mh_path,'data_mh1.txt')
    #上述两句表示文件调用并打开数据
    m = np.loadtxt(fname,delimiter=',',dtype=np.float,unpack=True)    
    print('m=',m)
    fname = os.path.join(_mh_path,'z_eff.txt')
    hz = np.loadtxt(fname,delimiter=',',dtype=np.float,unpack=True)
    #print('hz=',hz)
    #该处读入的红移为实际的红移数据
    #plt.plot(m[0,:],m[1,:],label='Doload-successfully')
    #plt.legend()
    #plt.show()
    lmw = len(m[:,0])
    lml = len(m[0,:])
    m_dex = np.linspace(11,15,lmw)
    # print('lmw=',lmw)
    # print('lml=',lml)
    # print('dex=',m_dex)
    return(m,m_dex,lmw,lml)
#dolad_data(m=True,hz=True)
#sction2:数据处理
##下面把deltasigma的计算转为Rp的函数关系，以方便调用
#part1:#下面这一部分为共动坐标情况（z=0）
def calcu_sigma(Rpc,m_,x):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
#a flat CDM model,omegak=0;omegalamd=0
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67*10**(-11)
    G = G_
  #  global ms_
  #  ms_ = 1.989*10**(30)
  #  ms = ms_  
    global c_
    c_ = 6
    ##调试结果c=3.25~3.75之间比较合适
    c = c_
    #下面开始计算
    #对导入的数据做单位转化转为:太阳质量/h
    m1 = h*m_*10**x
    L = len(Rpc)
    Rps = np.zeros(L,dtype=np.float)
    g_x = np.zeros(L,dtype=np.float)
    deltaSigma = np.zeros(L,dtype=np.float)
    Sigma = np.zeros(L,dtype=np.float)
    # rs =  np.zeros(LL,dtype=np.float)
    # r_200 =  np.zeros(LL,dtype=np.float)
    # rho_0 = np.zeros(LL,dtype=np.float)
    #单位修正因子
    #Qc = 0.1412#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    #Q200 = 0.7
    ##重新对单位修正因子订正如下：
    Qc = 3.084*10**(-2)/(1.989*h**2)#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    Q200 = 1#对r200的单位修正（此时对rs也做了修正）
    H = h*100
    #修正
    #E = np.sqrt((1+z)**3*omegam)
    #H = h*100*E
    ##修正可以选择两处，对于计算过程涉及的物理量修正，或者对于常数修正
    rhoc = Qc*(3*H**2)/(8*np.pi*G)
    rhom = 200*rhoc*omegam
    r_200 = Q200*(3*m1/(4*np.pi*200*rhoc))**(1/3)
    rs = r_200/c
    rho_0 = m1/((np.log(1+c)-c/(1+c))*4*np.pi*rs**3)
    for t in range(0,L):
        Rps[t] = Rpc[t]
        #引入中间函数
        f0 = Rps[t]/rs#这是考虑把R参数化的情况
        if Rps[t]<rs:
           f1 = np.arctanh(np.sqrt((1-f0)/(1+f0)))
           f2 = np.log(f0/2)
           f4 = f0**2*np.sqrt(1-f0**2)
           f5 = (f0**2-1)*np.sqrt(1-f0**2)
           g_x[t] = 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
           delta_segma = rs*rho_0*g_x[t]
          # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
          # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
           deltaSigma[t] = delta_segma*10**(-12)
           Sigma[t] = (2*rs*rho_0/(f0**2-1))*(1-2*f1/(np.sqrt(1-f0**2)))
        elif Rps[t]==rs:
             g_x[t] = 10/3+4*np.log(1/2)
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0/3
        else:
             f1 = np.arctan(np.sqrt((f0-1)/(f0+1)))
             f2 = np.log(f0/2)
             f4 = f0**2*np.sqrt(f0**2-1)
             f5 = (f0**2-1)**(3/2)
             g_x[t]= 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0*(1-2*f1/np.sqrt(f0**2-1))/(f0**2-1)
    return(Rpc,Sigma,deltaSigma)
#calcu_sigma(Rpc=True,m_=True,x=True)   
#part2:#下面这部分为考虑红移的情况
def calcu_sigmaz(Rpc,m_,x,z):
    global omega_m
    omega_m = 0.315
    omegam = omega_m
#a flat CDM model,omegak=0;omegalamd=0
    global h_
    h_ = 0.673
    h = h_
    global G_
    G_ = 6.67*10**(-11)
    G = G_
    #global ms_
    #ms_ = 1.989*10**(30)
    #ms = ms_  
    global c_
    c_ = 6
    ##调试结果c=3.25~3.75之间比较合适
    c = c_
    #下面开始计算
    #对导入的数据做单位转化转为:太阳质量/h
    m1 = h*m_*10**x
    L = len(Rpc)
    Rps = np.zeros(L,dtype=np.float)
    g_x = np.zeros(L,dtype=np.float)
    deltaSigma = np.zeros(L,dtype=np.float)
    Sigma = np.zeros(L,dtype=np.float)
    # rs =  np.zeros(LL,dtype=np.float)
    # r_200 =  np.zeros(LL,dtype=np.float)
    # rho_0 = np.zeros(LL,dtype=np.float)
    #单位修正因子
    #Qc = 0.1412#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    #Q200 = 0.7
    ##重新对单位修正因子订正如下：
    Qc = 3.084*10**(-2)/(1.989*h**2)#对rouc的单位修正（此时对roum、rou0的也完成了修正）
    Q200 = 1#对r200的单位修正（此时对rs也做了修正）
    H = h*100
    #修正
    #E = np.sqrt((1+z)**3*omegam)
    #H = h*100*E
    ##修正可以选择两处，对于计算过程涉及的物理量修正，或者对于常数修正
    rhoc = (Qc*(3*H**2)/(8*np.pi*G))/(1/(1+z))**3
    rhom = 200*rhoc*omegam
    r_200 = Q200*(3*m1/(4*np.pi*200*rhoc))**(1/3)
    rs = r_200/c
    rho_0 = m1/((np.log(1+c)-c/(1+c))*4*np.pi*rs**3)
    for t in range(0,L):
        Rps[t] = Rpc[t]
        #引入中间函数
        f0 = Rps[t]/rs#这是考虑把R参数化的情况
        if Rps[t]<rs:
           f1 = np.arctanh(np.sqrt((1-f0)/(1+f0)))
           f2 = np.log(f0/2)
           f4 = f0**2*np.sqrt(1-f0**2)
           f5 = (f0**2-1)*np.sqrt(1-f0**2)
           g_x[t] = 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
           delta_segma = rs*rho_0*g_x[t]
          # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
          # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
           deltaSigma[t] = delta_segma*10**(-12)
           Sigma[t] = (2*rs*rho_0/(f0**2-1))*(1-2*f1/(np.sqrt(1-f0**2)))
        elif Rps[t]==rs:
             g_x[t] = 10/3+4*np.log(1/2)
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0/3
        else:
             f1 = np.arctan(np.sqrt((f0-1)/(f0+1)))
             f2 = np.log(f0/2)
             f4 = f0**2*np.sqrt(f0**2-1)
             f5 = (f0**2-1)**(3/2)
             g_x[t]= 8*f1/f4+4*f2*f0**(-2)-2/(f0**2-1)+4*f1/f5
             delta_segma = rs*rho_0*g_x[t]
            # deltasigma[n,t] = delta_segma  #(原有计算,面积单位为Mpc^2)
            # 为了和观测对比，把单位转为Msun*h/pc^2,即把面积单位转为pc^2
             deltaSigma[t] = delta_segma*10**(-12)
             Sigma[t] = 2*rs*rho_0*(1-2*f1/np.sqrt(f0**2-1))/(f0**2-1)
    return(Rpc,Sigma,deltaSigma,rs)
#calcu_sigmaz(Rpc=True,m_=True,x=True,z=True)
#单独调用模块时需要对此句修改，引入直接的Rp和质量
def fig_f(inm,Rp,r,rs,r_200,rho_R,sigma,deltasigma):
#def fig_f(ff):
#inm,Rs,r,rs,r_200,rho_R,sigma,deltasigma=semula_tion(m_,x)
    lm = len(inm)
    plt.figure()
    plt.subplot(121)
    # delta1 = np.zeros(lm,dtype=np.float)
    # delta2 = np.zeros(lm,dtype=np.float)
    for k in range(0,lm):
        plt.loglog(r[k,:],rho_R[k])
        '''
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
        '''
    plt.xlabel(r'$r-Mpc/h$')
    plt.ylabel(r'$\rho(R)-M_\odot/h$')
    plt.grid()
    print('mh_inte=',inm)
    plt.subplot(122)
    #下面引入具体的deltasigma数据的观测值,和所给质量比较
    test_read_m16_ds()
    delta1 = np.zeros(lm,dtype=np.float)
    delta2 = np.zeros(lm,dtype=np.float)
    for k in range(0,lm):
        plt.loglog(Rp[k,:],deltasigma[k,:])
        delta1[k] = r_200[k]
        delta2[k] = rs[k]
        plt.axvline(delta1[k],ls='--',linewidth=0.5,color='red')
        plt.axvline(delta2[k],ls='--',linewidth=0.5,color='blue')
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')
    plt.tight_layout()
    plt.show()
    out=[inm,rho_R,deltasigma,sigma]
    return(out)
#fig_f(inm=True,Rs=True,r=True,rs=True,r_200=True,rho_R=True,sigma=True,deltasigma=True)
#fig_f(ff=True)
##定义一个五个质量最佳预测和观测的对比图象
def fig_ff(Rpc,ds_sim,lmw,r):
    test_read_m16_ds()
    for k in range(0,lmw):
        pa = r[k,1]
        plt.plot(Rpc[:],ds_sim[k,pa,:],'-*')
        nn = ds_sim[k,pa,:].size#求样本大小
        x_mean = np.mean(ds_sim[k,pa,:])#求算术平均值
        x_std = np.std(ds_sim[k,pa,:])#求标准偏差
        x_se = x_std/np.sqrt(nn)#求标准误差
        dof = nn-1#自由度计算
        alpha = 1.0-0.95
        conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
        np.sqrt(1.+1./nn)#设置置信区间
        plt.errorbar(Rpc[:],ds_sim[k,pa,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_ff(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_ff(f=True)
##定义一个质量最佳预测和观测的对比图象
def fig_fff1(Rpc,ds_sim,k):
    test_read_m16_ds()
    plt.plot(Rpc[:],ds_sim[k,:],'-*')
    nn = ds_sim[k,:].size#求样本大小
    x_mean = np.mean(ds_sim[k,:])#求算术平均值
    x_std = np.std(ds_sim[k,:])#求标准偏差
    x_se = x_std/np.sqrt(nn)#求标准误差
    dof = nn-1#自由度计算
    alpha = 1.0-0.95
    conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
    np.sqrt(1.+1./nn)#设置置信区间
    plt.errorbar(Rpc[:],ds_sim[k,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_fff1(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_fff1(f=True)
def fig_fff2(Rpc,ds_sim,k,t):
    test_read_m16_ds()
    plt.plot(Rpc[:],ds_sim[k,t,:],'-*')
    nn = ds_sim[k,t,:].size#求样本大小
    x_mean = np.mean(ds_sim[k,t,:])#求算术平均值
    x_std = np.std(ds_sim[k,t,:])#求标准偏差
    x_se = x_std/np.sqrt(nn)#求标准误差
    dof = nn-1#自由度计算
    alpha = 1.0-0.95
    conf_region = st.ppf(1-alpha/2.,dof)*x_std*\
    np.sqrt(1.+1./nn)#设置置信区间
    plt.errorbar(Rpc[:],ds_sim[k,t,:],yerr=x_std,fmt='-',linewidth=0.5)
    plt.xlabel(r'$R-Mpc/h$')
    plt.ylabel(r'$\Delta\Sigma-hM_\odot/{pc^2}$')   
    return()
#fig_fff2(Rpc=True,ds_sim=True,lmw=True,r=True)
#fig_fff2(f=True)
##
def fit_datar(y):
    h = 0.673
    #先找出观测值对应的rp
    rsa,dssa,ds_errsa = test_read_m16_ds()
    # print('ds=',dssa)
    # print('ds_err=',ds_errsa)
    #print('rp=',rsa)  
    a = np.shape(dssa)
    print('size a=',a)
    ##注意到观测数值选取的是物理坐标，需要靠里尺度银子修正，修正如下
    z_r = 0.230
    #z_r = 0.105
    a_r = 1/(1+z_r)
    ##此时对于预测的R也要做修正
    rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2])
    #rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=5]) 
    #rp = rsa[0,:]
    #该句直接对数组筛选，计算2Mpc以内（保证NFW模型适用的情况下）信号
    print(rp)
    b = len(rp)
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    #下面做两组数据的方差比较,注意不能对观测数据进行插值
    #比较方法，找出观测的sigma数值对应的Rp,再根据模型计算此时的模型数值Sigma（该步以完成）    
    ds_simr = np.zeros((lmw,lml,b),dtype=np.float)
    rr = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            m_ = m[k,t]
            x = m_dex[k]
            #Rpc = rp
            #Rpc,Sigma,deltaSigma = calcu_sigma(Rpc,m_,x)    
            #ds_simr[k,t,:] = deltaSigma
            #对比文献，加入尺度因子修正如下,把物理的距离转为共动的距离
            #计算模型在对应的投射距离上的预测信号取值
            Rpc = rp
            Rpc,Sigma,deltaSigma,rs = calcu_sigmaz(Rpc,m_,x,z_r)
            #对预测信号做相应的变化，把共动的转为物理的
            ds_simr[k,t,:] = deltaSigma
            rr[k,t] = rs
    yy = np.shape(ds_simr)
    #print('rr=',rr)#输出查看ds_sim的维度，即模型预测下的透镜信号    
    #比较观测的sigma和预测的Sigma,比较结果用fitr和fitb记录比较结果
    fitr = np.zeros((lmw,2),dtype=np.int)
    delta_r = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            d_r = 0
            for n in range(0,yy[2]):
                d_r = d_r+((ds_simr[k,t,n]-dssa[0,n])/ds_errsa[0,n])**2
            delta_r[k,t] = d_r
    #print(delta_r)
    #下面这段求每个质量级对应的方差最小值
    for k in range(0,lmw):
        #先把np的定义数组转为列表，再寻找索引号，tolist()表示内置函数调用
        aa = delta_r[k,:].tolist()
        xa = aa.index(min(delta_r[k,:]))
        fitr[k,0] = k
        fitr[k,1] = xa
    print(fitr)
    #下面做图比较几个最佳预测值与观测的对比情况
    Rpc = rp
    fig_ff(Rpc,ds_simr,lmw,fitr)
    plt.title('R-galaxy')
    plt.show() 
    #比较提取出的五个最小值的最小值为最后结果
    deltar = np.zeros(lmw,dtype=np.float)
    for k in range(0,lmw):
        deltar[k] = delta_r[fitr[k,0],fitr[k,1]]
    #下面用bestr和bestb记录最后比较结果
    bestfr = np.zeros(2,dtype=np.float)
    aa = deltar.tolist()
    xa = aa.index(min(deltar))
    bestfr = fitr[xa,:]  
    bestr = m[bestfr[0],bestfr[1]]
    best_mr = bestr
    dexr1 = m_dex[bestfr[0]]
    #下面两句是为与文献对比做的修正
    dexr = m_dex[bestfr[0]]+np.log10(bestr*h) 
    #作图对比最佳情况
    plt.figure()
    k = bestfr[0]
    t = bestfr[1]
    fig_fff2(rp,ds_simr,k,t)
    plt.title('Red')
    #plt.savefig('fit_r.png',dpi=600)
    plt.show()
    print('x^2=',min(deltar))
    print('mr=',best_mr) #此句质量计算正常
    print('positionr=',bestfr)
    print('dexr=',dexr1)
    print('co_dexr=',dexr)  
    print('corr_mr=',10**dexr)
    return rp,best_mr,dexr1,a,b,yy,dexr,a_r,delta_r 
#fit_datar(y=True)
def fit_datab(y):
    h = 0.673
    #先找出观测值对应的rp
    rsa,dssa,ds_errsa = test_read_m16_ds()
    # print('ds=',dssa)
    # print('ds_err=',ds_errsa)
    #print('rp=',rsa)  
    a = np.shape(dssa)
    print('size a=',a)
    ##注意到观测数值选取的是物理坐标，需要靠里尺度银子修正，修正如下
    z_b = 0.246
    #z_b = 0.124
    a_b = 1/(1+z_b)
    ##此时对于预测的R也要做修正
    rp = np.array([rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2]) 
   # rp = [rsa[0,k] for k in range(0,len(rsa[0,:])) if rsa[0,k]*h<=2]
   # rp = rsa[0,:]
    #该句直接对数组筛选，计算2Mpc以内（保证NFW模型适用的情况下）信号
    print(rp)
    b = len(rp)
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    #下面做两组数据的方差比较,注意不能对观测数据进行插值
    #比较方法，找出观测的sigma数值对应的Rp,再根据模型计算此时的模型数值Sigma（该步以完成）    
    ds_simb = np.zeros((lmw,lml,b),dtype=np.float)
    rb = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            m_ = m[k,t]
            x = m_dex[k]
            #Rpc = rp
            #Rpc,Sigma,deltaSigma = calcu_sigma(Rpc,m_,x,z_b)
            #ds_simb[k,t,:] = deltaSigma
            #对比文献，把物理距离转为共动距离
            Rpc = rp
            #计算模型在对应的投射距离上的预测信号取值
            Rpc,Sigma,deltaSigma,rs = calcu_sigmaz(Rpc,m_,x,z_b)
            #预测信号也做修正
            ds_simb[k,t,:] = deltaSigma
            rb[k,t] = rs
    yy = np.shape(ds_simb)
    #print('rb=',rb)#输出查看ds_sim的维度，即模型预测下的透镜信号    
    #比较观测的sigma和预测的Sigma,比较结果用fitr和fitb记录比较结果
    fitb = np.zeros((lmw,2),dtype=np.int)
    delta_b = np.zeros((lmw,lml),dtype=np.float)
    for k in range(0,lmw):
        for t in range(0,lml):
            d_b = 0
            for n in range(0,yy[2]):
                d_b = d_b+((ds_simb[k,t,n]-dssa[1,n])/ds_errsa[1,n])**2
            delta_b[k,t] = d_b
    #print(delta_b)
    #下面这段求每个质量级对应的方差最小值
    for k in range(0,lmw):
        #先把np的定义数组转为列表，再寻找索引号，tolist()表示内置函数调用
        bb = delta_b[k,:].tolist()
        xb = bb.index(min(delta_b[k,:]))
        fitb[k,0] = k
        fitb[k,1] = xb
    print(fitb)
    #下面做图比较几个最佳预测值与观测的对比情况
    Rpc = rp
    fig_ff(Rpc,ds_simb,lmw,fitb)
    plt.title('B-galaxy')
    plt.show() 
    #比较提取出的五个最小值的最小值为最后结果
    deltab = np.zeros(lmw,dtype=np.float)
    for k in range(0,lmw):
        deltab[k] = delta_b[fitb[k,0],fitb[k,1]]
    #下面用bestr和bestb记录最后比较结果
    bestfb = np.zeros(2,dtype=np.float)
    bb = deltab.tolist()
    xb = bb.index(min(deltab))
    bestfb = fitb[xb,:]    
    bestb = m[bestfb[0],bestfb[1]]
    best_mb = bestb
    dexb1 = m_dex[bestfb[0]]
    #下面是与文献对比做的修正
    dexb = m_dex[bestfb[0]]+np.log10(bestb*h)    
    #作图对比最佳情况
    plt.figure()
    k = bestfb[0]
    t = bestfb[1]
    fig_fff2(rp,ds_simb,k,t)
    plt.title('Blue')
    #plt.savefig('fit_b.png',dpi=600)
    plt.show()
    print('x^2=',min(deltab))    
    print('mb=',best_mb) #此句质量计算正常
    print('positionb=',bestfb)
    print('dexb=',dexb1)
    print('co_dexb=',dexb) 
    print('corr_mb=',10**dexb)
    return rp,best_mb,dexb1,a,b,yy,dexb,a_b,delta_b 
#fit_datab(y=True)
#下面做图比较最佳预测值与观测的对比情况
def fig_(T):
    m,m_dex,lmw,lml = dolad_data(m=True,hz=True)
    rp,best_mr,dexr1,a,b,yy,dexr,a_r,delta_r = fit_datar(y=True)
    rp,best_mb,dexb1,a,b,yy,dexb,a_b,delta_b = fit_datab(y=True)
    plt.subplot(121)
    plt.plot(m_dex,np.log10(delta_r))
    plt.title('R-galaxy')
    plt.xlabel(r'$log(\frac{M_h}{M_\odot})$')
    plt.ylabel(r'$log(\chi^2)$')
    plt.grid()
    plt.subplot(122)
    plt.plot(m_dex,np.log10(delta_b))
    plt.xlabel(r'$log(\frac{M_h}{M_\odot})$')
    plt.ylabel(r'$log(\chi^2)$')
    plt.title('B-galaxy')
    plt.grid()
    plt.tight_layout()
    plt.show()  
    print('mr=',dexr)
    print('mb=',dexb)
    return dexr,dexb
fig_(T=True)
'''
if __name__ == "__main__":
    dolad_data(m=True,hz=True)
    pass
'''