#!/usr/bin/env python
import sys, numpy, scipy.signal

class shape(object):
    px = []
    py = []
    hy = []
    h = None
    absorptionshape = False

    def __init__(self, stream):
        for line in stream:
            if line.startswith('#'):
                continue
            parts = line.strip().split('\t')
            if len(parts) < 2:
                continue
            self.px.append(float(parts[0]))
            self.py.append(float(parts[1]))
    
    def percentage(self):
        a = abs(min(self.py))
        d = max(self.py)
        A = min(a, d)
        D = max(a, d)
        a = abs(min(self.hy))
        d = max(self.hy)
        Ah = min(a, d)
        Dh = max(a, d)
        percent = 2*A/(A+D)
        p = ((2*Ah)/(Ah+Dh)-percent)/(1-percent) 
        if p < 0 or p > 1:
            percent = 2*Ah/(Ah+Dh)
            p = ((2*A)/(A+D)-percent)/(1-percent)
            self.absorptionshape = True
        return p

    def weight(self):
        self.h = scipy.signal.hilbert(self.py)
        self.hy = numpy.imag(self.h)

    def print_spectra(self):
        p = self.percentage()
        if self.absorptionshape:
            sys.stderr.write('shape: absorption\n')
        else:
            sys.stderr.write('shape: dispersion\n')
        sys.stderr.write('percentage %f\n' % p)
        b = (1-p)/(p*p+(1-p)*(1-p))
        a = p/(p*p+(1-p)*(1-p))
        if self.absorptionshape:
            t = a
            a = b
            b = t
        # chi' (dispersion)
        mvalx, mvaly = (0,0)
        for i in zip(self.px, self.py, self.hy):
            val = i[1]*a-i[2]*b
            if abs(val) > mvaly:
                mvaly = abs(val)
                mvalx = i[0]
            print i[0], val
        print '\n\n'
        dp = mvalx
        sys.stderr.write('dispersion peak position %f\n' % mvalx)
        # chi'' (absorption)
        mvalx = 0
        mvaly /= 10
        rbefore = False
        found = False
        for i in zip(self.px, self.py, self.hy):
            val = i[1]*b+i[2]*a
            if not found and abs(val) > mvaly:
                rbefore = True
                mvaly = val
                found = True
            if rbefore and abs(mvaly+val) < abs(mvaly):
                mvalx = i[0]
                rbefore = False
            print i[0], val
        sys.stderr.write('absorption peak position %f\n' % mvalx)
        sys.stderr.write('peak position (%f \pm %f)\n' % ((mvalx+dp)/2, abs(mvalx-dp)/2))
        

# read peak data from stdin
s = shape(sys.stdin)
s.weight()
s.print_spectra()