function avp=avPNumS(mu,s,T,a,b)
mu1=mu(a);
s1=s(a);
T1=T(a);
mu2=mu(b);
s2=s(b);
T2=T(b);
%2*(mu2-mu1)/(sqrt(s1)+sqrt(s2))
db=1/(0.001987191*T2)-1/(0.001987191*T1);
strt=mu1-3*sqrt(s1);
endd=mu2+3*sqrt(s2);
avp=0;
stp=(endd-strt)/500;
for u2=strt:stp:endd
    numr=1/4*2^(1/2)*pi^(1/2)*(1+erf(1/2*2^(1/2)*(-u2+mu1)/sqrt(s1)))/pi/sqrt(s2)*exp(-1/2*(mu2-u2)^2/s2);
    % check not out of range of arith
    tmp = -1/2*(-2*db*u2*s2+u2^2-2*u2*mu2+mu2^2+2*mu1*db*s2-s1*db^2*s2)/s2;
    if tmp < 200 
        numr=numr-1/4*2^(1/2)*pi^(1/2)*(-1+erf(1/2*2^(1/2)*(-u2+mu1-db*s1)/s1^(1/2)))*exp(tmp)/pi/s2^(1/2);
    end
    %
    if u2 ~= strt
        avp = avp + (numr+lastn)/2*stp;
    end
    lastn = numr;
end
