# ==== setting output and test-flag ====
sp          := Caro      # select to Cali, Caro or N
output      := 750_{$sp}
#output      := gnoesyCaliCaroNhsqc_{$sp}
test_ft     := n         # test-ft for phasing with XEASY
reduce_flag := n         # reduce size for test-ft by factor of 2
plane := n               # set plane := 12 for HH plane, or 13 for HN/HC
                         # set plane := n for 3D spectrum
lp := n
np := 512
ni := 140
ni2:= 48

# ==== reading of data ====
readvnmr3d fid $np $ni $ni2 byteswap            # byteswap is needed on Linux

if ('$plane' == '12') then
   dim 3
   reduce 1..1                          # select plane: 1..1, 2..2, etc.
   dim 1 2 3
   write integer tmp
   read integer tmp {$np}c {$ni}c
end if
if ('$plane' == '13') then
   dim 2
   reduce 1..1                          # select plane: 1..1, 2..2, etc.
   dim 1 3 2
   write integer tmp
   read integer tmp {$np}c {$ni2}c
end if


# ==== setting of parameters ====
#DIMENSION: 1 (1H)
si(1)      = 1024         # frequency domain size
aqm(1)    := RSHc             # acquisition mode (ST/TP, TPPI, RSH)
w0(1)      = 750.5138600      # frequency
delta(1)   = 0.000121200      # increment in time domain
ppmmax(1)  = 10.277        # maximum ppm for calibration

#DIMENSION: 2 (1H)
si(2)      = 256         # frequency domain size
aqm(2)    := RSHc             # acquisition mode (ST/TP, TPPI, RSH)
w0(2)      = 750.5138600      # frequency
delta(2)   = 0.000121200      # increment in time domain
ppmmax(2)  = 10.277        # maximum ppm for calibration
fpmult(2) := 0.5      # first point multiplier

#DIMENSION: 3 (13C ali, 13C aro, N)
si(3)      = 128         # frequency domain size
aqm(3)    := RSHc             # acquisition mode (ST/TP, TPPI, RSH)
fpmult(3) := 1.0      # first point multiplier
if ('$sp' == 'Cali') then
   w0(3)      = 188.7163574      # frequency
   delta(3)   = 0.000120413      # increment in time domain
   ppmmax(3)  = 92.940        # maximum ppm for calibration
else if ('$sp' == 'Caro') then
   w0(3)      = 188.7163574      # frequency
   delta(3)   = 0.000120413      # increment in time domain
   ppmmax(3)  = 224.960        # maximum ppm for calibration
else
   w0(3)      = 76.0489075      # frequency
   delta(3)   = 0.000454545      # increment in time domain
   ppmmax(3)  = 132.674        # maximum ppm for calibration
end if


if ('$plane' == '13') then
   si(2)      = $si(3)
   aqm(2)    := $aqm(3)
   w0(2)      = $w0(3)
   delta(2)   = $delta(3)
   ppmmax(2)  = $ppmmax(3)
end if

if ('$plane' == 'n') then
   set n_dim= 3
else
   set n_dim= 2
end if


window1 := "window cos2"
window2 := "window cos2"
if ('$plane' == 'n') then
   window3 := "window cos2"
end if

ph0_d(1) = 91.0+22.7
ph1_d(1) = -24.0
ph0_d(2) = 0.0
ph1_d(2) = 0.0
ph0_d(3) = 90
ph1_d(3) = -180
                                                                                                                                
if ('$plane' == '13') then
   ph0_d(2) = $ph0_d(3)
   ph1_d(2) = $ph1_d(3)
   fpmult(2) = fpmult(3)
end if


# ====================
# ==== processing ====
# ====================

dim 1
if ('$sp' == 'Cali') then
   suppress cos 128     # use mild water suppression for aliphatic
else
   suppress cos 32      # use aggressive water suppression
end if

# ==== ft ====
if ('$sp' == 'Cali') then
   strip_ft 1 1..$si(1)
else
   strip_ft 1 1..$si(1)/2        # FT 1H dimension and extract left half
end if

if ('$lp'!='y') then
   dim 2
   mul fpmult(2) 1
   strip_ft 2 1..$si(2)
   if ('$plane' == 'n') then
      dim 3
      mul fpmult(3) 1
      strip_ft 3 1..$si(3)
   end if
else if ('$plane' == 'n') then   # with linear prediction
   dim 2
   ft $si(2)
   dim 3
   predict lpsvd 16 24
   mul fpmult(3) 1
   strip_ft 3 1..$si(3)
   dim 2
   ift $si(2) 1 $ni
   predict lpsvd 16 70
   mul fpmult(2) 1
   strip_ft 2 1..$si(2)
else
   print "Linear prediction cannot be used for 2D planes!"
   print "Use either plane = 'n' lp = 'y' or plane = '12' lp = 'n'"
   exit
end if

# ==== baseline correction ====
if ('$test_ft'!='y') then
  dim 1
cflatt iterative rft 3 $si(1) 1 $ph1_d(1)
#  cflatt iterative polynom 5
  dim 2
cflatt iterative rft 3 $si(2) 1 $ph1_d(2)
#  cflatt iterative polynom 5
  if ('$plane' == 'n') then
     dim 3
cflatt iterative rft 3 $si(3) 1 $ph1_d(3)
#     cflatt iterative polynom 5
  end if
end if

if ('$plane' == 'n') then
   dim 1 2 3
else
   dim 1 2
end if

scale noise

# ==== write output ====
if ('$test_ft'.eq.'y') then
   if ('$plane' == 'n') then
      write easy16 $output.rrr r r r byteswap
      write easy16 $output.rri r r i byteswap         # 15N/13C imaginary part
      write easy16 $output.irr i r r byteswap         # 1H noe  imaginary part
      write easy16 $output.rir r i r byteswap         # 1H inept imaginary part
   else
      write easy16 {$output}.{$plane}.rr r r byteswap
      write easy16 {$output}.{$plane}.ir i r byteswap
      write easy16 {$output}.{$plane}.ri r i byteswap
   end if
else
   if ('$plane' == 'n') then
      write easy16 $output byteswap
      if ('$sp' == 'N') then
        system "echo 'Identifier for dimension w1 ... N15' >> $output.3D.param"
      else
        system "echo 'Identifier for dimension w1 ... C13' >> $output.3D.param"
      end if
      system "echo 'Identifier for dimension w2 ... H1'  >> $output.3D.param"
      system "echo 'Identifier for dimension w3 ... H1'  >> $output.3D.param"
   else
      write easy16 {$output}.{$plane} byteswap
   end if
end if