#! /usr/bin/python

# -c 3 -s 0-12-2,19,20,22-28,35,36  -o radftest_tail.S-490295028.mkcalib radftest_tail.S-490295028-patched

import sys

from getopt import getopt
opt,files = getopt(sys.argv[1:], "c:s:o:")
ch = 0
sa=range(16)
out = sys.stdout
for o,v in opt:

    if o=="-c":
        ch = int(v)
    if o=="-s":
        sa = []
        for s in v.split(","):
            ss = s.split("-")
            if len(ss)==1:
                sa.append(int(ss[0]))
            else:
                sa += list(range(*tuple([int(sss) for sss in ss])))
        print >>sys.stderr, "sa =", sa
                
    if o=="-o":
        out = open(v, "w")

if not files:
    f = sys.stdin
else:
    f = open(files[0])

d=[]
for l in f:
    d.append(int(l.split()[19-ch]))

from numarray import *

print >>sys.stderr, "d =", d

d = array(d)
m = d*[0.0]
for i in sa:
    m[i] = 1

a=(d-sum(m*d)/sum(m))*m
dd=d*[0.0]
dd[1:-1]=d[2:]-d[:-2]
dmax = max(d)
dmin = min(d)
d1 = dmin+0.2*(dmax-dmin)
d2 = dmin+0.8*(dmax-dmin)
mm=m*(d>d1)*(d<d2)
b=(dd-sum(mm*dd)/sum(mm))*mm

bl=(d<(dmin+0.02*(dmax-dmin)))*m
where(bl,sum(a*bl)/sum(bl),a,a)

a /= max(abs(min(a)), max(a))
b /= max(abs(min(b)), max(b))

class Interpolation:
    def __init__(self, d):
        self.d = d*1.0;
        y1 = self.d[:-2]
        y2 = self.d[1:-1]
        y3 = self.d[2:]
        self.a = (y3-2*y2+y1)/2
        self.b = (y3-y1)/2
        self.c = y2

    def y(self, x):
        r = self.d.copy()
        r[1:-1] = (self.a*x+self.b)*x+self.c
        return r
    
    def yy(self, x):
        r1 = (self.a*x+self.b)*x+self.c
        r = self.d.copy()
        if x>=0:
            r2 = (self.a*(x-1)+self.b)*(x-1)+self.c
            r[1:-2] = (r1[:-1]+r2[1:])/2
        else:
            r2 = (self.a*(x+1)+self.b)*(x+1)+self.c
            r[2:-1] = (r1[1:]+r2[:-1])/2
        return r

interpol = Interpolation(d)

N=64

for ix in range(-N,N):
    x = float(ix)/N
    r = interpol.yy(x)
    print >>out, x, sum(a*r), sum(b*r)

ii = [(float(ix)/N, interpol.yy(float(ix)/N)) for ix in range(-N/2,N/2)]

for i in range(len(d)):
    for x,r in ii:
        print >>out, i+x, r[i]

if files:
    out.close()

aa=around(a*4095)
bb=around(b*4095)

t0=int(max(arange(len(d))*m))
for i in range(len(d)):
    if m[t0-i]:
        print >>sys.stderr, "(%2d, %5d, %5d)," % (i, int(aa[t0-i]), int(bb[t0-i]))