#!/usr/local/bin/gawk -f
# 
#  
#   Use:
#
#     most parameters can be read from a .par file [belonging to .nv file], 
#     only d19 or dnul or dnulppm or dnulhz has to be given
#
#       unwatergate.awk -v d19=85us         spectrum.par spectrum.xpk > new.xpk
#       unwatergate.awk -v dnul=7.84174     spectrum.par spectrum.xpk > new.xpk
#       unwatergate.awk -v dnul=5882.35hz   spectrum.par spectrum.xpk > new.xpk
#       unwatergate.awk -v dnulhz=5882.35   spectrum.par spectrum.xpk > new.xpk  
#     #
#     unwatergate.awk -v dnul=7.84174 -v ref=4.774 spectrum.xpk
#
#     unwatergate.awk -v ref=4.731ppm -v dnulppm=4.861 spectrum.par spectrum.xpk
#     unwatergate.awk -v ref=4.731ppm -v dnul=4.861 spectrum.par spectrum.xpk
#
#
#   Purpose:
#
#@    undo the effect of watergate water suppression on peak intensities
#@    in nmrview xpk files
#
#   Command line options:
#   
#     -v dom = 1  [default=1]
#                       spectral domain on which to apply the correction 
#
#     -v d19 = val[s|ms|us|mus]      
#                       d19 bruker watergate delay in (milli/micro)seconds 
#                       can be e.g: ( all equivalent) 
#                           0.000085 0.000085s 0.085ms 85us 85mus  
#                       
#                       the nearest zero crossing will be 1/(d19*2) Hz
#
#      the parameters below are normally retrieved from a nmrview *.par file   
#
#     -v ref = val      water frequency in ppm
#
#     -v dnul = val[ppm|hz]  
#                       distance from water to nearest zero in 
#                       ppm[default] or Hz
#     -v dnulppm = val  distance from water to nearest zero in ppm
#     -v dnulhz = val   distance from water to nearest zero in ppm
#                       this should be ginven in combination with:
#
#
#     -v sf = val       spectrom frequency in MHz 
#     -v sw = val       sweep width in Hz
#
#
#     -v estimate=val   estimate nearest zero to be at val*sweep width 
#                       so with estimate=0.5 the nearest zero is at the 
#                       edge of the spectrum     
#
#
#
#     -v xpk=1          next (standard) input will be considered to be a xpk file 
#                       
#                       e.g. 
#                          cat f.xpk | unwatergate.awk -v d19=85us spec.par xpk=1 - > out.xpk
#                          cat f.xpk | unwatergate.awk -v d19=85us -v xpk=1 spec.par - > out.xpk
#
#     -v pedantic=1      be pedantic
#
#     -v mode=inv        undo the correction 
#

BEGIN {

    stderr="/dev/stderr"
    progname="unwatergate.awk"
    if ( ARGV[1]~/^(h|help)$/ || help!="" ) { 
        print ARGV[0]
        # because of the which command, help is only available if 
        # this script is in the PATH
        system("gawk '/^#/{print}/^BEGIN/{exit}' `which "progname"`")
        help=1
        exit
    }
    # print "ARGC",ARGC>stderr 
    # for (i in ARGV) print "ARGV",i,ARGV[i]>stderr 
    
    IGNORECASE=1
    setstr="    |"
    if (dom=="") { 
        dom=1
    }
    print "##  unwatergate.awk">stderr
    print "    | undoing the effect of watergate watersuppression">stderr
    print "    | using ppm-values in domain",dom,"for correction">stderr
    
    if (pi=="") pi = 3.141592653589793
    
    if ( ARGC==1 ) { 
        help++
        com  = "gawk '/^[^#]/{exit}1' `which unwatergate.awk`" 
        while ( (com | getline)>0) { 
            print > stderr
        } 
        exit 
    }
    if ( ARGC==1 || ( 1 in ARGV && ARGV[1]=="help") ) { 
        help++
        system("gawk '/^[^#]/{exit}1' `which unwatergate.awk`")
        exit 
    }

    
    
    if (dnul!="") { 
        if (dnul~/hz$/) dnulhz=dnul+0   
        else if (dnul~/ppm$/) dnulppm=dnul+0   
        else dnulppm = dnul+0
    }
    
    if (ref!="")     print setstr,"ref       :",( ref ),"ppm">stderr
    if (dnulhz!="")  print setstr,"dnulhz    :",( dnulhz ),"Hz">stderr
    if (dnulppm!="") print setstr,"dnul      :",( dnulppm ),"ppm">stderr
    if (sf!="")      print setstr,"sf        :",( sf ),"MHz">stderr
    if (sw!="")      print setstr,"sw        :",( sw ),"Hz">stderr
    
    
}

FILENAME ~ /[.]par$/ {
    parfileread++
    # only set them if they are not set before
    if ($1=="ref"   && $2 == dom && ref   == "") 
        print setstr,"ref[par]  :",( ref   = $3 ),"ppm">stderr
    if ($1=="label" && $2 == dom && label == "") 
        print setstr,"label[par]:",( label = $3 )>stderr
    if ($1=="sf"    && $2 == dom && sf    == "") 
        print setstr,"sf[par]   :",( sf    = $3 ),"MHz">stderr 
    if ($1=="sw"    && $2 == dom && sw    == "") 
        print setstr,"sw[par]   :",( sw    = $3 ),"Hz">stderr
    next 
}

!referenced {
    
    if (d19!="" && !dnulhzset++)  { 
        if (d19~/ms$/) d19=(d19+0)/1000
        else if (d19~/m?us$/) d19=(d19+0)/1000000
        else if (d19~/s$/) d19=(d19+0)/1
        if (d19!="") {
            print setstr,"d19       :",(d19*1000),"ms">stderr
            print setstr,"dnulhz    :","1/(2*"d19")=",(dnulhz=1/d19/2),"Hz">stderr # guestimate
        }
    }

    if (sf!="" && sw!="" && swppm=="")  
        print setstr,"swppm     :",sw"/"sf"=",(   swppm = sw/sf ),"Hz">stderr

    if (sf!="" && dnulhz!="" && dnulppm=="" ) 
        print setstr,"dnulppm   :",dnulhz"/"sf"=",(dnulppm = dnulhz/sf),"ppm">stderr 

    if (estimate!="")
        if (dnulppm=="" && swppm!="") 
            print setstr,"dnulppm[guess]  :",(dnulppm = swppm/(estimate*2)),"ppm">stderr # guestimate
    

    if (sf!="" && dnulhz=="" && dnulppm!="" ) 
        print "|"setstr,"dnulhz   :",dnulppm"*"sf"=",(dnulhz = dnulppm*sf),"ppm">stderr 

    if (ref!="" && dnulppm!="") { #print "referenced!">stderr
        referenced++
    }
}

FILENAME ~ /[.]xpk$/ && xpk==""  { xpk=1 }

xpk && FNR==2{ #label line
    if (label=="") print setstr,"label     :",(label=$dom)>stderr
    else if ($dom!=label) { 
        print "label in .par and .xpk files don't match ",label,"<>",$dom
    } 
}
#xpk { print FNR > stderr }

xpk && (FNR==4 || FNR==5) { #print "xpk5">stderr
    s=$0;gsub(/[{}]/," ",s)
    if (split(s,xxx)<dom ) {
        print "warning not enough dimensions in "FILENAME,"line:"FNR,$0
    }
    if (FNR==4) {
        if (sw=="" && !referenced) { 
        #if (sw=="" ) { 
            print setstr,"sw[xpk]   :",( sw=xxx[dom] ) >stderr
        } else { 
            if (sw!=xxx[dom] && pedantic) 
                print "warning: sw<>sw[xpk]   :",sw,"<>",xxx[dom]  >stderr
            #else print "       : sw==sw[xpk]   :",sw,"==",xxx[dom]  >stderr
        }
    }
    if (FNR==5) {
        if (sf=="" && !referenced) {
        #if (sf=="" ) { 
            print setstr,"sf[xpk]   :",( sf=xxx[dom] ) >stderr
        } else { 
            if (sf!=xxx[dom] && pedantic) 
                print "warning: sf<>sf[xpk]   :",sf,"<>",xxx[dom]  >stderr
            #else print "       : sf==sf[xpk]   :",sf,"==",xxx[dom]  >stderr
        }
    }
}

xpk && FNR==6{

    if (final) {
        print "      final values:" > stderr
        print "      ","label="label,"ref="ref,"sf="sf,"sw="sw,"swppm"swppm,"[/2]:"dnulppm > stderr
        print "      ","d19="d19*1000"ms ; dv[Hz]="dnulhz" ;dv[ppm]="dnulppm > stderr
    }
    if (dnulppm=="") { 
        print "dnulppm [dppm(H20-zero)] value not available" > stderr
        exit
    }
    if (ref=="") { 
        print "ref ppm value not available" > stderr
        exit
    }

    for (i=1;i<=NF;i++) if ($i=="int") intcol=i+1
    for (i=1;i<=NF;i++) if ($i=="vol") volcol=i+1
    for (i=1;i<=NF;i++) if ($i==label".P") ppmcol=i+1
    #print > stderr
}

xpk && FNR>6{
    # data lines
    gsub(/[{][ \t]+/,"{")
    gsub(/[ \t]+[}]/,"}")
    sepchar = "|"
    
    $0=toupper($0)
    while (match($0,/[{][^{}]*[ \t]+[^{}]*[}]/)) {
        lstr=substr($0,1,RSTART-1)
        mstr=substr($0,RSTART,RLENGTH)
        rstr=substr($0,RSTART+RLENGTH)
        if (verbose) print mstr
        gsub(/[ \t]+/,"|",mstr)
        $0=lstr mstr rstr
        if (verbose) print mstr
    }
    ppm=$ppmcol
    
    dppm = ppm - ref  # 

    adjfac = abs(sin(pi * (dppm)/dnulppm ) ) 
    
    # print ppm,adjfac | "sort -n|xmgr -pipe"

    # print intcol,volcol,$intcol,$volcol >stderr
    if (mode=="inv") {
        if ($volcol>0) newvol = sprintf("%.4f0",$volcol*adjfac)
        if ($intcol>0) newint = sprintf("%.4f0",$intcol*adjfac)
        sub($volcol" +"$intcol,newvol" "newint)
    } else { 
        if ($volcol>0) newvol = $volcol/adjfac
        if ($intcol>0) newint = $intcol/adjfac
        sub($volcol" +"$intcol,newvol" "newint)
    }    
    gsub("["sepchar"]"," ")


}

1

function abs(x) { return (x<0?-x:x) }