# Copyright (c) 2002 Peter Guntert. All rights reserved.
## 7MACROS: cflatt - baseline correction with standard base function sets
##
## Usage: cflatt flatt      <n> <tau>    <set> <m> [<ntotal>] [<nfirst>] [<phase>]
##        cflatt derivative <n> <tau>    <set> <m> [<ntotal>] [<nfirst>] [<phase>]
##        cflatt iterative [{<range>} -] <set> <m> [<ntotal>] [<nfirst>] [<phase>]
##        cflatt manual    [{<range>} -] <set> <m> [<ntotal>] [<nfirst>] [<phase>]
##
## "cflatt" ("convenient flatt") is a macro for baseline correction
## that provides a convenient interface to the "flatten" command.
## The parameters <range>, <n>, and <tau> have the same meaning as
## for the "flatten" command.
## 
## set       base function set
## --------- ----------------------------------------------------------
## cft/rft   trigonometric funtions that correspond to the first <m>
##           data points before complex/real Fourier transformation
## cftw/rftw same as "cft" or "rft" plus Lorentzian functions to account  
##           for contributions from the water line which is assumed to be 
##           located in the middle of the spectrum (before a possible
##           strip transform)
## polynom   polynomial of order <m>
##
## For the trigonometric methods it is assumed that the present data are a
## strip starting at position <nfirst> out of <ntotal> original data points 
## after Fourier transformation. The default is to assume that no strip 
## transform was performed, i.e. ntotal = n and nfirst = 1. <phase> 
## denotes the value (in degrees) of the first order phase correction,
## the default is no linear phase correction, i.e. phase = 0.0. The
## parameters <ntotal>, <nfirst> and<phase> must not be given if polynomial 
## base functions are used.
##
## See also: flatten

var par baseset norder ntotal nfirst phase fun i j t \
    t0 t1 t2 t3 t4 t5 t6 t7 t8 t9 t10 r c

syntax ** replace

do i 2 nparam
  if ('$p$i'.eq.'-') break
end do
if (i.gt.nparam) then
  if ('$p1'.eq.'flatt' .or. '$p1'.eq.'derivative') then
    i = min(3,nparam)
  else
    i = min(1,nparam)
  end if
end if
par := 
do j 1 i
  par := $par $p$j
end do

baseset := '$p$j'
norder = 3
ntotal = n
nfirst = 1
phase  = 0.0
do i j+1 nparam
  if (i.eq.j+1) norder = p$i
  if (i.eq.j+2) ntotal = p$i
  if (i.eq.j+3) nfirst = p$i
  if (i.eq.j+4) phase  = p$i
  if (i.gt.j+4 .or. (baseset.eq.'polynom' .and. i.gt.j+1)) \ 
    error "*** Error: Too many parameters."
end do

if (baseset.eq.'cft'  .or. baseset.eq.'rft' .or. \
    baseset.eq.'cftw' .or. baseset.eq.'rftw') then
  fun :=
  do i 0 norder-1
    t = i*pi/ntotal
    if (baseset(1:3).eq.'cft') t = 2*t
    t$i = t-phase*pi/180.0/ntotal
    if (abs(t$i).gt.1.0E-6) then 
      fun := $fun sin(t$i*k) cos(t$i*k)
    else
      fun := $fun 1.0
    end if
  end do
  if (baseset.eq.'cftw' .or. baseset.eq.'rftw' ) then
    j = ntotal/2-nfirst+2
    t$i = -phase*pi/180.0/ntotal
    r = -1.0/ntotal
    c = 2.0*pi/ntotal
#    fun := $fun real(exp(cmplx(0.0,t$i*k))/cmplx(1,$j-k))
#    fun := $fun aimag(exp(cmplx(0.0,t$i*k))/cmplx(1,$j-k)) 
    fun := $fun real(exp(cmplx(0.0,t$i*k))/(1-exp(cmplx(r,c*($j-k)))))
    fun := $fun aimag(exp(cmplx(0.0,t$i*k))/(1-exp(cmplx(r,c*($j-k)))))
  end if
else if (baseset.eq.'polynom') then
  fun := 1.0
  do i 1 norder
    fun := $fun (real(k)/n)**$i
  end do
else
  error "*** Error: Illegal baseset \"$baseset\"."
end if

flatten $par $if(replace,'replace',' ') $fun