function [f,dfdx,dfdy] = biMatCrackBranch(r,theta,e,alpha)
% Compute the branch functions spanning the near tip field for bimaterial
% cracks.
% Inputs:
%   (r,theta) : polar coordinates of points where the branch
%               functions are to be evaluated
%   alpha     : inclination of the crack tip segment w.r.t x axis
% Vinh Phu Nguyen, nvinhphu@gmail.com
% PhD at Delft University of Technology, The Netherlands
% April, 2012, Hue, Vietnam
% Adapted from MXFEM of Mathiew Pais, Florida Univ.

sr       = sqrt(r);
c        = 1/2/sr; 
st       = sin(theta);
ct       = cos(theta);
st2      = sin(theta/2);
ct2      = cos(theta/2);
s3t2     = sin(3*theta/2);
c3t2     = cos(3*theta/2);
lgr      = log(r);
coselogr = cos(e*lgr);
sinelogr = sin(e*lgr);
expeth   = exp(e*theta);
expmeth  = 1/expeth;%exp(-e*theta);
dd1      = e*lgr-theta;
dd2      = e*lgr+theta;
selogrmt = sin(dd1);
selogrpt = sin(dd2);
celogrmt = cos(dd1);
celogrpt = cos(dd2);

% 12 branch functions (checked)

a1 = sr*coselogr*expmeth;
a2 = sr*coselogr*expeth;
a4 = sr*sinelogr*expeth;
a5 = sr*sinelogr*expmeth;

f  = [a1*st2;          % Alpha 1 crack tip enrichment value
      a1*ct2;          % Alpha 2 crack tip enrichment value
      a2*st2;          % Alpha 3 crack tip enrichment value
      a2*ct2;          % Alpha 4 crack tip enrichment value
      a2*st2*st;       % Alpha 5 crack tip enrichment value
      a2*ct2*st;       % Alpha 6 crack tip enrichment value
      a5*st2;          % Alpha 7 crack tip enrichment value
      a5*ct2;          % Alpha 8 crack tip enrichment value
      a4*st2;          % Alpha 9 crack tip enrichment value
      a4*ct2;          % Alpha 10 crack tip enrichment value
      a4*st2*st;       % Alpha 11 crack tip enrichment value
      a4*ct2*st];

% Derivative of bimaterial crack tip enrichment functions 
% with respect to x1 (crack tip coordinate system) (checked)

px = c*[-expmeth*st2*(coselogr + 2*e*selogrmt);...
         expmeth*ct2*(coselogr - 2*e*selogrmt);...
        -expeth*st2*(coselogr + 2*e*selogrpt);...
         expeth*ct2*(coselogr - 2*e*selogrpt);...
        -expeth*st*(coselogr*s3t2 + 2*e*selogrpt*st2);...
        -expeth*st*(coselogr*c3t2 + 2*e*selogrpt*ct2);...
         expmeth*st2*(-sinelogr + 2*e*celogrmt);...
         expmeth*ct2*(sinelogr + 2*e*celogrmt);...
         expeth*st2*(-sinelogr + 2*e*celogrpt);...
         expeth*ct2*(sinelogr  + 2*e*celogrpt);...
         expeth*st*(-sinelogr*s3t2 + 2*e*celogrpt*st2);...
         expeth*st*(-sinelogr*c3t2 + 2*e*celogrpt*ct2)];

% Derivative of bimaterial crack tip enrichment functions 
% with respect to x2 (crack tip coordinate system)

py = c*[expmeth*(coselogr*ct2-2*e*celogrmt*st2);...
        expmeth*(coselogr*st2-2*e*celogrmt*ct2);...
        expeth*(coselogr*ct2+2*e*celogrpt*st2);...
        expeth*(coselogr*st2+2*e*celogrpt*ct2);...
        expeth*(coselogr*(st2+s3t2*ct)+2*e*celogrpt*st2*st);...
        expeth*(coselogr*(ct2+c3t2*ct)+2*e*celogrpt*ct2*st);...
        expmeth*(sinelogr*ct2-2*e*selogrmt*st2);...
        expmeth*(sinelogr*st2-2*e*selogrmt*ct2);...
        expeth*(sinelogr*ct2+2*e*selogrpt*st2);...
        expeth*(sinelogr*st2+2*e*selogrpt*ct2);...
        expeth*(sinelogr*(st2+s3t2*ct)+2*e*selogrpt*st2*st);...
        expeth*(sinelogr*(ct2+c3t2*ct)+2*e*selogrpt*ct2*st)];

CCS   = [cos(alpha) sin(alpha);-sin(alpha) cos(alpha)];
ct    = CCS(1,1); 
st    = CCS(1,2);

% Derivative of bimaterial crack tip enrichment functions with respect to X

dfdx = [px(1)*ct+py(1)*-st;...
        px(2)*ct+py(2)*-st;...
        px(3)*ct+py(3)*-st;...
        px(4)*ct+py(4)*-st;...
        px(5)*ct+py(5)*-st;...
        px(6)*ct+py(6)*-st;...
        px(7)*ct+py(7)*-st;...
        px(8)*ct+py(8)*-st;...
        px(9)*ct+py(9)*-st;...
        px(10)*ct+py(10)*-st;...
        px(11)*ct+py(11)*-st;...
        px(12)*ct+py(12)*-st];

% Derivative of bimaterial crack tip enrichment functions with respect to Y

dfdy = [px(1)*st+py(1)*ct;...
        px(2)*st+py(2)*ct;...
        px(3)*st+py(3)*ct;...
        px(4)*st+py(4)*ct;...
        px(5)*st+py(5)*ct;...
        px(6)*st+py(6)*ct;...
        px(7)*st+py(7)*ct;...
        px(8)*st+py(8)*ct;...  
        px(9)*st+py(9)*ct;...
        px(10)*st+py(10)*ct;...
        px(11)*st+py(11)*ct;...
        px(12)*st+py(12)*ct];
