@@ -0,0 +1,1960 @@
+# -*- coding: utf-8 -*-%0Aimport matplotlib.pyplot as plt%0Aimport numpy as np%0Aimport scipy.spatial as sp%0A#import scipy.integrate as integrate%0A%0Ag, bp , rp, age, bn, feh, m, ms = np.load('10MSTAR_POPULATION_GBRABFM.npy').T%0A%0Adef show_pop():%0A    fig, ax = plt.subplots()%0A    im = ax.scatter((bp-rp),g,c=age,s=1, alpha=1)%0A    fig.colorbar(im, ax=ax)%0A    ax.set_ylabel('G')%0A    ax.set_xlabel('Gaia B-R')%0A    ax.set_ylim(15.5,-10)%0A    ax.figure.set_size_inches(4, 7)%0A    fig.savefig('KerSynthfc.png',dpi=200)%0A    fig.savefig('KerSynthfc.pdf')%0A    fig.show1%0A    %0Azipnowork=np.vstack(((bp-rp),g)).T%0Azipnowork.shape%0Atree=sp.KDTree(zipnowork)%0A%0Adef treeGetAges(col,mag,radius,data=tree):%0A%0A    inds=data.query_ball_point((col,mag),radius)%0A    ages=age%5Binds%5D%0A    colres=bp%5Binds%5D-rp%5Binds%5D%0A    gres=g%5Binds%5D%0A%0A    return ages, colres, gres%0A%0Algage=np.arange(5,11.4, 0.1)%0A%0Adef make_hists(c,g,n=100,r=0.01):%0A    ages,colres,gres=treeGetAges(c, g, r)%0A    if len(ages)%3E(n*2):%0A        r=0.75*r%0A        ages,colres,gres=treeGetAges(c, g, r)%0A    if len(ages)%3Cn:%0A        r=1.2*r%0A        ages,colres,gres=treeGetAges(c, g, r)%0A    %0A    dists=np.sqrt((c-colres)**2+(g-gres)**2)%0A    ages=ages%5Bnp.argsort(dists)%5D%0A    ages=ages%5B:n%5D%0A    ags=np.round((ages-5)/0.1).astype(int)%0A    age_hists=np.zeros_like(lgage)%0A    for i in ags:%0A        age_hists%5Bi%5D+=1%0A    print(len(ages))%0A    return age_hists    %0A    %0Agaus  = np.exp(-(lgage-np.median(lgage))**2/ 0.1**2 /2)%0Agaus /= np.sum(gaus)%0Agausft = np.fft.rfft(np.fft.fftshift(gaus))%0A%0Adef g_kernal_den(c,g):%0A    age_h=make_hists(c, g)%0A    age_pdf=np.fft.irfft(np.fft.rfft(age_h)*gausft)%0A   # grated=integrate.simps(age_pdf)%0A    normed=age_pdf*10**lgage/sum(age_pdf*10**lgage)%0A    plt.plot(lgage,normed)%0A    return normed%0A%0Adef get_probage(c,g,A):%0A     if not(5%3CA%3C11.4):%0A        print('Err; log(age) out of 5 to 11.4 interval')%0A     pdf=g_kernal_den(c, g)%0A     iclosest=(np.abs(lgage-A)).argmin()%0A     return pdf%5Biclosest%5D%0A   %0A    %0A#if __name__=='__main__':%0A

Create honpy.py