!****************************************************************
      
      function bspder(x,c,y,n,m,l,cd,cd1,der,init)
!     Version: 02/12/99
      
!     calcule la derivee d'ordre der=d-1 au point x
!     y(l).le.x.le.y(l+1)
!     de la spline d'ordre m s(x):=somme(ci*Bi(x))i=1,n
!     Bi : i-ieme B-spline d'ordre m.le.6 (non teste),
!     les ci:=c(i) etant connus
!     cd1(n) : tableau de travail
!     cd(m,n) : coefficients calcules au premier appel (init=.false.)
!     d'apres schumaker algorithmes 5-10 et 5-11 p.196
      
!     Auteur: P. Morel, Departement J.D. Cassini, O.C.A., Observatoire de Nice
!     Version: 04/01/95
      
      implicit none
      
      integer n,m,l,der,i,j,d,mj
      
      real*8 c(1),cd(1),y(1),cd1(1),bspder,x,denom,schu58
      
      logical init
      
!------------S-----------A-----------V------------E--------------------
      save 
!----------------------------------------------------------------------
!     initialisation : calcul du tableau cd(j,i) j=1,m i=j,n
      
      if(.not. init)then
       do i=1,n
        cd(m*(i-1)+1)=c(i)
       enddo
       do j=2,m
        mj=m-j+1
        do i=n,j,-1
         denom=y(i+mj)-y(i)
         cd(m*(i-1)+j)=0.
         if(denom.ne.0.)cd(m*(i-1)+j)=float(mj)/denom*                  &
     &        (cd(m*(i-1)+j-1)-cd(m*(i-2)+j-1))
        enddo
       enddo
      endif
      
      d=der+1
      if(d .lt. 1 .or. d .gt. m)then
       write(*,*)'dans bspder d=',d,'   n est pas entre 1 et m=',m
       stop
      endif
      
      do i=1,n
       cd1(i)=cd(m*(i-1)+d)
      enddo
      
      bspder=schu58(x,y,m-d+1,l,cd1)
      
      return
      
      end