subroutine linecrossing(nlines,nsamp,maxnsamp,xdata,ydata,
     &params,xcross,ycross)
      implicit none
!fit nlines lines from nlines pairs of datasets and estimate the mean crossings of the lines
      integer nlines,nsamp(nlines),maxnsamp,i,j,k,n,m
      double precision xdata(nlines,maxnsamp),ydata(nlines,maxnsamp),
     &xcross,ycross,x(maxnsamp),y(maxnsamp),intercepts(nlines),
     &slopes(nlines),params(nlines,2)

      k=0
      do i=1,nlines
        n=0
        do j=1,nsamp(i)
          if(dabs(xdata(i,j)+9999.0d0).gt.1.0d-8.and.
     &dabs(ydata(i,j)+9999.0d0).gt.1.0d-8)then
            n=n+1
            x(n)=xdata(i,j)
            y(n)=ydata(i,j)
          endif
        enddo
        if(n.gt.1)then
          k=k+1
          call y_aPLUSbx(n,x,y,intercepts(k),slopes(k))
          params(i,1)=slopes(k)
          params(i,2)=intercepts(k)
        endif
      enddo
      if(k.gt.0)then
        call meancrossing(k,slopes,intercepts,xcross,ycross)
      else
        xcross=-9999.0d0
        ycross=-9999.0d0
      endif
      return
      end

      subroutine meancrossing(nlines,slopes,intercepts,xcross,ycross)
      implicit none
!calculate the average crossing point of the nlines lines. Parrell lines are excluded
      integer nlines,i,j,k
      double precision slopes(nlines),intercepts(nlines),xcross,ycross,
     &x(nlines),y(nlines),a1,b1,a2,b2
      k=0
      do i=1,nlines-1
        do j=i+1,nlines
          a1=slopes(i)
          b1=intercepts(i)
          a2=slopes(j)
          b2=intercepts(j)
          if(dabs(a1-a2).gt.1.0d-8)then
            k=k+1
            x(k)=(b2-b1)/(a1-a2)
            y(k)=a1*x(k)+b1
          endif
        enddo
      enddo
      if(k.gt.0)then
        xcross=0.0d0
        ycross=0.0d0
        do i=1,k
          xcross=xcross+x(i)
          ycross=ycross+y(i)
        enddo
        xcross=xcross/dble(k)
        ycross=ycross/dble(k)
      else
        xcross=-9999.0d0
        ycross=-9999.0d0
      endif
      return
      end