function out = ddBS( L0,p0,t0,e0,p,t,e )
  Nk = N( L0,p0,t0,e0 );
  Np = N( L0,p ,t ,e );
  out = ( Nk-Np ) / ( 1+Np*1e-6 );

% =========================================================================
function out = N( L, p, t, e )
  N0 = 287.604 + 3*1.6288/L^2 + 5*0.0136/L^4;
  out = N0 * p/760 * 1/(1+t/273.16) - 5.5e-2 *e /(1+t/273.16);


 