from rsf.proj import *

def read_chevron(target=None,source=None,env=None):
    inp = open(str(source[0]),'r')
    file = str(target[0])
    out = open(file,'w')
    n2=0
    for line in inp.readlines():
        n2=n2+1
        out.write(string.join(line.split()[2:],' ')+'\n')
    out.write("n1=3 n2=%d esize=0 in=%s data_format=ascii_float\n" %
              (n2,file))
    out.close()
    inp.close()

Chevron = Builder(action=Action(read_chevron),src_suffix='.5',suffix='.T')
env = Environment(BUILDERS={'Chevron':Chevron})

sets = ['seismic','well']
baseplate = -2782
def grey(title,clip='allpos=y clip=2000'):
    return '''
    grey %s crowd=0.88 transp=n yreverse=n
    title="%s" label1="X axis" label2="Y axis"
    ''' % (clip,title)

for data in sets:
    file = data[0]+'2000'
    Fetch(file+'.5','chevron')
    env.Chevron(file+'.T',file+'.5')
    Flow(file,file+'.T','dd form=native')
    Flow(data,file,
         '''
         window n1=1 f1=2 | math output="input-(%g)" |
         bin interp=1 nx=40 ny=40 xkey=0 ykey=1 head=$SOURCE
         xmin=2000 xmax=28000 ymin=-25000 ymax=12000
         ''' % baseplate)
    Plot(data,grey('Binned %s data' % data))
    Flow(data+'1',[data,file],
         'frame xyz=${SOURCES[1]} base=%g' % baseplate)
    Plot(data+'1',grey('Binned and framed %s data' % data))

Result('wellseis',sets,'SideBySideIso')

Flow('sx','seismic1','lapfill niter=500')
Plot('sx',grey('Extended seismic'))

Result('misseis',['seismic1','sx'],'SideBySideIso')

Flow('mask','well1','add mode=d $SOURCE')
Flow('vv',['well1','sx','mask'],
     '''
     add scale=1,-1 ${SOURCES[1]} |
     add mode=p ${SOURCES[2]}
     ''')
Plot('vv',grey('Wells - Seismic',clip='pclip=100'))

for case in [1,2]:
    diff = 'd%d' % case
    modl = 'm%d' % case
    cont = 'c%d' % case

    case = 2-case
    label = ['Laplacian','Gradient'][case]

    Flow(diff,'vv','lapfill niter=500 grad=%d' % case)
    Plot(diff,grey('%s Map - Seismic' % label,clip='pclip=100'))

    Flow(modl,[diff,'sx'],'add ${SOURCES[1]}')
    Plot(modl,grey('Map based on %s' % label))
    Plot(cont,modl,
         '''
         contour nc=50 c0=0 transp=n yreverse=n crowd=.88
         title="Map based on %s"
         ''' % label)

Result('diffdiff',['vv','d1','d2'],'SideBySideAniso')
Result('finalmap',['m2','c2'],'SideBySideAniso')

###########################################################################
End()