subroutine linuncertainty(nsamp,a1_mean,a1_sigma,b1_mean,
     &b1_sigma,y1_mean,y1_sigma,a2_mean,a2_sigma,b2_mean,b2_sigma,
     &y2_mean,y2_sigma,x1,x2,x1_mean,x1_sigma,x2_mean,x2_sigma)
      implicit none
!given:
!y1=a1*x1+b1*x2
!y2=a2*x1+b1*x2
!y1 ~ (y1_mean, y1_sigma)
!y2 ~ (y2_mean, y2_sigma)
!a1 ~ (a1_mean, a1_sigma)
!b1 ~ (b1_mean, b1_sigma)
!a2 ~ (a2_mean, a2_sigma)
!b2 ~ (b2_mean, b2_sigma)
!find:
!x1 ~ (x1_mean, x1_sigma)
!x2 ~ (x2_mean, x2_sigma)
      integer nsamp,i
      double precision a1_mean,a1_sigma,b1_mean,b1_sigma,y1_mean,
     &y1_sigma,a2_mean,a2_sigma,b2_mean,b2_sigma,y2_mean,y2_sigma,
     &x1(nsamp),x2(nsamp),x1_mean,x1_sigma,x2_mean,x2_sigma,gasdev2,
     &a1,b1,a2,b2,y1,y2
      x1_mean=0.0d0
      x2_mean=0.0d0
      do i=1,nsamp
        y1=gasdev2()*y1_sigma+y1_mean
        y2=gasdev2()*y2_sigma+y2_mean
        a1=gasdev2()*a1_sigma+a1_mean
        b1=gasdev2()*b1_sigma+b1_mean
        a2=gasdev2()*a2_sigma+a2_mean
        b2=gasdev2()*b2_sigma+b2_mean
        x1(i)=(b2*y1-b1*y2)/(a1*b2-a2*b1)
        x2(i)=(a1*y2-a2*y1)/(a1*b2-a2*b1)
        x1_mean=x1_mean+x1(i)
        x2_mean=x2_mean+x2(i)
      enddo
      x1_mean=x1_mean/dble(nsamp)
      x2_mean=x2_mean/dble(nsamp)
      x1_sigma=0.0d0
      x2_sigma=0.0d0
      do i=1,nsamp
        x1_sigma=x1_sigma+(x1(i)-x1_mean)*(x1(i)-x1_mean)
        x2_sigma=x2_sigma+(x2(i)-x2_mean)*(x2(i)-x2_mean)
      enddo
      x1_sigma=x1_sigma/dble(nsamp)
      x2_sigma=x2_sigma/dble(nsamp)
      return
      end