import os
import numpy as np

# This is table1 from 2019A&A...621A..86B
cs = np.genfromtxt('table1.txt')
dg = cs[:,0]
dists = cs[:,1:]


pairs_found = open('mock_pairs.csv')

with open('tmp_raw_sens.txt', 'w') as output_file:
    output_file.write("#sid, percentage\n")
    x_values = np.array([1,10,50,90,99])
    for i, line in enumerate(pairs_found):
        if i%1000000 == 0:
            print(i)
        # First line is the header
        if i == 0:
            t =[x.strip() for x in line.split(',')]
            print(t)
            continue

        # Extracting the values and converting to numbers
        # gmag2 needs to be fainter than gmag1
        t =[x.strip() for x in line.split(',')]
        #sid1 = np.int64(t[0])
        #sid2 = np.int64(t[1])
        gmag1 = float(t[2])
        gmag2 = float(t[3])
        dist = float(t[4])
        deltagmag = gmag2-gmag1
        # Find where the deltagmag is closest to the contras sensitivity table entry
        cut = np.where(np.abs(dg-deltagmag) == np.min(np.abs(dg-deltagmag)))
        # simply use the distance limits in the contrast sensitivity table of that line
        func = dists[cut[0]][0]
        #print(i, deltagmag,dist,func)
        # if distance is bigger than the 99% distance of the cs.txt than 100% visibility, jump to next source
        if func[-1] < dist:
            continue
        # if distance is smaller than 1% distance, visibility is 0%
        elif func[0] > dist:
            p_value = 0
        # if distance is somewhere in between 1% and 99% distance, then interpolate linearly and assign the respective probability
        else:
            p_value = np.interp(dist, func, x_values)
        #print('#############', p_value)
        # write source_id of the fainter star and of its visibility
        output_file.write("%d, %.0f\n" %(np.int64(t[1]),p_value))

pairs_found.close()

# Now erase the double entries
x = np.genfromtxt('tmp_raw_sens.txt',
  delimiter=',',
  dtype=[('sid', np.int64),
  ('percentage', np.int16)],
  skip_header=1)
# Sort first by source_id and then by visibility probability
x = np.sort(x, kind = 'mergesort', order = ['sid', 'percentage'])


# sort worked and if np.unique is applied only the first value is retained,
# therefore if the same star is blended multiple times, the worst percentage of
# visibility will be shown
u, indices = np.unique(x['sid'], return_index=True)
print(len(x))
x = x[indices]
print('after cutting multiple entries len: ', len(x))

np.savetxt('../data/detectability.csv',
  x,
  fmt='%d, %d',
  header = 'sid, detectability_percentage',
  comments = '#')

os.unlink("tmp_raw_sens.txt")