from rsf.proj import *
from math import *
from rsf.recipes.beg import server as private

#-----read data from segy file
Fetch('ln472_old.sgy','guochang',private)

Flow('old','ln472_old.sgy',
     '''
     segyread bfile=/dev/null hfile=/dev/null read=data |  
     put n1=751 n2=201 o1=0 d1=0.002 02=1 d2=0.01 label1="Time" unit1="s" 
     label2="X" unit2="km" | window min1=0
     ''')

ns=80


def ref(trace):
    out = 'ref%d' % trace
    Flow(out+'.asc',None,
         '''
         echo %d 0 %d 1 n1=4 in=$TARGET data_format=ascii_float
         ''' % (trace,trace))
    Plot(out,out+'.asc',
         '''
         dd form=native type=complex | 
         graph min1=0 max1=201 min2=0 title="" wantaxis=n scalebar=y pad=n plotfat=4 dash=2
         ''')
    return out



Plot('old','old','grey title=" " scalebar=y label2="Distance" label1=Time unit1=s ')
Result('old-1',['old',ref(ns)],'Overlay')


Result('spectra','old','spectra all=y | scale axis=1 | graph pad=n label2="Scaled amplitude" wanttitle=n unit2="" plotfat=4')

############################### LTFT

Flow('ltft','old',
     '''
     timefreq rect=8 verb=n nw=377 dw=0.665779 niter=100
     ''')

Result('tf-l','ltft',
       '''
       window n3=1 f3=%d max2=180 | smooth rect2=5 | scale axis=2 | 
       grey color=j allpos=y scalebar=y wanttitle=n screenratio=1.5 
       ''' % ns)
### Freq slice
Result('l1','ltft',
       '''
       window n2=1 min2=15 | scale axis=2 | 
       sfgrey color=j allpos=y scalebar=y title="15HZ" label2="Distance" unit2=km clip=1

       ''')

Result('l2','ltft',
       '''
       window n2=1 min2=31 | scale axis=2 | 
       sfgrey color=j allpos=y scalebar=y title="31HZ" label2="Distance" unit2=km clip=1

       ''')

Result('l3','ltft',
       '''
       window n2=1 min2=70 | scale axis=2 | 
       sfgrey color=j allpos=y scalebar=y title="70HZ" label2="Distance" unit2=km clip=1

       ''')

# Average frequency
Flow('num-l','ltft','window max2=110 | math output="input*x2" | stack axis=2 ')

Flow('num-l-1','num-l','window f1=200 | smooth rect1=1 rect2=1')
Flow('num-l-2','num-l','window n1=200 ')
Flow('num-l-a','num-l-2 num-l-1','cat ${SOURCES[1]} axis=1 o=0 d=0.002')

Flow('den-l','ltft','window max2=110 | stack axis=2')
Flow('lcf-l','num-l den-l','divn den=${SOURCES[1]} rect1=1  rect2=1')

Result('lcf-l','lcf-l',
       '''
        grey allpos=y scalebar=y  clip=86  barunit=Hz  color=j title="" label2="Distance" label1=Time unit1=s
       ''')




############## ST
Flow('st','old','st  | math output="abs(input)" | real ')
Result('tf-s','st',
       '''
       window n3=1 f3=%d max2=180| scale axis=2 | 
       grey color=j allpos=y scalebar=y wanttitle=n screenratio=1.5
       ''' % ns)

# Freq Slice
Result('s1','st',
       '''
       window n2=1 min2=15 | scale axis=2 | 
       sfgrey color=j allpos=y scalebar=y title="15HZ" label2="Distance" unit2=km clip=1

       ''')

Result('s2','st',
       '''
       window n2=1 min2=31 | scale axis=2 | 
       sfgrey color=j allpos=y scalebar=y title="31HZ" label2="Distance" unit2=km clip=1

       ''')

Result('s3','st',
       '''
       window n2=1 min2=70 | scale axis=2 | 
       grey color=j allpos=y scalebar=y title="70HZ" label2="Distance" unit2=km clip=1
       ''')




End()