function f=dancoff(a,b)
% Dancoff factor calculation in a square lattice of cylindrical pins
% Reference: J. Bussac and P. Reuss, Traite de neutronique, Hermann (1978).
% f=dancoff(a,b)
% a: rod radius / mean free path
% b: lattice pitch / rod radius
  amesh=0:0.2:2.0 ;
  bmesh=[2, 2.5, 3, 4, 5] ;
  if b == 2
    dan=[1., 0.903, 0.825, 0.761, 0.708, 0.664, 0.626, 0.594, 0.566, 0.542, 0.521] ;
    f=interp1(amesh, dan, a, 'cubic') ;
  elseif b == 2.5
    dan=[1., 0.710, 0.533, 0.413, 0.328, 0.264, 0.216, 0.178, 0.148, 0.124, 0.105] ;
    f=interp1(amesh, dan, a, 'cubic') ;
  elseif b == 3
    dan=[1., 0.552, 0.346, 0.232, 0.161, 0.115, 0.083, 0.062, 0.046, 0.034, 0.026] ;
    f=interp1(amesh, dan, a, 'cubic') ;
  elseif b == 4
    dan=[1., 0.328, 0.152, 0.079, 0.044, 0.025, 0.015, 0.009, 0.005, 0.004, 0.003] ;
    f=interp1(amesh, dan, a, 'cubic') ;
  elseif b == 5
    dan=[1., 0.197, 0.071, 0.030, 0.014, 0.006, 0.003, 0.002, 0.002, 0.001, 0.000] ;
    f=interp1(amesh, dan, a, 'cubic') ;
  else
    error('invalid value of b')
  end