<jupyter_start><jupyter_text>Code to determine iterative distance estimates for stars in the anticentre  As used by "Gaia Early Data Release 3: The Galactic anticentre" Table of results used for that paper available via zenodo<jupyter_code># input file - taken from a request to the Gaia archive like that in the 
# appendix of "Gaia Early Data Release 3: The Galactic anticentre"
filename = 'EDR3_input_data.csv'
# File for tabulated density prior
prior_filename = 'Density_Prior_' + filename
# Output files (csv and hdf)
outname = 'EDR3_anticentre_distances'


# Parallax zeropoint correction
correctiontype = 'const17'
zeropoint_filename = '' # Only needed if tabulauted from e.g. the Lindegren correction

# Range in l,b (note that symmetry of field in b is assumed)
lmin,lmax = 170,190
bmax = 10
# file containing parallax, parallax_error, l, b
df = pd.read_csv(filename)
#df = pd.read_hdf('EDR3_input_data.h5','table')
zpcode = ''

if correctiontype == 'zeropoint_file' :
# file containing zero-point values
    df_z = pd.read_csv(zeropoint_filename)
    df = pd.merge(df_z,df,on='source_id')
    df['parallax_corrected'] = df['parallax'] - df['zpt']
    df_z =  []
    zpcode = '_zp_file'
elif correctiontype == 'none' :
    df['parallax_corrected'] = df['parallax']
    zpcode = ''
elif correctiontype == 'const17' :
    df['parallax_corrected'] = df['parallax'] + 0.017
    zpcode = '_zp17'
else :
    print('No correction given, assuming none')
    df['parallax_corrected'] = df['parallax']
    zpcode = ''
# Produce prior file
subprocess.run(["./Prior_calc_anticentre.exe", 
                prior_filename,'%f' % lmin,'%f' % lmax, '%f' % bmax])
# and read it back in
df_prior = pd.read_csv(prior_filename)
# radial range over which this is tabulated
rrange = np.linspace(0,10,51)<jupyter_output><empty_output><jupyter_text>For $\varpi/\sigma_\varpi<3$ distance estimates become much less useful, so we only use  values with $\varpi/\sigma_\varpi\geq3$. To enable easier comparison between samples, we use the quoted parallax, rather than the corrected one.<jupyter_code>df_good = df[df.parallax>=3.*df.parallax_error]
df_good<jupyter_output><empty_output><jupyter_text>Put into format for reading by the C++ code<jupyter_code>df_good.to_csv(filename+zpcode+'.input',sep=' ',
               columns=['source_id', 'l','b','parallax_corrected','parallax_error'], 
              index=False, header=False)<jupyter_output><empty_output><jupyter_text>Run distance code with estimated selection function. Then iterate to improve estimated selection functionSee appendix C of "Gaia Early Data Release 3: The Galactic anticentre" for details<jupyter_code>nIterations = 5
runDistanceCode=True

plt.gcf().set_size_inches(15,10)
plt.subplot(2,3,1)

# Initial estimate
hist,bin_edges = np.histogram(1./(df_good['parallax_corrected']),rrange,density=True)
rcentres = 0.5*( bin_edges[:-1]+bin_edges[1:] )

# Estimate comes from the range 1-3 kpc
maskdist = (rcentres< 3.) & (rcentres>1.)
color = 0
plt.plot(rcentres,hist/df_prior.integrated_prior.values,
         color='C%d' % color)
slope, intercept, _, _, _ = stats.linregress(rcentres[maskdist],
                                            np.log(hist[maskdist]
                                                   /df_prior.integrated_prior[maskdist]))
# Slope = 1/scalelength for selection function
plt.plot(rcentres,np.exp(intercept+slope*rcentres),linestyle=':',color='C%d' % color)


plt.yscale('log')
plt.title('slope = %.4f' % slope)

for i in range (nIterations) : 
    plt.subplot(2,3,i+2)
    outname = filename+zpcode+'.out%d' % (i+1)
    
    if runDistanceCode :
        Runcode = subprocess.run(["./Distance_calc_anticentre.exe", 
                                  filename+zpcode+'.input',
                                  outname,
                                 '%f' % slope])
    print ('.',end='')
    
    df_tmp_out = pd.read_csv(outname)
    df_tmp = pd.merge(df_good,df_tmp_out,on='source_id')
    parcuts = [5,4,3]
    rrange = np.linspace(0,10,51)

    color = 0
    for pc in parcuts :
        testmask = df_tmp['parallax_corrected']/df_tmp['parallax_error'] > pc
        hist,bin_edges = np.histogram((df_tmp['distance'])[testmask],rrange,density=True)
        rcentres = 0.5*( bin_edges[:-1]+bin_edges[1:] )
        
        maskdist = (rcentres< 3.) & (rcentres>1.)
        plt.plot(rcentres,hist/df_prior.integrated_prior.values,
                label='%d' % pc,color='C%d' % color)
        slope, intercept,  _, _, _  = stats.linregress(rcentres[maskdist],
                                                        np.log(hist[maskdist]
                                                            /df_prior.integrated_prior[maskdist]))
        plt.plot(rcentres,np.exp(intercept+slope*rcentres),linestyle=':',color='C%d' % color)
        color+=1
        #print(pc,slope)
    plt.yscale('log')
    plt.legend()
    plt.title('slope = %.4f' % slope)

plt.show()
df_write = pd.read_csv(filename+zpcode+'.out5')


df_write.to_csv(outname+zpcode+'.csv',float_format='%g',index=False)
df_write.to_hdf(outname+zpcode+'.h5','table')

df_write<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13122083981337482, 'fraction_numerical': 0.014774494556765163, 'mean_word_length': 3.4011976047904193, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '206994', 'license': 'MIT', 'n_tokens_mistral': 1843, 'n_tokens_neox': 1740, 'n_words': 362, 'path': 'PrepareForPureDistanceCode-Final.ipynb', 'repo_name': 'PaulMcMillan-Astro/AnticentreBayesianDistances', 'size': '5144'}