function b=loads2(ni,x,p,d,Dd,c,r,f,g,bctype);
disp('loads')
h=x(2:ni+1)-x(1:ni);
%
%setting constant in streamline diffusion test quantity: v+chv'
%
csd=zeros(size(c));
if p==2; hh=(1:2:2*ni-1)=h; hh(2:2:2*ni)=h; end
if max(hh-d./(abs(c)+eps*ones(size(c))))>0; 
csd=c./(abs(c)+eps*ones(size(c)))/2; 
end
%
if p==1;
%
b(1)=f(1)*(h(1)/2+csd(1));
for j=2:ni-1;
b(j-1)=b(j-1)+f(j)*(h(j)/2-csd(j));
b(j)=f(j)*(h(j)/2+csd(j));
end
b(ni-1)=b(ni-1)+f(ni)*(h(ni)/2-csd(ni));
%
if bctype(1)==0
b(1)=b(1)-(-d(1)/h(1)-c(1)/2+r(1)*h(1)/6)*g(1);
b(1)=b(1)-csd(1)*(Dd(1)-c(1)+r(1)*h(1)/2)*g(1);
b(2)=b(2)-
b(2)=b(2)-
else
b=[g(1),b];
end
if bctype(2)==0
ndf=2*ni-1+bctype(1);
b(ndf-1)=b(ndf-1)-
b(ndf-1)=b(ndf-1)-
b(ndf)=b(ndf)-(-d(ni)/h(ni)+c(ni)/2+r(ni)*h(ni)/6)*g(2);
b(ndf)=b(ndf)-csd(ni)*(Dd(ni)-c(ni)-r(ni)*h(ni)/2)*g(2);
else
b=[b,g(2)];
end
%
elseif p==2;
%
gauss1=1/2-1/sqrt(3)/2;
gauss2=1/2+1/sqrt(3)/2;
v=4*gauss1*(1-gauss1);
dv=4-8*gauss1;
w1=(1-gauss1)*(1-2*gauss1);
w2=(1-gauss2)*(1-2*gauss2);
dw1=4*gauss1-3;
dw2=4*gauss2-3;
b(1)=f(1)*v*h(1)/2+
b(1)=b(1)+f(2)*v*h(1)/2+
b(2)=f(1)*w2*h(1)/2+
b(2)=b(2)+f(2)*w1*h(1)/2+
for j=2:ni-1;
b(2*j-2)=b(2*j-2)+f(2*j-1)*w1*h(j)/2+
b(2*j-2)=b(2*j-2)+f(2*j)*w2*h(j)/2+
b(2*j-1)=f(2*j-1)*v*h(j)/2+
b(2*j-1)=b(2*j-1)+f(2*j)*v*h(j)/2+
b(2*j)=f(2*j-1)*w2*h(j)/2+
b(2*j)=b(2*j)+f(2*j)*w1*h(j)/2+
end
b(2*ni-2)=b(2*ni-2)+f(2*ni-1)*w1*h(ni)/2+
b(2*ni-2)=b(2*ni-2)+f(2*ni)*w2*h(ni)/2+
b(2*ni-1)=f(2*ni-1)*v*h(ni)/2+
b(2*ni-1)=b(2*ni-1)+f(2*ni)*v*h(ni)/2+
%
if bctype(1)==0
b(1)=b(1)-(-d(1)/h(1)-c(1)/2+r(1)*h(1)/6)*g(1);
b(1)=b(1)-csd(1)*(Dd(1)-c(1)+r(1)*h(1)/2)*g(1);
b(2)=b(2)-
b(2)=b(2)-
else
b=[g(1),b];
end
if bctype(2)==0
ndf=2*ni-1+bctype(1);
b(ndf-1)=b(ndf-1)-
b(ndf-1)=b(ndf-1)-
b(ndf)=b(ndf)-(-d(ni)/h(ni)+c(ni)/2+r(ni)*h(ni)/6)*g(2);
b(ndf)=b(ndf)-csd(ni)*(Dd(ni)-c(ni)-r(ni)*h(ni)/2)*g(2);
else
b=[b,g(2)];
end
end %of if p branch
