from rsf.proj import *
import math

###### Draw the Box ######
min1,max1=0.3,0.7
min2,max2=0.85,1.1
Flow('box.asc',None,'''echo %s n1=2 n2=5 data_format=ascii_float
    in=$TARGET'''%' '.join(map(str,(min1,min2,max1,min2,
                                    max1,max2,min1,max2,
                                    min1,min2))))
Plot('box','box.asc','''dd form=native type=complex |
    window | graph transp=y yreverse=y min1=0 max1=1.2 min2=0 max2=1.28
    wanttitle=n plotfat=5 plotcol=1 wantaxis=n''')

#####################################
###### Extract Fault Attribute ######
#####################################
Flow('sigmoid',None,'''sigmoid d1=.0032 n1=400 d2=.0032 n2=400 taper=n |
    smooth rect1=5 diff1=1 | smooth rect1=5''')
Plot('sigmoid','''window max1=1.2 |
    grey title="2D Synthetic" label1=Time label2=Distance''')
Result('sigmoid','sigmoid box','Overlay')

###### Dip Estimation ######
Flow('dip','sigmoid','dip2 rect1=10 rect2=10 order=4')
Result('syndip','dip','''window max1=1.2 |
    grey color=j label2="Distance" title=Dip scalebar=y''')

###### Plane Wave Destruction ######
Flow('pwd','sigmoid dip','pwd2 dip=${SOURCES[1]} order=3')
Flow('pws','sigmoid dip','pwsmooth dip=${SOURCES[1]} ns=3 order=3')

###### Sobel Filter ######
Flow('sobel1','pwd pws','''math pwd=${SOURCES[0]} pws=${SOURCES[1]}
    output="pwd*pwd+pws*pws"''')
Result('sobel1','grey title="Old Sobel" allpos=y')

Flow('dip12','dip','transp plane=12')
Flow('sobel2','pwd dip12',
     '''transp plane=12 |
     pwsmooth dip=${SOURCES[1]} ns=3 order=3 |
     transp plane=12''')

Flow('sobel','sobel1 sobel2','''math s1=${SOURCES[0]} s2=${SOURCES[1]}
    output="s1*s1+s2*s2"''')
Result('synsobel','sobel','''window max1=1.2 |
    grey allpos=y label2="Distance" title="Sobel"''')

Flow('dip2','sobel','odip rect1=10 rect2=10 | transp')
Result('dip2','transp | grey color=j')
Flow('smooth','sobel dip2','''transp |
    pwsmooth dip=${SOURCES[1]} ns=5 order=3 | transp |
    smooth rect1=5 rect2=5''')
Result('smooth','''window max1=1.2 |
    grey allpos=y label2="Distance" title="Smoothed Sobel"''')


##########################################
###### Well Log Predictive Painting ######
##########################################
amp=20
wellLoc = [30,200,370]
scale=10.0
r0=100.0

for well in range(3):
    log = 'log%d'%well
    Flow(log,'sigmoid','window n2=1 f2=%d'%wellLoc[well])
    Plot(log,'''scale axis=1 | math output="input*%g+%g" | window max1=1.2 |
        graph transp=y yreverse=y min2=0 max2=400 min1=0 max1=1.2
        wherexlabel=top wheretitle=bottom wanttitle=n wantaxis=n
        plotfat=3'''%(amp,wellLoc[well]))

    paint = 'paint%d'%well
    Flow(paint,['dip',log],'''pwpaint order=3 seed=${SOURCES[1]} eps=0.1
        i0=%d'''%wellLoc[well])
    Plot(paint,'''window max1=1.2 |
        grey title="Predictive Painting from Log"''')
    Result(paint,[paint,log,'box'],'Overlay')

    rbfold = 'rbfold%d'%well
    Flow(rbfold,'smooth','''faultrbf1d useinput=n
        i0=%d scale=%g r0=%g'''%(wellLoc[well],scale,r0))

    old = 'old%d'%well
    Flow(old,[paint,rbfold],'add mode=m ${SOURCES[1:2]}')

    rbf = 'rbf%d'%well
    Flow(rbf,'smooth','''faultrbf1d useinput=y
        i0=%d scale=%g r0=%g'''%(wellLoc[well],scale,r0))
    Result(rbf,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
        title="RBF of log %d"'''%(well+1))

    weighted = 'weighted%d'%well
    Flow(weighted,[paint,rbf],'add mode=m ${SOURCES[1:2]}')

    seed = 'seed%d'%well
    Flow(seed,log,'math output="input*0.0"')
    dist = 'dist%d'%well
    Flow(dist,['dip',seed,'smooth'],'''distpaint order=3 eps=0.1 i0=%d
        seed=${SOURCES[1]} flt=${SOURCES[2]} faultscale=5 |
        clip2 upper=100'''%wellLoc[well])
    rbfnew = 'rbfnew%d'%well
    Flow(rbfnew,dist,'''math output="input*300" |
        math output="1./(1.+(input/%g)*(input/%g))"'''%(r0,r0))
    Result(rbfnew,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
        title="New RBF of log %d"'''%(well+1))

    rbfweird = 'rbfweird%d'%well
    Flow(rbfweird,['dip',seed,'smooth'],'''distpaint order=3 eps=0.1 i0=%d
        seed=${SOURCES[1]} flt=${SOURCES[2]} faultscale=1 |
        math output="input*300" |
        math output="1./(1.+(input/%g)*(input/%g))"
        '''%(wellLoc[well],r0,r0))
    Result(rbfweird,'''window max1=1.2 | grey color=j scalebar=y bias=0.5
        title="RBF of log %d"'''%(well+1))

    weightednew = 'weightednew%d'%well
    Flow(weightednew,[paint,rbfnew],'add mode=m ${SOURCES[1:2]}')

Flow('rbfoldsum','rbfold0 rbfold1 rbfold2','add ${SOURCES[1:3]}')
Flow('oldsum','old0 old1 old2','add ${SOURCES[1:3]}')
Flow('interpold','oldsum rbfoldsum','add mode=d ${SOURCES[1:2]}')
Plot('interpold','''window max1=1.2 |
    grey title="RBF interpolation without fault"''')
Result('interpold','interpold log0 log1 log2 box','Overlay')

Flow('errold','sigmoid interpold',
    'add scale=1,-1 ${SOURCES[1:2]} | math output="abs(input)"')
Result('errold','''window max1=1.2 |
    grey scalebar=y bias=0.0015 clip=0.001 maxval=0.003
    title="Error of interpolation without fault"''')

Flow('rbfsum','rbf0 rbf1 rbf2','add ${SOURCES[1:3]}')
Flow('weightedsum','weighted0 weighted1 weighted2','add ${SOURCES[1:3]}')
Flow('interp','weightedsum rbfsum','add mode=d ${SOURCES[1:2]}')
Plot('interp','''window max1=1.2 |
    grey title="RBF interpolation with fault"''')
Result('interp','interp log0 log1 log2 box','Overlay')

Flow('err','sigmoid interp',
    'add scale=1,-1 ${SOURCES[1:2]} | math output="abs(input)"')
Result('err','''window max1=1.2 |
    grey scalebar=y bias=0.0015 clip=0.001 maxval=0.003
    title="Error of interpolation with Sobel"''')

Flow('rbfnewsum','rbfnew0 rbfnew1 rbfnew2','add ${SOURCES[1:3]}')
Flow('weightednewsum','weightednew0 weightednew1 weightednew2',
    'add ${SOURCES[1:3]}')
Flow('interpnew','weightednewsum rbfnewsum','add mode=d ${SOURCES[1:2]}')
Plot('interpnew','''window max1=1.2 |
    grey title="New RBF interpolation with fault"''')
Result('interpnew','interpnew log0 log1 log2 box','Overlay')

End()