from rsf.proj import *
from math import *


#----------------------## parameters of synthetic signal

nt=501
dt = 0.008
nf = 257
ot=0

wf = 2*pi

def grey(title,bias):
    return '''
    transp | grey title="%s" bias=%g color=j screenratio=0.8 scalebar=y 
    barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
    label2=Time unit2=s label1=Frequency unit1=Hz 
    ''' % (title,bias)
###########################################################
# EXAMPLE 0000---cos+cos+cos+spike
###########################################################
nlevel0=0
nf0=247
df0=0.244141
Flow('s1-0',None,
     '''
     math n1=%d o1=0 d1=%g n2=1 output="cos(%g*(10*x1))+cos(%g*(20*x1))+cos(%g*(30*x1))" | 
     put label1=Time unit1=s label2=Amplitude 
     ''' % (nt,dt,wf,wf,wf))
Flow('add-0',None,'spike n1=%d o1=0 d1=%g n2=1 nsp=2 mag=10,10 k1=%d,%d' % (nt,dt,nt/2,nt/2+40))
Flow('s-0','s1-0 add-0','add ${SOURCES[1]} | noise seed=2012 var=%g' % (nlevel0))

Result('s-0','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n')

Flow('proj-0','s-0','timefreq rect=15 | window max2=60')

Result('proj-0',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y
       ''')
Flow('proj-0-30','s-0','timefreq rect=30 | window max2=60')

Result('proj-0-30',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y
       ''')


Flow('st-0','s-0','st  | math output="abs(input)" | real |  window max2=60')

Result('st-0',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y
       ''')



######################################
#EXAMPLE 1111 chirp signal: cos(g(f)*t+a)+cos(g(f)*t+b)
#####################################
nlevel1=0
nf1=257
df1=0.488281
dt1 = 0.004

Flow('s-1',None,
    '''
    math n1=%d o1=0 d1=%g n2=1 output="cos(%g*(2.5*x1*x1+8)*x1)+cos(%g*(5*x1*x1+25)*x1)"
    ''' % (nt, dt1, wf, wf))
Result('s-1','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n ')

Flow('proj-1','s-1','timefreq rect=17')

Result('proj-1',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=30 grid=y
       ''')


Flow('st-1','s-1','st  | math output="abs(input)" | real')

Result('st-1',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=30 grid=y
       ''')



















##########################################################
#EXAMPLE 2222: synthetic trace by ricker wavelet
##########################################################
nlevel2=0
dt2=0.002
nf2=500


################
Flow('s2-1',None,
     '''
     spike n1=%d d1=%g o1=0 nsp=1 k1=%g mag=1 | ricker1 frequency=60 | window max1=1.2
     ''' % (1000,dt2,0.602*nt))
Flow('s2-4',None,
     '''
     spike n1=%d d1=%g o1=0 nsp=1 k1=%g mag=1 | ricker1 frequency=15 | window max1=1.2
     ''' % (1000,dt2,0.602*nt))
Flow('s2-2',None,
     '''
     spike n1=%d d1=%g o1=0 nsp=2 k1=%g,%g mag=1,1 | ricker1 frequency=50 | window max1=1.2
     ''' % (1000,dt2,1.02*nt,1.045*nt))
Flow('s2-3',None,
     '''
     spike n1=%d d1=%g o1=0 nsp=1 k1=%g mag=1 | ricker1 frequency=30 | window max1=1.2
     ''' % (1000,dt2,0.202*nt))


Flow('s-2','s2-1 s2-2 s2-3 s2-4','add ${SOURCES[1]} ${SOURCES[2]} ${SOURCES[3]}| noise seed=12012 var=5e-5 | math output="input/2" ')



Result('s-2','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n pad=n')





Flow('proj-2','s-2','timefreq  rect=6 niter=100 | smooth rect2=5 rect1=2')

Result('proj-2',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz  n2tic=20 grid=y g2num=50 max1=150
       ''')





#---------S Tran

Flow('st-2','s-2','st  | math output="abs(input)" | real ')

Result('st-2',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y g2num=50 max1=150
       ''')




#######################################################
## EXAMPLE 66666666: ricker with frequency f=at+b
########################################################



Flow('rickers',None,
     '''
     math n1=2001 n2=2001 o1=0 o2=0 d1=0.001 d2=0.001 
     output="(1-2*(3.1415926*(25*x2*x2+15)*(x1-1))*(3.1415926*(25*x2*x2+15)*(x1-1)))*exp(-(3.1415926*(25*x2*x2+15)*(x1-1))*(3.1415926*(25*x2*x2+15)*(x1-1)))" | scale axis=1 
     ''')


Flow('spec-r','rickers','spectra')


Plot('spec-r',
       '''
       window max1=200 | scale axis=1 |
       grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y g2num=50
       ''')



Flow('num-r','spec-r',' window  | math output="input*x1" | stack axis=1')
Flow('den-r','spec-r','window   | stack axis=1')
Flow('lcf-r','num-r den-r','divn rect1=1 den=${SOURCES[1]}')





nf6=600
nt6=2001
dt6=0.001
nlevel6=0
Flow('ref-6',None,
     '''
     math n1=%d d1=%g o1=%g output=0 | noise seed=22012 var=1 
     ''' % (nt6,dt6,ot))
Flow('freq-6','ref-6','math output="25*x1*x1+15"')
Flow('phase-6','ref-6','math output=0')
Flow('s-6','ref-6 freq-6 phase-6',
     '''
     ricker2 tfreq=${SOURCES[1]} tphase=${SOURCES[2]} norm=n hiborder=100  | 
     noise seed=32012 var=%g  | agc rect1=100 
     ''' % nlevel6)
Result('ref-6','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n')
Result('s-6','graph pad=n screenratio=0.8 crowd1=0.75 label2=Amplitude label1=Time unit1=s crowd2=0.3 wanttitle=n')



Flow('proj-6','s-6','timefreq rect=30 | window max2=201')

Result('proj-6',
       '''
       transp | scale axis=2 | grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y g2num=50 max1=200
       ''')


#---test local frequency 
Flow('numer-6','proj-6',' math output="input*x2" | stack axis=2')
Flow('denom-6','proj-6','stack axis=2')

Flow('lcf-6','numer-6 denom-6','math x=${SOURCES[1]} output="input/(x+0.01)"')
Flow('theo-6','lcf-6','math output="25*x1*x1+15" ')


Flow('flcf-6','s-6','iphase rect1=90 order=10 hertz=y')


######################

Flow('st-6','s-6','st  | math output="abs(input)" | real | window max2=201')

Result('st-6',
       '''
       transp | scale axis=2|
       grey wanttitle=n color=j screenratio=0.8 scalebar=y bartype=v
       barwidth=0.2 crowd1=0.75  crowd2=0.3 wherexlabel=b allpos=y
       label2=Time unit2=s label1=Frequency unit1=Hz n2tic=20 grid=y g2num=50 
       ''')


Flow('numer-s-6','st-6',' window max2=150 | math output="input*x2" | stack axis=2')
Flow('denom-s-6','st-6','window max2=150 | stack axis=2')
Flow('lcf-s-6','numer-s-6 denom-s-6','divn rect1=1 den=${SOURCES[1]}')
Result('lcf-6','lcf-6 lcf-s-6 lcf-r',
       '''
        cat axis=2 ${SOURCES[1]} ${SOURCES[2]}  |
        graph pad=n screenratio=0.8 yreverse=y min2=0 max2=200 n2tic=30 
        crowd1=0.75 label2=Frequency unit2=Hz label1=Time unit1=s crowd2=0.3 wanttitle=n 
        n2tic=20 grid=y g2num=50 plotcol=2,4,7 plotfat=3 dash=0,5,7 
       ''')



############ test for dominant Freq and Ava

Plot('lcf-r',
       '''
       graph wanttitle=n wantaxis=n screenratio=0.8 scalebar=y yreverse=y pad=n
       barwidth=0.2 crowd1=0.75  crowd2=0.3  n2tic=20 grid=y g2num=50 max2=200 min2=0 plotfat=2 plotcol=7 
       ''')
Plot('theo-6',
       '''
       graph wanttitle=n wantaxis=n screenratio=0.8 scalebar=y yreverse=y pad=n
       barwidth=0.2 crowd1=0.75  crowd2=0.3  n2tic=20 grid=y g2num=50 max2=200 min2=0 plotfat=2 plotcol=0 
       ''')

Result('spec','spec-r lcf-r theo-6','Overlay')




End()