from rsf.proj import *
import math

## CHIRP SIGNAL
###############

Flow('chirp',None,
     'math  n1=1024 d1=0.004 output="cos(%g*x1*(1+0.5*x1))" ' % (math.pi*10))
Flow('chirp2',None,
     '''
     math  n1=1424 d1=0.004 o1=-0.8 output="cos(%g*x1*(1+0.5*x1))"
     ''' % (math.pi*10))

Plot('chirp',
     '''
     graph title="(a) Chirp Signal" label1="Time (s)" label2=Amplitude
     labelsz=10 titlesz=12
     ''')

for case in (1,2):
    ic = 'ic%d' % case
    Plot(ic,'chirp2',
         '''     
         iphase order=200 rect1=%d hertz=y |
         window min1=0 n1=1024 |
         graph title="(a) Chirp %s Frequency" min2=0 max2=30 plotcol=5
         label1="Time (s)" label2="Frequency (Hz)"
         labelsz=10 titlesz=12
         ''' % ((1,"Instantaneous"),(10,"Local"))[case-1])
    Result(ic,['chirp',ic],'OverUnderAniso')

## SYNTHETIC SIGNAL
###################

Fetch('ai1.rsf','attr')

Flow('ai','ai1.rsf','dd form=native')

Flow('refl','ai1','ai2refl')
Flow('ricker',None,'spike n1=250 d1=1e-04 k1=125 | ricker1')
Flow('conv','refl','ricker1')

plots = {'ai1':'Acoustic Impedance',
         'refl':'Reflectivity',
         'ricker': 'Wavelet',
         'conv': '(b) Synthetic Trace'}
for plot in plots.keys():
    Plot(plot,
         '''
         graph title="%s" label1=Time unit1=s label2=Amplitude
         labelsz=10 titlesz=12
         ''' % plots[plot])
    Plot('spec'+plot,plot,
         '''
         spectra |
         graph title="%s Spectrum" wantaxis1=n
         labelsz=10 titlesz=12
         ''' % plots[plot])
Result('inv1','ai1 refl ricker conv','TwoRows')
Result('inv2','specai1 specrefl specricker specconv','TwoRows')

for case in (1,2):
    sf = 'sf%d' % case
    Flow(sf,'conv','iphase order=10 rect1=%d hertz=y complex=%d' % ((1,500)[case-1],case-1))
    Plot(sf,
         '''
         scale dscale=0.05 |
         graph title="(b) Synthetic %s Frequency"
         label1=Time unit1=s min2=-5 max2=100
         label2="Frequency (Hz)" plotcol=5
         labelsz=10 titlesz=12
         ''' % ('Instantaneous','Local')[case-1])
    Result(sf,['conv',sf],'OverUnderAniso')


## REAL SIGNAL
##############

Fetch('pptrace.rsf','attr')
Flow('pptrace1','pptrace.rsf','dd form=native | window max1=5')

Plot('seis','pptrace1',
     '''
     graph title="(c) Seismic Trace" label1="Time (s)" label2=Amplitude
     labelsz=10 titlesz=12 
     ''')

for case in (1,2):
    seis = 'seis%d' % case
    Plot(seis,'pptrace1',
         '''
         iphase rect1=%d order=10 hertz=y complex=%d |
         graph title="(c) Seismic %s Frequency"
         min2=-20 max2=100 plotcol=5
         label1="Time (s)" label2="Frequency (Hz)"
         labelsz=10 titlesz=12
         ''' % ((1,0,"Instantaneous"),(60,1,"Local"))[case-1])
    Result(seis,['seis',seis],'OverUnderAniso')

## COLLECTING TOGETHER
######################

Result('sign','chirp conv seis','OverUnderAniso',vppen='txscale=1.2')
Result('inst','ic1 sf1 seis1','OverUnderAniso',vppen='txscale=1.2')
Result('locl','ic2 sf2 seis2','OverUnderAniso',vppen='txscale=1.2')

End()