est+ from est_noise

#!/bin/sh
#  sample script to run est_noise6ac
#  test data is CARH north from Tom Herring's 'combined' SCIGN solutions


rm -f seed.dat
echo 82789427 > seed.dat     # create a random number for doing hypothesis testing a la Langbein [2004] Figure 4

rm -f junk.in
cat > junk.in <otr                  # format style of inputting specifications (using julian day style; year, month and day could be input, too)
1                   #  Number of data sets  (almost always 1)
1995 1 2005 1       # time span to be considered
y   # rate             #  Estimate the secular rate; almost alway 'y' for yes
0   # num of rate change; If you think there are any rate changes in data, this should be the number of rate changes;   if non-zero, then you'll need to supply the time of the rate change
1  # num of periods;   I'm solving for a 1-year periodicity in the data; this is NOT a band-pass filter
365.25         #  period in days
1              #  Number of offsets in data
2003 357.0     # time of offset due to San Simeon EQ (M6+)
0             # number of exponential exp(-t/tau) or Log (1 + t/tau) or Log (tau + t ) terms for characterizing post-seismic responce. This seems to work okay; outputs searchs for optimal tau on file fort.89
otr      #  format style of input data
CARH.n 
0    #  Number of time-series of auxillary data --- This is useful for analyzing Dilatometer data -- this would be the presure time series
1   # mininum sampling in days  (usually 1 day for CGPS)
n   # This is usually n for 'no'.  With 2-color EDM stuff, sometimes I'd have more than one measurement on a single day and I'd want to combine these data
n  # n for 'no' unless you want to do simulations for hypothesis testing (Figures 2-4 in Langbein, 2004)
0   #  decimation parmeter; choice of 0, 1, 2, or 3.  Use 0 for 'best' results as it uses all of the available data; the other choices toss-out observations but are useful for 'quick and dirty' results-- more exploritory in nature
1 float    #white noise  Initial guess of white noise is 1 mm and you are asking the program to estimate the white noise amplitude
1 float   #power law amplitude; initial guess of PL amplitude is 1 and you are asking the program to estimate a better amplitude
1 float    # power law index; initial guess of PL index is 1 (flicker) and you are asking the program to estimate a better index
0 fix   #Gauss Markov;  fixing the GM time constant to 0; do not make any estimates
.5 2.   # Pass-band limits of the band-pass filtered noise  0.5 cycles/year to 2 cycles/year
1      #  Number of poles in the BP filter; (goes 1, 2, 3, 4, 5, and 6)
0 fix  # Amplitude of BP noise; in this example, it is set to 0 and not estimated
2 fix   # power law index;  Index for a second power law; it is set to 2 (random walk) and not updated in this example
0.0 fix  # power law amplitude; amplitude for a second power law; it is set to 0 and not updated in this example
0  #  Additional white noise to add to data (data.in); usually set to 0 but for cases such as strain or creep data where the power-law index is > 2.5, it becomes critical to add some white noise to the data in order to keep the data covariance from becoming singular.
EOF

time ~/proglib/Strain/Progs/est_noise6ac < junk.in

exit
EOF