from rsf.proj import *

# Download input trace

trace = 'trace_il1190_xl1155.trace'

Fetch(trace,'data',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master/1410_Phase')

# Convert format to RSF

Flow('trace',trace,
     '''
     echo n1=1503 in=$SOURCE data_format=ascii_float 
     o1=0 d1=0.004 label1=Time unit1=s |
     dd form=native | window f1=2
     ''')

# Apply Hilbert transform

Flow('hilb','trace','envelope hilb=y')

# Compute envelope and phase

Flow('envp','trace hilb','math h=${SOURCES[1]} output="sqrt(input*input+h*h)" ')
Flow('envm','envp','scale dscale=-1')

Flow('phase','trace hilb','math h=${SOURCES[1]} output="h&input" ')

Plot('trace','trace hilb envp envm',
     '''
     cat ${SOURCES[1:4]} axis=2 |
     window f1=500 n1=201 | 
     graph title="Trace Amplitudes" dash=0,1,0,0 plotcol=5,5,4,4
     ''')

Plot('phase',
     '''
     window f1=500 n1=201 | 
     graph title="Instantaneous Phase" plotcol=3 label2=Degree unit2=radians
     ''')

Result('trace','trace phase','OverUnderAniso')

# Repeat on 2-D section

sgy = 'penobscot_xl1155.sgy'

Fetch(sgy,'data',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master/1410_Phase')

Flow('data tdata',sgy,'segyread tfile=${TARGETS[1]}')
Flow('hdata','data','envelope hilb=y')
Flow('edata','data hdata','math h=${SOURCES[1]} output="sqrt(input*input+h*h)" ')
Flow('pdata','data hdata','math h=${SOURCES[1]} output="h&input" ')

Plot('data','window min1=2 max1=3 | grey title="Original" scalebar=y')
Plot('hdata','window min1=2 max1=3 | grey title="Hilbert" scalebar=y')
Plot('edata','window min1=2 max1=3 | grey title="Envelope" scalebar=y color=j allpos=y')
Plot('pdata','window min1=2 max1=3 | grey title="Phase" scalebar=y color=j')

Result('data','data hdata edata pdata','TwoRows')

# Plot envelope peaks

Plot('env','trace envp envm',
     '''
     cat ${SOURCES[1:3]} axis=2 |
     window f1=500 n1=201 | 
     graph title="Envelope" plotcol=5,4,4 grid=y 
     min1=2 max1=2.8 min2=-7200 max2=7200
     ''')

Flow('peaks','envp','max1')
Flow('meaks','peaks','math output="conj(input)" ')

Plot('peaks','peaks meaks',
     '''
     cat axis=1 ${SOURCES[1]} | 
     graph wanttitle=n wantaxis=n symbol=o plotcol=3 symbolsz=6
     min1=2 max1=2.8 min2=-7200 max2=7200
     ''')

Plot('peaks2','peaks meaks',
     '''
     cat axis=2 ${SOURCES[1]} | transp | 
     graph wanttitle=n wantaxis=n plotcol=3 
     min1=2 max1=2.8 min2=-7200 max2=7200
     ''')

Plot('envpeaks','env peaks peaks2','Overlay')

Flow('xpeaks','peaks','real')
Flow('phasepeaks','phase xpeaks','inttest1 coord=${SOURCES[1]} interp=spline nw=4')
Flow('cpeaks','xpeaks','rtoc')

Plot('ppeaks','xpeaks phasepeaks cpeaks',
     '''
     cmplx ${SOURCES[1]} | cat axis=2 ${SOURCES[2]} | transp |
     graph wanttitle=n wantaxis=n plotcol=3 
     min1=2 max1=2.8 min2=-3.5 max2=3.5
     ''')
Plot('ppeaks2','xpeaks phasepeaks',
     '''
     cmplx ${SOURCES[1]} |
     graph title=Phase symbol=o plotcol=3 symbolsz=6 grid=y
     min1=2 max1=2.8 min2=-3.5 max2=3.5
     ''')

# ideal phase

from math import pi

Flow('iphase0','phasepeaks',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="(1-input)*p" 
     ''' % (-pi/2,pi/2))
Flow('iphase1','iphase0',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="%g*input+(1-input)*p"
     ''' % (-3*pi/2,-pi/2,-pi))
Flow('iphase2','iphase1',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="%g*input+(1-input)*p"
     ''' % (pi/2,3*pi/2,pi))

Plot('ipeaks','xpeaks iphase2 cpeaks',
     '''
     cmplx ${SOURCES[1]} | cat axis=2 ${SOURCES[2]} | transp |
     graph wanttitle=n wantaxis=n plotcol=2 plotfat=3
     min1=2 max1=2.8 min2=-3.5 max2=3.5
     ''')
Plot('ipeaks2','xpeaks iphase2',
     '''
     cmplx ${SOURCES[1]} |
     graph title=Phase symbol=o plotcol=2 symbolsz=8 grid=y
     min1=2 max1=2.8 min2=-3.5 max2=3.5 plotfat=3
     ''')

Plot('phasepeaks','ipeaks ipeaks2 ppeaks ppeaks2','Overlay')

Result('phasepeaks','envpeaks phasepeaks','OverUnderAniso')

# Compute phase error for seismic section

Flow('dpeaks','edata','max1 | real')
Flow('ppeaks','pdata dpeaks','inttest1 coord=${SOURCES[1]} interp=spline nw=4 same=n')

Flow('ipeaks0','ppeaks',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="(1-input)*p" 
     ''' % (-pi/2,pi/2))
Flow('ipeaks1','ipeaks0',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="%g*input+(1-input)*p"
     ''' % (-3*pi/2,-pi/2,-pi))
Flow('ipeaks2','ipeaks1',
     '''
     mask min=%g max=%g | dd type=float | math p=$SOURCE output="%g*input+(1-input)*p"
     ''' % (pi/2,3*pi/2,pi))

Flow('epeaks bar','ipeaks2 ppeaks',
     'math p=${SOURCES[1]} output="abs(p-input)" | byte allpos=y bar=${TARGETS[1]} pclip=100')

Result('epeaks','dpeaks epeaks',
       '''
       rtoc | math output="input+I*x2" |
       graph symbol='.' min1=2 max1=3 min2=0 max2=600 
       title="Absolute Phase Error at Envelope Peaks" 
       transp=y yreverse=y depth=${SOURCES[1]} 
       color=j scalebar=y 
       ''')

# Zero-phase correction?

Flow('skew','trace',
     'pad beg1=250 end1=250 | localskew rect=250 a0=-90 da=3 na=121 | window n1=1501 f1=250')

Plot('skew',
     '''
     grey color=j allpos=y title="Inverse Local Skewness" 
     label2="Phase Rotation" Xunit2="\^o\_" 
     ''')

# Pick & Slice

Flow('pick','skew','pick vel0=-45 rect1=50')

Plot('pick','graph plotcol=7 plotfat=5 transp=y yreverse=y min2=-90 max2=270 wanttitle=n pad=n wantaxis=n')

Result('skew','skew pick','Overlay')

# Phase rotation

Flow('trace90','trace pick','phaserot phase=${SOURCES[1]}')
Flow('trace0','trace90','envelope order=100 hilb=y phase=-90')

Result('trace2','trace trace0',
       '''
       cat axis=2 ${SOURCES[1]} |
       dots labels="Input:Output" yreverse=y
       ''')


End()