from rsf.proj import *
from math import *
import math

## --------------------- model with constant Q -----------------------
## synthetic data with Q attenuation
Flow('asignal',None,
     ''' 
     modtraceq n1=1001 d1=0.001 o1=0. nc=6 fm=50
     at=100,200,400,600,700,900 q=99999,60,60,60,60,60
     ''')
Result('asignal',
       '''
       graph pad=n screenratio=0.8 crowd1=0.3 crowd2=0.75 
       label2=Amplitude unit2= label1=Time unit1=s title= 
       yreverse=y transp=y wherexlabel=t parallel2=n n2tic=20
       labelfat=3 plotfat=4 plotcol=7 font=2 axisfat=3 
       ''')

## --------------------- estimate Q using LTFT -----------------------
## LTFT
Flow('attf','asignal','ltft rect=9')
Flow('mattf','attf','math output="abs(input)" |real|window min1=0.')
Plot('mattf',
     '''
     scale axis=2|     
     grey wanttitle=n max2=150 min2=0 color=j screenratio=0.8 scalebar=n
     barwidth=0.2 crowd1=0.3 crowd2=0.75 wherexlabel=t allpos=y clip=0.8
     n1tic=20 grid=y g1num=50 label2="Frequency" label1="Time" unit1=s
     labelfat=3 font=2 axisfat=3 parallel2=n n2tic=20
     ''')

## determine bandwidth using peak frequency
dw = 0.976562    ## need to be modified
Flow('pf','mattf',
     '''
     transp | findmax1 verb=n | dd type=float | smooth rect1=1 |
     math output="(input-1)*%g*2.5" | dd type=float | window min1=0. | put o1=0.
     ''' % dw)

## local centroid frequency
Flow('bw cf var2','attf pf',
     '''          
     window max2=150 | 
     lcf trange=${SOURCES[1]} avef=${TARGETS[1]} var2=${TARGETS[2]} rect1=60
     ''')
Plot('cf',
       '''
       graph max2=150 min2=0 pad=n screenratio=0.8 crowd1=0.3
       crowd2=0.75 wanttitle=n yreverse=y transp=y wherexlabel=t dash=0
       label2="" unit2="" g1num=50 plotfat=6 grid=n
       labelfat=0 plotcol=7 font=2 parallel2=n n2tic=20
       ''')
Result('fp',' mattf cf','Overlay')

## equivalent Q value
Flow('repos',None,'math n1=1 o1=0. d1=1 output="100"|dd type=int')
Flow('qts','cf var2 repos',
     ''' 
     lcfseq var2=${SOURCES[1]} repos=${SOURCES[2]}
     ''')

## true equivalent Q value
Flow('trueq','qts repos','convert0eq repos=${SOURCES[1]}')

## theoretical Q-value
Flow('tq teq',None,
     ''' 
     theoreqiq teq=${TARGETS[1]} n1=1001 d1=0.001 o1=0. nc=6 
     at=100,200,400,600,700,900 q=99999,60,60,60,60,60 |
     dd type=float
     ''')

## Q estimation using the LCFS method and theoretical Q
Result('qq','tq qts',
       '''
       cat axis=2 ${SOURCES[1:2]}  |
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=150 o1num=0 d1num=50 n1tic=3
       ''')

## inverse Q filtering
Flow('asf','asignal','fft1 | spray axis=2 d=0.001 n=1001 | transp')
Flow('cmtf','asf trueq',
     '''
     invqfilt eqt=${SOURCES[1]} gim=25
     ''')
Flow('cmsig','cmtf','transp | fft1 inv=y | window n1=1 f1=0')
Result('cmsig',
       '''
       graph pad=n screenratio=0.8 crowd1=0.3 crowd2=0.75
       label2=Amplitude unit2= label1=Time unit1=s 
       title= yreverse=y transp=y wherexlabel=t plotfat=4
       labelfat=3 plotcol=7 font=2 axisfat=3 parallel2=n n2tic=20 
       ''')

## --------------------- estimate Q using ST -----------------------
## ST
Flow('attf1','asignal','st')
Flow('mattf1','attf1','math output="abs(input)" |real|window min1=0.')
Plot('mattf1',
     '''
     scale axis=2|     
     grey wanttitle=n max2=150 min2=0 color=j screenratio=0.8 scalebar=n
     barwidth=0.2 crowd1=0.3 crowd2=0.75 wherexlabel=t allpos=y clip=0.8
     n1tic=20 grid=y g1num=50 label2="Frequency" label1="Time" unit1=s
     labelfat=3 font=2 axisfat=3 parallel2=n n2tic=20
     ''')

## determining bandwidth using peak frequency
dw1 = 0.999001    ## need to be modified
Flow('pf1','mattf1',
     '''
     transp | findmax1 verb=n | dd type=float | smooth rect1=1 |
     math output="(input-1)*%g*2.5"|dd type=float|window min1=0.|put o1=0.
     ''' % dw1)

## local centroid frequency
Flow('bw1 cf1 var21','attf1 pf1',
     '''          
     window max2=150 | 
     lcf trange=${SOURCES[1]} avef=${TARGETS[1]} var2=${TARGETS[2]} rect1=65
     ''')
Plot('cf1',
       '''
       graph max2=150 min2=0 pad=n screenratio=0.8 crowd1=0.3
       crowd2=0.75 wanttitle=n yreverse=y transp=y wherexlabel=t dash=0
       label2="" unit2="" g1num=50 plotfat=6 grid=n
       labelfat=0 plotcol=7 font=2 parallel2=n n2tic=20
       ''')
Result('fp1',' mattf1 cf1','Overlay')

## equivalent Q value
Flow('qts1','cf1 var21 repos',
     ''' 
     lcfseq var2=${SOURCES[1]} repos=${SOURCES[2]}
     ''')

## true equivalent Q value
Flow('trueq1','qts1 repos','convert0eq repos=${SOURCES[1]}')

## Q estimated using the LCFS method and theoretical Q
Result('qq1','tq qts1',
       '''
       cat axis=2 ${SOURCES[1:2]}  |
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=150 o1num=0 d1num=50 n1tic=3
       ''')

## inverse Q filtering
Flow('cmtf1','asf trueq1',
     '''
     invqfilt eqt=${SOURCES[1]} gim=25
     ''')
Flow('cmsig1','cmtf1','transp|fft1 inv=y|window n1=1 f1=0')
Result('cmsig1',
       '''
       graph pad=n screenratio=0.8 crowd1=0.3 crowd2=0.75
       label2=Amplitude unit2= label1=Time unit1=s 
       title= yreverse=y transp=y wherexlabel=t plotfat=4
       labelfat=3 plotcol=7 font=2 axisfat=3 parallel2=n n2tic=20 
       ''')

## ----------------------------- SR1 -------------------------------
## LTFT
Flow('attf2','asignal','ltft rect=10')
Flow('mattf2','attf2','math output="abs(input)" |real|window min1=0.')
## pick spectrum
Flow('sp1','mattf2','window n1=150')
Flow('sp2','mattf2','window f1=150 n1=150')
Flow('sp3','mattf2','window f1=300 n1=200')
Flow('sp4','mattf2','window f1=500 n1=150')
Flow('sp5','mattf2','window f1=650 n1=150')
Flow('sp6','mattf2','window f1=800 n1=200')

Flow('f1-1','sp1','window n1=1 f1=101')
Flow('f2-1','sp2','window n1=1 f1=52')
Flow('f3-1','sp3','window n1=1 f1=104')
Flow('f4-1','sp4','window n1=1 f1=106')
Flow('f5-1','sp5','window n1=1 f1=57')
Flow('f6-1','sp6','window n1=1 f1=110')

## spectrum ratio
Flow('ldiv1-1','f2-1 f1-1',
     '''
     math den=${SOURCES[1]} output="input/den"|window max1=200|
     math output="-log(abs(input))"
     ''')
Flow('ldiv2-1','f3-1 f2-1',
     '''
     math den=${SOURCES[1]} output="input/den"|window max1=200|
     math output="-log(abs(input))"
     ''')
Flow('ldiv3-1','f4-1 f3-1',
     '''
     math den=${SOURCES[1]} output="input/den"|window max1=200|
     math output="-log(abs(input))"
     ''')
Flow('ldiv4-1','f5-1 f4-1',
     '''
     math den=${SOURCES[1]} output="input/den"|window max1=200|
     math output="-log(abs(input))"
     ''')
Flow('ldiv5-1','f6-1 f5-1',
     '''
     math den=${SOURCES[1]} output="input/den"|window max1=200|
     math output="-log(abs(input))"
     ''')

## sr-q
Flow('srq1-1','ldiv1-1',
     '''
     window min1=30 max1=80 | lineslope |
     window n1=1 | math output="%g*0.101/input" |
     spray o=0.1 n=100 | transp
     ''' %(math.pi))
Flow('srq2-1','ldiv2-1' ,
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.202/input" |
     spray o=0.2 n=200 | transp
     ''' %(math.pi))
Flow('srq3-1','ldiv3-1',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.202/input" |
     spray o=0.4 n=200 | transp
     ''' %(math.pi))
Flow('srq4-1','ldiv4-1',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.101/input" |
     spray o=0.6 n=100 | transp
     ''' %(math.pi))
Flow('srq5-1','ldiv5-1',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.203/input" |
     spray o=0.7 n=301 | transp
     ''' %(math.pi))
Flow('srq-1','srq1-1 srq2-1 srq3-1 srq4-1 srq5-1',
     '''
     cat axis=1 ${SOURCES[1:5]}|put o1=0.1 d1=0.001
     ''')
Result('difsrq-1','tq srq-1',
       '''
       window min1=0.1|cat axis=2 ${SOURCES[1:2]}|
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=155 o1num=0 d1num=50 n1tic=3     
       ''')

## ----------------------------- SR2 -------------------------------
## pick spectra
Flow('f1-2','sp1','window n1=1 f1=90')
Flow('f2-2','sp2','window n1=1 f1=40')
Flow('f3-2','sp3','window n1=1 f1=110')
Flow('f4-2','sp4','window n1=1 f1=100')
Flow('f5-2','sp5','window n1=1 f1=60')
Flow('f6-2','sp6','window n1=1 f1=110')

## spectrum ratio
Flow('ldiv1-2','f2-2 f1-2',
     '''
     math den=${SOURCES[1]} output="input/den" | window max1=200 |
     math output="-log(abs(input))"
     ''')
Flow('ldiv2-2','f3-2 f2-2',
     '''
     math den=${SOURCES[1]} output="input/den" | window max1=200 |
     math output="-log(abs(input))"
     ''')
Flow('ldiv3-2','f4-2 f3-2',
     '''
     math den=${SOURCES[1]} output="input/den" | window max1=200 |
     math output="-log(abs(input))"
     ''')
Flow('ldiv4-2','f5-2 f4-2',
     '''
     math den=${SOURCES[1]} output="input/den" | window max1=200 |
     math output="-log(abs(input))"
     ''')
Flow('ldiv5-2','f6-2 f5-2',
     '''
     math den=${SOURCES[1]} output="input/den" | window max1=200 |
     math output="-log(abs(input))"
     ''')

## sr-q
Flow('srq1-2','ldiv1-2',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.1/input" |
     spray o=0.1 n=100 | transp
     ''' %(math.pi))
Flow('srq2-2','ldiv2-2',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.22/input" |
     spray o=0.2 n=200 | transp
     ''' %(math.pi))
Flow('srq3-2','ldiv3-2',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.19/input" |
     spray o=0.4 n=200 | transp
     ''' %(math.pi))
Flow('srq4-2','ldiv4-2',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.11/input" |
     spray o=0.6 n=100 | transp
     ''' %(math.pi))
Flow('srq5-2','ldiv5-2',
     '''
     window min1=40 max1=80 | lineslope |
     window n1=1 | math output="%g*0.2/input" |
     spray o=0.7 n=301 | transp
     ''' %(math.pi))
Flow('srq-2','srq1-2 srq2-2 srq3-2 srq4-2 srq5-2',
     '''
     cat axis=1 ${SOURCES[1:5]}|put o1=0.1 d1=0.001
     ''')
Result('difsrq-2','tq srq-2',
       '''
       window min1=0.1|cat axis=2 ${SOURCES[1:2]}|
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=155 o1num=0 d1num=50 n1tic=3     
       ''')

## ---------------------------- CFS1 -------------------------------
## centroid frequency
Flow('zbw zcf zvar2','attf pf',
     '''          
     window max2=150 |
     lcf trange=${SOURCES[1]} avef=${TARGETS[1]} var2=${TARGETS[2]} rect1=1
     ''')
Flow('cf1-1','zcf','window n1=1 f1=101')
Flow('cf2-1','zcf','window n1=1 f1=202')
Flow('cf3-1','zcf','window n1=1 f1=404')
Flow('cf4-1','zcf','window n1=1 f1=606')
Flow('cf5-1','zcf','window n1=1 f1=707')
Flow('cf6-1','zcf','window n1=1 f1=910')
Flow('var1-1','zvar2','window n1=1 f1=101')
Flow('var2-1','zvar2','window n1=1 f1=202')
Flow('var3-1','zvar2','window n1=1 f1=404')
Flow('var4-1','zvar2','window n1=1 f1=606')
Flow('var5-1','zvar2','window n1=1 f1=707')
Flow('var6-1','zvar2','window n1=1 f1=910')

## cfsq
Flow('cq1-1','cf2-1 cf1-1 var1-1',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.101*r1/(r-input)"|
     spray axis=1 n=100 d=0.001 o=0.1
     '''%(math.pi))
Flow('cq2-1','cf3-1 cf2-1 var2-1',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.201*r1/(r-input)"|
     spray axis=1 n=200 d=0.001 o=0.2
     '''%(math.pi))
Flow('cq3-1','cf4-1 cf3-1 var3-1',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.202*r1/(r-input)"|
     spray axis=1 n=200 d=0.001 o=0.4
     '''%(math.pi))
Flow('cq4-1','cf5-1 cf4-1 var4-1',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.101*r1/(r-input)"|
     spray axis=1 n=100 d=0.001 o=0.6
     '''%(math.pi))
Flow('cq5-1','cf6-1 cf5-1 var5-1',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.203*r1/(r-input)"|
     spray axis=1 n=301 d=0.001 o=0.7
     '''%(math.pi))
Flow('cfsq-1','cq1-1 cq2-1 cq3-1 cq4-1 cq5-1',
     '''
     cat axis=1 ${SOURCES[1:5]}|put o1=0.1 d1=0.001
     ''')
Result('difcfsq-1','tq cfsq-1',
       '''
       window min1=0.1|cat axis=2 ${SOURCES[1:2]}|
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=155 o1num=0 d1num=50 n1tic=3     
       ''')

## ---------------------------- CFS2 -------------------------------
## centroid frequency
Flow('cf1-2','zcf','window n1=1 f1=100')
Flow('cf2-2','zcf','window n1=1 f1=200')
Flow('cf3-2','zcf','window n1=1 f1=400')
Flow('cf4-2','zcf','window n1=1 f1=600')
Flow('cf5-2','zcf','window n1=1 f1=700')
Flow('cf6-2','zcf','window n1=1 f1=900')
Flow('var1-2','zvar2','window n1=1 f1=100')
Flow('var2-2','zvar2','window n1=1 f1=200')
Flow('var3-2','zvar2','window n1=1 f1=400')
Flow('var4-2','zvar2','window n1=1 f1=600')
Flow('var5-2','zvar2','window n1=1 f1=700')
Flow('var6-2','zvar2','window n1=1 f1=900')

## cfsq
Flow('cq1-2','cf2-2 cf1-2 var1-2',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.1*r1/(r-input)"|
     spray axis=1 n=100 d=0.001 o=0.1
     '''%(math.pi))
Flow('cq2-2','cf3-2 cf2-2 var2-2',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.2*r1/(r-input)"|
     spray axis=1 n=200 d=0.001 o=0.2
     '''%(math.pi))
Flow('cq3-2','cf4-2 cf3-2 var3-2',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.2*r1/(r-input)"|
     spray axis=1 n=200 d=0.001 o=0.4
     '''%(math.pi))
Flow('cq4-2','cf5-2 cf4-2 var4-2',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.1*r1/(r-input)"|
     spray axis=1 n=100 d=0.001 o=0.6
     '''%(math.pi))
Flow('cq5-2','cf6-2 cf5-2 var5-2',
     ''' 
     math r=${SOURCES[1]} r1=${SOURCES[2]} output="%g*0.2*r1/(r-input)"|
     spray axis=1 n=301 d=0.001 o=0.7
     '''%(math.pi))
Flow('cfsq-2','cq1-2 cq2-2 cq3-2 cq4-2 cq5-2',
     '''
     cat axis=1 ${SOURCES[1:5]}|put o1=0.1 d1=0.001
     ''')
Result('difcfsq-2','tq cfsq-2',
       '''
       window min1=0.1|cat axis=2 ${SOURCES[1:2]}|
       graph plotcol=7,5 dash=3,0 pad=n screenratio=0.8 yreverse=y  
       crowd1=0.3 crowd2=0.75 wanttitle=n transp=y
       label2=Q unit2= label1=Time unit1=s 
       wherexlabel=t plotfat=4 font=2 parallel2=n labelfat=3 
       n2tic=20 min1=0.1 min2=0 max2=155 o1num=0 d1num=50 n1tic=3     
       ''')

End()