c	INCLUDE '../SOURCE95/mod_splines95.f'

c*********************************************************************

	PROGRAM test_sbsp_dis
	
c	test pour sbsp_dis, version F95 de sbsp_dis

c	Auteur: P. Morel, Departement J.D. Cassini, O.C.A.
	
c--------------------------------------------------------------------
	
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, bsp_dis
      
	IMPLICIT NONE
	
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: f, fd
	REAL (kind=dp), ALLOCATABLE, DIMENSION(:) :: dfdx, fx, x, xt
	
	REAL (kind=sp), ALLOCATABLE, DIMENSION(:) :: xd, yd, yc
	REAL (kind=sp) :: pas, xmax, xmin, ymax, ymin
	INTEGER, ALLOCATABLE, DIMENSION(:) :: id
	INTEGER :: i, knot, l=1, m, n, nd, ndes, nx
	
	CHARACTER (len=3) :: device='?'
	
c----------------------------------------------------------------
	
2000	FORMAT(8es10.3)

	nx=10 ; nd=2 ; n=2 ; m=4 ; nx=MAX(nx,3*m+10)

 	ALLOCATE(xt(nx+nd+m),x(nx),fd(n,nd),f(n,nx+nd),
 	1 dfdx(n),fx(n),id(0:nd+1),xd(nx),yd(nx+nd))
 
	id=-100 ; id(1)=m+3 ; id(2)=nx-(m+3) ; f=-100.d0

	DO i=1,nx
	 x(i)=i ; f(1,i)=x(i)**2 ; f(2,i)=7.d0
	ENDDO

	DO i=id(1)+1,nx			!premiere discontinuite
	 f(1,i)=f(1,i)+10.d0 ; f(2,i)=f(2,i)+12.d0
	ENDDO
	fd(1,1)=f(1,id(1))+10.d0 ; fd(2,1)=f(2,id(1))+12.d0	!a droite

	DO i=id(2)+1,nx			!seconde discontinuite
	 f(1,i)=f(1,i)+15.d0 ; f(2,i)=f(2,i)+20.d0
	ENDDO
	fd(1,2)=f(1,id(2))+15.d0 ; fd(2,2)=f(2,id(2))+20.d0	!a droite

c	PRINT*,'nd,id/les f/les fd',nd,id ; WRITE(*,2000)x
c	DO l=1,n
c	 WRITE(*,2000)f(l,1:nx)
c	ENDDO	!i
c	DO l=1,n
c	 WRITE(*,2000)fd(l,1:nd)
c	ENDDO	!i
c	PAUSE'données'
	
c	PRINT*,'pour le device'
c	PRINT*,'sous window entrer: /WS, /NULL ou ?'
c	PRINT*,'sous UNIX/LINUX entrer: /XW, /NULL ou ?'
	device='/xw'
	CALL pgbegin(0,device,1,1)
	xmin=0. ; xmax=MAXVAL(x)*1.05 ; ymax=MAXVAL(f)*1.05 ; ymin=0.
	CALL pgsls(1) ; CALL pgenv(xmin,xmax,ymin,ymax,0,0)
	CALL pglabel('X','Y','bsp_dis')	

	xd=x ; yd(:)=f(1,:)	
	WRITE(*,2000)xd ; WRITE(*,2000)yd
	CALL pgpoint(nx,xd,yd,2) !; PAUSE'f(1,nx)'

	yd(:)=f(2,:)
	CALL pgpoint(nx,xd,yd,3) !; WRITE(*,2000)yd ; PAUSE'f(2,nx)'	
		
	CALL bsp_dis(n,x,f,nd,id,fd,nx,m,xt,knot)
	
	ndes=1000
	DEALLOCATE(xd,yd) ; ALLOCATE(xd(ndes),yd(ndes),yc(ndes))
	pas=(x(nx)-x(1))/float(ndes-1)
	DO i=1,ndes
	 xd(i)=x(1)+pas*(i-1)
	 CALL bsp1dn(n,f,x,xt,nx,m,knot,.TRUE.,DBLE(xd(i)),l,fx,dfdx)
	 yd(i)=fx(1) ; yc(i)=fx(2)
c	 WRITE(*,2000)xd(i),yd(i),yc(i)
	ENDDO
	CALL pgline(ndes,xd,yd) ; CALL pgsci(3) ; CALL pgline(ndes,xd,yc)

	CALL pgend	

	STOP

	END PROGRAM test_sbsp_dis