# Time-warping with randomly sampled inputs for inversion
# We do 2D by cutting from 3D to expedite revision process
#########################################################

from rsf.proj import *
import math

## Plot font and screen ratio, width, height for uniform ppt figs.
p1=0.7
p2=0.7
sr=1.0
sw=7.0
sh=10.0
swc=8.0
llw=5.0
srw=1.0


## Number of time slices to do inversion on (= number of events)
numtime=1


#####
##Synthetic CMP Parameters
nx=151
delx=0.02*2
ox=-3
ny=nx
dely=delx
oy=ox
nt=501
dt=0.004
fw=10.0
####

###############################################################################################################
####   Prepare Geometry Header Files

## Absolute offset
Flow('offset','spikes','window n1=1 | math output="sqrt(x1*x1+x2*x2)" ',local=1)
#Result('offset','grey title=Offset allpos=y color=j scalebar=y',local=1)

## Azimuth
Flow('azimuth','spikes','window n1=1 | math output="%g*(x2&x1)" ' % (180/math.pi),local=1)
#Result('azimuth','grey title=Azimuth color=j scalebar=y',local=1)

## Header File
Flow('head','offset azimuth',
     'cat axis=3 ${SOURCES[1]} | transp plane=23 | transp plane=12 | put n3=1 n2=10201',local=1)

####################################################################################################################
## Single cmp with W(t)

####### Parameter values#####################################################
####### Generalized moveout approximation ###################################
# Dry Green Shale
# {c11 -> 15.06379696`, c33 -> 10.837264`, c13 -> 1.638114075744646`, c55 -> 3.125824`}
# {v -> 2.463501288062909`, vs -> 1.7679988687779187`, vp -> 3.292005467796188`, \[Eta] -> 0.741077659341173`}
# {W -> 0.164776, A -> -0.0804846, B -> 0.751616, C -> 0.00440688}

## Wx
Wx0=0.164776
WxE=0.164776
dWx=(WxE-Wx0)/500.0
Flow('begWx',None,'spike n1=126 | scale dscale=%g' % Wx0)
Flow('midWx',None,'spike n1=250 mag=%g d1=%g | math output="input+x1"' % (Wx0,dWx))
Flow('endWx',None,'spike n1=125 | scale dscale=%g' % WxE)
Flow('Wxoft','begWx midWx endWx','cat ${SOURCES[1]} ${SOURCES[2]} axis=1')
## Wy
Flow('Wyoft','Wxoft','math output="input" ')
## Wxy
Flow('Wxyoft','Wxoft','math output="0" ')

## A
Flow('A1oft','Wxoft','math output="-0.0804846" ')
Flow('A2oft','Wxoft','math output="0" ')
Flow('A3oft','Wxoft','math output="-0.01" ') # dummy
Flow('A4oft','Wxoft','math output="0" ')
Flow('A5oft','Wxoft','math output="-0.0804846" ')

## B
Flow('B1oft','Wxoft','math output="0.751616" ')
Flow('B2oft','Wxoft','math output="0" ')
Flow('B3oft','Wxoft','math output="0.751616" ')

## C
Flow('C1oft','Wxoft','math output="0.00440688" ')
Flow('C2oft','Wxoft','math output="0" ')
Flow('C3oft','Wxoft','math output="0.001" ') # dummy
Flow('C4oft','Wxoft','math output="0" ')
Flow('C5oft','Wxoft','math output="0.00440688" ')

## Values of parameters in time
Flow('Woft','Wxoft Wxyoft Wyoft A1oft A2oft A3oft A4oft A5oft B1oft B2oft B3oft C1oft C2oft C3oft C4oft C5oft','cat ${SOURCES[1:16]} axis=2')

## Create spikes and inverse nmo 3D CMP gathers ########################################################
it1=251
itplot=300


Flow('spikes',None,
     '''
     spike n1=%d nsp=1 k1=%d | ricker1 frequency=%g |
     spray axis=2 n=%d o=%g d=%g |
     spray axis=3 n=%d o=%g d=%g 
     ''' % (nt,it1,fw,nx,ox,delx,ny,oy,dely))

Flow('cmp','spikes Woft','inmo3gma velocity=${SOURCES[1]} | put label2="x" label3="y" ')
Flow('cmp2d','cmp','window n3=1 min3=0')
#Result('cmpcube','cmp',
#     '''
#     byte gainpanel=all |
#     grey3 point1=%g point2=%g title="d(t,x,y)" label2="x (km)" label3="y #(km)" flat=n 
#     frame1=%d frame2=75 frame3=75 screenratio=%g screenwd=%g screenht=%g
#     labelfat=4 titlefat=4
#     ''' % (p1,p2,itplot,sr,swc,sh),local=1)



###############################################################################################################


################################################################################################################
##### Measure Local Slopes of CMP

## Measure smoothed slope field.  Used for PNMO.
tsmooth=10
xsmooth=10

##  Smooth x-slope
Flow('pxsmooth','cmp2d',
     '''
     window squeeze=n |
     dip2 rect1=%d rect2=%d order=3
       ''' % (tsmooth,xsmooth))

Plot('cmp2d',
     '''
     grey title="d(t,x)" label2="x (km)" 
     screenratio=%g screenwd=%g screenht=%g
     labelfat=4 titlefat=4
     ''' % (sr,swc,sh),local=1)

## Inline(px) slopes
#Result('pxsmooth2d','pxsmooth',
#     '''
#     grey title="Slope (px)" label2="x (km)" color=j
#     screenratio=%g screenwd=%g screenht=%g
#     labelfat=4 titlefat=4
#     ''' % (sr,swc,sh),local=1)
Flow('timecube','cmp2d','math output="x1"')

################################################################################################################
#### Predictive flattening

# Time shifts via predictive painting
Flow('pick','pxsmooth',
     '''
     pwpaint i0=75 verb=y |
     clip2 lower=0 upper=4
     ''')

Flow('t0.asc',None,'echo %g n1=1 data_format=ascii_float in=$TARGET' % (it1))
Flow('t0','t0.asc','dd form=native | math output="(input-1)*0.004" ')
Plot('pick','pick t0',
     '''
     contour wanttitle=n wantaxis=n plotfat=5
     screenratio=%g screenwd=%g screenht=%g
     cfile=${SOURCES[1]} flat=y 
     ''' % (sr,swc,sh),local=1)
Result('pick2d','cmp2d pick','Overlay')


# Flattening via time warping both data and time cubes # T0(T) -> T(T0) ############################################################

Flow('flat','cmp2d pick','iwarp warp=${SOURCES[1]} eps=1')
Flow('warpedtime','timecube pick','iwarp warp=${SOURCES[1]} eps=1')
Flow('pickedtime','pick','iwarp warp=$SOURCE eps=1')

#Result('flat2d','flat',
#     '''
#     grey title="d(t\_0\^,x)" label2="x (km)"
#     screenratio=%g screenwd=%g screenht=%g
#     ''' % (sr,sw,sh),local=1)
#Result('warpedtime2d','warpedtime',
#     '''
#     grey title="t(t\_0\^,x,y)" label2="x (km)"
#     screenratio=%g screenwd=%g screenht=%g color=j
#     ''' % (sr,sw,sh),local=1)

# t^2(t0,x,y) - t0^2(t0,x,y) ###########################################################
Flow('moveout','timecube warpedtime','math time=${SOURCES[0]} warpt=${SOURCES[1]} output="warpt^2-time^2" | put d4="0" |clip2 lower=0')
#Result('moveout2d','moveout',
#     '''
#     window n3=1 min3=0|
#     grey title="t\^2\_-t\_0\^\^2\_" label2="x (km)"
#     screenratio=%g screenwd=%g screenht=%g color=j
#     ''' % (sr,sw,sh),local=1)

###########################################################################################################################################
###########################################   Effective W Inversion by warping

# Extract only the timewarp data at the event (using local max to distinguish)
Flow('count','cmp2d',
    '''
    window n2=1 f2=75 | envelope | 
    max1 | real | window n1=%d | put d1=1 o1=0 | sort ascmode=y
    ''' %(numtime))

# Compute the SD from flattened event 
Flow('est','flat',
    '''
    envelope | smooth rect1=5| max1 | real | window n1=%d | sort ascmode=y | put n1=%d n2=%d
    ''' %(numtime,nx,numtime))

# Compute the SD from pick traveltime accuracy (sd very small)
# Flow('countn.par','count','math output="input/0.004" | disfil number=n format="f1=%g"')
# Flow('est','pickedtime countn.par','window n1=1 par=${SOURCES[1]}')

Flow('wantedtime','count','spray axis=2 n=%d | sort ascmode=y | put n1=%d n2=%d' %(nx,nx,numtime))
Flow('sd','est wantedtime',
        '''
        math in1=${SOURCES[0]} in2=${SOURCES[1]} output="(in2-in1)^2/%d" |
        stack axis=1 norm=n | stack axis=1 norm=n | math output="sqrt(input)" | math output="1"
        ''' %(nx*ny-1))

# Compute limitx
Flow('limitx','count','math output="1.25" ') #offset/depth ratio ~ 1
Flow('maxx','count','math output="3" ')

# Interpolate in the timewarp volume
Flow('timecubepicked',['timecube','count'], 'put d4=0| inttest1 coord=${SOURCES[1]} interp=lag nw=2')
Flow('warpedtimepicked',['warpedtime','count'], 'inttest1 coord=${SOURCES[1]} interp=lag nw=2')

# Window transpose t,x,y -> x,y,t
Flow('timewarpsq','warpedtimepicked',' math output="input^2"| transp plane=13 | transp plane=12')
# Find t0^2
Flow('t0sqcube','timecubepicked',' math output="input^2" | transp plane=13 | transp plane=12 ')


# Range of parameters coefficients
def arr2str(array,sep=' '):
    return sep.join(map(str,array))

sW1=0.1;lW1=0.3;
sA1=-0.1;lA1=0.0;
sB1=0.5;lB1=1.0
sC1=0.0;lC1=0.006;
rangepar = (sW1,lW1,sA1,lA1,sB1,lB1,sC1,lC1)

Flow('range.asc',None,
     '''
         echo %s
     n1=%d o1=0 d1=1
     data_format=ascii_float in=$TARGET     
     ''' % (arr2str(rangepar),8) )
Flow('rangecoeff','range.asc','dd form=native')

# MCMC workflow ###################################################################
ntime=numtime
nmodel=20000
Flow('timewarpsqcut','timewarpsq','window n2=%d f3=0' %ntime)
Flow('t0sqcubecut','t0sqcube','window n2=%d f3=0' %ntime)
Flow('offsetx','cmp2d','math output=x2 | window n1=%d' %ntime)
Flow('rms','timewarpsqcut','stack axis=1 rms=y')
rms = 1.37393 # RMS of timewarpsqcut from sfattr

list = ['0.1','0.5','1','2','5','10','25','50']
for i in range(8):
        inum=float(list[i])
        istr=str(i)
        # First run for W
        Flow('histshort'+istr,['timewarpsqcut','t0sqcubecut','rangecoeff','limitx','offsetx','rms'],
                ''' 
                nmomcmc seed=62007 nmodel=%d saveiter=100 prior=n sdperc=%f datrms=${SOURCES[5]}
                t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
                ''' %(nmodel,inum))

        # Manually get mu and sigma for W from histshort
        Flow('mu'+istr,'histshort'+istr,'window n1=1 | stack axis=1')
        Flow('muext'+istr,'mu'+istr,'spray axis=1 n=%d' %nmodel)
        Flow('sigma'+istr,'histshort'+istr+' muext'+istr,
                '''
                window n1=1 | math mean=${SOURCES[1]} output="(input-mean)^2/%d" | 
                stack axis=1 norm=n | math output="sqrt(input)"
                ''' %(nmodel-1))
        Flow('rangecoeffcut'+istr,'rangecoeff','window n1=6 f1=2')
        Flow('rangecoeffnew'+istr,'mu'+istr+' sigma'+istr+' rangecoeffcut'+istr,'cat axis=1 ${SOURCES[1]} ${SOURCES[2]} ')

        # Second run for A B C
        Flow('hist'+istr,['timewarpsqcut','t0sqcubecut','rangecoeffnew'+istr,'maxx','offsetx','rms'],
                ''' 
                nmomcmc seed=62007 nmodel=%d saveiter=100 prior=n sdperc=%f wgauss=y datrms=${SOURCES[5]}
                t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
                ''' %(nmodel,inum))

        # Second run for A B C with enforce eta (only the range of A is used)
        Flow('histeta'+istr,['timewarpsqcut','t0sqcubecut','rangecoeffnew'+istr,'maxx','offsetx','rms'],
                ''' 
                nmomcmc seed=62007 nmodel=%d saveiter=100 prior=n sdperc=%f wgauss=y eta=y datrms=${SOURCES[5]}
                t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
                ''' %(nmodel,inum))


# Transdimensional workflow testing with artificially added noise ###################
# Test passed Hooray !!!
#####################################################################################
ssign=0.001; lsign=60;
trangeparnoise = (sW1,lW1,sA1,lA1,sB1,lB1,sC1,lC1,ssign,lsign)
Flow('trangenoise.asc',None,
     '''
         echo %s
     n1=%d o1=0 d1=1
     data_format=ascii_float in=$TARGET     
     ''' % (arr2str(trangeparnoise),10) )
Flow('trangenoisecoeff','trangenoise.asc','dd form=native')

prog = Program('Mnmomcmctrans.c')
sc = str(prog[0])

list = ['0.1','0.5','1','2','5','10','25','50']
for i in range(8):
        inum=float(list[i])
        istr=str(i)
        Flow('timewarpsqcutnoise'+istr,'timewarpsqcut',
                '''
                noise var=%g seed=2009 
                ''' %((inum/100*rms)*(inum/100*rms)))

        # First run for W
        Flow('thistshortnoise'+istr,['timewarpsqcutnoise'+istr,'t0sqcubecut','trangenoisecoeff','limitx','offsetx','rms',sc],
                ''' 
                ./${SOURCES[6]} seed=1990 nmodel=%d saveiter=500 prior=n datrms=${SOURCES[5]}
                t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
                ''' %(nmodel))
        #Result('thistnoiseW'+istr,'thistshortnoise'+istr,
        #               '''
        #       window n1=1 | histogram n1=101 o1=0.1 d1=0.002 |
        #               dd type=float | graph
        #       ''')
        #Result('thistnoisesig'+istr,'thistshortnoise'+istr,
        #       '''
        #       window n1=1 f1=-1 | histogram n1=101 o1=0 d1=0.1 |
        #       dd type=float | graph
        #       ''')            


# Transdimensional workflow ###################################################################
ssig=0.01; lsig=5;
trangepar = (sW1,lW1,sA1,lA1,sB1,lB1,sC1,lC1,ssig,lsig)
Flow('trange.asc',None,
     '''
         echo %s
     n1=%d o1=0 d1=1
     data_format=ascii_float in=$TARGET     
     ''' % (arr2str(trangepar),10) )
Flow('trangecoeff','trange.asc','dd form=native')

# First run for W
Flow('thistshort',['timewarpsqcut','t0sqcubecut','trangecoeff','limitx','offsetx','rms',sc],
        ''' 
        ./${SOURCES[6]} seed=62007 nmodel=%d saveiter=100 prior=n datrms=${SOURCES[5]}
        t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
        ''' %(nmodel))
# Manually get mu and sigma for W from histshort
Flow('tmu','thistshort','window n1=1 | stack axis=1')
Flow('tmuext','tmu','spray axis=1 n=%d' %nmodel)
Flow('tsigma','thistshort tmuext',
        '''
        window n1=1 | math mean=${SOURCES[1]} output="(input-mean)^2/%d" | 
        stack axis=1 norm=n | math output="sqrt(input)"
        ''' %(nmodel-1))
Flow('trangecoeffcut','trangecoeff tmu tsigma','window n1=8 f1=2')
Flow('trangecoeffnew','tmu tsigma trangecoeffcut','cat axis=1 ${SOURCES[1]} ${SOURCES[2]} ')

#Second run for A B C
Flow('thist',['timewarpsqcut','t0sqcubecut','trangecoeffnew','maxx','offsetx','rms',sc],
        ''' 
        ./${SOURCES[6]} seed=120 nmodel=%d saveiter=500 prior=n wgauss=y datrms=${SOURCES[5]}
        t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
        ''' %(nmodel))

#Second run for A B C with different seed
Flow('thist2',['timewarpsqcut','t0sqcubecut','trangecoeffnew','maxx','offsetx','rms',sc],
        ''' 
        ./${SOURCES[6]} seed=12000 nmodel=%d saveiter=500 prior=n wgauss=y datrms=${SOURCES[5]}
        t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
        ''' %(nmodel))

#Third run for A B C with enforce eta (only the range of A is used)
Flow('thisteta',['timewarpsqcut','t0sqcubecut','trangecoeffnew','maxx','offsetx','rms',sc],
        ''' 
        ./${SOURCES[6]} seed=120 nmodel=%d saveiter=500 prior=n wgauss=y datrms=${SOURCES[5]} eta=y
        t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
        ''' %(nmodel))





#Result('thistW','thist',
#       '''
#       window n1=1 | histogram n1=101 o1=0.1 d1=0.002 |
#       dd type=float | graph
#       ''')
#Result('thistA','thist',
#       '''
#       window n1=1 f1=1| histogram n1=101 o1=-0.1 d1=0.001 |
#       dd type=float | graph
#       ''')    
#Result('thistB','thist',
#       '''
#       window n1=1 f1=2| histogram n1=101 o1=0.5 d1=0.005 |
#       dd type=float | graph
#       ''')
#Result('thistC','thist',
#       '''
#       window n1=1 f1=3| histogram  n1=101 o1=0.0 d1=0.0001 |
#       dd type=float | graph
#       ''')    
#Result('thistsig','thist',
#       '''
#       window n1=1 f1=-1 | histogram n1=101 o1=-5 d1=0.1 |
#       dd type=float | graph
#       ''')    

#       #Second run for A B C with enforce eta (only the range of A is used)
#       Flow('thisteta',['timewarpsqcut','t0sqcubecut','trangecoeff','maxx','offsetx',sc],
#               ''' 
#               ./${SOURCES[5]} seed=62007 nmodel=%d saveiter=100 prior=n wgauss=n eta=y
#               t0sq=${SOURCES[1]} rangecoef=${SOURCES[2]} limitx=${SOURCES[3]} offsetx=${SOURCES[4]}
#               ''' %nmodel)




# Flow('tW','thist','window n1=1 f1=0')
# Flow('tA','thist','window n1=1 f1=1')
# Flow('tB','thist','window n1=1 f1=2')
# Flow('tC','thist','window n1=1 f1=3')
# 
# Flow('tWeta','thisteta','window n1=1 f1=0')
# Flow('tAeta','thisteta','window n1=1 f1=1')
# Flow('tBeta','thisteta','window n1=1 f1=2')
# Flow('tCeta','thisteta','window n1=1 f1=3')


# <histshort.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &
# <hist.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &
# <histeta.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &

# <histshort.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# <hist.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# <histeta.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# 
# <histshort.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# <hist.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# <histeta.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# 
# <histshort.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &
# <hist.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &
# <histeta.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &



# <thistshort.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &
# <thist.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &
# <thisteta.rsf sfwindow n1=1 | sfhistogram n1=101 o1=0.1 d1=0.002 | sfdd type=float | sfgraph | sfpen &

# <thistshort.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# <thist.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=1| sfhistogram n1=101 o1=-0.1 d1=0.001 | sfdd type=float | sfgraph | sfpen &
# 
# <thistshort.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# <thist.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=2 | sfhistogram n1=101 o1=0.5 d1=0.005 | sfdd type=float | sfgraph | sfpen &
# 
# <thistshort.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &
# <thist.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=3 | sfhistogram n1=101 o1=0.0 d1=0.0001 | sfdd type=float | sfgraph | sfpen &

# <thistshort.rsf sfwindow n1=1 f1=4 | sfhistogram n1=101 o1=0.0 d1=0.05 | sfdd type=float | sfgraph | sfpen &
# <thist.rsf sfwindow n1=1 f1=4 | sfhistogram n1=101 o1=0.0 d1=0.05 | sfdd type=float | sfgraph | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=4 | sfhistogram n1=101 o1=0.0 d1=0.05 | sfdd type=float | sfgraph | sfpen &



# <thist.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tW.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.1 d2=0.002 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thist.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tA.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=-0.1 d2=0.001 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thist.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tB.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.5 d2=0.005 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thist.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tC.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.0 d2=0.0001 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &


# <thisteta.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tWeta.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.1 d2=0.002 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tAeta.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=-0.1 d2=0.001 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tBeta.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.5 d2=0.005 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &
# <thisteta.rsf sfwindow n1=1 f1=4 | sfhist2 inp2=tCeta.rsf n1=101 o1=0.0 d1=0.05 n2=101 o2=0.0 d2=0.0001 | sfdd type=float | sfgrey color=j scalebar=y | sfpen &




End()