# Smoothing the optimum initial parameters of Gaussian beams
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
# Deleting old files
  del.pl: "r0m1.out"    
  del.pl: "y0m1.out"    
  del.pl: "r0mnew.out"    
  del.pl: "y0mnew.out"    
  del.pl: "ry2.out"    

# Grid dimensions
  N1=26 N2=48 N3=36  D1=1. D2=1. D3=1.  O1=1. O2=1. O3=1.

# Smoothing the parameter R0S - 1st iteration
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 
  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)
  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=2.3    # 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
  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)
  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
  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:
  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: