subroutine addq(lquad,nquad,xyquad,ipoin,coeg)
      implicit double precision (a-h,o-z)
      dimension lquad(7,*),xyquad(3,*),coeg(2,*),ipnt(4)
*
*     add to quadtree
*
*     lquad(1-4,iquad)   the points  (lquad(7,iquad).gt.0)
*                        new quads   (lquad(7,iquad).lt.0)
*     lquad(5,iquad)     position from which iquad came
*     lquad(6,iquad)     quad from which iquad came
*     lquad(7,iquad)     >0 number of points,=0 empty,<0 full
*
      x=coeg(1,ipoin)
      y=coeg(2,ipoin)
*
*     find quad
*
      iquad=1
10    continue
      call findp(lquad,x,y,iquad,xyquad)
*
*     to be put in quad no. iquad
*
      num=lquad(7,iquad)+1
      if(num.eq.5) goto 40
      goto 9000
40    continue
*
*     quad full, divide
*
      lquad(7,iquad)=-1
      do 50 i=1,4
      lquad(5,nquad+i)=i
      lquad(6,nquad+i)=iquad
*
      a=-i*i+5.d0*i-5.d0
      b=-2.d0*i*i*i/3.d0+5.d0*i*i-31.d0*i/3.d0+5.d0
      xyquad(1,nquad+i)=xyquad(1,iquad)+a*xyquad(3,iquad)/4.d0
      xyquad(2,nquad+i)=xyquad(2,iquad)+b*xyquad(3,iquad)/4.d0
      xyquad(3,nquad+i)=xyquad(3,iquad)/2.d0
50    continue
*
*     redistribute points
*
      do 75 ip=1,4
      ipnt(ip)=lquad(ip,iquad)
75    lquad(ip,iquad)=nquad+ip
*
      do 80 ip=1,4
      iact=ipnt(ip)
      xm=coeg(1,iact)
      ym=coeg(2,iact)
      itemp=iquad
      call findp(lquad,xm,ym,itemp,xyquad)
      iquad2=itemp
      num=lquad(7,itemp)+1
      lquad(num,itemp)=iact
      lquad(7,itemp)=num
80    continue
      nquad=nquad+4
      goto 10
9000  continue
      lquad(num,iquad)=ipoin
      lquad(7,iquad)=num
      return
      end
      subroutine downq(xm,ym,dm,inr)
      implicit double precision (a-h,o-z)
*
*     in father, out son inr
*     (xm,ym) midpoint, dm sidelenght
*
      a=-inr*inr+5.d0*inr-5.d0
      b=-2.d0*inr*inr*inr/3.d0+5.d0*inr*inr-31.d0*inr/3.d0+5.d0
      xm=xm+a*dm/4.d0
      ym=ym+b*dm/4.d0
      dm=dm/2.d0
      return
      end
      subroutine findp(lquad,x,y,iquad,xyquad)
      implicit double precision (a-h,o-z)
      dimension lquad(7,*),xyquad(3,*)
*
*     find iquad in which (x,y) is situated
*
10    continue
      if(lquad(7,iquad).ge.0) goto 9000
*
*     quad divided
*
      xm=xyquad(1,iquad)
      ym=xyquad(2,iquad)
      dm=xyquad(3,iquad)
*
      if(x.le.xm) then
      it=4
      if(y.le.ym) it=1
      else
      it=3
      if(y.le.ym) it=2
      endif
      iquad=lquad(it,iquad)
*
      goto 10
9000  continue
      return
      end