from rsf.proj import *

# download well-log data
Fetch('L-30.las','1406_Make_a_synthetic',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master')

# Convert to RSF
Flow('L-30','L-30.las','las2rsf $SOURCE $TARGET',stdin=0,stdout=-1)

# Examine with "< L-30.rsf sfheaderattr segy=n desc=y"

# Mask bad values
Flow('mask1','L-30','headermath segy=n output=RHOB | mask min=0')
Flow('mask2','L-30','headermath segy=n output=DT   | mask min=0')
Flow('mask','mask1 mask2','mul ${SOURCES[1]}')

scale = dict(RHOB=1000,DT=3.28084,DEPTH=1/3.28084) # convert to SI units
threshold = dict(RHOB=100,DT=10)

for case in ('DT','RHOB','DEPTH'):
     # Extract values
     Flow(case,'L-30 mask',
          '''
          headermath segy=n output=%s | 
          headerwindow mask=${SOURCES[1]} | window |
          scale dscale=%g | 
          put o1=3058.5 d1=0.5 label1=Depth unit1=ft
          ''' % (case,scale[case]))

     # Clip spikes
     if case != 'DEPTH':
          despike = case+'-despike'
          Flow(despike,case,'despike wide=13')

          spikemask = case+'-spikemask'
          Flow(spikemask,[despike,case],
               '''
               math orig=${SOURCES[1]} output="abs(input-orig)" |
               mask min=%g
               ''' % threshold[case])

          Flow(despike+'2',[spikemask,despike,case],
               '''
               dd type=float |
               math orig=${SOURCES[2]} filt=${SOURCES[1]}
               output="input*filt+(1-input)*orig"
               ''')
          Result(case,[case,despike+'2'],
                 '''
                 cat axis=2 ${SOURCES[1]} |
                 graph title=%s
                 ''' % case)

# Compute time-to-depth relationship

log_start_time = 0.410553230106 # (s)

# windowed samples - 3837

Flow('dt1','L-30',
     '''
     headermath segy=n output=DT | 
     window n2=3837 | scale dscale=%g 
     ''' % scale['DT'])
Flow('dt2','dt1','window f1=21')
Flow('dt0','dt1 dt2','window n1=21 | math output=507.2 | cat axis=1 ${SOURCES[1]}')
Flow('deltat','dt0','stack axis=1 norm=n | scale dscale=%g' % (2*0.1524/1e6))

delta_time = 0.5596 # (s)

Flow('tdr','DT-despike','causint | math output="%g+(%g)*input" ' % (log_start_time+delta_time,2*0.1524/1e6))

Result('tdr','graph title="Twoway Time" yreverse=y transp=y label2=Time unit2=s')

# Compute acoustic impedance

Flow('z','RHOB slowness','math rho=${SOURCES[0]} slow=${SOURCES[1]} output="1e6*rho/slow" ')

# Compute reflection coefficient series

Flow('rc','z','ai2refl')

# Plots

Flow('dz','DEPTH','igrad')
Flow('slowness','DT-despike dz','div ${SOURCES[1]} | clip2 lower=300')

Plot('slowness','DEPTH slowness',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=function unit1=m label2="P-wave slowness" unit2="\F10 m\F3 m/s" labelsz=12 titlesz=15
     ''')

Plot('twt','DEPTH tdr',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=integral unit1=m label2="two-way time" unit2="s" labelsz=12 titlesz=15
     ''')

Plot('impedance','DEPTH z',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=impedance unit1=m label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')

Plot('reflectivity','DEPTH rc',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=reflectivity unit1=m label2= unit2= labelsz=12 titlesz=15
     ''')

Plot('four','slowness twt impedance reflectivity','SideBySideAniso')

# add labels for well tops

welltops = {
    'Abenaki': 3404.311,
    'Base O-Marker': 2469.207,
    'Dawson Canyon': 984.504,
    'L Baccaro': 3964.534,
    'L Missisauga': 3190.646,
    'Logan Canyon': 1136.904,
    'Mid Baccaro': 3485.083,
    'U Missisauga': 2251.253,
    'Wyandot': 867.156
    }

a = (8.95-1.3)/(750-4250)
b = 1.3-4250*a
plots=[]
for top in welltops.keys():
    name = top.replace(' ','-').lower()
    Plot(name,None,'box font=2 x0=12.5 y0=%g label="%s" size=0.15 xt=0 yt=0' % (a*welltops[top]+b,top))
    plots.append(name)
Plot('welltops',plots,'Overlay')

Result('four','four welltops','Overlay')

Flow('tops.txt',None,
     'echo %s in=$TARGET n1=%d data_format=ascii_float' % (' '.join(map(str,welltops.values())),len(welltops)))
Flow('tops','tops.txt','dd form=native | sort')
Plot('tops','spray axis=1 n=2 o=0 d=1 | graph min2=750 max2=4250 pad=n yreverse=y wantaxis=n wanttitle=n')

Flow('ttops','tdr DEPTH tops',
     'iwarp warp=${SOURCES[1]} o1=932 d1=1 n1=3307 | inttest1 nw=4 interp=spline coord=${SOURCES[2]}')
Plot('ttops','spray axis=1 n=2 o=0 d=1 | graph min2=0.9 max2=2.9 pad=n yreverse=y wantaxis=n wanttitle=n')

# Convert to two-way-traveltime

Flow('zt','z tdr','cut f1=-20 | iwarp warp=${SOURCES[1]} o1=0 d1=0.004 n1=750')
Flow('rt','rc tdr','cut f1=-20 | iwarp warp=${SOURCES[1]} o1=0 d1=0.004 n1=750')

Plot('z','DEPTH z',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 plotfat=3 plotcol=7
     title=impedance unit1=m label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')
Plot('zw','z tops','Overlay')

Plot('zt',
     '''
     window min1=0.9 max1=2.9 |
     graph transp=y yreverse=y plotfat=3 plotcol=7
     title=impedance label1=Time unit1=s pad=n 
     label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')
Plot('ztw','zt ttops','Overlay')

Plot('rt',
     '''
     window min1=0.9 max1=2.9 |
     graph transp=y yreverse=y plotfat=3 plotcol=7
     title=reflectivity label1=Time unit1=s pad=n 
     label2= unit2= labelsz=12 titlesz=15
     ''')
Plot('rtw','rt ttops','Overlay')

# Create synthetic

Flow('synth','rt','ricker1 frequency=25')
Plot('synth',
     '''
     window min1=0.9 max1=2.9 | 
     wiggle poly=y yreverse=y title=synthetic transp=y yreverse=y
     label1=Time unit1=s labelsz=12 titlesz=15 clip=0.002 grid=n plotcol=7
     ''')
Plot('synthw','synth ttops','Overlay')

Result('synth','zw ztw rtw synthw','SideBySideAniso')

# Get seismic data

Fetch('PenobXL_1155.txt','1406_Make_a_synthetic',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master')

Flow('seismic','PenobXL_1155.txt',
     '''
     echo in=$SOURCE data_format=ascii_float n1=1003 n2=401 
     label2=Trace label1=Time unit1=s o1=0 d1=0.004 |
     dd form=native | reverse which=2 | window n1=750
     ''',stdin=0)

Plot('seismic','grey title=Seismic')
Flow('xpos',None,'spike n1=1 mag=77')
Plot('trace','synth xpos',
     '''
     wiggle transp=y yreverse=y wanttitle=n wantaxis=n plotcol=4
     xpos=${SOURCES[1]} xmin=400 xmax=0 poly=y zplot=10
     ''')
Result('seismic','seismic trace','Overlay')

# improve registration

End()