**************************************************************************************************
*
* This macro goes through the waveforms of a single event and computes the source parameters. Although
* the computation is performed for a number of different parameters (Mw, fc, Dsigma), focus is on Mw.
* Corner frequency and stress parameter are unstable, and their automatic computation will need a 
* substantial amount of work to be done on this macro. 
* Computation can be done completely automatically; the manual processing option is kept for special cases,
* and only involves the proper cutting time window for S-waves. Just mark t5 and t6 to comprise the wavetrain
* to be used for the spectral analysis.
* Local magnitude is used in order to properly window the S-waves, and also in order to apply a hi-pass
* filter to the current waveforms before the spectral analysis is performed. This is done in order to
* eliminate the microseismic noise in those cases when the noise spectral level at low frequency would be higher
* than the spectral displacement plateau. Corner frequencies are computed on each waveform by using Andrews (1986)
* formula (Earthquake Source Mechanics, S. Das, J. Boatwright and C.H. Scholz Editors, A.G.U. Monograph 
* 37, 259-267).
*
* Luca Malagnini, January 4, 2002
*
* The code needs the program ttimes and the script comp_times_new.sh, that are to be located in the same 
* SAC macro directory, and it is used by running the following shell script, that needs to be located in the 
* Data directory:
*
****************************************************************************************************
*#!/bin/sh
*MAIN_DIR=`pwd`
*cat<<eof>get_average.f
*        sum=0.
*        n=0
*1       read(5,*,end=999,err=999)amag
*        n=n+1
*        sum=sum+amag
*        goto 1
*999     continue
*        ave=sum/n
*        write(6,*)ave
*        end
*eof
*make get_average
*rm Moment_Magnitudes
*ls -d [0-9][0-9][0-9] > list_dir
*ls -d [0-9][0-9][0-9][0-9] >> list_dir
*while read dir
*do
*        cd ${MAIN_DIR}/${dir}
*        sac2000<<EOF
*m comp_t0_mac
*m comp_ml_mac
*Wood-Anderson
*m prova_friuli
*no
*EOF
*MW=`cat tmp_mw|${MAIN_DIR}/get_average`
*echo " "
*echo Event ${dir}>>${MAIN_DIR}/Moment_Magnitudes
*cat Mw_output>>${MAIN_DIR}/Moment_Magnitudes
*echo Average Moment Magnitude   Mw= ${MW}>>${MAIN_DIR}/Moment_Magnitudes
*echo " "
*echo "------------">>${MAIN_DIR}/Moment_Magnitudes
*done<list_dir
*
*
*************** Here is the script comp_times_new.sh ******************************
*
*#!/bin/sh
*cat<<EOF>make_ptimes_stations.awk
*{
*        print \$1 ,  \$5 
*}
*EOF
*cat<<EOF>get_ptimes.awk
*BEGIN { FS=" "
*        RS="\n"
*}
*NF == 9  && ( \$2 ~ /1/ && \$3 ~ /Pg/ || \$2 ~ /1/ && \$3 ~ /Pn/ ||  \$2 ~ /1/ && \$3 ~ /Pb/ )  {
*        print "Delta " "\t" \$1 "\t" \$3 "\t" \$4 "\t" \$7
*}
*EOF
*
*cat<<EOF>get_stimes.awk
*BEGIN { FS=" "
*        RS="\n"
*}
*
*NF == 9  && ( \$2 ~ /1/ && \$3 ~ /Sg/ || \$2 ~ /1/ && \$3 ~ /Sn/ ||  \$2 ~ /1/ && \$3 ~ /Sb/ )  {
*        print "Delta " "\t" \$1 "\t" \$3 "\t" \$4 "\t" \$7
*}
*EOF
*read depth<depth.out
*read gcarc<gcarc.out
*cat > ttimes.in <<EOF
*S
*
*${depth}
*${gcarc}
*-1
*-1
*EOF
*ttimes<ttimes.in>ttimes.out
*awk -f get_stimes.awk<ttimes.out>stimes.out
*paste statname.out stimes.out |awk -f make_ptimes_stations.awk>stime.out
*read statname<statname.out
*read stat stime<stime.out
*sac2000<<eof
* echo on
* r $1.r $1.t
* evaluate to t0 ( &1,o + ${stime} ) 
* ch t0 %t0%
* wh over
* exit
*eof
*
*
**************************************************************************************************
plabel 1 off
plabel 2 off
plabel 3 off
plabel 4 off
plabel 5 off
plabel 6 off
plabel 7 off
cut off
setbb manual (REPLY "Type yes for manual processing, no for automatic...")
*setbb woodanderson (REPLY "Choose which way to compute local magnitude (Wood-Anderson - Richter formula; Udias - Udias formula) ...")
*
setbb woodanderson Wood-Anderson
*setbb woodanderson Udias

echo on
*
* Remove old versions of files that will be re-created later...
*
sc rm tmp_* *_tmp.? Mw_output
*
* Remove the previous S-wave marker and substitute it with a theoretical one... Then write out a new file
*
* At the same time, check the local magnitude given by FVGSN and hipass filter waveforms accordingly...
*
do file wild ???.r ????.r
 message 'file:  $file '
 setbb statname (BEFORE ".r" $file )
 r %statname%.r %statname%.t
 rtr
 rmean
 setbb t0_r &1,t0
 setbb t0_t &2,t0
 ch lcalda true
 ch lovrok true
 sc rm ttimes.in ttimes.out gcarc.out depth.out stimes.out stime.out
 evaluate to depth ( &1,evdp / 1000 )
 TRANSCRIPT CREATE FILE gcarc.out CONTENTS OUTPUT
 message " &1,gcarc "
 TRANSCRIPT CLOSE
 TRANSCRIPT CREATE FILE depth.out CONTENTS OUTPUT
 message " %depth% "
 TRANSCRIPT CLOSE
 TRANSCRIPT CREATE FILE statname.out CONTENTS OUTPUT
 message " %statname% "
 TRANSCRIPT CLOSE
 w %statname%_tmp.r %statname%_tmp.t
 if %t0_r% eq UNDEFINED then
  if %t0_t% eq UNDEFINED then 
   sc sh /home2/RAMI1/Data/comp_times_new.sh %statname%_tmp
  else
   r %statname%_tmp.r 
   ch t0 %t0_t%
   w over
  endif
 endif
enddo
*
*
* Now make sure that the h-comp. have the same number of points...
*
do file wild *_tmp.r
 message 'file:  $file '
 setbb statname (BEFORE ".r" $file )
 r %statname%.r %statname%.t
 ch lcalda true
 ch lovrok true
 synch
 wh over
*
 setbb c0 &1,b
*
 if &2,b gt %c0 then
 setbb c0 &2,b
 endif
*
 setbb c1 &1,e
*
 if &2,e lt %c1 then
 setbb c1 &2,e
 endif
*
 cut %c0 %c1
 r %statname%.r %statname%.t
 cut off
 synch
 w over
enddo
sc rm -rf tmp*.wa
do file wild *_tmp.r
 message 'file:  $file '
 setbb statname (BEFORE ".r" $file )
 r %statname%.r %statname%.t
 setbb dist &1,dist&
 setbb t0 &1,t0&
 rtr
 rmean
 taper
 ch lcalda true
 ch lovrok true
********************************
 if %woodanderson% eq Udias
*
* Compute local magnitude from Udias, p. 286: Ml = log(A) + 2.56*log(dist) - 1.67 
*
  int
  hp bu n 4 p 2 c .1
  lp bu n 4 p 2 c 20.
*
* Needs ground motion (displacement) in micrometers...
*
  mul 1.e+06
  w tmp1.wa tmp2.wa
  evaluate to cut1 ( %t0% - 3.0 )
  evaluate to cut2 ( %t0% + 40.0 )
  cut %cut1% %cut2%
  r tmp1.wa tmp2.wa
  rtr
  rmean
  cut off
  qdp off
*  width 1
*  win 3 x .01 .4 y .10 .40
*  xdiv nice
*  xdiv p off
*  beginwindow 3
*  xlim off
*  title "Waveforms for the Automatic Computation of Local Magnitude @(Udias formula@)" l t s l
*  ylabel "Corrected Ground Displacement @(micrometers@)" l l s m
*  xlabel "Time @(s@)" l b s m
*  plabel 1 'Station: %statname% ' position .22 .80 s m
*  p1
*  plabel 1 off
*  title off 
*  xlabel off 
*  ylabel off
*
* Now compute distance correction for Ml...
*
  evaluate to dist_corr ( alog10 ( %dist% ) )
  evaluate to dist_corr ( 2.56 * %dist_corr% )
*
* Now compute max value for each component and average them out...
*
  evaluate to maxmin ( abs ( &1,depmin ) )
  evaluate to maxmax ( abs ( &1,depmax ) )
*
  if %maxmin% gt %maxmax%
   setbb peak_val %maxmin%
  else
   setbb peak_val %maxmax%
  endif
  evaluate to loc_mag_1 ( alog10 ( %peak_val% ) )
  evaluate to loc_mag_1 ( %loc_mag_1% - 1.67 )
*
* Now insert distance correction...
*
  evaluate to loc_mag_1 ( %loc_mag_1% + %dist_corr% )
*
  evaluate to maxmin ( abs ( &2,depmin ) )
  evaluate to maxmax ( abs ( &2,depmax ) )
 *
  if %maxmin% gt %maxmax%
   setbb peak_val %maxmin%
  else
   setbb peak_val %maxmax%
  endif
  evaluate to loc_mag_2 ( alog10 ( %peak_val% ) )
  evaluate to loc_mag_2 ( %loc_mag_2% - 1.67 )
*
* Now insert distance correction...
*
  evaluate to loc_mag_2 ( %loc_mag_2% + %dist_corr% )
*
  evaluate to loc_mag ( ( %loc_mag_1% + %loc_mag_2% ) / 2 )
 endif
*
*
*
 if %woodanderson% eq Wood-Anderson
*
* Compute local magnitude from classic Ml table...
*
* r_tab_low: table of 121 distance values...
* r_tab_hi: table of 121 distance values...
* rloga0_tab_low: table of 121 log values, from Richter...
* rloga0_tab_low_hi: tab_lowle of 121 log values, from Richter...
*
  setbb r_tab_low                        "0 5 10 15 20 25 30 35 40 45 50"
  setbb r_tab_low            "%r_tab_low% 55 60 65 70 75 80 85 90 95 100"
  setbb r_tab_low   "%r_tab_low% 105 110 115 120 125 130 135 140 145 150"
  setbb r_tab_low   "%r_tab_low% 155 160 165 170 175 180 185 190 195 200"
  setbb r_tab_low   "%r_tab_low% 205 210 215 220 225 230 235 240 245 250"
  setbb r_tab_low   "%r_tab_low% 255 260 265 270 275 280 285 290 295 300"
  setbb r_tab_low   "%r_tab_low% 305 310 315 320 325 330 335 340 345 350"
  setbb r_tab_low   "%r_tab_low% 355 360 365 370 375 380 385 390 395 400"
  setbb r_tab_low   "%r_tab_low% 405 410 415 420 425 430 435 440 445 450"
  setbb r_tab_low   "%r_tab_low% 455 460 465 470 475 480 485 490 495 500"
  setbb r_tab_low   "%r_tab_low% 505 510 515 520 525 530 535 540 545 550"
  setbb r_tab_low   "%r_tab_low% 555 560 565 570 575 580 585 590 595 junk"
*
  setbb r_tab_hi                         "5 10 15 20 25 30 35 40 45 50"
  setbb r_tab_hi            "%r_tab_hi% 55 60 65 70 75 80 85 90 95 100"
  setbb r_tab_hi   "%r_tab_hi% 105 110 115 120 125 130 135 140 145 150"
  setbb r_tab_hi   "%r_tab_hi% 155 160 165 170 175 180 185 190 195 200"
  setbb r_tab_hi   "%r_tab_hi% 205 210 215 220 225 230 235 240 245 250"
  setbb r_tab_hi   "%r_tab_hi% 255 260 265 270 275 280 285 290 295 300"
  setbb r_tab_hi   "%r_tab_hi% 305 310 315 320 325 330 335 340 345 350"
  setbb r_tab_hi   "%r_tab_hi% 355 360 365 370 375 380 385 390 395 400"
  setbb r_tab_hi   "%r_tab_hi% 405 410 415 420 425 430 435 440 445 450"
  setbb r_tab_hi   "%r_tab_hi% 455 460 465 470 475 480 485 490 495 500"
  setbb r_tab_hi   "%r_tab_hi% 505 510 515 520 525 530 535 540 545 550"
  setbb r_tab_hi   "%r_tab_hi% 555 560 565 570 575 580 585 590 595 600 junk"
*
  setbb rloga0_tab_low             "-1.3 -1.40 -1.50 -1.60 -1.70 -1.90 -2.10 -2.32 -2.43 -2.54 -2.63"
  setbb rloga0_tab_low "%rloga0_tab_low% -2.70 -2.77 -2.79 -2.83 -2.87 -2.90 -2.94 -2.96 -2.98 -3.00"
  setbb rloga0_tab_low "%rloga0_tab_low% -3.30 -3.32 -3.35 -3.38 -3.40 -3.43 -3.45 -3.47 -3.50 -3.53"
  setbb rloga0_tab_low "%rloga0_tab_low% -3.56 -3.59 -3.62 -3.65 -3.68 -3.70 -3.72 -3.74 -3.77 -3.79"
  setbb rloga0_tab_low "%rloga0_tab_low% -3.81 -3.83 -3.85 -3.88 -3.92 -3.94 -3.97 -3.98 -4.00 -4.02"
  setbb rloga0_tab_low "%rloga0_tab_low% -3.03 -3.08 -3.10 -3.12 -3.15 -3.19 -3.21 -3.23 -3.28 -3.29"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.05 -4.08 -4.10 -4.12 -4.15 -4.17 -4.20 -4.22 -4.24 -4.26"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.28 -4.30 -4.32 -4.34 -4.36 -4.38 -4.40 -4.42 -4.44 -4.46"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.48 -4.50 -4.51 -4.52 -4.54 -4.56 -4.57 -4.59 -4.61 -4.62"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.63 -4.64 -4.66 -4.68 -4.69 -4.70 -4.71 -4.72 -4.73 -4.74"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.75 -4.76 -4.77 -4.78 -4.79 -4.80 -4.81 -4.82 -4.83 -4.84"
  setbb rloga0_tab_low "%rloga0_tab_low% -4.85 -4.86 -4.87 -4.88 -4.89 -4.90 -4.91 -4.92 -4.93  junk"
*
  setbb rloga0_tab_hi                 "-1.40 -1.50 -1.60 -1.70 -1.90 -2.10 -2.32 -2.43 -2.54 -2.63"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -2.70 -2.77 -2.79 -2.83 -2.87 -2.90 -2.94 -2.96 -2.98 -3.00"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -3.30 -3.32 -3.35 -3.38 -3.40 -3.43 -3.45 -3.47 -3.50 -3.53"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -3.56 -3.59 -3.62 -3.65 -3.68 -3.70 -3.72 -3.74 -3.77 -3.79"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -3.81 -3.83 -3.85 -3.88 -3.92 -3.94 -3.97 -3.98 -4.00 -4.02"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -3.03 -3.08 -3.10 -3.12 -3.15 -3.19 -3.21 -3.23 -3.28 -3.29"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.05 -4.08 -4.10 -4.12 -4.15 -4.17 -4.20 -4.22 -4.24 -4.26"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.28 -4.30 -4.32 -4.34 -4.36 -4.38 -4.40 -4.42 -4.44 -4.46"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.48 -4.50 -4.51 -4.52 -4.54 -4.56 -4.57 -4.59 -4.61 -4.62"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.63 -4.64 -4.66 -4.68 -4.69 -4.70 -4.71 -4.72 -4.73 -4.74"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.75 -4.76 -4.77 -4.78 -4.79 -4.80 -4.81 -4.82 -4.83 -4.84"
  setbb rloga0_tab_hi "%rloga0_tab_hi% -4.85 -4.86 -4.87 -4.88 -4.89 -4.90 -4.91 -4.92 -4.93 -4.94 junk"
*
 setbb delta_local %dist%
 do inum from 1 to 120 by 1
   if %delta_local% gt 600
     message "distance greater than 600: not a LOCAL magnitude " 
     setbb alogazero 999
     break
   endif
   setbb r1 "(before ' ' "%r_tab_low%")"
   setbb r2 "(before ' ' "%r_tab_hi%")"
   getbb r1
   getbb r2
   setbb ralog1 "(before ' ' "%rloga0_tab_low%")"
   setbb ralog2 "(before ' ' "%rloga0_tab_hi%")"
   getbb ralog1
   getbb ralog2
   setbb r_tab_low "(after ' ' "%r_tab_low%)"
   setbb r_tab_hi "(after ' ' "%r_tab_hi%)"
   setbb rloga0_tab_low "(after ' ' "%rloga0_tab_low%)"
   setbb rloga0_tab_hi "(after ' ' "%rloga0_tab_hi%)"
*
   if %delta_local% ge %r1% 
    if %delta_local% lt %r2%
*                          0 1 2 3                      3   3                     3 2   2             2 1            0
     evaluate to alogazero ( ( ( ( %delta_local% - %r1% ) * ( %ralog2% - %ralog1% ) ) / ( %r2% - %r1% ) ) + %ralog1% )
     break
    endif
   endif
  enddo
*
* evaluate to alogazero ( %delta_local% -r(i-1))*(rloga0(i)-rloga0(i-1))/(r(i)-r(i-1)) + rloga0(i-1)
*
  rtr
  transfer from vel to wa freq .05 .1  15.0 20.0
*
* Remember that since I start from m/s, now the W/A time history is in m 
* (displacement), and needs to be converted in mm...
*
  if %alogazero% ne 999 
   mul 1.e+03
   w tmp1.wa tmp2.wa
   evaluate to cut1 ( %t0% - 3.0 )
   evaluate to cut2 ( %t0% + 60.0 )
   cut %cut1% %cut2%
   r tmp1.wa tmp2.wa
   rtr
   rmean
   cut off
   qdp off
*   width 1
*   win 3 x .01 .4 y .10 .40
*   xdiv nice
*   xdiv p off
*   beginwindow 3
*   xlim off
*   title "Waveforms for the Automatic Computation of Local Magnitude @(Synthetic Wood-Anderson@)" l t s l
*   ylabel "Displacement on Wood-Anderson @(mm@)" l l s m
*   xlabel "Time @(s@)" l b s m
*   fileid off
*   plabel 1 'Station: %statname% ' position .22 .80 s m
*   p1
*   plabel 1 off
*   title off 
*   xlabel off 
*   ylabel off
*
* Now compute max value for each component and average them out...
*
   evaluate to maxmin ( abs ( &1,depmin ) )
   evaluate to maxmax ( abs ( &1,depmax ) )
*
   if %maxmin% gt %maxmax%
    setbb peak_val %maxmin%
   else
    setbb peak_val %maxmax%
   endif
*
   evaluate to loc_mag_1 ( alog10 ( %peak_val% ) )
   evaluate to loc_mag_1 ( %loc_mag_1% - %alogazero% )
*
   evaluate to maxmin ( abs ( &2,depmin ) )
   evaluate to maxmax ( abs ( &2,depmax ) )
*
   if %maxmin% gt %maxmax%
    setbb peak_val %maxmin%
   else
    setbb peak_val %maxmax%
   endif
*
   evaluate to loc_mag_2 ( alog10 ( %peak_val% ) )
   evaluate to loc_mag_2 ( %loc_mag_2% - %alogazero% )
*
   evaluate to loc_mag ( ( %loc_mag_1% + %loc_mag_2% ) / 2 )
  endif
  if %alogazero% eq 999 
   setbb loc_mag UNDEFINED
  endif
 endif
*
* Ifs on magnitude computation finished here...
*
*******************************
 r %statname%.r %statname%.t
 ch lcalda true
 ch lovrok true
 rtr

* For a general situation, when the local magnitude is not available from the origin network,
* do setbb magn &1,user8&. Remember that (in this case) local magnitude is written in field user4 
*
 ch user8 %loc_mag%  
 wh over
 setbb magn &1,user4&
 setbb ml %magn%
 if %magn% eq UNDEFINED
  setbb magn &1,user8&
 endif
 if %magn% ge 5
  message " no hipass filter large events... " 
 endif
 if %magn% ge 4
  if %magn% lt 5 
   hp bu n 2 p 2 c .2 
  endif
 endif
 if %magn% ge 3
  if %magn% lt 4 
   hp bu n 2 p 2 c .3 
  endif
 endif
 if %magn% ge 2
  if %magn% lt 3 
   hp bu n 2 p 2 c 0.5
  endif
 endif
 if %magn% lt 2 
  hp bu n 2 p 2 c 1.5
 endif
 w over
enddo
*
* Now starts the spectral analysis for the magnitude computation...
*
do file wild *_tmp.r
 message 'Processing file $file ...'
 setbb statname (BEFORE ".r" $file ) 
 sc echo %statname% > statname.out
 setbb inputfile_t %statname%.t
 setbb inputfile_r %statname%.r
*
* get the header information from the SAC file
*
 r %inputfile_t% %inputfile_r%
 setbb xmin &1,o
 evaluate to xmax ( &1,t0 + 30 )
 xlim %xmin% %xmax% 
 qdp off
 width 1
* win 1 x .01 .4 y .5 .99
* win 2 x .41 .99 y .5 .99
* xdiv nice
* xdiv p off
*** qdp term 1000
* beginwindow 1
* fileid off
* title "Station: %statname% - Radial @(bottom@) and Transverse @(top@) Components" l t s l
* ylabel "Corrected Ground Velocity" l l s m
* xlabel "Time @(s@)" l b s m
* p1
* xlim off
* ylim off
* title off
************************************************************************
*			HERE I SET THE VARIABLES I'LL NEED LATER
 setbb dist &1,dist&
 if %dist% eq UNDEFINED
  message 'HEADER FIELD dist UNDEFINED FOR FILE %inputfile% '
  break
 endif
 setbb baz &1,baz&
 setbb depth &1,evdp
*
* evdp in m already; dist in km in the original header field...
*
 evaluate to dist ( &1,dist * 1000 )
 evaluate to dist2 ( %dist% ** 2)
 evaluate to dist ( &1,dist )
 evaluate to h2 ( %depth% ** 2 )
 evaluate to hypo_dist ( %dist2% + %h2% )
 evaluate to hypo_dist ( sqrt ( %hypo_dist% ) )
 setbb gcarc &1,gcarc&
 setbb evla &1,evla&
 setbb evlo &1,evlo&
 setbb begin &1,b&
 setbb t0 &1,t0&
 setbb magn &1,user4&
 setbb newmagn &1,user8&
 if %magn% eq UNDEFINED
*  setbb ml %magn%
  setbb ml UNAVAILABLE
  setbb magn %newmagn%
 endif
 setbb gcarc &1,gcarc
 evaluate to depth ( &1,evdp / 1000 )
 setbb a &1,a&
 setbb o &1,o&
 setbb vel1 4.0
 setbb vel2 1.0
 if %dist% ge 200.0
  setbb vel1 4.2
  setbb vel2 1.0
 endif
 if %dist% ge 300.0
  setbb vel1 4.5
  setbb vel2 1.0
 endif
*
* Here is the 2-Hz duration function
*
setbb dist_1 20.000 ; setbb dur_1 4.405 
setbb dist_2 30.000 ; setbb dur_2 4.844 
setbb dist_3 40.000 ; setbb dur_3 5.284 
setbb dist_4 50.000 ; setbb dur_4 5.724 
setbb dist_5 60.000 ; setbb dur_5 6.163 
setbb dist_6 70.000 ; setbb dur_6 6.603 
setbb dist_7 80.000 ; setbb dur_7 7.134 
setbb dist_8 90.000 ; setbb dur_8 7.666 
setbb dist_9 100.000 ; setbb dur_9 8.197 
setbb dist_10 110.000 ; setbb dur_10 8.729 
setbb dist_11 120.000 ; setbb dur_11 10.367 
setbb dist_12 140.000 ; setbb dur_12 12.005 
setbb dist_13 180.000 ; setbb dur_13 15.281 
setbb dist_14 200.000 ; setbb dur_14 16.918 

********************************************
 if %hypo_dist% ge %dist_14%
  evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_13% ) / ( %dist_14% - %dist_13% ) ) * %dur_14% ) 
 endif
 if %hypo_dist% ge %dist_13%
  if %hypo_dist% lt %dist_14% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_13% ) / ( %dist_14% - %dist_13% ) ) * %dur_13% ) 
  endif
 endif
 if %hypo_dist% ge %dist_12%
  if %hypo_dist% lt %dist_13% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_12% ) / ( %dist_13% - %dist_12% ) ) * %dur_12% ) 
  endif
 endif
 if %hypo_dist% ge %dist_11%
  if %hypo_dist% lt %dist_12% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_11% ) / ( %dist_12% - %dist_11% ) ) * %dur_11% ) 
  endif
 endif
 if %hypo_dist% ge %dist_10%
  if %hypo_dist% lt %dist_11% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_10% ) / ( %dist_11% - %dist_10% ) ) * %dur_10% ) 
  endif
 endif
 if %hypo_dist% ge %dist_9%
  if %hypo_dist% lt %dist_10% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_9% ) / ( %dist_10% - %dist_9% ) ) * %dur_9% ) 
  endif
 endif
 if %hypo_dist% ge %dist_8%
  if %hypo_dist% lt %dist_9% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_8% ) / ( %dist_9% - %dist_8% ) ) * %dur_8% ) 
  endif
 endif
 if %hypo_dist% ge %dist_7%
  if %hypo_dist% lt %dist_8% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_7% ) / ( %dist_8% - %dist_7% ) ) * %dur_7% ) 
  endif
 endif
 if %hypo_dist% ge %dist_6%
  if %hypo_dist% lt %dist_7% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_6% ) / ( %dist_7% - %dist_6% ) ) * %dur_6% ) 
  endif
 endif
 if %hypo_dist% ge %dist_5%
  if %hypo_dist% lt %dist_6% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_5% ) / ( %dist_6% - %dist_5% ) ) * %dur_5% ) 
  endif
 endif
 if %hypo_dist% ge %dist_4%
  if %hypo_dist% lt %dist_5% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_4% ) / ( %dist_5% - %dist_4% ) ) * %dur_4% ) 
  endif
 endif
 if %hypo_dist% ge %dist_3%
  if %hypo_dist% lt %dist_4% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_3% ) / ( %dist_4% - %dist_3% ) ) * %dur_3% ) 
  endif
 endif
 if %hypo_dist% ge %dist_2%
  if %hypo_dist% lt %dist_3% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_2% ) / ( %dist_3% - %dist_2% ) ) * %dur_2% ) 
  endif
 endif
 if %hypo_dist% ge %dist_1%
  if %hypo_dist% lt %dist_2% 
   evaluate to t_dist_corr ( ( ( %hypo_dist% - %dist_1% ) / ( %dist_2% - %dist_1% ) ) * %dur_1% ) 
  endif
 endif
 if %hypo_dist% lt %dist_1%
  evaluate to t_dist_corr ( ( %hypo_dist% / %dist_1% ) * %dur_1% )
 endif
 evaluate to t_dist_corr ( %t_dist_corr% * 0.5 )
*
* Duration is taken as a function of magnitude, too...
*
 if %magn% ge 5
  evaluate to cut1 ( %begin% + %t0% - 2 )
  evaluate to cut2 ( %begin% + %t0% + 12 + %t_dist_corr% )
 endif
 if %magn% ge 4
  if %magn% lt 5 
   evaluate to cut1 ( %begin% + %t0% - 2 )
   evaluate to cut2 ( %begin% + %t0% + 9 + %t_dist_corr% )
  endif
 endif
 if %magn% ge 3
  if %magn% lt 4 
   evaluate to cut1 ( %begin% + %t0% - 2 )
   evaluate to cut2 ( %begin% + %t0% + 7 + %t_dist_corr% )
  endif
 endif
 if %magn% ge 2
  if %magn% lt 3 
   evaluate to cut1 ( %begin% + %t0% - 2 )
   evaluate to cut2 ( %begin% + %t0% + 4 + %t_dist_corr% )
  endif
 endif
 if %magn% lt 2 
  evaluate to cut1 ( %begin% + %t0% - 2 )
  evaluate to cut2 ( %begin% + %t0% + 3 + %t_dist_corr% )
 endif
*
*
*
 evaluate to time_window ( %cut2% - %cut1% )
 evaluate to noise1 (add %begin% 1.0)
 evaluate to noise2 (add %noise1% 6.0)
*
* cut the direct wave & compute spectrum
*
 if %manual% EQ no
  r %inputfile_t% %inputfile_r%
  setbb end1 &1,e
  if %cut1% gt %end1%
   message " start cut greater than file end... break loop for files %inputfile_t% %inputfile_r% "
   break
  else
   cut %cut1% %cut2%
   r %inputfile_t% %inputfile_r%
   cut off
   setbb peak_amp_1 &1,depmax&
   setbb peak_amp_2 &2,depmax&
   rtr
   taper 
   fft
   wsp amph %inputfile_t% %inputfile_r%
  endif
 elseif %manual% EQ yes
  r %inputfile_r% 
  ppk
  setbb cut1 &1,t5&
  setbb cut2 &1,t6&
  cut %cut1% %cut2%
  r %inputfile_r% %inputfile_t%
  cut off
  setbb peak_amp_1 &1,depmax&
  setbb peak_amp_2 &2,depmax&
  rtr
  taper
  fft
  wsp amph %inputfile_r% %inputfile_t%
 endif
*
* cut out a noise sample & get the mean amplitude, "noise_amp"
*
 cut %noise1% %noise2%
 evaluate to window_length ( %noise2% - %noise1% )
 r %inputfile_r% %inputfile_t%
 cut off
 sqr
 int
 div %window_length%
 sqrt
*
* RMS is the max of the integrated a**2(t) dt /T 
*
 setbb noise_amp_1 &1,depmax&
 setbb noise_amp_2 &2,depmax&
*
* Define the signal-to-noise ratio, "SNR"
*
 evaluate to SNR_R (%peak_amp_1% / %noise_amp_1%)
 evaluate to SNR_T (%peak_amp_2% / %noise_amp_2%)
*
* if SNR is larger than 1.0, then plot; check both SNR_R & SNR_T...
*
* RELAX the S/N condition... Needs to be more sophisticated
*
if %SNR_R% ge 0.1
  if %SNR_T% ge 0.1
   r %inputfile_r%.am
   ch iftype IXY
   smooth mean h 1
   sqr
   w a
   r %inputfile_t%.am
   ch iftype IXY
   smooth mean h 1
   sqr
   addf a
   sqrt
   w %statname%.am
*
* Apply correction...
*
   r %statname%.am
   setbb npoints &1,npts
   setbb delta &1,delta
*
* Correct for crustal attenuation and kappa...
*
********  BAY AREA PARAMETERS  *******
*
   setbb kappa 0.1
   setbb Q0 160
   setbb eta 0.45
   setbb pi 3.14159
   evaluate to const ( %Q0% * 3500 )
   evaluate to const ( %pi% / %const% )
   evaluate to const ( %const% * %hypo_dist% )
   evaluate to power ( 1.0 - %eta% )
   fg line 1 0 DELTA %delta% NPTS %npoints% BEGIN %delta%
*   sqrt
   log10
   mul %power%
   exp10
   mul %const% 
   w tmp_qcomp
*
* Here is the contribution of kappa...
*
   evaluate to const ( %pi% * %kappa% )
   fg line %const% 0 DELTA %delta% NPTS %npoints% BEGIN %delta%
*
* Put everything together...
*
   w tmp_kcomp
   r tmp_kcomp
   addf tmp_qcomp
   exp
   w attenuation 
   r %statname%.am
*
* Correct for a complex geometrical spreading... 
*
* The following constants are specific for the region...
*
* Remember hypo_dist in m...
*
   setbb crossover_dist_1 500
   setbb eta1 -1.0
   setbb crossover_dist_2 50000
   setbb eta2  0.0
   setbb crossover_dist_3 70000
   setbb eta3  0.0
   setbb crossover_dist_4 90000
   setbb eta4 -0.6
   setbb crossover_dist_5 100000
   setbb eta5 -0.6
*
   if %hypo_dist% gt %crossover_dist_1%
    if %hypo_dist% gt %crossover_dist_5%
     evaluate to corr_1 ( %crossover_dist_2% )
     evaluate to corr_1 ( %corr_1% ** %eta1% )
     evaluate to corr_2 ( %crossover_dist_3% / %crossover_dist_2% )
     evaluate to corr_2 ( %corr_2% ** %eta2% )
     evaluate to corr_3 ( %crossover_dist_4% / %crossover_dist_3% )
     evaluate to corr_3 ( %corr_3% ** %eta3% )
     evaluate to corr_4 ( %crossover_dist_5% / %crossover_dist_4% )
     evaluate to corr_4 ( %corr_4% ** %eta4% )
     evaluate to corr_5 ( %hypo_dist% / %crossover_dist_5% )
     evaluate to corr_5 ( %corr_5% ** %eta5% )
     evaluate to dist_correction ( %corr_1% * %corr_2% * %corr_3% * %corr_4% * %corr_5% )
    endif
*
    if %hypo_dist% gt %crossover_dist_4%
     evaluate to corr_1 ( %crossover_dist_2% )
     evaluate to corr_1 ( %corr_1% ** %eta1% )
     evaluate to corr_2 ( %crossover_dist_3% / %crossover_dist_2% )
     evaluate to corr_2 ( %corr_2% ** %eta2% )
     evaluate to corr_3 ( %crossover_dist_4% / %crossover_dist_3% )
     evaluate to corr_3 ( %corr_3% ** %eta3% )
     evaluate to corr_4 ( %hypo_dist% / %crossover_dist_4% )
     evaluate to corr_4 ( %corr_4% ** %eta4% )
     evaluate to dist_correction ( %corr_1% * %corr_2% * %corr_3% * %corr_4% )
    endif
*
    if %hypo_dist% gt %crossover_dist_3%
     if %hypo_dist% le %crossover_dist_4%
      evaluate to corr_1 ( %crossover_dist_2% )
      evaluate to corr_1 ( %corr_1% ** %eta1% )
      evaluate to corr_2 ( %crossover_dist_3% / %crossover_dist_2% )
      evaluate to corr_2 ( %corr_2% ** %eta2% )
      evaluate to corr_3 ( %hypo_dist% / %crossover_dist_3% )
      evaluate to corr_3 ( %corr_3% ** %eta3% )
      evaluate to dist_correction ( %corr_1% * %corr_2% * %corr_3% )
     endif
    endif
*
    if %hypo_dist% gt %crossover_dist_2%
     if %hypo_dist% le %crossover_dist_3%
      evaluate to corr_1 ( %crossover_dist_2% )
      evaluate to corr_1 ( %corr_1% ** %eta1% )
      evaluate to corr_2 ( %hypo_dist% / %crossover_dist_2% )
      evaluate to corr_2 ( %corr_2% ** %eta2% )
      evaluate to dist_correction ( %corr_1% * %corr_2% )
     endif
    endif
*
    if %hypo_dist% le %crossover_dist_2%
     evaluate to corr_1 ( %hypo_dist% )
     evaluate to corr_1 ( %corr_1% ** %eta1% )
     evaluate to dist_correction ( %corr_1% )
    endif
   endif
   div %dist_correction%
*
* I multiply the spectra for the 'attenuation' function to remove
* the effects of anelastic attenuation...
*
   mulf attenuation
   w corrected_v_%statname%.am
*
* Compute corner frequency by using Andrews' (1986) relationship
*
   r corrected_v_%statname%.am
   setbb bb &1,delta&
   cut %bb% 20
   r corrected_v_%statname%.am
   cut off
   setbb n &1,npts&
   sqr
   int
   setbb v2 &1,depmax
   evaluate to v2 ( 2.0 * ( %v2% - &1,depmin ) )
   fg line 1 0 delta %bb% npts %n% begin %bb%
   mul 2
   mul 3.14159 
   w twopif
   cut %bb% 20
   r corrected_v_%statname%.am
   cut off
   divf twopif
   sqr
   int
   setbb d2 &1,depmax
   evaluate to d2 ( 2.0 * ( %d2% - &1,depmin ) )
*
* Transform v2 and d2 in cgs units...
*
   evaluate to v2 ( %v2% * 1.0e+04 )
   evaluate to d2 ( %d2% * 1.0e+04 )
*
   evaluate to f0 ( %v2% / %d2% )
   evaluate to f0 ( sqrt ( %f0% ) )
   evaluate to f0 ( %f0% / %pi% )
   evaluate to f0 ( %f0% / 2 )
   setbb corner_freq %f0%
*
   evaluate to min_freq ( %corner_freq% / 2.0 )
   evaluate to max_freq ( %corner_freq% * 1.0 )
*
* Get RMS at low frequency...
*
   r corrected_v_%statname%.am
   setbb bb &1,delta&
   cut %bb% 30
   r corrected_v_%statname%.am
   cut off
   setbb bb &1,delta&
   setbb d &1,delta&
   setbb n &1,npts&
   fg line 1 0 delta %bb% npts %n% begin %bb%
   mul 2
   mul 3.14159
   w twopif
   cut %bb% 30
   r corrected_v_%statname%.am
   cut off
   divf twopif
   w corrected_d_%statname%.am
   cut %min_freq% %max_freq%
   evaluate to freq_window ( %max_freq% - %min_freq% )
   r corrected_d_%statname%.am
   cut off
   sqr
   int
   div %freq_window%
   sqrt
   setbb omega0 &1,depmax
*----
* Parameters to be used in the Brune's spectral model
*----
*       ff = partition of energy
*       V  = free surface factor
*       df = frequency increment
*       fmin,fmax = min and max frequencies in the spectrum; no Hanks' fmax
*       beta = shear-wave velocity (crustal average) (cm/s)
*       rho  = density (crustal average) (g/cm^3)
*       Rtp  = average radiation pattern
*----
*   setbb ff 0.707
   setbb ff 0.8
   setbb V 2.0
   setbb beta 3.5e+05
   setbb rho 2.8
   setbb Rtp 0.8
   setbb f0 1.0
   evaluate to beta3 ( %beta% ** 3 )
   setbb pi 3.14159
   evaluate to omega0dynecm ( %omega0% * 1.0e+02 ) 
*
* Needs to translate in cm the distance correction (it was m)....
*
   evaluate to omega0dynecm ( %omega0dynecm% * 1.0e+02 )
*
   evaluate to M0 ( %omega0dynecm% * 4 * %pi% * %rho% * %beta3% ) 
*   evaluate to den ( 0.55 * 2 * 0.707 )
*
* Need to remove 1/sqrt(2) from den, since I took the 
* geometric mean of the observed horizontal spectra...
*
   evaluate to den ( %Rtp% * 2 )
   evaluate to M0 ( %M0% / %den% )
   evaluate to MW ( ALOG10 ( %M0% ) )
   evaluate to MW ( %MW% / 1.5  )
   evaluate to MW ( %MW% - 10.73 )
*
* Evaluate stress drop by using Andrews (1986) relationships... (Ds_2 is the Brune stress drop)
*
* Ds_2 = (k'/(4*pi**3 k")) * (I_v**5/I_d**3)**(1/4)   Here I_v and I_d are indicated by %v2% and by %d2%, respectively
* k'=(7/16)*(2*pi/(2.34*beta))**3 ; k"= Rtp*ff/(4*pi*beta**3)
*
* Snoke (1987) observed that a more stable estimate of the Brune stress drop may be obtained by using the following expression:
* (Snoke, J.A. (1987). Stable determination of (Brune) stress drops, Bull. Seism. Soc. Am., 77, 530-538.)
* Ds_3=(k'/(2*pi**3*k"**2))*(I_v/M0), that is ~= 3.4 Ds_a   (Ds_a=apparent stress drop)...
* Iv must be translated to cm^2/s^2 (it was computed in  m^2/s^2)!!!!!
*
*                          0 1            1 2 3          3   3               3 2      1 0
   evaluate to kappa_prime ( ( 7 / 16 ) * ( ( ( 2 * %pi% ) / ( 2.34 * %beta% ) ) ** 3 ) )
*                            0 1             1   1                            1 0
   evaluate to kappa_d_prime ( ( %Rtp% * %V% ) / ( 4 * %pi% * %rho% * %beta3% ) )
**
** Andrews (DSIGMA in dyne/cm^2)...
**                     
*   evaluate to DSIGMA ( %kappa_prime% / ( 4 * ( %pi% **3 ) * %kappa_d_prime%  ) )
*   evaluate to aexp ( 5 / 4 )
*   evaluate to bexp ( 3 / 4 )
*   evaluate to aaa ( ( ( ( %v2% * 1.0e+04 ) ** %aexp% ) / ( ( %d2% * 1.0e+04 )  ** %bexp% ) ) )
*   evaluate to DSIGMA ( %DSIGMA% * %aaa% )
*
* Snoke (DSIGMA in dyne/cm^2)...
*
*                     0                 1     2          2   2                      2 1 0  
   evaluate to DSIGMA ( %kappa_prime% / ( 2 * ( %pi% ** 3 ) * ( %kappa_d_prime% ** 2 ) ) )  
*                     0            1             1 0
   evaluate to DSIGMA ( %DSIGMA% * ( ( %v2% * 1.0e+04 )  / %M0% ) )  
*
*
* Transform DSIGMA in MPa... (1 dyna/cm^2 = 10^-5 N /10^-4 m^2 = 10^-1 N/m^2 = 10^-1 Pa = 10^-7 MPa)
*
  evaluate to DSIGMA ( %DSIGMA% * 1.e-07 )
*
* Apparent stress = Dsigma/3.4 
*
  evaluate to app_stress ( %DSIGMA% / 3.4 ) 
*
* Compute M0 in N-m ...
*
   evaluate to M0 ( %M0% * 1.0e-07 )
*******************************************
***
*** DSIGMA, computed by using Brune's model formula; in DSIGMA, beta in km/s...
***
**   setbb beta 3.5
**   evaluate to DSIGMA ( %corner_freq% / 4.9e+06 )
**   evaluate to DSIGMA ( %DSIGMA% / %beta% )
**   evaluate to DSIGMA ( %DSIGMA% ** 3 )
**   evaluate to DSIGMA ( %DSIGMA% * %M0% )
*****************
***
*** DSIGMA in bars, now translate in MPa...
***
**   evaluate to DSIGMA ( %DSIGMA% / 10 )
**
*******************************************
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' moment magnitude %MW% ' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_mw
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' low-frequency spectral amplitude %omega0% ' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_omega0
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message " %M0% "
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_m0
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' corner frequency %corner_freq% ' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_fc
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' stress drop %DSIGMA% ' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_sigma
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' local magnitude %newmagn%  @(%woodanderson%@)' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> tmp_newm
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' network local magnitude %ml%  @(given@)' 
   TRANSCRIPT CLOSE
   sc cat ./tmp > tmp_ml
   sc rm ./tmp
   TRANSCRIPT CREATE FILE tmp CONTENTS OUTPUT
   message ' %statname%: fc= %corner_freq% M0= %M0% Mw= %MW% Dsigma= %DSIGMA% Ml= %ml% Automatic Ml= %newmagn% @(%woodanderson%@)' 
   TRANSCRIPT CLOSE
   sc cat ./tmp >> Mw_output
   sc rm ./tmp
*
* omega0 is the RMS spectral value at low-freq
*
   evaluate to npoints ( %corner_freq% / %delta% )
   fg line 0 %omega0% DELTA %delta% BEGIN %min_freq%
   ch iftype IAMPH
   wsp am omega0_line
*   w omega0_line
   ch t1 %corner_freq%
*   plabel 1 'Corner Freq=%corner_freq% Hz ' position .22 .80 s m
*   plabel 2 'M0=%M0% N-m ' position .22 .75 s m
*   plabel 3 'Automatic Moment Magnitude=%MW% ' position .22 .50 s m
*   plabel 4 'Network Local Magnitude= %ml% @(given@) ' position .22 .45 s m
*   plabel 5 'Automatic Local Magnitude= %newmagn% @(%woodanderson%@)' position .22 .40 s m
*   plabel 6 'Brune Stress Drop: Dsigma=%DSIGMA% MPa ' position .22 .35 s m
*   plabel 7 'Time Window =%time_window% s ' position .22 .30 s m
*   plabel 8 'Starting Time =%cut1% s ' position .22 .25 s m
*   title " Station: %statname% - Corrected Displacement Spectrum " l t s l
*   xlabel "Frequency @(Hz@)" l b s m
*   ylabel "Displacement Spectrum @(m * s@)" l l s m
*   fileid off
*   width 1
*   win 1 x .01 .5 y .5 .99
*   xdiv nice
*   xdiv p off
**   qdp term 1000
*   beginwindow 2
   cut .005 50
   r corrected_d_%statname%.am
   cut off
   evaluate to ymax ( &1,depmax * 10.0 )
   setbb ymin ( &1,depmin )
   ylim %ymin% %ymax%
   ch iftype IAMPH
   xlim .01 30
*   loglog
*   beginframe
**   message " cut %min_freq% %max_freq% "
*   psp am
   cut %min_freq% %max_freq%
   r omega0_line.am
   cut off
   ch iftype IAMPH
*   width 2
*   psp am
*   plabel 1 off
*   plabel 2 off
*   plabel 3 off
*   plabel 4 off
*   plabel 5 off
*   plabel 6 off
*   plabel 7 off
*   plabel 8 off
*   title off
*   xlabel off
*   ylabel off
*   endframe
   ylim off 
   pause p 100
  else 
   message 'low S/N ratio... check data'
   sc rm -f %inputfile_t%.am %inputfile_r%.am
  endif
 else
  message 'low S/N ratio... check data'
  sc rm -f %inputfile_t%.am %inputfile_r%.am
 endif
 linlin
*
* Cleanup...
*
enddo
sc rm *_tmp.? core *.am *.ph *.wa 
echo off





*exit
*
* HERE ARE JUST EXAMPLES...
*
*fg line 0 0
*w zero
*r zero
*setbb count 0
*do file wild *.am
* addf $file
* evaluate to count (add %count% 1.0)
*enddo
*div %count%
*w average_hspectra.am
*r average_hspectra.am
*
*
*
*
*rtr
*setbb temp_slope %rtr_slp%
*r temp_GF
*rtr
*setbb gf_slope %rtr_slp%
**
*evaluate to diff_slope (abs (sub %gf_slope% %temp_slope%))
*evaluate to ratio_slope (%diff_slope% / (abs %gf_slope%))
*if %step% ge 1.0
*  if %ratio_slope% gt 0.45
*    getbb diff_slope
*    getbb temp_slope
*    getbb step
**    evaluate to end_window (sub %wayahead% %window_size%)
*    setbb end_window %ahead2%
*    setbb message 'Stopped by Slope %ratio_slope%' 
*    break
*  endif
*endif
*sc rm *_tmp.?