MINICERN/mathlib/gen/e/dspps2.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:26  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       FUNCTION DSPPS2(KX,KY,MX,MY,NDERX,NDERY,X,Y,TX,TY,C,NDIMC,W,NERR)
  11
  12 #include "gen/imp64.inc"
  13       DIMENSION TX(*),TY(*),C(NDIMC,*),W(*),BX(27),BY(27)
  14       CHARACTER NAME*(*)
  15       CHARACTER*80 ERRTXT
  16       PARAMETER (NAME = 'DSPPS2')
  17
  18 ************************************************************************
  19 *   NORBAS, VERSION: 15.03.1993
  20 ************************************************************************
  21 *
  22 *   DSPPS2 COMPUTES FUNCTION VALUES, VALUES OF DERIVATIVES, AND THE
  23 *   VALUE OF INTEGRAL, RESPECTIVELY, OF A TWO-DIMENSIONAL POLYNOMIAL
  24 *   SPLINE  S(X,Y)  IN REPRESENTATION OF NORMALIZED TWO-DIMENSIONAL
  25 *   B-SPLINES  B(I,J)(X,Y)
  26 *
  27 *     S(X,Y) = SUMME(I=1,...,MX-KX-1)
  28 *              SUMME(J=1,...,MY-KY-1) C(I,J) * B(I,J)(X,Y) .
  29 *
  30 *   THE TWO-DIMENSIONAL B-SPLINES  B(I,J)(X,Y)  ARE THE PRODUCT OF TWO
  31 *   ONE-DIMENSIONAL B-SPLINES  BX , BY
  32 *          B(I,J)(X,Y) = BX(I,KX)(X) * BY(J,KY)(Y)
  33 *   OF DEGREE  KX AND KY ( 0 <= KX , KY <= 25 )   WITH INDICES  I , J
  34 *   ( 1 <= I <= MX-KX-1 , 1 <= J <= MY-KY-1 ) OVER TWO SETS OF SPLINE-
  35 *   KNOTS
  36 *       TX(1),TX(2),...,TX(MX)   ( MX >= 2*KX+2 )
  37 *       TY(1),TY(2),...,TY(MY)   ( MY >= 2*KY+2 ) ,
  38 *   RESPECTIVELY.
  39 *   FOR FURTHER DETAILS TO THE ONE-DIMENSIONAL NORMALIZED B-SPLINES SEE
  40 *   THE COMMENTS TO  DSPNB1.
  41 *
  42 *   C(I,J) (I=1,...,MX-KX-NDERX-1 , J=1,...,MY-KY-NDERY-1) MUST CONTAIN
  43 *   THE (I,J)-TH C   COEFFICIENT OF THE POLYNOMIAL SPLINE  S(X,Y)  OR
  44 *   OF ONE OF ITS PARTIAL DERIVATIVE, REPRESENTED BY NORMALIZED
  45 *   TWO-DIMENSIONAL B-SPLINES OF DEGREE  (KX-NDERX) AND (KY-NDERY),
  46 *   RESPECTIVELY.
  47 *
  48 *   FOR TRANSFORMATION THE COEFFICIENTS OF THE POLYNOMIAL SPLINE  S(X,Y)
  49 *   TO THE CORRESPONDING COEFFICIENTS OF THE  NDERX-TH  AND  NDERY-TH
  50 *   PARTIAL DERIVATIVE OF  S(X,Y)  THE ROUTINE  DSPCD2  MAY BE USED.
  51 *
  52 *   ESPECIALLY FOR COMPUTING THE COEFFICIENTS  C(I,J)  OF THE TWO-
  53 *   DIMENSIONAL POLYNOMIAL VARIATION DIMINISHING SPLINE APPREOXIMATION
  54 *   OF A USER SUPPLIED FUNCTION  Z = F(X,Y)  THE ROUTINE  DSPVD2  MAY BE
  55 *   USED.
  56 *
  57 *   PARAMETERS:
  58 *
  59 *   KX    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN X-DIRECTION
  60 *         OVER THE SET OF KNOTS  TX.
  61 *   KY    (INTEGER) DEGREE OF ONE-DIMENSIONAL B-SPLINES IN Y-DIRECTION
  62 *         OVER THE SET OF KNOTS  TY.
  63 *   MX    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN X-DIRECTION.
  64 *   MY    (INTEGER) NUMBER OF KNOTS FOR THE B-SPLINES IN Y-DIRECTION.
  65 *   NDERX (INTEGER) ON ENTRY, NDERX MUST CONTAIN AN INTEGER VALUE >= -1.
  66 *         = -1: DSPPS2 COMPUTES THE INTEGRAL OF  S(TAU,Y)  OVER THE
  67 *               RANGE  TAU <= X.
  68 *         =  0: DSPPS2 COMPUTES THE FUNCTION VALUE  S(X,Y)  FOR
  69 *               FOR THE SPECIFIED VALUES OF  X,Y.
  70 *         >= 1: DSPNB2 COMPUTES THE VALUE OF THE NDERX-TH PARTIAL
  71 *               DERIVATIVE OF  S(X,Y)  WITH RESPECT TO X
  72 *               FOR THE SPECIFIED VALUES OF  X,Y.
  73 *               (IF  NDERX > KX  ZERO RETURNS).
  74 *   NDERY (INTEGER) ON ENTRY, NDERY MUST CONTAIN AN INTEGER VALUE >= -1.
  75 *         THE MEANING OF  NDERY  IS THE SAME AS THAT OF THE PARAMETER
  76 *         NDERX WITH RESPECT TO Y-DIRECTION INSTEAD OF X-DIRECTION.
  77 *   NDIMC (INTEGER) DECLARED FIRST DIMENSION OF ARRAY  C  IN THE
  78 *         CALLING PROGRAM, WITH  NDIMC >= MX-KX-NDERX-1 .
  79 *   X     (DOUBLE PRECISION) ON ENTRY, X MUST CONTAIN THE VALUE OF THE
  80 *         INDEPENDENT VARIABLE X OF  S(X,Y)
  81 *   Y     (DOUBLE PRECISION) ON ENTRY, Y MUST CONTAIN THE VALUE OF THE
  82 *         INDEPENDENT VARIABLE Y OF  S(X,Y)
  83 *   TX    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MX CONTAINING THE
  84 *         KNOTS IN X-DIRECTION, ON ENTRY.
  85 *   TY    (DOUBLE PRECISION) ARRAY OF AT LEAST ORDER MY CONTAINING THE
  86 *         KNOTS IN Y-DIRECTION, ON ENTRY.
  87 *   C     (DOUBLE PRECISION) ARRAY OF ORDER (NDIMC, >= MY-KY-NDERY-1),
  88 *         CONTAINING THE COEFFICIENTS OF THE TWO-DIMENSIONAL B-SPLINE
  89 *         REPRESENTATION OF  S(X,Y) , ON ENTRY.
  90 *   W     (DOUBLE PRECISION) WORKING ARRAY OF AT LEAST ORDER  MY-KY-1.
  91 *   NERR  (INTEGER) ERROR INDICATOR. ON EXIT:
  92 *         = 0: NO ERROR DETECTED
  93 *         = 1: AT LEAST ONE OF THE CONSTANTS KX , KY , MX , MY ,
  94 *              NDERX , NDERY  IS ILLEGAL
  95 *
  96 *   USAGE:
  97 *
  98 *   THE FUNCTION-CALL
  99 *         Z = DSPPS2(KX,KY,MX,MY,NDERX,NDERY,X,Y,TX,TY,C,NDIMC,W,NERR)
 100 *   RETURNS
 101 *       - THE VALUE OF THE INTEGRAL (NDERX =-1  .AND.  NDERY =-1) OR
 102 *       - THE FUNCTION VALUE        (NDERX = 0  .AND.  NDERY = 0) OR
 103 *       - THE VALUE OF THE NDERX-TH AND NDERY-TH PARTIAL DERIVATIVE
 104 *         (NDERX >= 0  .AND.  NDERY >= 0  .AND.  NDERX + NDERY > 0)
 105 *   OF THE POLYNOMIAL SPLINE  S(X,Y)  AT (X,Y).
 106 *
 107 *   ERROR MESSAGES:
 108 *
 109 *   IF ONE OF THE FOLLOWING RELATION IS SATISFIED BY THE CHOSEN INPUT-
 110 *   PARAMETERS THE PROGRAM RETURNS, AND AN ERROR MESSAGE IS PRINTED:
 111 *     KX < 0       OR    KX > 25      OR    KY < 0   OR   KY > 25  OR
 112 *     MX < 2*KX+2  OR    MY < 2*KY+2  OR
 113 *     NDERX < -1   OR    NDERY < -1   .
 114 *
 115 ************************************************************************
 116
 117       PARAMETER (Z0 = 0 , Z1 = 1)
 118
 119       NERR=1
 120       IF(KX .LT. 0 .OR. KX .GT. 25) THEN
 121        WRITE(ERRTXT,101) 'KX',KX
 122        CALL MTLPRT(NAME,'E210.1',ERRTXT)
 123       ELSEIF(KY .LT. 0 .OR. KY .GT. 25) THEN
 124        WRITE(ERRTXT,101) 'KY',KY
 125        CALL MTLPRT(NAME,'E210.1',ERRTXT)
 126       ELSEIF(MX .LT. 2*KX+2) THEN
 127        WRITE(ERRTXT,101) 'MX',MX
 128        CALL MTLPRT(NAME,'E210.2',ERRTXT)
 129       ELSEIF(MY .LT. 2*KY+2) THEN
 130        WRITE(ERRTXT,101) 'MY',MY
 131        CALL MTLPRT(NAME,'E210.2',ERRTXT)
 132       ELSEIF(NDERX .LT. -1) THEN
 133        WRITE(ERRTXT,101) 'NDERX',NDERX
 134        CALL MTLPRT(NAME,'E210.5',ERRTXT)
 135       ELSEIF(NDERY .LT. -1) THEN
 136        WRITE(ERRTXT,101) 'NDERY',NDERY
 137        CALL MTLPRT(NAME,'E210.5',ERRTXT)
 138       ELSEIF(NDERX .EQ. -1  .AND.  NDERY .NE. -1  .OR.
 139      +       NDERX .NE. -1  .AND.  NDERY .EQ. -1      ) THEN
 140        WRITE(ERRTXT,102) 'NDERX',NDERX,'NDERY',NDERY
 141        CALL MTLPRT(NAME,'E210.6',ERRTXT)
 142       ELSE
 143
 144        NERR=0
 145        IF(NDERX .EQ. -1  .AND.  NDERY .EQ. -1) THEN
 146         DO 60 JJ=1,MY-KY-1
 147         IF(X .GE. TX(MX-KX)) THEN
 148          R=Z0
 149          DO 10 I=1,MX-KX-1
 150    10    R=R+C(I,JJ)*(TX(I+KX+1)-TX(I))
 151          R=R/(KX+1)
 152         ELSE
 153          KK=LKKSPL(X,TX(KX+1),MX-2*KX-1)+KX
 154          R=Z0
 155          DO 20 I=1,KK-KX-2
 156    20    R=R+C(I,JJ)*(TX(I+KX+1)-TX(I))
 157          R=R/(KX+1)
 158          IF(KX .EQ. 0) THEN
 159           K1=MAX(1,KK-1)
 160           R=R+(X-TX(K1))*C(K1,JJ)
 161          ELSE
 162           DO 50 I=MAX(1,KK-KX-1),KK-1
 163           CALL DVSET(KX+1,Z0,BX(1),BX(2))
 164           BX(KK-I)=1/(TX(KK)-TX(KK-1))
 165           DO 30 L=1,KX
 166           DO 30 J=MAX(1,KK-I-L),MIN(KX+1-L,KK-I)
 167           DIF=TX(I+J+L)-TX(I+J-1)
 168           B0=Z0
 169           IF(DIF .NE. 0)
 170      +    B0=((X-TX(I+J-1))*BX(J)+(TX(I+J+L)-X)*BX(J+1))/DIF
 171    30     BX(J)=B0
 172           S=Z0
 173           DO 40 L=1,KK-I
 174    40     S=S+(X-TX(I+L-1))*BX(L)
 175           S=S*(TX(I+KX+1)-TX(I))/(KX+1)
 176    50     R=R+C(I,JJ)*S
 177          ENDIF
 178         ENDIF
 179    60   W(JJ)=R
 180         IF(Y .GE. TY(MY-KY)) THEN
 181          R=Z0
 182          DO 70 I=1,MY-KY-1
 183    70    R=R+W(I)*(TY(I+KY+1)-TY(I))
 184          R=R/(KY+1)
 185         ELSE
 186          KK=LKKSPL(Y,TY(KY+1),MY-2*KY-1)+KY
 187          R=Z0
 188          DO 80 I=1,KK-KY-2
 189    80    R=R+W(I)*(TY(I+KY+1)-TY(I))
 190          R=R/(KY+1)
 191          IF(KY .EQ. 0) THEN
 192           K1=MAX(1,KK-1)
 193           R=R+(Y-TY(K1))*W(K1)
 194          ELSE
 195           DO 110 I=MAX(1,KK-KY-1),KK-1
 196           CALL DVSET(KY+1,Z0,BY(1),BY(2))
 197           BY(KK-I)=1/(TY(KK)-TY(KK-1))
 198           DO 90 L=1,KY
 199           DO 90 J=MAX(1,KK-I-L),MIN(KY+1-L,KK-I)
 200           DIF=TY(I+J+L)-TY(I+J-1)
 201           B0=Z0
 202           IF(DIF .NE. 0)
 203      +    B0=((Y-TY(I+J-1))*BY(J)+(TY(I+J+L)-Y)*BY(J+1))/DIF
 204    90     BY(J)=B0
 205           S=Z0
 206           DO 100 L=1,KK-I
 207   100     S=S+(Y-TY(I+L-1))*BY(L)
 208           S=S*(TY(I+KY+1)-TY(I))/(KY+1)
 209   110     R=R+W(I)*S
 210          ENDIF
 211         ENDIF
 212         DSPPS2=R
 213         RETURN
 214        ENDIF
 215
 216        DSPPS2=Z0
 217        IF(X  .LT. TX(KX+1)  .OR.  X  .GT. TX(MX-KX) .OR.
 218      +    Y  .LT. TY(KY+1)  .OR.  Y  .GT. TY(MY-KY) .OR.
 219      +    KX .LT. NDERX     .OR.  KY .LT. NDERY        ) RETURN
 220
 221        KKX=LKKSPL(X,TX(KX+1),MX-2*KX-1)+KX
 222        KKY=LKKSPL(Y,TY(KY+1),MY-2*KY-1)+KY
 223        E1=X-TX(KKX-1)
 224        BX(1)=Z1
 225        DO 120 J=2,KX-NDERX+1
 226   120  BX(J)=E1*BX(J-1)/(TX(KKX-2+J)-TX(KKX-1))
 227        E2=TX(KKX)-X
 228        DO 130 J=1,KX-NDERX
 229        E3=X-TX(KKX-1-J)
 230        BX(1)=E2*BX(1)/(TX(KKX)-TX(KKX-J))
 231        DO 130 L=2,KX-NDERX+1-J
 232   130  BX(L)=E3*BX(L-1)/(TX(KKX-2+L)-TX(KKX-1-J))+
 233      +       (TX(KKX-1+L)-X)*BX(L)/(TX(KKX-1+L)-TX(KKX-J))
 234        E1=Y-TY(KKY-1)
 235        BY(1)=Z1
 236        DO 140 J=2,KY-NDERY+1
 237   140  BY(J)=E1*BY(J-1)/(TY(KKY-2+J)-TY(KKY-1))
 238        E2=TY(KKY)-Y
 239        DO 150 J=1,KY-NDERY
 240        E3=Y-TY(KKY-1-J)
 241        BY(1)=E2*BY(1)/(TY(KKY)-TY(KKY-J))
 242        DO 150 L=2,KY-NDERY+1-J
 243   150  BY(L)=E3*BY(L-1)/(TY(KKY-2+L)-TY(KKY-1-J))+
 244      +       (TY(KKY-1+L)-Y)*BY(L)/(TY(KKY-1+L)-TY(KKY-J))
 245        R=Z0
 246        DO 160 I=1,KX-NDERX+1
 247        DO 160 J=1,KY-NDERY+1
 248   160  R=R+C(KKX-2-KX+I,KKY-2-KY+J)*BX(I)*BY(J)
 249        DSPPS2=R
 250       ENDIF
 251       RETURN
 252
 253   101 FORMAT(1X,A5,' =',I6,'   NOT IN RANGE')
 254   102 FORMAT(1X,A5,' =',I6,A7,' =',I6,'   INCONSISTENT')
 255       END
 256
 257
 258