# Calculation and smoothing of the optimum initial 
# parameters of Gaussian beams
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# see  Zacek, K. (2006): Optimization of the Shape of Gaussian Beams. 
#      Stud. Geophys. Geod., 50, 346-366.
#
# This example shows calculation for shot 77.
# Programs matmul and invpts need MRAM=12000000 in file "ram.inc".

# Input files required
  chk.pl: "data/mar/" "mgb-crt.dat"
  chk.pl: "data/mar/" "mgp-wri.dat"
  chk.pl: "data/mar/" "mgb-m0.dat"
  chk.pl: "data/mar/" "mar-mod.dat"
  chk.pl: "data/mar/" "mgb-opt1.h"
  chk.pl: "data/mar/" "gb-ry.cal"
  chk.pl: "data/mar/" "gb-s.cal"
  chk.pl: "data/mar/" "gb-yd.cal"
  chk.pl: "data/mar/" "gb-tau.cal"
  chk.pl: "data/mar/" "mgbstep.ini"
  chk.pl: "data/mar/" "gb-sob11.dat"
  chk.pl: "data/mar/" "mgp-step.h"
  chk.pl: "data/mar/" "mgb-grd.h"
  chk.pl: "crt/"      "codep.dat"
  chk.pl: "crt/"      "writall.dat"
  chk.pl: "forms/"    "echo.pl"
  chk.pl: "forms/"    "append.pl"
  chk.pl: "forms/"    "div.cal"

# Calculation of the optimum initial parameters of Gaussian beams
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  echo.pl: "ISRC=77" "> gbo.tmp"
  append.pl: "gbo.tmp" "mgb-opt1.h"
  go.pl: "gbo.tmp" "gbo.out"

# Smoothing the optimum initial parameters of Gaussian beams
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Smoothing the parameter R0S - 1st iteration
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=0. O2=0. O3=0.
  CAL='gb-ry.cal' GRD1='r0s.out' GRD2='y0s.out'  GRD3='ry2r.out'
  grdcal:
  CAL='div.cal'   GRD1='r0s.out' GRD2='ry2r.out' GRD3='r0s1.out'
  grdcal:
  CAL='div.cal'   GRD1='y0s.out' GRD2='ry2r.out' GRD3='y0s1.out'
  grdcal:
  CAL='gb-s.cal'  GRD1='sigmar.out' GRD2='y0s.out' 
                  GRD3='ry2r.out'   GRD4='sigmar1.out'
  grdcal:
  
# Initial and updated models
  MODEL='mgb-m0.dat'  MODOUT='r0mod.out'  NEGPAR=1

# Calculating matrices for inversion (minimizing first derivatives)
  M1='m1.out'  M2='m2.out'  MODL2=' '  MODSOB='modsob.out'
  SOBOLEV='gb-sob11.dat'  SOBW01=1
  invsoft:

# Unit given standard deviation: ERRMUL=SQRT(N1*N2*N3)
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=1. O2=1. O3=1.
  CAL='gb-tau.cal' GRD1='tau.out' GRD2='tauint.out' 
  DTAU=0.200    # DTAU is 2xSTORE
  grdcal:
  GRD='r0s1.out'     KOLUMN=4 
  GRD1='sigmar1.out' KOLUMN1=5 
  GRD2='tauint.out'  KOLUMN2=1
  PTS='r0pts.out'
  grdpts:
  GRD=' ' 
  KOLUMN=4 KOLERR=5
  INDFUN=1  MPAR=1  POWERM=1  ERRMUL=211.962 
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  invpts:

# Matrix operations
  MATIN1='dm1.out'                       MATOUT='dm2.out'
  MATFUN='inv'
  matfun:
  MATIN1='dm2.out'   MATIN2='gm1.out'    MATOUT='dm2gm1.out'
  SYMMETRY=      MATT1=  MATT2=1
  matmul:
  MATIN1='gm1.out'   MATIN2='dm2gm1.out' MATOUT='sm1.out'
  SYMMETRY='sym' MATT1=  MATT2=
  matmul:
  MATIN1='sm1.out'   MATIN2='modsob.out' MATOUT='sm2.out'
  COEF1=1.  COEF2=400.   # range of Im(N0) approximately one order
 #COEF1=1.  COEF2=41.    # range of Im(N0) approximately two orders
  matlin:
  MATIN1='sm2.out'                       MATOUT='sm3.out'
  matinv:
  MATIN1='dm2.out'   MATIN2='gm2.out'    MATOUT='gm3.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='gm1.out'   MATIN2='gm3.out'    MATOUT='gm4.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='sm3.out'   MATIN2='gm4.out'    MATOUT='gm5.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'
  modmod:
#
# Calculating gridded values of the parameter R0S
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=1. O2=1. O3=1.
  MODEL='r0mod.out'  MPAR=1  ICB=' '  VEL='r0m1.out'
  grid:

# Smoothing the parameter Y0S - 1st iteration
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Data to be fit
  CAL='gb-yd.cal' 
  GRD1='y0s1.out' GRD2='r0m1.out' GRD3='r0s1.out' GRD4='y0d1.out'
  grdcal:
  
# Initial and updated models
  MODEL='mgb-m0.dat'  MODOUT='y0mod.out'  NEGPAR=1

# Calculating matrices for inversion (minimizing first derivatives)
  M1='m1.out'  M2='m2.out'  MODL2=' '  MODSOB='modsob.out'
  SOBOLEV='gb-sob11.dat'  SOBW01=1
  invsoft:

# Unit given standard deviation: ERRMUL=SQRT(N1*N2*N3)
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=1. O2=1. O3=1.
  CAL='gb-tau.cal' GRD1='tau.out' GRD2='tauint.out' 
  DTAU=0.200    # DTAU is 2xSTORE
  grdcal:
  GRD='y0d1.out'     KOLUMN=4  
  GRD1='sigmar1.out' KOLUMN1=5 
  GRD2='tauint.out'  KOLUMN2=1
  PTS='y0pts.out'
  grdpts:
  GRD=' ' 
  KOLUMN=4 KOLERR=5
  INDFUN=1  MPAR=1  POWERM=1  ERRMUL=211.962 
  GM1='gm1.out'  GM2='gm2.out'  GM3=' '  DM1='dm1.out'
  invpts:

# Matrix operations
  MATIN1='dm1.out'                       MATOUT='dm2.out'
  MATFUN='inv'
  matfun:
  MATIN1='dm2.out'   MATIN2='gm1.out'    MATOUT='dm2gm1.out'
  SYMMETRY=      MATT1=  MATT2=1
  matmul:
  MATIN1='gm1.out'   MATIN2='dm2gm1.out' MATOUT='sm1.out'
  SYMMETRY='sym' MATT1=  MATT2=
  matmul:
  MATIN1='sm1.out'   MATIN2='modsob.out' MATOUT='sm2.out'
  COEF1=1.  COEF2=400.   # range of Im(N0) approximately one order
 #COEF1=1.  COEF2=41.    # range of Im(N0) approximately two orders
  matlin:
  MATIN1='sm2.out'                       MATOUT='sm3.out'
  matinv:
  MATIN1='dm2.out'   MATIN2='gm2.out'    MATOUT='gm3.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='gm1.out'   MATIN2='gm3.out'    MATOUT='gm4.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:
  MATIN1='sm3.out'   MATIN2='gm4.out'    MATOUT='gm5.out'
  SYMMETRY=' '   MATT1=  MATT2=
  matmul:

# Updating the model
  M1='m1.out'  MODNEW='gm5.out'
  modmod:

# Calculating gridded values of the parameter Y0S
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=1. O2=1. O3=1.
  MODEL='y0mod.out'  MPAR=1  ICB=' '  VEL='y0m1.out'
  grid:

# 
  CAL='gb-ry.cal' GRD1='r0m1.out' GRD2='y0m1.out' GRD3='ry2.out'
  grdcal:
  CAL='div.cal'   GRD1='r0m1.out' GRD2='ry2.out'  GRD3='r0m.out'
  grdcal:
  CAL='div.cal'   GRD1='y0m1.out' GRD2='ry2.out'  GRD3='y0m.out'
  grdcal:

# Plotting optimized parameters of Gaussian beams
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  GRD='r0s.out'      KOLUMN=4 
  GRD1='sigmar.out'  KOLUMN1=5 
  GRD2='tauint.out'  KOLUMN2=1
  PTS='r0pts1.out'
  grdpts:
  GRD='y0s.out'      KOLUMN=4 
  GRD1='sigmar.out'  KOLUMN1=5 
  GRD2='tauint.out'  KOLUMN2=1
  PTS='y0pts1.out'
  grdpts:
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1. O1=0. O2=0. O3=0.
  CREF=0.166667  VMIN=-99999999999.  VMAX=0.5E-06
  R=0.  G=0.  B=0.

# Before smoothing and interpolation
  VREF=0.  VCIRC=
  PTS='r0pts1.out'  PS='r0s001.ps'
  ptsps:
  
  VREF=0.19E-08  VCIRC=0.6E-06
  PTS='y0pts1.out'  PS='y0s001.ps'
  ptsps:
  
# After smoothing and interpolation
  VREF=0.  VCIRC=
  GRD='r0m.out'  PS='r0m.ps'
  grdps:
  VREF=0.19E-08  VCIRC=0.6E-06
  GRD='y0m.out'  PS='y0m.ps'
  grdps:

# Calculation of discretization steps for decomposition and migration 
  echo.pl: "ISRC=77" "> mgp-step.tmp"
  append.pl: "mgp-step.tmp" "mgp-step.h"
  go.pl: "mgp-step.tmp" "mgp-step.out"