#!/bin/sh
#
# calculates chelator binding constant using  monte-carlo type fitting procedure
# finds optimal fit by randomly choosing values of the fitted parameters within an
# n-dimensional cube. the cube shrinks until the error sum improves less than 
# mindiff

if test $# = 0
        then
                echo "usage:   kjell_ator inputfile kd_est tries denominator"
                exit
fi

# the inputfile can look as follows:
# cchel 28.8	chelator concentration in uM
# c_0 1		initial ion concentration in uM 
# v_0 2500	volume in uL
# 0 5		added vol in uL and concentration in mM (for each of the following rows)
# 0.5379	absorbance
# 5 2.35        added vol in uL and concentration in mM (for each of the following rows)
# 0.4514        absorbance
# 0.3678        absorbance
# 0.2915        absorbance
#   .
#   .
#   .
# kd_est is a starting value for kd in M
# tries is the number of random values used 
# denominator is the shrinkage factor of the cube
# anders malmendal 1997-2001

nawk 'BEGIN{
  k1_0=1/'$2'
  fact_0=1.0		#concentration correction factor
  tries='$3'
  denom='$4'
  mindiff=0.0001
  k1_range=0.1
  fact_range=0.1
  aapo_range=0.1        #range for allowed values of aapo
  aholo_range=0.1       #range for allowed values of aholo
  pi=3.1415926359
  srand()
  minsum=100000
  print "round\tkd\t\tcorrection\taapo\t\taholo\t\terrsum"
}

$1~/cchel/ {cchel_0=$2/1e6}
$1~/v_0/ {v_0=$2/1e6}
$1~/c_0/ {c_0=$2/1e6}

$1~/^[0-9.]/ && NF==1 {
a[++n]=$1
ionvol[n]=ionvol_
ioncons[n]=ioncons_
}

$1~/^[0-9.]/ && NF==2 {
  ionvol_=$1/1e6
  ioncons_=$2/1e3
}

END {
  nmax=n
  ctot=c_0
  vtot_=v_0
  nion=c_0*v_0
  for (n=1;n<=nmax;n++) {
    vtot_=vtot_+ionvol[n]
    vtot[n]=vtot_
    if (v_0!=0) korr[n] = vtot[n]/v_0
    if (korr[n]!=0)cchel[n]=cchel_0/korr[n]
    nion=nion+ioncons[n]*ionvol[n]
    iontot[n]=nion/vtot_
  }
  aapo_0=a[1]
  aholo_0=a[nmax]*vtot[nmax]/vtot[1]
  while (sqrt((1-oldsum/minsum)^2)>mindiff || oldsum==0) {
    oldsum=minsum
    if (++round!=1) {
      k1_best=k1_0
      fact_best=fact_0
      aapo_best=aapo_0
      aholo_best=aholo_0 
    }
    for (i=1;i<=tries;i++) {
      sum=0
      k1=k1_0^(1-k1_range/2+k1_range*rand())
      kdchel=1/k1
      fact=fact_0*(1-fact_range/2+fact_range*rand())
      aapo=aapo_0*(1-aapo_range/2+aapo_range*rand())
      aholo=aholo_0*(1-aapo_0/aholo_0*aholo_range/2+aapo_0/aholo_0*aholo_range*rand())
      y = 1e-12
      for (n=1;n<=nmax;n++) {
        f1=(cchel[n]*fact+iontot[n]+kdchel)/2
        f2=cchel[n]*fact*iontot[n] 
        cchelion=f1-sqrt(f1^2-f2)
        acalc=(aapo-(aapo-aholo)*cchelion/(cchel[n]*fact))/korr[n]
        sum=sum+(a[n]-acalc)^2
      }
      if (sum<minsum) {
        minsum=sum
        k1_0=sqrt(k1*k1_0)
        fact_0=(fact+fact_0)/2
        aapo_0=(aapo+aapo_0)/2
        aholo_0=(aholo+aholo_0)/2
        k1_best=k1
        fact_best=fact
        aapo_best=aapo
        aholo_best=aholo
#      print i,1/k1,fact,aapo,aholo,sum
      }
    }
    if ( round%1==0 ) {
      print round "\t" 1/k1_best "\t" fact_best "\t" aapo_best "\t" aholo_best "\t" minsum
    }
    k1_0=k1_best
    fact_0=fact_best
    aapo_0=aapo_best
    aholo_0=aholo_best
    k1_range=k1_range/denom
    fact_range=fact_range/denom
    aapo_range=aapo_range/denom
    aholo_range=aholo_range/denom
  }
}' $1 > $5