#########################
#########################
# Need to account for limit in input period
#########################
#########################

# Baseline M67 long script -- NO crowding
# New script copied from quest - want to take p and ecc from each population (all, obs, rec) and put them into separate file
# Doing this so we don't have to run analyse each time
# Can write separate script for p-ecc plots
# Quest paths in this version of script

import pandas as pd
import numpy as np
import os
from astropy.coordinates import SkyCoord
from astropy import units, constants
from astropy.modeling import models, fitting
import scipy.stats
from scipy.integrate import quad

#for Quest
import matplotlib
matplotlib.use('Agg')

doIndividualPlots = True

from matplotlib import pyplot as plt
def file_len(fname):
	i = 0
	with open(fname) as f:
		for i, l in enumerate(f):
			pass
	return i + 1

def getPhs(sigma, m1=1*units.solMass, m2=1*units.solMass, m3=0.5*units.solMass):
	Phs = np.pi*constants.G/np.sqrt(2.)*(m1*m2/m3)**(3./2.)*(m1 + m2)**(-0.5)*sigma**(-3.)
	return Phs.decompose().to(units.day)

#similar to field, but limiting by the hard-soft boundary
def fitRagfb():
	x = [0.05, 0.1, 1, 8, 15]  #estimates of midpoints in bins, and using this: https://sites.uni.edu/morgans/astro/course/Notes/section2/spectralmasses.html
	y = [0.20, 0.35, 0.50, 0.70, 0.75]
	init = models.PowerLaw1D(amplitude=0.5, x_0=1, alpha=-1.)
	fitter = fitting.LevMarLSQFitter()
	fit = fitter(init, x, y)

	return fit

def RagNormal(x, cdf = False):
	mean = 5.03
	std = 2.28
	if (cdf):
		return scipy.stats.norm.cdf(x,mean,std)

	return scipy.stats.norm.pdf(x,mean,std)


def saveHist(histAll, histObs, histRec, bin_edges, xtitle, fname, filters = ['u_', 'g_', 'r_', 'i_', 'z_', 'y_','all']):
	c1 = '#5687A6' #Dali Blue (Andrew's AAS Poster)
	c2 = '#A62B1F' #Dai Red 
	c3 = '#BF8A26' #Dali Beige
	fig,ax1 = plt.subplots(figsize=(8,6), sharex=True)#can change to include cdf with ax1, ax2

	histAll = np.insert(histAll,0,0)
	histObs = np.insert(histObs,0,0)
	for f in filters:
		histRec[f] = np.insert(histRec[f],0,0)

	#PDF
	ax1.step(bin_edges, histAll/np.sum(histAll), color=c1)
	ax1.step(bin_edges, histObs/np.sum(histObs), color=c2)
	for f in filters:
		lw = 1
		if (f == 'all'):
			lw = 0.5
	ax1.step(bin_edges, histRec[f]/np.sum(histRec[f]), color=c3, linewidth=lw)
	ax1.set_ylabel('PDF')
	ax1.set_yscale('log')
	ax1.set_title('Globular Clusters - Baseline', fontsize = 16)
	ax1.set_xlabel(xtitle)
	#CDF
	#cdfAll = []
	#cdfObs = []
	#cdfRec = dict()
	#for f in filters:
#		cdfRec[f] = []

#	for i in range(len(histAll)):
#		cdfAll.append(np.sum(histAll[:i])/np.sum(histAll))
#	for i in range(len(histObs)):
#		cdfObs.append(np.sum(histObs[:i])/np.sum(histObs))
#	for f in filters:
	#	for i in range(len(histRec[f])):
	#		cdfRec[f].append(np.sum(histRec[f][:i])/np.sum(histRec[f]))
	#ax2.step(bin_edges, cdfAll, color=c1)
	#ax2.step(bin_edges, cdfObs, color=c2)
	#for f in filters:
	#	lw = 1
	#	if (f == 'all'):
	#		lw = 0.5
	#	ax2.step(bin_edges, cdfRec[f], color=c3, linewidth=lw)
	#ax2.set_ylabel('CDF')

	#ax2.set_xlabel(xtitle)
	fig.subplots_adjust(hspace=0)
	fig.savefig('./plots/'  + fname+'.pdf',format='pdf', bbox_inches = 'tight')

	#write to a text file
	with open('./eblsst_files/' + fname+'.csv','w') as fl:
		outline = 'binEdges,histAll,histObs'
		for f in filters:
			outline += ','+f+'histRec'
		outline += '\n'
		fl.write(outline)
		for i in range(len(bin_edges)):
			outline = str(bin_edges[i])+','+str(histAll[i])+','+str(histObs[i])
			for f in filters:
				outline += ','+str(histRec[f][i])
			outline += '\n'
			fl.write(outline)

if __name__ == "__main__":

	filters = ['u_', 'g_', 'r_', 'i_', 'z_', 'y_', 'all']

	#get the Raghavan binary fraction fit
	fbFit= fitRagfb()
	print(fbFit)
		
	#to normalize
	intAll, err = quad(RagNormal, -20, 20)
	intCut, err = quad(RagNormal, -20, np.log10(365*10.))
	intNorm = intCut/intAll

	#cutoff in percent error for "recovered"
	Pcut = 0.1

	#assumed mean stellar mass
	mMean = 0.5

	#minimum number of lines to consider in file
	Nlim = 3

	if (doIndividualPlots):
		fmass, axmass = plt.subplots()
		fqrat, axqrat = plt.subplots()
		fecc, axecc = plt.subplots()
		flper, axlper = plt.subplots()
		fdist, axdist = plt.subplots()
		fmag, axmag = plt.subplots()
		frad, axrad = plt.subplots()

	#bins for all the histograms
	Nbins = 25
	mbins = np.arange(0,10, 0.1, dtype='float')
	qbins = np.arange(0,1, 0.1, dtype='float')
	ebins = np.arange(0, 1.05, 0.05, dtype='float')
	lpbins = np.arange(-2, 10, 0.5, dtype='float')
	dbins = np.arange(0, 40, 1, dtype='float')
	magbins = np.arange(11, 25, 1, dtype='float')
	rbins = np.arange(0, 100, 0.2, dtype='float')

	#blanks for the histograms
	#All
	m1hAll = np.zeros_like(mbins)[1:]
	qhAll = np.zeros_like(qbins)[1:]
	ehAll = np.zeros_like(ebins)[1:]
	lphAll = np.zeros_like(lpbins)[1:]
	dhAll = np.zeros_like(dbins)[1:]
	maghAll = np.zeros_like(magbins)[1:]
	rhAll = np.zeros_like(rbins)[1:]
	#Observable
	m1hObs = np.zeros_like(mbins)[1:]
	qhObs = np.zeros_like(qbins)[1:]
	ehObs = np.zeros_like(ebins)[1:]
	lphObs = np.zeros_like(lpbins)[1:]
	dhObs = np.zeros_like(dbins)[1:]
	maghObs = np.zeros_like(magbins)[1:]
	rhObs = np.zeros_like(rbins)[1:]
	#Recovered
	m1hRec = dict()
	qhRec = dict()
	ehRec = dict()
	lphRec = dict()
	dhRec = dict()
	maghRec = dict()
	rhRec = dict()
	for f in filters:
		m1hRec[f] = np.zeros_like(mbins)[1:]
		qhRec[f] = np.zeros_like(qbins)[1:]
		ehRec[f] = np.zeros_like(ebins)[1:]
		lphRec[f] = np.zeros_like(lpbins)[1:]
		dhRec[f] = np.zeros_like(dbins)[1:]
		maghRec[f] = np.zeros_like(magbins)[1:]
		rhRec[f] = np.zeros_like(rbins)[1:]

	RA = []
	Dec = []
	recFrac = []
	recN = []
	rawN = []
	obsN = []
	fileN = []
	fileObsN = []
	fileRecN = []

	allNPrsa = []
	obsNPrsa = []
	recNPrsa = []

	# Lists for period and eccentricity for Andrew's circularization plots
	eccAll = []
	eccObs = []
	eccRec = []

	pAll = []
	pObs = []
	pRec = []
	# Using prsa dataframes for these lists because of period cutoff at 1000 days

	# Dataframes to write to files later; 3 files for each sub-population - append everything to these
	peccAll = pd.DataFrame(columns = ['e', 'p'])
	peccObs = pd.DataFrame(columns = ['e', 'p'])
	peccRec = pd.DataFrame(columns = ['e', 'p'])

	#Read in all the data and make the histograms
	d = "./input_files/"
	files = os.listdir(d)
	IDs = []
	for i, f in enumerate(files):
		print(round(i/len(files),4), f)
		fl = file_len(d+f)
		if (fl >= 4):
			#read in the header
			header = pd.read_csv(d+f, nrows=1)
	######################
	#NEED TO ACCOUNT FOR THE BINARY FRACTION when combining histograms
	#####################
			Nmult = header['clusterMass'][0]/mMean
			#Nmult = 1.

			RA.append(header['OpSimRA'])
			Dec.append(header['OpSimDec'])

			#read in rest of the file
			data = pd.read_csv(d+f, header = 2).fillna(-999)
			rF = 0.
			rN = 0.
			Nrec = 0.
			Nobs = 0.
			raN = 0.
			obN = 0.
			fiN = 0.
			fioN = 0.
			firN = 0.
			NallPrsa = 0.
			NobsPrsa = 0.
			NrecPrsa = 0.
			Nall = len(data.index)/intNorm ###is this correct? (and the only place I need to normalize?)
			prsa = data.loc[(data['appMagMean_r'] <= 19.5) & (data['appMagMean_r'] > 15.8) & (data['p'] < 1000) & (data['p'] > 0.5)]

			# Appending for Andrew
			eccAll.append(prsa['e'].values)
			pAll.append(prsa['p'].values)

			NallPrsa = len(prsa.index)
			if (Nall >= Nlim):
				#create histograms
				#All
				m1hAll0, m1b = np.histogram(data["m1"], bins=mbins)
				qhAll0, qb = np.histogram(data["m2"]/data["m1"], bins=qbins)
				ehAll0, eb = np.histogram(data["e"], bins=ebins)
				lphAll0, lpb = np.histogram(np.ma.log10(data["p"].values).filled(-999), bins=lpbins)
				dhAll0, db = np.histogram(data["d"], bins=dbins)
				maghAll0, magb = np.histogram(data["appMagMean_r"], bins=magbins)
				rhAll0, rb = np.histogram(data["r2"]/data["r1"], bins=rbins)

				if (doIndividualPlots):
					axmass.step(m1b[0:-1], m1hAll0/np.sum(m1hAll0), color='black', alpha=0.1)
					axqrat.step(qb[0:-1], qhAll0/np.sum(qhAll0), color='black', alpha=0.1)
					axecc.step(eb[0:-1], ehAll0/np.sum(ehAll0), color='black', alpha=0.1)
					axlper.step(lpb[0:-1], lphAll0/np.sum(lphAll0), color='black', alpha=0.1)
					axdist.step(db[0:-1], dhAll0/np.sum(dhAll0), color='black', alpha=0.1)
					axmag.step(magb[0:-1], maghAll0/np.sum(maghAll0), color='black', alpha=0.1)
					axrad.step(rb[0:-1], rhAll0/np.sum(rhAll0), color='black', alpha=0.1)

				#account for the binary fraction, as a function of mass
				dm1 = np.diff(m1b)
				m1val = m1b[:-1] + dm1/2.
				fb = np.sum(m1hAll0/len(data.index)*fbFit(m1val))
				#account for the hard-soft boundary
				Phs = getPhs(header['clusterVdisp'].iloc[0]*units.km/units.s).to(units.day).value
				fb *= RagNormal(np.log10(Phs), cdf = True)
				print("fb, Phs = ", fb, Phs)
				Nmult *= fb

							
				m1hAll += m1hAll0/Nall*Nmult
				qhAll += qhAll0/Nall*Nmult
				ehAll += ehAll0/Nall*Nmult
				lphAll += lphAll0/Nall*Nmult
				dhAll += dhAll0/Nall*Nmult
				maghAll += maghAll0/Nall*Nmult
				rhAll += rhAll0/Nall*Nmult

				#Obs
				obs = data.loc[data['LSM_PERIOD'] != -999]
				Nobs = len(obs.index)
				prsaObs = data.loc[(data['appMagMean_r'] <= 19.5) & (data['appMagMean_r'] > 15.8) & (data['p'] < 1000) & (data['p'] >0.5) & (data['LSM_PERIOD'] != -999)]
				NobsPrsa = len(prsaObs.index)

				# Appending for Andrew's files
				eccObs.append(prsaObs['e'].values)
				pObs.append(prsaObs['p'].values)

				if (Nobs >= Nlim):
					m1hObs0, m1b = np.histogram(obs["m1"], bins=mbins)
					qhObs0, qb = np.histogram(obs["m2"]/obs["m1"], bins=qbins)
					ehObs0, eb = np.histogram(obs["e"], bins=ebins)
					lphObs0, lpb = np.histogram(np.ma.log10(obs["p"].values).filled(-999), bins=lpbins)
					dhObs0, db = np.histogram(obs["d"], bins=dbins)
					maghObs0, magb = np.histogram(obs["appMagMean_r"], bins=magbins)
					rhObs0, rb = np.histogram(obs["r2"]/obs["r1"], bins=rbins)
					m1hObs += m1hObs0/Nall*Nmult
					qhObs += qhObs0/Nall*Nmult
					ehObs += ehObs0/Nall*Nmult
					lphObs += lphObs0/Nall*Nmult
					dhObs += dhObs0/Nall*Nmult
					maghObs += maghObs0/Nall*Nmult
					rhObs += rhObs0/Nall*Nmult

					#Rec
					recCombined = pd.DataFrame()
					prsaRecCombined = pd.DataFrame()
					for filt in filters:
						key = filt+'LSS_PERIOD'
						if (filt == 'all'):
							key = 'LSM_PERIOD'
						fullP = abs(data[key] - data['p'])/data['p']
						halfP = abs(data[key] - 0.5*data['p'])/(0.5*data['p'])
						twiceP = abs(data[key] - 2.*data['p'])/(2.*data['p'])
						rec = data.loc[(data[key] != -999) & ( (fullP < Pcut) | (halfP < Pcut) | (twiceP < Pcut))]
						prsaRec = data.loc[(data['appMagMean_r'] <= 19.5) & (data['appMagMean_r'] >15.8) & (data['p'] < 1000) & (data['p'] >0.5) & (data['LSM_PERIOD'] != -999) & ( (fullP < Pcut) | (halfP < Pcut) | (twiceP < Pcut))]
						Nrec = len(rec.index)

						#I'd like to account for all filters here to have more accurate numbers
						recCombined = recCombined.append(rec)
						prsaRecCombined = prsaRecCombined.append(prsaRec)

						# Going to use prsaRecCombined for ecc-p plots to account for all filters
						eccRec.append(prsaRec['e'].values)
						pRec.append(prsaRec['p'].values)


						if (filt == 'all'):
							recCombined.drop_duplicates(inplace=True)
							prsaRecCombined.drop_duplicates(inplace=True)

						if (Nrec >= Nlim):
							m1hRec0, m1b = np.histogram(rec["m1"], bins=mbins)
							qhRec0, qb = np.histogram(rec["m2"]/rec["m1"], bins=qbins)
							ehRec0, eb = np.histogram(rec["e"], bins=ebins)
							lphRec0, lpb = np.histogram(np.ma.log10(rec["p"].values).filled(-999), bins=lpbins)
							dhRec0, db = np.histogram(rec["d"], bins=dbins)
							maghRec0, magb = np.histogram(rec["appMagMean_r"], bins=magbins)
							rhRec0, rb = np.histogram(rec["r2"]/rec["r1"], bins=rbins)
							m1hRec[filt] += m1hRec0/Nall*Nmult
							qhRec[filt] += qhRec0/Nall*Nmult
							ehRec[filt] += ehRec0/Nall*Nmult
							lphRec[filt] += lphRec0/Nall*Nmult
							dhRec[filt] += dhRec0/Nall*Nmult
							maghRec[filt] += maghRec0/Nall*Nmult
							rhRec[filt] += rhRec0/Nall*Nmult

							#for the mollweide
							if (filt == 'all'):
								Nrec = len(recCombined.index)
								rF = Nrec/Nall
								rN = Nrec/Nall*Nmult
								raN = Nmult
								obN = Nobs/Nall*Nmult
								fiN = Nall
								fioN = Nobs
								firN = Nrec

								NrecPrsa = len(prsaRecCombined.index)
								NrecPrsa = NrecPrsa/Nall*Nmult
								NobsPrsa = NobsPrsa/Nall*Nmult
								NallPrsa = NallPrsa/Nall*Nmult		



			recFrac.append(rF)
			recN.append(rN)
			rawN.append(raN)
			obsN.append(obN)
			fileN.append(fiN)
			fileObsN.append(fioN)
			fileRecN.append(firN)
			allNPrsa.append(NallPrsa)
			obsNPrsa.append(NobsPrsa)
			recNPrsa.append(NrecPrsa)
			#print(np.sum(lphRec), np.sum(recN), np.sum(lphRec)/np.sum(recN), np.sum(lphRec0), Nrec, np.sum(lphRec0)/Nrec, np.sum(lphObs), np.sum(obsN), np.sum(lphObs)/np.sum(obsN))

	# Concatenating p and ecc lists
	eccAll = np.concatenate(eccAll)
	eccObs = np.concatenate(eccObs)
	eccRec = np.concatenate(eccRec)

	pAll = np.concatenate(pAll)
	pObs = np.concatenate(pObs)
	pRec = np.concatenate(pRec)


	# print('Ecc lists:', eccAll, eccObs, eccRec)
	# print('P lists:', pAll, pObs, pRec)
	# Appending lists with all the p/ecc values to our dataframes
	# All dataframe
	peccAll['e'] = eccAll
	peccAll['p'] = pAll

	# Observable dataframe
	peccObs['e'] = eccObs
	peccObs['p'] = pObs

	# Recovered dataframe
	peccRec['e'] = eccRec
	peccRec['p'] = pRec

	# print('Final Dataframes:', peccAll, peccObs, peccRec)
	# print(peccRec.columns)
	
	# 3 letter code corresponds to scenario (OC/GC, baseline/colossus, crowding/no crowding)
	peccAll.to_csv('./pecc/all-M67BN-ecc-p.csv', header = ['e', 'p'])
	peccObs.to_csv('./pecc/obs-M67BN-ecc-p.csv', header = ['e', 'p'])
	peccRec.to_csv('./pecc/rec-M67BN-ecc-p.csv', header = ['e', 'p'])

	#plot and save the histograms
	saveHist(m1hAll, m1hObs, m1hRec, m1b, 'm1 (Msolar)', 'EBLSST_m1hist')
	saveHist(qhAll, qhObs, qhRec, qb, 'q (m2/m1)', 'EBLSST_qhist')
	saveHist(ehAll, ehObs, ehRec, eb, 'e', 'EBLSST_ehist')
	saveHist(lphAll, lphObs, lphRec, lpb, 'log(P [days])', 'EBLSST_lphist')
	saveHist(dhAll, dhObs, dhRec, db, 'd (kpc)', 'EBLSST_dhist')
	saveHist(maghAll, maghObs, maghRec, magb, 'mag', 'EBLSST_maghist')
	saveHist(rhAll, rhObs, rhRec, rb, 'r2/r1', 'EBLSST_rhist')


	#make the mollweide
	coords = SkyCoord(RA, Dec, unit=(units.degree, units.degree),frame='icrs')	
	lGal = coords.galactic.l.wrap_at(180.*units.degree).degree
	bGal = coords.galactic.b.wrap_at(180.*units.degree).degree
	RAwrap = coords.ra.wrap_at(180.*units.degree).degree
	Decwrap = coords.dec.wrap_at(180.*units.degree).degree

	f, ax = plt.subplots(subplot_kw={'projection': "mollweide"}, figsize=(8,5))
	ax.grid(True)
	#ax.set_xlabel(r"$l$",fontsize=16)
	#ax.set_ylabel(r"$b$",fontsize=16)
	#mlw = ax.scatter(lGal.ravel()*np.pi/180., bGal.ravel()*np.pi/180., c=np.log10(np.array(recFrac)*100.), cmap='viridis_r', s = 4)
	ax.set_xlabel("RA",fontsize=16)
	ax.set_ylabel("Dec",fontsize=16)
	mlw = ax.scatter(np.array(RAwrap).ravel()*np.pi/180., np.array(Decwrap).ravel()*np.pi/180., c=np.array(recFrac)*100., cmap='viridis_r', s = 4)
	cbar = f.colorbar(mlw, shrink=0.7)
	cbar.set_label(r'% recovered')
	f.savefig('./plots/'  + 'mollweide_pct.pdf',format='pdf', bbox_inches = 'tight')

	f, ax = plt.subplots(subplot_kw={'projection': "mollweide"}, figsize=(8,5))
	ax.grid(True)
	#ax.set_xlabel(r"$l$",fontsize=16)
	#ax.set_ylabel(r"$b$",fontsize=16)
	#mlw = ax.scatter(lGal.ravel()*np.pi/180., bGal.ravel()*np.pi/180., c=np.log10(np.array(recN)), cmap='viridis_r', s = 4)
	ax.set_xlabel("RA",fontsize=16)
	ax.set_ylabel("Dec",fontsize=16)
	mlw = ax.scatter(np.array(RAwrap).ravel()*np.pi/180., np.array(Decwrap).ravel()*np.pi/180., c=np.log10(np.array(recN)), cmap='viridis_r', s = 4)
	cbar = f.colorbar(mlw, shrink=0.7)
	cbar.set_label(r'log10(N) recovered')
	f.savefig('./plots/'  + 'mollweide_N.pdf',format='pdf', bbox_inches = 'tight')

	if (doIndividualPlots):
		fmass.savefig('./plots/'  + 'massPDFall.pdf',format='pdf', bbox_inches = 'tight')
		fqrat.savefig('./plots/'  + 'qPDFall.pdf',format='pdf', bbox_inches = 'tight')
		fecc.savefig('./plots/'  + 'eccPDFall.pdf',format='pdf', bbox_inches = 'tight')
		flper.savefig('./plots/'  + 'lperPDFall.pdf',format='pdf', bbox_inches = 'tight')
		fdist.savefig('./plots/'  + 'distPDFall.pdf',format='pdf', bbox_inches = 'tight')
		fmag.savefig('./plots/'  + 'magPDFall.pdf',format='pdf', bbox_inches = 'tight')
		frad.savefig('./plots/'  + 'radPDFall.pdf',format='pdf', bbox_inches = 'tight')

	print("###################")
	print("number of binaries in input files (raw, log):",np.sum(fileN), np.log10(np.sum(fileN)))
	print("number of binaries in tested with gatspy (raw, log):",np.sum(fileObsN), np.log10(np.sum(fileObsN)))
	print("number of binaries in recovered with gatspy (raw, log):",np.sum(fileRecN), np.log10(np.sum(fileRecN)))
	print("recovered/observable*100 with gatspy:",np.sum(fileRecN)/np.sum(fileObsN)*100.)
	print("###################")
	print("total in sample (raw, log):",np.sum(rawN), np.log10(np.sum(rawN)))
	print("total observable (raw, log):",np.sum(obsN), np.log10(np.sum(obsN)))
	print("total recovered (raw, log):",np.sum(recN), np.log10(np.sum(recN)))
	print("recovered/observable*100:",np.sum(recN)/np.sum(obsN)*100.)
	print("###################")
	print("total in Prsa 15.8<r<19.5 P<1000d sample (raw, log):",np.sum(allNPrsa), np.log10(np.sum(allNPrsa)))
	print("total observable in Prsa 15.8<r<19.5 P<1000d sample (raw, log):",np.sum(obsNPrsa), np.log10(np.sum(obsNPrsa)))
	print("total recovered in Prsa 15.8<r<19.5 P<1000d sample (raw, log):",np.sum(recNPrsa), np.log10(np.sum(recNPrsa)))
	print("Prsa 15.8<r<19.5 P<1000d rec/obs*100:",np.sum(recNPrsa)/np.sum(obsNPrsa)*100.)



[((24, 0, 24, 21), 'matplotlib.use', 'matplotlib.use', ({(24, 15, 24, 20): '"""Agg"""'}, {}), "('Agg')", False, 'import matplotlib\n'), ((44, 8, 44, 58), 'astropy.modeling.models.PowerLaw1D', 'models.PowerLaw1D', (), '', False, 'from astropy.modeling import models, fitting\n'), ((45, 10, 45, 35), 'astropy.modeling.fitting.LevMarLSQFitter', 'fitting.LevMarLSQFitter', ({}, {}), '()', False, 'from astropy.modeling import models, fitting\n'), ((63, 11, 63, 51), 'matplotlib.pyplot.subplots', 'plt.subplots', (), '', True, 'from matplotlib import pyplot as plt\n'), ((65, 11, 65, 33), 'numpy.insert', 'np.insert', ({(65, 21, 65, 28): 'histAll', (65, 29, 65, 30): '0', (65, 31, 65, 32): '0'}, {}), '(histAll, 0, 0)', True, 'import numpy as np\n'), ((66, 11, 66, 33), 'numpy.insert', 'np.insert', ({(66, 21, 66, 28): 'histObs', (66, 29, 66, 30): '0', (66, 31, 66, 32): '0'}, {}), '(histObs, 0, 0)', True, 'import numpy as np\n'), ((132, 15, 132, 39), 'scipy.integrate.quad', 'quad', ({(132, 20, 132, 29): 'RagNormal', (132, 31, 132, 34): '-20', (132, 36, 132, 38): '20'}, {}), '(RagNormal, -20, 20)', False, 'from scipy.integrate import quad\n'), ((156, 9, 156, 44), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((157, 9, 157, 43), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((158, 9, 158, 48), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((159, 10, 159, 47), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((160, 9, 160, 43), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((161, 11, 161, 46), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((162, 9, 162, 46), 'numpy.arange', 'np.arange', (), '', True, 'import numpy as np\n'), ((223, 11, 223, 45), 'pandas.DataFrame', 'pd.DataFrame', (), '', True, 'import pandas as pd\n'), ((224, 11, 224, 45), 'pandas.DataFrame', 'pd.DataFrame', (), '', True, 'import pandas as pd\n'), ((225, 11, 225, 45), 'pandas.DataFrame', 'pd.DataFrame', (), '', True, 'import pandas as pd\n'), ((229, 9, 229, 22), 'os.listdir', 'os.listdir', ({(229, 20, 229, 21): 'd'}, {}), '(d)', False, 'import os\n'), ((407, 10, 407, 32), 'numpy.concatenate', 'np.concatenate', ({(407, 25, 407, 31): 'eccAll'}, {}), '(eccAll)', True, 'import numpy as np\n'), ((408, 10, 408, 32), 'numpy.concatenate', 'np.concatenate', ({(408, 25, 408, 31): 'eccObs'}, {}), '(eccObs)', True, 'import numpy as np\n'), ((409, 10, 409, 32), 'numpy.concatenate', 'np.concatenate', ({(409, 25, 409, 31): 'eccRec'}, {}), '(eccRec)', True, 'import numpy as np\n'), ((411, 8, 411, 28), 'numpy.concatenate', 'np.concatenate', ({(411, 23, 411, 27): 'pAll'}, {}), '(pAll)', True, 'import numpy as np\n'), ((412, 8, 412, 28), 'numpy.concatenate', 'np.concatenate', ({(412, 23, 412, 27): 'pObs'}, {}), '(pObs)', True, 'import numpy as np\n'), ((413, 8, 413, 28), 'numpy.concatenate', 'np.concatenate', ({(413, 23, 413, 27): 'pRec'}, {}), '(pRec)', True, 'import numpy as np\n'), ((450, 10, 450, 75), 'astropy.coordinates.SkyCoord', 'SkyCoord', (), '', False, 'from astropy.coordinates import SkyCoord\n'), ((456, 9, 456, 76), 'matplotlib.pyplot.subplots', 'plt.subplots', (), '', True, 'from matplotlib import pyplot as plt\n'), ((468, 9, 468, 76), 'matplotlib.pyplot.subplots', 'plt.subplots', (), '', True, 'from matplotlib import pyplot as plt\n'), ((68, 15, 68, 40), 'numpy.insert', 'np.insert', ({(68, 25, 68, 35): 'histRec[f]', (68, 36, 68, 37): '0', (68, 38, 68, 39): '0'}, {}), '(histRec[f], 0, 0)', True, 'import numpy as np\n'), ((133, 36, 133, 53), 'numpy.log10', 'np.log10', ({(133, 45, 133, 52): '365 * 10.0'}, {}), '(365 * 10.0)', True, 'import numpy as np\n'), ((146, 18, 146, 32), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((147, 18, 147, 32), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((148, 16, 148, 30), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((149, 18, 149, 32), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((150, 18, 150, 32), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((151, 16, 151, 30), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((152, 16, 152, 30), 'matplotlib.pyplot.subplots', 'plt.subplots', ({}, {}), '()', True, 'from matplotlib import pyplot as plt\n'), ((166, 10, 166, 30), 'numpy.zeros_like', 'np.zeros_like', ({(166, 24, 166, 29): 'mbins'}, {}), '(mbins)', True, 'import numpy as np\n'), ((167, 9, 167, 29), 'numpy.zeros_like', 'np.zeros_like', ({(167, 23, 167, 28): 'qbins'}, {}), '(qbins)', True, 'import numpy as np\n'), ((168, 9, 168, 29), 'numpy.zeros_like', 'np.zeros_like', ({(168, 23, 168, 28): 'ebins'}, {}), '(ebins)', True, 'import numpy as np\n'), ((169, 10, 169, 31), 'numpy.zeros_like', 'np.zeros_like', ({(169, 24, 169, 30): 'lpbins'}, {}), '(lpbins)', True, 'import numpy as np\n'), ((170, 9, 170, 29), 'numpy.zeros_like', 'np.zeros_like', ({(170, 23, 170, 28): 'dbins'}, {}), '(dbins)', True, 'import numpy as np\n'), ((171, 11, 171, 33), 'numpy.zeros_like', 'np.zeros_like', ({(171, 25, 171, 32): 'magbins'}, {}), '(magbins)', True, 'import numpy as np\n'), ((172, 9, 172, 29), 'numpy.zeros_like', 'np.zeros_like', ({(172, 23, 172, 28): 'rbins'}, {}), '(rbins)', True, 'import numpy as np\n'), ((174, 10, 174, 30), 'numpy.zeros_like', 'np.zeros_like', ({(174, 24, 174, 29): 'mbins'}, {}), '(mbins)', True, 'import numpy as np\n'), ((175, 9, 175, 29), 'numpy.zeros_like', 'np.zeros_like', ({(175, 23, 175, 28): 'qbins'}, {}), '(qbins)', True, 'import numpy as np\n'), ((176, 9, 176, 29), 'numpy.zeros_like', 'np.zeros_like', ({(176, 23, 176, 28): 'ebins'}, {}), '(ebins)', True, 'import numpy as np\n'), ((177, 10, 177, 31), 'numpy.zeros_like', 'np.zeros_like', ({(177, 24, 177, 30): 'lpbins'}, {}), '(lpbins)', True, 'import numpy as np\n'), ((178, 9, 178, 29), 'numpy.zeros_like', 'np.zeros_like', ({(178, 23, 178, 28): 'dbins'}, {}), '(dbins)', True, 'import numpy as np\n'), ((179, 11, 179, 33), 'numpy.zeros_like', 'np.zeros_like', ({(179, 25, 179, 32): 'magbins'}, {}), '(magbins)', True, 'import numpy as np\n'), ((180, 9, 180, 29), 'numpy.zeros_like', 'np.zeros_like', ({(180, 23, 180, 28): 'rbins'}, {}), '(rbins)', True, 'import numpy as np\n'), ((490, 55, 490, 68), 'numpy.sum', 'np.sum', ({(490, 62, 490, 67): 'fileN'}, {}), '(fileN)', True, 'import numpy as np\n'), ((491, 62, 491, 78), 'numpy.sum', 'np.sum', ({(491, 69, 491, 77): 'fileObsN'}, {}), '(fileObsN)', True, 'import numpy as np\n'), ((492, 65, 492, 81), 'numpy.sum', 'np.sum', ({(492, 72, 492, 80): 'fileRecN'}, {}), '(fileRecN)', True, 'import numpy as np\n'), ((495, 37, 495, 49), 'numpy.sum', 'np.sum', ({(495, 44, 495, 48): 'rawN'}, {}), '(rawN)', True, 'import numpy as np\n'), ((496, 38, 496, 50), 'numpy.sum', 'np.sum', ({(496, 45, 496, 49): 'obsN'}, {}), '(obsN)', True, 'import numpy as np\n'), ((497, 37, 497, 49), 'numpy.sum', 'np.sum', ({(497, 44, 497, 48): 'recN'}, {}), '(recN)', True, 'import numpy as np\n'), ((500, 62, 500, 78), 'numpy.sum', 'np.sum', ({(500, 69, 500, 77): 'allNPrsa'}, {}), '(allNPrsa)', True, 'import numpy as np\n'), ((501, 73, 501, 89), 'numpy.sum', 'np.sum', ({(501, 80, 501, 88): 'obsNPrsa'}, {}), '(obsNPrsa)', True, 'import numpy as np\n'), ((502, 72, 502, 88), 'numpy.sum', 'np.sum', ({(502, 79, 502, 87): 'recNPrsa'}, {}), '(recNPrsa)', True, 'import numpy as np\n'), ((71, 29, 71, 44), 'numpy.sum', 'np.sum', ({(71, 36, 71, 43): 'histAll'}, {}), '(histAll)', True, 'import numpy as np\n'), ((72, 29, 72, 44), 'numpy.sum', 'np.sum', ({(72, 36, 72, 43): 'histObs'}, {}), '(histObs)', True, 'import numpy as np\n'), ((77, 32, 77, 50), 'numpy.sum', 'np.sum', ({(77, 39, 77, 49): 'histRec[f]'}, {}), '(histRec[f])', True, 'import numpy as np\n'), ((190, 14, 190, 34), 'numpy.zeros_like', 'np.zeros_like', ({(190, 28, 190, 33): 'mbins'}, {}), '(mbins)', True, 'import numpy as np\n'), ((191, 13, 191, 33), 'numpy.zeros_like', 'np.zeros_like', ({(191, 27, 191, 32): 'qbins'}, {}), '(qbins)', True, 'import numpy as np\n'), ((192, 13, 192, 33), 'numpy.zeros_like', 'np.zeros_like', ({(192, 27, 192, 32): 'ebins'}, {}), '(ebins)', True, 'import numpy as np\n'), ((193, 14, 193, 35), 'numpy.zeros_like', 'np.zeros_like', ({(193, 28, 193, 34): 'lpbins'}, {}), '(lpbins)', True, 'import numpy as np\n'), ((194, 13, 194, 33), 'numpy.zeros_like', 'np.zeros_like', ({(194, 27, 194, 32): 'dbins'}, {}), '(dbins)', True, 'import numpy as np\n'), ((195, 15, 195, 37), 'numpy.zeros_like', 'np.zeros_like', ({(195, 29, 195, 36): 'magbins'}, {}), '(magbins)', True, 'import numpy as np\n'), ((196, 13, 196, 33), 'numpy.zeros_like', 'np.zeros_like', ({(196, 27, 196, 32): 'rbins'}, {}), '(rbins)', True, 'import numpy as np\n'), ((236, 12, 236, 37), 'pandas.read_csv', 'pd.read_csv', (), '', True, 'import pandas as pd\n'), ((490, 79, 490, 92), 'numpy.sum', 'np.sum', ({(490, 86, 490, 91): 'fileN'}, {}), '(fileN)', True, 'import numpy as np\n'), ((491, 89, 491, 105), 'numpy.sum', 'np.sum', ({(491, 96, 491, 104): 'fileObsN'}, {}), '(fileObsN)', True, 'import numpy as np\n'), ((492, 92, 492, 108), 'numpy.sum', 'np.sum', ({(492, 99, 492, 107): 'fileRecN'}, {}), '(fileRecN)', True, 'import numpy as np\n'), ((495, 60, 495, 72), 'numpy.sum', 'np.sum', ({(495, 67, 495, 71): 'rawN'}, {}), '(rawN)', True, 'import numpy as np\n'), ((496, 61, 496, 73), 'numpy.sum', 'np.sum', ({(496, 68, 496, 72): 'obsN'}, {}), '(obsN)', True, 'import numpy as np\n'), ((497, 60, 497, 72), 'numpy.sum', 'np.sum', ({(497, 67, 497, 71): 'recN'}, {}), '(recN)', True, 'import numpy as np\n'), ((500, 89, 500, 105), 'numpy.sum', 'np.sum', ({(500, 96, 500, 104): 'allNPrsa'}, {}), '(allNPrsa)', True, 'import numpy as np\n'), ((501, 100, 501, 116), 'numpy.sum', 'np.sum', ({(501, 107, 501, 115): 'obsNPrsa'}, {}), '(obsNPrsa)', True, 'import numpy as np\n'), ((502, 99, 502, 115), 'numpy.sum', 'np.sum', ({(502, 106, 502, 114): 'recNPrsa'}, {}), '(recNPrsa)', True, 'import numpy as np\n'), ((271, 19, 271, 55), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((272, 17, 272, 64), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((273, 17, 273, 52), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((275, 17, 275, 52), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((276, 21, 276, 69), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((277, 17, 277, 64), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((289, 10, 289, 22), 'numpy.diff', 'np.diff', ({(289, 18, 289, 21): 'm1b'}, {}), '(m1b)', True, 'import numpy as np\n'), ((463, 95, 463, 112), 'numpy.array', 'np.array', ({(463, 104, 463, 111): 'recFrac'}, {}), '(recFrac)', True, 'import numpy as np\n'), ((475, 104, 475, 118), 'numpy.array', 'np.array', ({(475, 113, 475, 117): 'recN'}, {}), '(recN)', True, 'import numpy as np\n'), ((493, 47, 493, 63), 'numpy.sum', 'np.sum', ({(493, 54, 493, 62): 'fileRecN'}, {}), '(fileRecN)', True, 'import numpy as np\n'), ((493, 64, 493, 80), 'numpy.sum', 'np.sum', ({(493, 71, 493, 79): 'fileObsN'}, {}), '(fileObsN)', True, 'import numpy as np\n'), ((498, 35, 498, 47), 'numpy.sum', 'np.sum', ({(498, 42, 498, 46): 'recN'}, {}), '(recN)', True, 'import numpy as np\n'), ((498, 48, 498, 60), 'numpy.sum', 'np.sum', ({(498, 55, 498, 59): 'obsN'}, {}), '(obsN)', True, 'import numpy as np\n'), ((503, 47, 503, 63), 'numpy.sum', 'np.sum', ({(503, 54, 503, 62): 'recNPrsa'}, {}), '(recNPrsa)', True, 'import numpy as np\n'), ((503, 64, 503, 80), 'numpy.sum', 'np.sum', ({(503, 71, 503, 79): 'obsNPrsa'}, {}), '(obsNPrsa)', True, 'import numpy as np\n'), ((37, 25, 37, 36), 'numpy.sqrt', 'np.sqrt', ({(37, 33, 37, 35): '(2.0)'}, {}), '(2.0)', True, 'import numpy as np\n'), ((247, 10, 247, 38), 'pandas.read_csv', 'pd.read_csv', (), '', True, 'import pandas as pd\n'), ((294, 20, 294, 33), 'numpy.log10', 'np.log10', ({(294, 29, 294, 32): 'Phs'}, {}), '(Phs)', True, 'import numpy as np\n'), ((318, 20, 318, 55), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((319, 18, 319, 63), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((320, 18, 320, 52), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((322, 18, 322, 52), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((323, 22, 323, 69), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((324, 18, 324, 63), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((334, 19, 334, 33), 'pandas.DataFrame', 'pd.DataFrame', ({}, {}), '()', True, 'import pandas as pd\n'), ((335, 23, 335, 37), 'pandas.DataFrame', 'pd.DataFrame', ({}, {}), '()', True, 'import pandas as pd\n'), ((463, 18, 463, 34), 'numpy.array', 'np.array', ({(463, 27, 463, 33): 'RAwrap'}, {}), '(RAwrap)', True, 'import numpy as np\n'), ((463, 55, 463, 72), 'numpy.array', 'np.array', ({(463, 64, 463, 71): 'Decwrap'}, {}), '(Decwrap)', True, 'import numpy as np\n'), ((475, 18, 475, 34), 'numpy.array', 'np.array', ({(475, 27, 475, 33): 'RAwrap'}, {}), '(RAwrap)', True, 'import numpy as np\n'), ((475, 55, 475, 72), 'numpy.array', 'np.array', ({(475, 64, 475, 71): 'Decwrap'}, {}), '(Decwrap)', True, 'import numpy as np\n'), ((274, 32, 274, 61), 'numpy.ma.log10', 'np.ma.log10', ({(274, 44, 274, 60): "data['p'].values"}, {}), "(data['p'].values)", True, 'import numpy as np\n'), ((280, 36, 280, 51), 'numpy.sum', 'np.sum', ({(280, 43, 280, 50): 'm1hAll0'}, {}), '(m1hAll0)', True, 'import numpy as np\n'), ((281, 34, 281, 48), 'numpy.sum', 'np.sum', ({(281, 41, 281, 47): 'qhAll0'}, {}), '(qhAll0)', True, 'import numpy as np\n'), ((282, 33, 282, 47), 'numpy.sum', 'np.sum', ({(282, 40, 282, 46): 'ehAll0'}, {}), '(ehAll0)', True, 'import numpy as np\n'), ((283, 36, 283, 51), 'numpy.sum', 'np.sum', ({(283, 43, 283, 50): 'lphAll0'}, {}), '(lphAll0)', True, 'import numpy as np\n'), ((284, 34, 284, 48), 'numpy.sum', 'np.sum', ({(284, 41, 284, 47): 'dhAll0'}, {}), '(dhAll0)', True, 'import numpy as np\n'), ((285, 37, 285, 53), 'numpy.sum', 'np.sum', ({(285, 44, 285, 52): 'maghAll0'}, {}), '(maghAll0)', True, 'import numpy as np\n'), ((286, 33, 286, 47), 'numpy.sum', 'np.sum', ({(286, 40, 286, 46): 'rhAll0'}, {}), '(rhAll0)', True, 'import numpy as np\n'), ((361, 22, 361, 57), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((362, 20, 362, 65), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((363, 20, 363, 54), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((365, 20, 365, 54), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((366, 24, 366, 71), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((367, 20, 367, 65), 'numpy.histogram', 'np.histogram', (), '', True, 'import numpy as np\n'), ((321, 33, 321, 61), 'numpy.ma.log10', 'np.ma.log10', ({(321, 45, 321, 60): "obs['p'].values"}, {}), "(obs['p'].values)", True, 'import numpy as np\n'), ((364, 35, 364, 63), 'numpy.ma.log10', 'np.ma.log10', ({(364, 47, 364, 62): "rec['p'].values"}, {}), "(rec['p'].values)", True, 'import numpy as np\n')]