source: SH_SHM/trunk/util/extract_restfil.csh @ 92

Revision 16, 3.1 KB checked in by marcus, 15 years ago (diff)

r1 | svn | 2007-12-13 11:10:29 +0100 (Do, 13 Dez 2007) | 2 lines

Initial import

  • Property svn:executable set to *
Line 
1#! /bin/csh
2#
3# file extract_restfil.csh
4#      ===================
5#
6# version 1, 7-Mar-2006
7#
8# Extract restitution filters from SEED file.
9# K. Stammler, 7-Mar-2006
10
11if  ("$4" == "")  then
12        echo "Usage: $0 <seedfile> <station> <year> <doy>"
13        exit
14endif
15
16# get parameters
17set seedfile=$1
18set station=$2
19set year=$3
20set doy=$4
21
22if  (! -e $seedfile)  then
23        echo "$0 : input file $seedfile not found.  Abort."
24        exit
25endif
26
27# set constants
28set per1=1000.0
29set per2=0.05
30set tmpdir=$HOME/xresp
31set tmpfile=fil_$$.000
32set chanlist=chanlist.000
33set rdseed=rdseed4.6
34set evalresp=evalresp
35
36set res=`which $rdseed`
37if  ("`echo $res | grep 'not found'`" != "")  then
38        echo "$0 : no program rdseed ($rdseed) found.  Abort."
39        exit
40endif
41set res=`which $evalresp`
42if  ("`echo $res | grep 'not found'`" != "")  then
43        echo "$0 : no program evalresp ($evalresp) found.  Abort."
44        exit
45endif
46
47if  (-e $tmpdir)  \rm -rf $tmpdir
48mkdir $tmpdir
49cd $tmpdir
50
51if  (-e $chanlist)  \rm $chanlist
52$SEED_PROG/seedquickdump -locid $seedfile | awk '{print $4,$7,$10}' \
53        | grep -i "^$station-" | sort -u | awk '{print $1,1.0/$2,$3}' >$chanlist
54
55set frq1=`echo "scale=3; 1.0 / $per1" | bc`
56set frq2=`echo "scale=3; 1.0 / $per2" | bc`
57echo "frq range $frq1 $frq2"
58
59set upstat=`$DPROG/stringop upper $station`
60
61# get response files of station
62$rdseed <<END >/dev/null
63$seedfile
64
65
66R
67$station
68
69
70
71quit
72END
73
74$evalresp $upstat '*' $year $doy $frq1 $frq2 1000 -f $PWD -u dis
75
76#set echo
77set cnt=1
78while  (1 > 0)
79        set line=`sed -n "$cnt"p $chanlist | sed 's/-/ /g'`
80        if  ("$line" == "")  break
81        set chan=`$DPROG/stringop upper $line[2]$line[3]`
82        set nyq=`$DPROG/floatop $line[4] div 2.0`
83        if  ($#line >= 5)  then
84                set locid=`$DPROG/stringop upper $line[5]`
85                set f="*.$upstat.$locid.$chan"
86        else
87                set locid="__"
88                set f="*.$upstat.*.$chan"
89        endif
90        set resp=`ls AMP.$f`
91        set phase=`ls PHASE.$f`
92        if  ($#resp > 1)  then
93                echo "more than one $f"
94        else if  ($#resp == 1 && $#phase == 1)  then
95                set rnum=`wc -l $resp`
96                set rnum=$rnum[1]
97                set pnum=`wc -l $phase`
98                set pnum=$pnum[1]
99                if  ($rnum != $pnum)  then
100                        echo "different number of lines in AMP and PHASE file ($resp)"
101                else
102                        set fname=`$DPROG/stringop upper ${station}:${chan}:$locid.FLT`
103                        # remove all filter entries above 92% of the nyquist frequency
104                        # otherwise the FIR filter stages create very large numbers
105                        # in the inversion filter towards the Nyquist frequencies
106                        set maxfrq=`$DPROG/floatop $nyq mul 0.92`
107                        if  (-e $tmpfile)  \rm $tmpfile
108                        paste $resp $phase | awk '{print $1,1.0e9/$2,-1.0*$4}' \
109                                | $DPROG/colrange 1 0 $maxfrq >$tmpfile
110                        set linenum=`wc -l $tmpfile`
111                        set linenum=$linenum[1]
112                        if  (-e $fname)  \rm $fname
113                        # two more line for frequencies 0 and $nyq
114                        @ linenum = $linenum + 2
115                        touch $fname
116                        echo "1357913578" >>$fname
117                        echo "2"          >>$fname
118                        echo "$linenum"   >>$fname
119                        # the first line with zeroes is needed to prevent out of range errors
120                        # But the filter function for periods larger than $per1 are wrong!
121                        echo "0.0 0.0 0.0">>$fname
122                        cat $tmpfile      >>$fname
123                        \rm $tmpfile
124                        echo "$nyq 0.0 0.0">>$fname
125                endif
126        endif
127        @ cnt ++
128end
129
130cd $HOME
131#\rm -rf $tmpdir
Note: See TracBrowser for help on using the repository browser.