subroutine curvecrossing(ncurves,nsamp,maxnsamp,xdata,ydata,
     &params,xcross,ycross)
      implicit none
!fit ncurves curves from ncurves pairs of datasets and estimate the mean crossings of the curves
      integer ncurves,nsamp(ncurves),maxnsamp,i,j,k,n,m,ndim,
     &iderivative,INFO
      double precision xdata(ncurves,maxnsamp),ydata(ncurves,maxnsamp),
     &xcross,ycross,x(maxnsamp),y(maxnsamp),a(ncurves),b(ncurves),
     &c(ncurves),d(ncurves),beta(10),betamin(10),betamax(10),
     &weitx(maxnsamp),weity(maxnsamp),xmin(maxnsamp),xmax(maxnsamp),
     &y_pred(maxnsamp),x_pred(maxnsamp),sumsquare,
     &params(ncurves,4)

      k=0
      do i=1,ncurves
        do j=1,4
          params(i,j)=-9999.0d0
        enddo
        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
          beta(1)=30.0d0
          betamin(1)=0.0d0
          betamax(1)=200.0d0
          beta(2)=-4.0d0
          betamin(2)=-100.0d0
          betamax(2)=0.0d0
          beta(3)=8.0d0
          betamin(3)=0.0d0
          betamax(3)=200.0d0
          beta(4)=-0.5d0
          betamin(4)=-10.0d0
          betamax(4)=0.0d0
          do j=1,n
            weitx(j)=1.0d0
            weity(j)=1.0d0
            xmin(j)=0.0d0
            xmax(j)=x(j)+20.0d0
          enddo
          ndim=4
          iderivative=1
          call GenericRegres(n,1,y,1,x,weity,weitx,ndim,beta,betamin,
     &betamax,xmin,xmax,iderivative,0,y_pred,x_pred,sumsquare)
          a(k)=beta(1)
          b(k)=beta(2)
          c(k)=beta(3)
          d(k)=beta(4)
          do j=1,ndim
            params(i,j)=beta(j)
          enddo
        endif
      enddo
      if(k.gt.0)then
        call curmeancrossing(k,a,b,c,d,xcross,ycross)
      else
        xcross=-9999.0d0
        ycross=-9999.0d0
      endif
      return
      end

      subroutine curmeancrossing(ncurves,a,b,c,d,xcross,ycross)
      implicit none
!calculate the average crossing point of the ncurves curves.
      integer ncurves,i,j,k
      double precision a(ncurves),b(ncurves),c(ncurves),d(ncurves),
     &xcross,ycross,x(ncurves),y(ncurves),a1,b1,c1,d1,a2,b2,c2,d2,
     &p,q,u,root1,root2,small,term,x1_root1,x2_root1,x1_root2,x2_root2,
     &x_root1,x_root2
      parameter(small=1.0d-7)
      k=0
      do i=1,ncurves-1
        do j=i+1,ncurves
          a1=a(i)
          b1=b(i)
          c1=c(i)
          d1=d(i)
          a2=a(j)
          b2=b(j)
          c2=c(j)
          d2=d(j)
          term=dabs(a1-a2)+dabs(b1-b2)+dabs(c1-c2)
          if(term.lt.small)goto 100
          p=a1-a2+d1-d2
          q=a1*(c2+b1)-a2*(c1+b2)+(d1-d2)*(c2+c1)
          u=a1*b1*c2-a2*b2*c1+(d1-d2)*c1*c2
          call quadraticroots(p,q,u,root1,root2)
          if(root1.gt.0.0d0.or.root2.gt.0.0d0)then
            k=k+1
            if(root1.gt.0.0d0.and.root2.gt.0.0d0)then
              x(k)=root1
            else
              if(root1.gt.0.0d0)then
                x(k)=root1
              else
                x(k)=root2
              endif
            endif
            y(k)=a1*(x(k)+b1)/(x(k)+c1)+d1
          else
            if(dabs(b1-c1).gt.small)then
              term=a2*(b2-c2)/(a1*(b1-c1))
              if(term.ge.0.0d0)then
                term=dsqrt(term)
                root1=((d1+a1)*term+d2+a2)/(term+1.0d0)
                if(dabs(term-1.0d0).gt.small)then
                  root2=((d1+a1)*term-d2-a2)/(term-1.0d0)
                else
                  root2=-9999.0d0
                endif
                x1_root1=-c1+a1*(b1-c1)/(root1-d1-a1)
                x2_root1=-c2+a2*(b2-c2)/(root1-d2-a2)
                if(dabs(root2+9999.0d0).gt.small)then
                  x1_root2=-c1+a1*(b1-c1)/(root2-d1-a1)
                  x2_root2=-c2+a2*(b2-c2)/(root2-d2-a2)
                else
                  x1_root2=-9999.0d0
                  x2_root2=-9999.0d0
                endif
                if(x1_root1.gt.0.0d0.and.x2_root1.gt.0.0d0)then
                  x_root1=(x1_root1+x2_root1)/2.0d0
                else
                  x_root1=-9999.0d0
                endif
                if(x1_root2.gt.0.0d0.and.x2_root2.gt.0.0d0)then
                  x_root2=(x1_root2+x2_root2)/2.0d0
                else
                  x_root2=-9999.0d0
                endif
                if(x_root1.gt.0.0d0.or.x_root2.gt.0.0d0)then
                  k=k+1
                  if(x_root1.gt.0.0d0.and.x_root2.gt.0.0d0)then
                    x(k)=x_root2
                    y(k)=root2
                  else
                    if(x_root1.gt.0.0d0)then
                      x(k)=x_root1
                      y(k)=root1
                    else
                      x(k)=x_root2
                      y(k)=root2
                    endif
                  endif
                endif
              endif
            endif
          endif
100       continue
        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