*
* $Id: dspcd2.F,v 1.1.1.1 1996/04/01 15:02:25 mclareni Exp $
*
* $Log: dspcd2.F,v $
* Revision 1.1.1.1  1996/04/01 15:02:25  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DSPCD2(KX,KY,MX,MY,NDERX,NDERY,TX,TY,C,NDIMC,D,NERR)

#include "gen/imp64.inc"
      DIMENSION TX(*),TY(*),C(NDIMC,*),D(NDIMC,*)
      CHARACTER NAME*(*)
      CHARACTER*80 ERRTXT
      PARAMETER (NAME = 'DSPCD2')

************************************************************************
*   NORBAS, VERSION: 15.03.1993
************************************************************************
*
*   DSPCD2 COMPUTES FROM GIVEN COEFFICIENTS
*          C(I,J)   (I=1,...,MX-KX-1 , J=1,...,MY-KY-1)
*   OF A TWO-DIMENSIONAL POLYNOMIAL SPLINE  S(X,Y) IN REPRESENTATION OF
*   NORMALIZED TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)
*
*     S(X,Y) = SUMME(I=1,...,MX-KX-1)
*              SUMME(J=1,...,MY-KY-1) C(I,J) * B(I,J)(X,Y)
*
*   THE CORRESPONDING COEFFICIENTS
*       D(I,J)   (I=1,...,MX-KX-NDERX-1 , J=1,...,MY-KY-NDERY-1)
*   OF THE  NDERX-TH , NDERY-TH  PARTIAL DERIVATIVE OF  S(X,Y).
*
*   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
*   ONE-DIMENSIONAL B-SPLINES BX , BY
*          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
*   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
*   ( 1 <= I <= MX-KX-1 , 1 <= J <= MY-KY-1 ) OVER TWO SETS OF SPLINE-
*   KNOTS
*       TX(1),TX(2),...,TX(MX)   ( MX >= 2*KX+2 )
*       TY(1),TY(2),...,TY(MY)   ( MY >= 2*KY+2 ) ,
*   RESPECTIVELY.
*   FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
*   THE COMMENTS TO  DSPNB1.
*
*   PARAMETERS:
*
*   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
*         OVER THE SET OF KNOTS  TX.
*   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
*         OVER THE SET OF KNOTS  TY.
*   MX    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN X-DIRECTION.
*   MY    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN Y-DIRECTION.
*   NDERX (INTEGER) ORDER OF PARTIAL DERIVATIVE IN X-DIRECTION.
*   NDERY (INTEGER) ORDER OF PARTIAL DERIVATIVE IN Y-DIRECTION.
*   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAYS  C AND D  IN THE
*         CALLING PROGRAM, WITH  NDIMC >= MX-KX-1 .
*   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX CONTAINING THE
*         KNOTS IN X-DIRECTION, ON ENTRY.
*   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY CONTAINING THE
*         KNOTS IN Y-DIRECTION, ON ENTRY.
*   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-1).
*         ON ENTRY C(I,J) MUST CONTAIN THE (I,J)-TH COEFFICIENT OF THE
*         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF  S(X,Y) .
*   D     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMD, >= MY-KY-1).
*         ON EXIT  D(I,J)  CONTAINS THE (I,J)-TH COEFFICIENT OF THE
*         TWO-DIMENSIONAL B-SPLINE REPRESENTATION OF THE NDERX-TH,
*         NDERY-TH PARTIAL DERIVATIVE OF  S(X,Y).
*   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
*         = 0: NO ERROR DETECTED
*         = 1: AT LEAST ONE OF THE CONSTANTS KX , KY , MX , MY , NDERX ,
*              NDERY  IS ILLEGAL
*
*   ERROR MESSAGES:
*
*   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
*   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
*     KX < 0       OR    KX > 25      OR    KY < 0   OR   KY > 25   OR
*     MX < 2*KX+2  OR    MY < 2*KY+2  OR
*     NDERX < 0    OR    NDERY < 0    OR
*     NDERX < 1    AND   NDERY < 1    OR
*     NDERX > KX   OR    NDERY > KY .
*
************************************************************************

      PARAMETER (Z0 = 0)

      NERR=1
      IF(KX .LT. 0 .OR. KX .GT. 25) THEN
       WRITE(ERRTXT,101) 'KX',KX
       CALL MTLPRT(NAME,'E210.1',ERRTXT)
      ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
       WRITE(ERRTXT,101) 'KY',KY
       CALL MTLPRT(NAME,'E210.1',ERRTXT)
      ELSEIF(MX .LT. 2*KX+2) THEN
       WRITE(ERRTXT,101) 'MX',MX
       CALL MTLPRT(NAME,'E210.2',ERRTXT)
      ELSEIF(MY .LT. 2*KY+2) THEN
       WRITE(ERRTXT,101) 'MY',MY
       CALL MTLPRT(NAME,'E210.2',ERRTXT)
      ELSEIF(NDERX .LT. 0 .OR. NDERX .GT. KX) THEN
       WRITE(ERRTXT,101) 'NDERX',NDERX
       CALL MTLPRT(NAME,'E210.5',ERRTXT)
      ELSEIF(NDERY .LT. 0 .OR. NDERY .GT. KY) THEN
       WRITE(ERRTXT,101) 'NDERY',NDERY
       CALL MTLPRT(NAME,'E210.5',ERRTXT)
      ELSEIF(NDERX .LT. 1 .AND. NDERY .LT. 1)THEN
       WRITE(ERRTXT,102) 'NDERX',NDERX,'NDERY',NDERY
       CALL MTLPRT(NAME,'E210.6',ERRTXT)
      ELSE

      NERR=0
      CALL DMCPY(MX-KX-1,MY-KY-1,C(1,1),C(1,2),C(2,1),
     +                           D(1,1),D(1,2),D(2,1))
      IF(NDERX .GT. 0) THEN
       DO 10 J = 1,MY-KY-1
       DO 10 L = 1,NDERX
       A=KX-L+1
       DO 10 I = 1,MX-KX-1-L
       DIF=TX(I+KX+1)-TX(I+L)
       D0=Z0
       IF(DIF .NE. Z0) D0=A*(D(I+1,J)-D(I,J))/DIF
   10  D(I,J)=D0
      ENDIF

      IF(NDERY .GT. 0) THEN
       DO 20 I = 1,MX-KX-1
       DO 20 L = 1,NDERY
       A=KY-L+1
       DO 20 J = 1,MY-KY-1-L
       DIF=TY(J+KY+1)-TY(J+L)
       D0=Z0
       IF(DIF .NE. Z0) D0=A*(D(I,J+1)-D(I,J))/DIF
   20  D(I,J)=D0
      ENDIF
      ENDIF

      RETURN

  101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
  102 FORMAT(1X,A5,' =',I6,A7,' =',I6,'   INCONSISTENT')
      END