subroutine planarfit(nsamp,x,y,z,a,b,d,fatbeta)
!fit a plane in the form of ax+by+cz+d=0 where c is set to 1.
!Because we set c=1, we cannot fit for vertical planes in parallel to the 
!x-z or y-z planes. But perfectly flat planes in parallel to the x-y plane can be included.
!no initial guess for a, b and d are needed
      implicit none
      integer nsamp,i
      double precision x(nsamp),y(nsamp),z(nsamp),
     &  a,b,d,fatbeta,p1,p2,p3,q1,q2,q3,xmean,ymean,
     &  zmean,x2mean,y2mean,xymean,xzmean,yzmean
!--------------------------------------------------------
      xmean=0.0d0
      ymean=0.0d0
      zmean=0.0d0
      x2mean=0.0d0
      y2mean=0.0d0
      xymean=0.0d0
      xzmean=0.0d0
      yzmean=0.0d0
      do i=1,nsamp
        xmean=xmean+x(i)
        ymean=ymean+y(i)
        zmean=zmean+z(i)
        x2mean=x2mean+x(i)*x(i)
        y2mean=y2mean+y(i)*y(i)
        xymean=xymean+x(i)*y(i)
        xzmean=xzmean+x(i)*z(i)
        yzmean=yzmean+y(i)*z(i)
      enddo
      xmean=xmean/dble(nsamp)
      ymean=ymean/dble(nsamp)
      zmean=zmean/dble(nsamp)
      x2mean=x2mean/dble(nsamp)
      y2mean=y2mean/dble(nsamp)
      xymean=xymean/dble(nsamp)
      xzmean=xzmean/dble(nsamp)
      yzmean=yzmean/dble(nsamp)
      p1=x2mean-xmean*xmean
      p2=xymean-xmean*ymean
      p3=xmean*zmean-xzmean
      q1=xymean-xmean*ymean
      q2=y2mean-ymean*ymean
      q3=ymean*zmean-yzmean
      call linearsys_dim2(p1,p2,p3,q1,q2,q3,a,b)
      d=-(a*xmean+b*ymean+zmean)
      fatbeta=0.0d0
      do i=1,nsamp
        fatbeta=fatbeta+(a*x(i)+b*y(i)+z(i)+d)*
     &    (a*x(i)+b*y(i)+z(i)+d)
      enddo
      return
      end