<jupyter_start><jupyter_text>Getting the catalog, and combining it with the subimg_data cat to get the NMS per cluster.<jupyter_code>Cat=Table.read('M33_catalog_with_Quick_BestFit_CMD_Estimates.fits')

subimg_data=Table.read('phatm33_subimg-meta.fits')
subimg_data.rename_column('SUBIMG', 'IMG')
subimg_data.rename_column('FILENAME', 'SUBIMG_old')
subimg_data.rename_column('RA', 'IMG_RA')
subimg_data.rename_column('DEC', 'IMG_DEC')

si_lst=list(subimg_data['SUBIMG_old'])
for i in range(len(si_lst)):
    si_lst[i]= si_lst[i].rstrip()
    
new_col=Column(si_lst, name='SUBIMG')
subimg_data.remove_column('SUBIMG_old')
subimg_data.add_column(new_col, index=4)

Cat_w_nms=join(Cat, subimg_data)

Cat_w_Est=Cat_w_nms[np.where(Cat_w_nms['Quick_BestFit_LogAge'] > 0)]

#Trim to exclude clusters younger than 7.0 and older than 8.5

Trimmed=Cat_w_Est[np.where((Cat_w_Est['Quick_BestFit_LogAge'] > 7.0 ) & (Cat_w_Est['Quick_BestFit_LogAge'] < 8.5))]
#Only using the sample in my age range

masses=np.array(Trimmed['Quick_BestFit_LogMass'])
ages=np.array(Trimmed['Quick_BestFit_LogAge'])
nmses=np.array(np.log10(Trimmed['NMS']))
#Setting the Mlim to a fixed value which is the median 50% compleness for the full sample
med_mlim=2.9677
med_mlim

fixed_mlim_masses=masses[np.where(masses>med_mlim)]
use_med_mlim=10**med_mlim

def pobs(M, mlim):
    k=6.3665
    
    y=(1.+ exp(-k*(M-mlim)))**(-1)
    
    return y

def lnobs_like(M, mlim):
    k=6.3665
    
    return -np.log(1.+ exp(-k*(M-mlim)))
        
def Shecter_Z(M, mlim, alpha, M_c):
    x = M/M_c
    k=6.3665
    pobs= 1./(1.+ exp((-k)*(np.log10(M)-mlim)))
    
    return (x**alpha) * exp(-x) * pobs
        
    
def lnlike(theta, M, mlim):
    alpha, M_c = theta
    lin_M_c= 10.**M_c
    
    lin_M= 10**M
    x= lin_M/lin_M_c
    
    
    ln_pobs=lnobs_like(M, np.log10(mlim))
        
    norm, err = quad(Shecter_Z, mlim, 1.e7, args=(np.log10(mlim), alpha, lin_M_c))
    
    
    lnlike = np.sum((-x) + alpha*np.log(x) + ln_pobs - np.log(norm))    
    return lnlike
    
def lnprior(theta):
    alpha, M_c = theta
    if -3 <= alpha <= -1 and 3 <= M_c <= 8:
        return 0.0
    return -np.inf

def lnprob(theta, M, mlim):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, M, mlim)

nll = lambda *args: -lnprob(*args)
starting_point=np.array([-2., 4.])

fixed_mlim_result=opt.minimize(nll, x0=starting_point, args=(fixed_mlim_masses, use_med_mlim))
fixed_mlim_result['x']<jupyter_output><empty_output>