Processing SCG data


Contents:

Human effects
Pointing project



# HOW.TO

# make downloaded bz2-files unique (delete co-existing uncompressed files)
set cmd=echo
foreach file ( *.054 )
if ( -r $file.bz2 ) then
   $cmd rm -f $file
endif
end
#
# inspect and, if/when o.k., then run once more, but now after "set cmd="
#__________________________________________________________________________
#
# HOW.TO
# ... get 1-min gravity, baro, and predicted tides into one MC-file
#
rm -f d/MC0906-1m.ts
cat MON_o054/GW09??00.GGP |\
 awk '/^2009/{printf "%s %s %11.7f\n",$1,$2,$3}' |\
   tslist - -gi4,2i2,i3,i2,t16,f12.0 -k2 -I -O:`label GRAV,VAL` d/MC0906-1m.ts
cat MON_o054/GW09??00.GGP |\
 awk '/^2009/{printf "%s %s %11.4f\n",$1,$2,$5}' |\
   tslist - -gi4,2i2,i3,i2,t16,f12.0 -k2 -I -O:`label BARO,VAL` d/MC0906-1m.ts
run_urtip d/MC0906-1m.ts
tslist o/pt.tsf -gI4,5i3,d18.0 -k3 -I -O:`label GRAV,PRED` d/MC0906-1m.ts
#
# See script ./myfo for a recipe how to generate de-tided and de-baroed gravity
#
# See the new script resids e.g.
#
resids -xp -d 10 2009 06 15 2009 10 31
#
# creates 10-min data  d/MC0906-10m.ts and a plot in plot/MC0906-10m.ps
#__________________________________________________________________________
#
# ... get 1s-data gravity into a daily, labelled ts-file
#
foreach day ( 8 9 10 )
tsf2ts-decim -r $day,$day -W -c 5 -d 1 -O GRAV,VAL d RAW_o054/ A2 0909
end
#
#or
tsf2ts-decim -w -r 8,10 -W -c 5 -d 1 -O GRAV,VAL d RAW_o054/ A2 0909
#
# ... make an RMS-plot of this segment, not admitting but extreme shocks
setenv SHOCKSIZE 2.
tslist d/A2090908-1s.ts -Ermsmode.tse,DRMS -I -SO776. -o d/A2090908-1s-rms.ts
setenv TSPLOT_FILE ../d/A2090908-1s-rms.ts
setenv TSPLOT_KALY 0.05
set legtxt='1-min RMS'
set PSOUT=A2090908-1s-rms.ps
cd plot
source hans-tsplot.env
# BUT: edit hans-tsplot.env first and take away the -D option
#
#__________________________________________________________________________
#
# ... get 10min-data from a month of RAW-A2 files:
 tsf2ts-decim -w -r 1,31 -c 11,19 -o d -d 600,0 -k RAW_o054/ A2 0908
#
# Explanation:
# -w        : create new collecting ascii files (if not, we would append).
# -r 1,31   : day range 1 to 31 (it's August)
# -c 11,19  : channels 11 and 19, the tilt controllers
# -o d      : output to subdir d
# -d 600,0  : decimate according to file deci.tse label 600
#             The overlap 4*int(600/2)
# -k        : keep all intermediate ascii files
# RAW_o054/ : input directory
# A2        : A2 first in file name
# 0908      : year month
#
# This job could abort if the data of September 01 is not available.
# To create the binary files of what has been collected, do
tsf2ts-decim -t -c 11,19 -o d -k RAW_o054/ A2 0908
#____________________________________________________________________________
#
# Create an empirical residual
# ( after urtap ~/Ttide/SCG/urtap-tides.ins on all 10-min data,
# and ~/TD/expfitm.ins )
#
# here after tsf2ts-decim -w -W -c 1,3 -r 1,53 -d 300 -n 1200 RAW_o054/ G1 0910
# getting 7 weeks of 300-s data
ts2mc-resid -PP -drift exp.tse,EXP -shocks newshocks.tse,1.0 \
    -o o/G1_1_091001-300s.r.ts
    d/G1_1_091001-300s.ts d/G1_3_091001-300s.ts o/pt.ts
#____________________________________________________________________________
#
# ... get 1s data undecimated, one day with a tail of 1200 next day.
tsf2ts-decim -w -r 12,12 -c 11 -d 1,0 -n 1200 -o d RAW_o054/ A2 0908
#
# This is for channel 11 only. -d 1,0 will not invoke the tsfedit subroutine.
#____________________________________________________________________________
#
# ... append a data set
cd d
setenv APPFILE G1_1_090813.ts
tslist gw10907-10m.ts -Eappend.tse,APPEND -I -qbe -C3 -O:GRAV_________VAL G1_13_090701-10m.ts
setenv APPFILE G1_3_090813.ts
tslist gw30907-10m.ts -Eappend.tse,APPEND -I -qbe -C3 -O:BARO_________VAL G1_13_090701-10m.ts
#
#____________________________________________________________________________
#
# A simple oneliner to get mode 1-min data

xchan-units -x 15 MON_o054/GW090900.TSF |\
 tslist - -gi4,5i3,f15.0 -k3 -rs60 -C3 -I -qbe -o o/rt.ts

#__________________________________________________________________________
#
# Daily power spectrum:
daily-psp -lf3600 2009 09 03
#
# creates a.o. a file o/brt.ts of residuals for that day
#
# The final reduction for a drift model, deletion of shocks, and
# scaling with calibration factor can be done
tslist o/brt.ts -Eexp.tse,EXP -SO-776. -o o/dbrt.ts
#
# ...make a spectrum plot
set df=`tslist o/brt.ts -I -Ermsmode.tse,PDG -Ni5 -n1-1 -Ff15.5 | awk '/FNy/{print 1000*$7}'`

#tslist o/rt.ts -Ermsmode.tse,PDG -Ni5 -n1-1 -Ff15.5 |\
# fgrep -v '>' | awk '{print $1*0.002596053997923157,$2}' > ! o/modesp.dat

tslist o/brt.ts -Ermsmode.tse,PDG -Ni5 -n1-1 -Ff15.5 |\
 fgrep -v '>' | awk -v df=$df '{print $1*df,$2}' > ! o/modesp.dat

psxy o/modesp.dat -R0/2/-10/10 \
     -Ba1f0.1:"f [mHz]":/a10f1:"A [dB]"::."Mode spectrum 2009-OCT-01":WeSn \
     -JX9/4 -W1/0 > ! plot/modesp.ps

# the strange factor is "df" and it is printed by tslist; look for
# <TsfEdi>>> Spectrum: FNy, df, dt:
# Multiply with 1000 to get mHz
#
# ... make a plot of sliding signal RMS, a kind of VU-meter
 
tslist o/rt.ts -Ermsmode.tse,RMS -I -C3 -qbe -o o/moderms.ts
tslist o/moderms.ts -J -Bc2009,9,29 | fgrep -v '>' |\
 psxy -R55103/55108/-0.5/10 \
      -Ba1f0.0833333333:"MJD":/a2f0.5:RMS::."Mode channel RMS":WeSn \
      -JX9/4 -W1/0 > moderms.ps

# where rmsmode.tse is:
# TSF EDIT RMS
# MOVRMS 10 10 From 2009 09 29 0 0 0 0 To 2009 10 04 0 0 0 0
# END
#
# TSF EDIT PDG
# PERIODOGRAM  From 2009 09 30 13 0 0 0 To 2009 10 04 23 59 59 0
# END

#______________________________________________________________________
#
# Own mode filter starting with 1-s data, 15-s downsample and tsplot
#
# Use script ./myfo
#
# The filter ( rmsmode.tse, MYFO ) is not well tuned
#
#TSF EDIT MYFO
#FILTER D:WD KAIS 240 1.5 0.012 0.1 1.0 0.02
#FILTER 0 1 ' '
#1.0 -1.0
#END
#
# using this for the earthquake 2009-AUG-10:
#
myfo -rm 2009 08 10 13
#
# Plot the time series:
setenv TSPLOT_TSLIST '-BHc2009,08,10,20,00,00,00
setenv TSPLOT_TSLIST '-BHc2009,08,10,20,00,00,00'
setenv PSOUT FO090810.ps
setenv TSPLOT_FILE ../o/fo.ts
# adjust the limits and ticks inside plot/hans-tsplot.env
# setenv TSPLOT_YR -0.002/0.002
# setenv TSPLOT_YTIX a0.001f0.0002
#
cd plot
source hans-tsplot.env
#
cd ~/TD
tslist o/fo.ts -Ermsmode.tse,PDG -N -n-1'*'0.003968253968253968 -S1000. -FF12.3 |\
 fgrep -v '>' > ! o/fosp090810.dat
#
# You can get the factor (in mHz !) with

set df=`tslist o/fo.ts -Ermsmode.tse,PDG -FF12.3 -I | awk '/Spectrum: FNy/{printf "%-15.8E\n", 1000*$7}'`

tslist o/fo.ts -Ermsmode.tse,PDG -N -n-1'*'${df}+0 -S1000. -FF12.3 |\
 fgrep -v '>' > ! o/fosp090810.dat

#
awk '($1>0){print $1, $2/2/sin(3.1415926535*$1/2/33.3333333)}' o/fosp090810.dat |\
 psxy -R0/2.5/0/1000 -Ba1f0.1:"f [mHz]":/a200f100:"A"::."Mode spectrum 2009-AUG-10":WeSn \
  -JX9/4 -W1/0 > ! plot/fosp090810.ps
#
# The awk above restores the diff-filter's attenuation.
#_________________________________________________________________________________________
#
# ... make a daily suite of power spectra from tide+baro residuals
#
# check parameters and files in sasm03-daily.ins
# check that run_urtip is o.k.
# check that t/urtap.tsd is available
# i.e. make a 1-day test taking the first day:
daily-psp -C o/collect_sp.xyz 09 06 13
# CAPITAL -C !
#
# If the results are o.k. you can issue
set db=`JDC -d -m -fi5 2009 06 14`
set de=`JDC -d -m -fi5 2009 09 30`
foreach d ( `fromto $db $de` )
daily-psp -lf3600 -y-10,40 -r0.0002,0.005 -dB -c o/collect_sp.xyz `JDC -j -m -y2k $d`
end
# lower-case -c !
#
# There is a super-script
run-daily-psp-vs-time -px 2009 06 15 2009 10 10
# -px: do both data preparation and plot
#
# The power spectra are intermittently perturbed by earthquakes.
# Make a list of DELete records for tsfedit using
setenv SHOCKSIZE 0.02
echo "TSF EDIT SHOCKS" > ! newshocks.tse
foreach file ( d/G109*-1s.ts )
tslist $file -Ermsmode.tse,CDRMS -C3 -I
end
echo END >> newshocks.tse
# eventually update the long list shocks.tse
#
# TSF EDIT DRMS
# FILTER 0 1 ' '
# 1.0 -1.0
# MOVRMS 60 120
# OPEN 41 A newshocks.tse
# WDR ${SHOCKSIZE} U=41
# END
#
# This procedure is carried out by
#
#    detect-shocks #shocksize
#______________________________________________________________________________
#

# Here's a show of human induced disturbance:

#10 visitors (CC Tscherning + students + HGS
# standing in a semicircle around the SCG, 1 minute up and one minute kneeling.
# Day and time was 2009-11-11 12:20-12:45 UT
#
# overwrite ascii, day Range 11,11, decimate 1, overlap n=0, column 1, output-dir d/
# from RAW_o054 file type G1, yymm=0911
#
tsf2ts-decim -w -r 11,11 -d 1 -n 0 -c 1 -o d RAW_o054/ G1 0911
tslist d/G1091111-1s.ts -Elpf.tse,LPF100 -BHc2009,11,11,12,20 -U2009,11,11,12,45 -S-77.6 -C3 -I -o t.ts
cd plot
set TSPLOT_FILE=../t.ts
set PSOUT=t.ps
set legtxt="Gravity LoPassFiltered"
source hans-tsplot.env
ps2png -m -rr t.ps
#
# The resulting PNG/t.png was renamed to CCT-On-Tour.png
#
# with a few modifications: y-label uGal, title
# The low-pass filter in lpf.tse was
# TSF EDIT LPF100
# FILTER D:WD KAIS 100 2.1 0.0 0.05 0.5 0.0
# /FILTHR S 30 20 10 5 1
# END
#_______________________________________________________________________________
#
# A similar case was on the day before, Lantmteriet visiting.
# That day the tide was strongly going up, so we subtract the predicted
# tide using the empirical model.
#
tsf2ts-decim -w -r 10,10 -d 1 -n 0 -c 1 -o d RAW_o054/ G1 0911

rm t.mc # !!!!!

tslist d/G1091110-1s.ts -BHc2009,11,10,15,00 -U2009,11,10,16,45 -C3 -I -O:`label SCG,VAL` t.mc
run_urtip t.mc
tslist o/pt.tsf -gi4,5i3,t20,d18.0 -k3 -Ff10.3 -C3 -qbe -rs1.0000 -O:`label SCG,PRED` t.mc -I
tslist t.mc -D -L'S|V' -ML1,1,-:'S|P' -Elpf.tse,LPF100 -SO-776. -I -o t.ts
#
#_______________________________________________________________________________
#
# A tslist option sequence for getting a residual from an mc-file:
set pc=`awk '/^ BARO /{print $4}' t/urtap.prl`
# -SO-776. -L'G|V' -ML'1,1,-:G|P' -ML$pc',1,-:B|V]Rwd'

resids -xpo -d 60,0 ~/BOOS/Ringhals-0911.ts
# creates GRAV-TIDE-BARO residual:  d/r0911-60m.ts
#
tslist ~/BOOS/Ringhals-0911.ts -O:`label RITGG,VAL` d/MC0911-60m.ts
# adds tide gauge to MC-file
#
tslist d/MC0911-60m.ts -L'B|V' -ML1,1,+:'RITGG|V' -F2D12.4 -D -I -o d/tggap.ts
# creates tidegauge + airpressure (factor 1 [hPa/cm] is slightly wrong)
setenv TSPLOT_YR -30/50
setenv TSPLOT_FILE "../d/r0911-60m.ts ../d/tggap.ts"
set PSOUT=grav-tgg.ps
set TSPLOT_TSLIST = ( -SO1. "-SO1. -Ul4" )
set legtxt = ( SCG-TIDE-BARO TGG+BARO )
cd plot
source hans-kattegat-tsplot.env


#-------------------------------------------------------------------------
#

# The pointing project

cd astro
rm pointing.mc
awk -F, -f logsrc.awk *.log | eq2ho.out 11.9 57.6 - |\
 tslist - -gi5,2i3,i5,2i3,t39,f10.0 -k3 -C3 -rs30 -00 -Esh.tse,SH -O:`label POINTING,COS` pointing.mc
awk -F, -f logsrc.awk *.log | eq2ho.out 11.9 57.6 - |\
 tslist - -gi5,2i3,i5,2i3,t49,f10.0 -k3 -C3 -rs30 -00 -Esh.tse,SH -O:`label POINTING,SIN` pointing.mc
# careful with the -0 option, Eugene! In this case it really rounds to the right
# second. Without -00 the first time is 2009 06 23  08 57 28
# With if the first time is 2009 06 23  08 57 30
#
cd ..
tsf2ts-decim -w -W -d 30,30 -c 1 -d 22,131 -O GRAV,VAL o RAW_o054/G1 0906
tsf2ts-decim -w -W -d 30,30 -c 3 -d 22,131 -O BARO,VAL o RAW_o054/G1 0906
run_urtip o/MC090622-30s.ts
#
# exec the suggested tslist command o/pt.tsf -> -I -O:`label GRAV,PRED`  o/MC090622-30s.ts
#
# trim the data sets to common origin for astro/pointing.mc and o/MC090622-30s.ts -> o/gdr4pointing.mc (DEFER! *)
# specifically,
set trimdate=2009,06,23,08,57,30
set pc=`awk '/^ BARO /{print $4}' t/urtap.prl`
#
# Without warranty
# tslist  o/MC090622-30s.ts -I -L'G|V' -ML'1,1,-:G|P' -ML$pc',1,-:B|V]Rwd' -Bc$trimdate -00.025 \
#   -Eshocks.tse,SHOCKS -o tmp.ts
# This -00.025 is special; it moves a 40 sec origin to the nearest 30 sec.
# tslist tmp.ts -I -S-776. -Eexp.tse,EXPM -O:`label GRAV,GDR` o/gdr4pointing.mc
# DEFER *) now, copy tslist astro/pointing.mc  -> o/gdr4pointing.mc
#
# Now, the first column in o/gdr4pointing.mc is residual gravity, drift corrected, earthquake-edited,
# and scaled to nm/s^2. We want to remove DC-levels in each experiment
#
tslist o/gdr4pointing.mc -L'P|C' -C3 -F1p,d12.4 -MV | \
   awk 'BEGIN{print "TSF EDIT REMDC"} \
        \!/>/{if(q){print s,"To",substr($0,1,20),"00"};q=0} \
        /VALID/{s= "REMDC   From "substr($0,1,20)" 00";q=1} \
        END{print "END"}' \
 >! astro/remdc.tse
tslist o/gdr4pointing.mc -L'G|G' -C3 -F1p,d12.4 -Eastro/remdc.tse,REMDC -I -o temp.ts
tslist temp.ts -I -O:`label GRAV,GDRBDC` o/gdr4pointing.mc
#
# The gaps between VLBI experiments are removed as follows (time information losing its meaning then)
set olabs=`tslist o/gdr4pointing.mc -L | awk '/LocGMc/{i++;printf "%s|",i":"substr($3,7)}'`
tslist o/gdr4pointing.mc -L'G|G' -L'P|C' -L'P|S' -MSYNCON -F1p,3d12.4 -C3 -O$olabs astro/gdrc4pointing.mc
#
# Column labelled  G|GDRBDC     is  the time series on the r.h.s. of the lsq-problem,
#                  P|C and P|S  are the l.h.s. signals
# Good luck!