# Author: G. Hennenfent
#         Seismic Laboratory for Imaging and Modeling
#         Department of Earch & Ocean Sciences
#         The University of British Columbia
#         
# Date  : July, 07

from rsfproj import *

n1 = 2**13
d1 = 1./(n1-1)
clip = 1.8
lsz1 = 12
lsz2 = 12
ysc = 1.2

graphx = '''
sfput o1=-.5 d1=%(d1)f | sfwindow min1=0 |
sfgraph label1="Normalized wavenumber" label2="Amplitude" unit1="cycle" unit2=""
title="" labelsz=%(lsz1)g crowd2=.7
'''%vars()

# generate sparse DCT vector
Flow('x',None,
     'sfspike n1=%(n1)i n2=1 k1=750,250,5000 mag=20,50,7 nsp=3'%vars() )

Flow('spec specREG specIRREG model REG IRREG','x script.py',
     '''
     ./${SOURCES[1]} specREG=${TARGETS[1]} specIRREG=${TARGETS[2]}
     model=${TARGETS[3]} REG=${TARGETS[4]} IRREG=${TARGETS[5]}
     ''')
plots = '''
sfwindow f1=98 n1=50 | sfput d1=1 | sfdots dots=1 connect=0
radius=.1 clip=%(clip)f unit1="sample" labelsz=%(lsz2)g
yyscale=%(ysc)g
'''%vars()

Plot('REG',plots)
Plot('IRREG',plots)
Result('IRREGpres','IRREG',
       '''sfwindow f1=98 n1=50 |
       sfput d1=1 |
       sfdots dots=2 connect=0 radius=.1
       clip=%(clip)f
       labelsz=%(lsz2)g
       yyscale=%(ysc)g
       unit1="sample"
       '''%vars() )
Plot('mod','model',plots)
Plot('model',
     'sfwindow f1=98 n1=50 | sfput d1=1 | sfdots dots=0 strings=n radius=.1 clip=%(clip)f labelsz=%(lsz2)g yyscale=%(ysc)g unit1="sample"'%vars() )

Result('model','mod model','Overlay')
Result('REG','REG model','Overlay')
Result('IRREG','IRREG model','Overlay')

Result('spec',graphx)
Plot('specREG',graphx+' min2=0 max2=8')
Plot('specIRREG',graphx+' min2=0 max2=8')

Flow('thr',None,'sfmath n1=%(n1)i output=1'%vars() )
Plot('thr',graphx+' min2=0 max2=8 plotcol=5 dash=1')
Result('specIRREG','specIRREG thr','Overlay')
Result('specREG','specREG thr','Overlay')

Flow('maiIRREG','specIRREG spec','sfmath x=${SOURCES[1]} output="input-.34*x" ')
Result('maiIRREG',graphx+' min2=0 max2=8')

Flow('maiREG','specREG spec','sfmath x=${SOURCES[1]} output="input-.34*x" ')
Result('maiREG',graphx+' min2=0 max2=8')

End()