# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import scipy.spatial as sp
#import scipy.integrate as integrate

g, bp , rp, age, bn, feh, m, ms = np.load('10MSTAR_POPULATION_GBRABFM.npy').T

def show_pop():
    fig, ax = plt.subplots()
    im = ax.scatter((bp-rp),g,c=age,s=1, alpha=1)
    fig.colorbar(im, ax=ax)
    ax.set_ylabel('G')
    ax.set_xlabel('Gaia B-R')
    ax.set_ylim(15.5,-10)
    ax.figure.set_size_inches(4, 7)
    fig.savefig('KerSynthfc.png',dpi=200)
    fig.savefig('KerSynthfc.pdf')
    fig.show1
    
zipnowork=np.vstack(((bp-rp),g)).T
zipnowork.shape
tree=sp.KDTree(zipnowork)

def treeGetAges(col,mag,radius,data=tree):

    inds=data.query_ball_point((col,mag),radius)
    ages=age[inds]
    colres=bp[inds]-rp[inds]
    gres=g[inds]

    return ages, colres, gres

lgage=np.arange(5,11.4, 0.1)

def make_hists(c,g,n=100,r=0.01):
    ages,colres,gres=treeGetAges(c, g, r)
    if len(ages)>(n*2):
        r=0.75*r
        ages,colres,gres=treeGetAges(c, g, r)
    if len(ages)<n:
        r=1.2*r
        ages,colres,gres=treeGetAges(c, g, r)
    
    dists=np.sqrt((c-colres)**2+(g-gres)**2)
    ages=ages[np.argsort(dists)]
    ages=ages[:n]
    ags=np.round((ages-5)/0.1).astype(int)
    age_hists=np.zeros_like(lgage)
    for i in ags:
        age_hists[i]+=1
    print(len(ages))
    return age_hists    
    
gaus  = np.exp(-(lgage-np.median(lgage))**2/ 0.1**2 /2)
gaus /= np.sum(gaus)
gausft = np.fft.rfft(np.fft.fftshift(gaus))

def g_kernal_den(c,g):
    age_h=make_hists(c, g)
    age_pdf=np.fft.irfft(np.fft.rfft(age_h)*gausft)
   # grated=integrate.simps(age_pdf)
    normed=age_pdf*10**lgage/sum(age_pdf*10**lgage)
    plt.plot(lgage,normed)
    return normed

def get_probage(c,g,A):
     if not(5<A<11.4):
        print('Err; log(age) out of 5 to 11.4 interval')
     pdf=g_kernal_den(c, g)
     iclosest=(np.abs(lgage-A)).argmin()
     return pdf[iclosest]
   
    
#if __name__=='__main__':

Create honpy.py