function [gamma, gamma_deriv] = build_gamma(ti,varargin)
% function [gamma, gamma_deriv] = build_gamma(p, x, Th, n)
%
% Create the gamma function to approach the the amount of water that would
% produce the same beam-hardening effect as the given amount of bone traversed

%
% in
%  ti    Bone projection
% 
% options
%  'type'   'exp' exponential, 'pol' JS approximation or 'aprox_pol'
%  approximated polynomial for iterative BH reconstruction method
%  'A','B'  Parameters to aproximate the function
%
%
% out
% gamma     Gamma function
% Gamma_deriv   Derivative gamma function
%
% Written to be used with sir_bh_os.m
%
% Jeff Fessler, University of Michigan
% Monica Abella and Cristobal Martinez, Universidad Carlos III de Madrid
arg.type='aprox_pol';
arg.alpha=0.3;
arg.beta1=4;
arg.beta2=2;
arg.A=4;
arg.B=2;

arg = vararg_pair(arg, varargin);

switch arg.type
case 'exp'
    alpha=arg.alpha;
    beta1=arg.beta1;
    beta2=arg.beta2;
    aux_gamma=alpha*exp(-beta1*ti)+(1-alpha)*exp(-beta2*ti);
    gamma = -log(aux_gamma);
    gamma_deriv = (alpha*beta1*exp(-beta1*ti)+(1-alpha)*beta2*exp(-beta2*ti))./aux_gamma;

case 'pol'
    A=arg.A;
    B=arg.B;
    gamma=A*ti-B*ti.^2;
    gamma_deriv=A- 2*B*ti;

case 'aprox_pol'
    A=arg.A;
    B=arg.B;
    gamma=A*ti./(1+(B/A)*ti);
    gamma_deriv=A./(1+(B/A)*ti).^2;
otherwise
    fail('unknown type %s', arg.type);
end