#!/usr/local/bin/gawk -f
#  
#  Calculate N and C chemical shift references 
#   on base of the zero-ppm frequency of hydrogen 
#   using the following gamma-ratio's:
#  
#  
#
#    15N/1H = 0.101329118 (DSS) [1]
#    13C/1H = 0.251449530 (DSS) [1]
#  
#@  determining calibration parameters from bruker acqu file
#@  use: calib.awk temp=298 acqu
#@  [1] Wishart, D.S; Sykes, B.D. (1995) J. Biomol. NMR., 6, 135-140
#@      1H, 13C and 15N chemical shift referencing in biomolecular NMR
#

BEGIN { 
    stderr="/dev/stderr"
    calibout="HNC_calib.out"
    if(!t0) t0=298
    if(!temp) temp=t0
    progname="calib.awk"
    if ( ARGC==1 || ARGV[1]~/^(h|help)$/ || help!="" ) { 
        print "## calib.awk ## determining calibration parameters" 
        print "#  use: calib.awk temp=298 acqu" 
        print ARGV[0]
        #print "  or  calib.awk [hit return]" > stderr
        # 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
    }
}
END  {  if (help) exit }

FNR==1 && NF==0 { 
    exit
}

FILENAME ~ /(^|[/])acqu/ && /^##[$]SFO[123]=/ { 
    if ($1=="##$SFO1=") fH=$2
    else if ($1=="##$SFO2=") fC=$2
    else if ($1=="##$SFO3=") fN=$2
    else { next }  
    strf=strf sprintf("%-27s",$0)
}

FILENAME ~ /(^|[/])acqu/ && /^##[$]BF[123]=/ { 
#    file=FILENAME
    if ($1=="##$BF1=") bH=$2
    else if ($1=="##$BF2=") bC=$2
    else if ($1=="##$BF3=") bN=$2
    else { next }  
    strb=strb sprintf("%-27s",$0)
}

FILENAME ~ /(^|[/])acqu/ && /^##[$]O[123]=/ { 
    if ($1=="##$O1=") oH=$2
    else if ($1=="##$O2=") oC=$2
    else if ($1=="##$O3=") oN=$2
    else { next }  
    stro=stro sprintf("%-27s",$0)
}

END {
    if (strf || strb || stro ) { 
        print strf
        print " "strb
        print "  "stro
    }

    while (!fH || !fC || !fN ) {
        if (!fH) { 
            printf "give proton carrier freq:   "
            getline fH
        } 
        if (!fN) { 
            printf "give nitrogen carrier freq: "
            getline fN
        } 
        if (!fC) { 
            printf "give carbon carrier freq:   "
            getline fC 
        }
    }

    rHN=0.101329118
    rHC=0.251449530
    brHN=0.101329120
    brHC=0.251450201
    if (!bH) {
        printf "# getting bruker ratio from memory\n",FILENAME
    } else {
        brHN=bN/bH
        brHC=bC/bH
        printf "# getting ratio's from file:   %s\n",FILENAME
    }
    printf "# bruker   r(N/H)=%-12.9f  r(C/H)=%-12.9f\n",brHN,brHC
    printf "# DSS:     r(N/H)=%-12.9f  r(C/H)=%-12.9f\n",rHN,rHC
    print  "\n# [1] Wishart, D.S; Sykes, B.D. (1995) "
    print  "#     1H, 13C and 15N chemical shift referencing"
    print  "#     in biomolecular NMR. JBNMR., 6, 135-140"
 
    if (bruker) {
        rHN=bN/bH
        rHC=bC/bH
    }


    if(fC && fN && fH && !found) {
        if (fH<fC) { tmp=fC;fC=fH;fH=tmp }
        if (fH<fN) { tmp=fN;fN=fH;fH=tmp }
        if (fC<fN) { tmp=fN;fN=fC;fC=tmp }
    }

    if (verbose) { 
        printf "%12.5f %12.5f %12.5f %12.5f\n",fH,bH,oH,bH+oH/1000000
        printf "%12.5f %12.5f %12.5f %12.5f\n",fN,bN,oN,bN+oN/1000000
        printf "%12.5f %12.5f %12.5f %12.5f\n",fC,bC,oC,bC+oC/1000000
    }

    printf "# T = %-15s  %43s\n",temp" K","\\"
    calcppms(temp,brHN,brHC)
    obsform="  -xOBS %14.4f   -yOBS %14.4f   -zOBS %14.4f   \\\n"
    carform="  -xCAR %14.4f   -yCAR %14.4f   -zCAR %14.4f   \\\n"
    
    if(refH) {
        printf "# bruker convention  %45s\n","\\"
        printf obsform,fH,fN,fC
        printf carform,refH,refN,refC
    }
    refH=""

    calcppms(temp,rHN,rHC)
    if(refH) { 
        printf "# DSS: (Wishart/Sykes;JBNMR 6,135-140) %27s\n","\\"
        printf obsform,fH,fN,fC
        printf carform,refH,refN,refC
    }
    print ""
}

function calcppms(temp,rHN,rHC) {
    refH = 7.83 - (temp/96.9)
    if (nowat && bruker) refH=( ( fH-bH)/ (fH/1000000) )
    refN = (fN * 1000000 - ( 1000000 * fH - fH * refH ) * rHN ) / fN
    refC = (fC * 1000000 - ( 1000000  -  refH ) * fH * rHC ) / fC
}