C ALGORITHM 626 COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 4, DEC., 1984,
C DEC., 1984, P. 473.
C PROGRAM TPDEM1
C
C T R I ANGLE C ONTOUR P LOTTING
C
C DEMONSTRATION OF THE TRICP SUBROUTINES
C (MODE=0 AND MODE=2)
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
C ONLY (STANDARD)-CALCOMP PLOT COMMANDS ARE USED
C
C ND= NUMBER OF DATA POINTS
C
DIMENSION WK(150),IWK(930)
C WK(5*ND),IWK(31*ND)
C
DIMENSION XD(30),YD(30),ZD(30),C(20),X(3),Y(3)
DATA ND/30/
C
C EXAMPLE 1. IS TAKEN FROM ALGORITHM 526 BY H.AKIMA
C (FIGURE 1E OF THE CORRESPONDING PAPER)
C CONTOUR LINE WITH VALUE C=15. IS PLOTTED IN ADDITION
C
DATA XD
1 /11.16, 24.20, 19.85, 10.35, 19.72, 0.00,
2 20.87, 19.99, 10.28, 4.51, 0.00, 16.70,
3 6.08, 25.00, 14.90, 0.00, 9.66, 5.22,
4 11.77, 15.10, 25.00, 25.00, 14.59, 15.20,
5 5.23, 2.14, 0.51, 25.00, 21.67, 3.31/
DATA YD
1 / 1.24, 16.23, 10.72, 4.11, 1.39, 20.00,
2 20.00, 4.62, 15.16, 20.00, 4.48, 19.65,
3 4.58, 11.87, 3.12, 0.00, 20.00, 14.66,
4 10.47, 17.19, 3.87, 0.00, 8.71, 0.00,
5 10.72, 15.03, 8.37, 20.00, 14.36, 0.13/
DATA ZD
1 /22.15, 2.83, 7.97, 22.33, 16.83, 34.60,
2 5.74, 14.72, 21.59, 15.61, 61.77, 6.31,
3 35.74, 4.40, 21.70, 58.20, 4.73, 40.36,
4 13.62, 12.57, 8.74, 12.00, 14.81, 21.60,
5 26.50, 53.10, 49.43, 0.60, 5.52, 44.08/
C
DATA C /5.,10.,15.,20.,30.,40.,50.,60./
DATA NC /8/
C
CMSC= 1.
C SET CMSC TO 2.54 WHEN USING AN INCH-CALIBRATED PLOTTER
C FOR ADJUSTING THE HEIGHT OF SYMBOLS
C
C SHIFT XD,YD-COORDINATES
DO 5 I=1,ND
YD(I)= YD(I) + 1.
XD(I)= XD(I) + 1.
5 CONTINUE
C INITIALIZE PLOTTING
CALL PLOTS
C GET CPU-TIME
CALL SECOND (T1)
C EXAMPLE 1.
C PLOT CONTOUR LINES (MODE=0)
CALL TRICP (XD,YD,ZD,ND,C,NC,WK,IWK,0)
CALL SECOND (T2)
C MARK DATA POINTS
DO 15 I=1,ND
15 CALL SYMBOL (XD(I),YD(I),0.18/CMSC,1,0.,-1)
C CLOSE PLOT
CALL PLOT (0.,0.,999)
C
C EXAMPLE 2.
C
C SAME TRIANGLES, BUT DIFFERENT ZD-VALUES (MODE=2)
C COMPUTING TIME ABOUT 5*(TIME OF FIRST EXAMPLE)
DO 20 I=1,ND
20 ZD(I)= 5 + MOD(I,2)*10
C DEFINE NEW C-VALUES
NC= 20
DO 30 I=1,NC
30 C(I)= I
CALL PLOTS
CALL SECOND (T3)
C PLOT CONTOUR LINES (MODE=2)
CALL TRICP (XD,YD,ZD,ND,C,NC,WK,IWK,2)
CALL SECOND (T4)
C
C PLOT TRIANGLES
NT= IWK(1)
C = NUMBER OF TRIANGLES
DO 130 I=1,NT
C LOAD VERTICES
DO 120 J=1,3
N= 3*I+2-J
NH= IWK(N)
X(J)= XD(NH)
Y(J)= YD(NH)
120 CONTINUE
C PLOT TRIANGLE SIDES
CALL PLOT (X(1),Y(1),3)
CALL PLOT (X(2),Y(2),2)
CALL PLOT (X(3),Y(3),2)
CALL PLOT (X(1),Y(1),2)
C PLOT TRIANGLE NUMBERS
X0=(X(1) + X(2) + X(3))/3. - 0.3/CMSC
Y0=(Y(1) + Y(2) + Y(3))/3. - 0.15/CMSC
RI= I
CALL NUMBER (X0,Y0,0.3/CMSC,RI,0.,-1)
130 CONTINUE
C ANNOTATE DATA POINTS
DO 200 I=1,ND
AI= I
CALL NUMBER(XD(I)-0.3/CMSC,YD(I)+0.1/CMSC,0.2/CMSC,AI,0.,-1)
CALL NUMBER(XD(I)+0.1/CMSC,YD(I)-0.1/CMSC,0.15/CMSC,ZD(I),0.,1)
200 CONTINUE
C
C CLOSE PLOT
CALL PLOT (0.,0.,999)
C
C WRITE TO OUTPUT CPU-TIME USED FOR THE CALLS TO TRICP
T12= T2-T1
T34= T4-T3
WRITE (*,9001) T12,T34
9001 FORMAT ('0CPU-TIME FOR EX.1',F10.2,/' CPU-TIME FOR EX.2',F10.2)
C
STOP
END
SUBROUTINE TRICP (XD,YD,ZD,ND,C,NC,WK,IWK,MODE)
C
C T R I ANGLE C ONTOUR P LOTTING
C
C USER INTERFACE (VERSION 0.1)
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
DIMENSION XD(ND),YD(ND),ZD(ND),C(NC),WK(*),IWK(*)
C
C XD,YD,ZD COORDINATES OF DATA POINTS
C XD AND YD MUST BE GIVEN IN CM.
C FOR USING INCH, SET INSTALLATION PARAMETER
C CMSCAL= 2.54 .
C ND NUMBER OF DATA POINTS
C C CONTOUR LEVELS
C NC NUMBER OF CONTOUR LEVELS
C WK REAL WORK ARRAY, DIMENSION 5*ND
C
C ON EXIT, WK CONTAINS THE PARTIAL DERIVATIVES
C (ZX(I),ZY(I),ZXX(I),ZXY(I),ZYY(I),I=1,ND)
C
C IWK INTEGER WORK ARRAY. FOR MODE=0, THE DIMENSION IS
C 31*ND .
C FOR OTHER MODES, SEE BELOW.
C
C ON EXIT, IWK CONTAINS
C A) IWK(1)= NT= NUMBER OF TRIANGLES
C B) IWK(2...3*NT+1) POINT NUMBERS FOR THE TRIANGLES.
C THE FIRST 3 NUMBERS DETERMINE THE VERTICES OF THE
C FIRST TRIANGLE, THE NEXT 3 FOR THE SECOND AND SO ON.
C THE NUMBERS COORESPOND TO THE INDICES OF THE
C XD,YD,ZD ARRAYS. THEY ARE ARRANGED COUNTER-
C CLOCKWISE WITHIN A TRIANGLE.
C
C C) IWK(3*NT+2 ... 3*NT+NCP*ND+1)
C A LIST OF NCP*ND POINT NUMBERS, REPRESENTING
C NCP POINTS FOR EACH DATA POINT THAT ARE
C THE CLOSEST TO THIS POINT. THE FIRST NCP
C NUMBERS ARE FOR THE FIRST DATA POINT, THE NEXT
C NCP FOR THE SECOND AND SO ON. THESE NUMBERS
C WERE USED FOR THE COMPUTATION OF THE PARTIAL
C DERIVATIVES.
C NCP IS AN INSTALLATION PARAMETER AND WILL BE SET
C TO 4.
C
C MODE MODE OF USAGE
C 0, NORMAL MODE, SEE ABOVE.
C 1, NO TRIANGULATION REQUESTED.
C IWK MUST CONTAIN THE INFORMATION ABOUT
C THE TRIANGLES AS DESCRIBED UNDER A) AND B)
C ON ENTRY.
C THESE LOCATIONS OF IWK WILL NOT BE CHANGED
C AND IN ADDITION, THE INFORMATION DESCRIBED
C UNDER C) WILL BE AVAILABLE ON EXIT.
C 2, NO TRIANGULATION AND NO DETERMINATION OF THE
C NCP CLOSEST POINTS FOR EACH DATA POINT.
C IWK MUST CONTAIN THE INFORMATION AS DESCRIBED
C UNDER A), B), AND C) AS INPUT INFORMATION.
C THE CONTENTS OF IWK WILL NOT BE CHANGED.
C 3, NO TRIANGULATION AND NO COMPUTATION OF THE
C PARTIAL DERIVATIVES.
C IWK MUST CONTAIN THE INFORMATION A) AND B)
C AND
C WK MUST CONTAIN THE PARTIAL DERIVATIVES
C (ZX(I),ZY(I),ZXX(I),ZXY(I),ZYY(I),I=1,ND)
C ON ENTRY.
C THE CONTENTS OF WK AND IWK WILL NOT BE CHANGED.
C THIS MODE IS ESPECIALLY USEFUL WHEN TRICP IS
C CALLED AGAIN AFTER A PREVIOUS CALL. FOR INSTANCE,
C ONLY THE CONTOUR LEVELS MAY HAVE CHANGED.
C WHEN DESIGNING A SURFACE, IT MAY BE APPROPRIATE
C TO CHANGE THE XD,YD,ZD PARAMETERS AND THE PARTIAL
C DERIVATIVES INTERACTIVELY AND TO CALL TRICP AGAIN
C USING THIS MODE.
C
C FOR MODES 1,2, AND 3 THE DIMENSION OF IWK CAN
C BE REDUCED TO
C 3*NT + NCP*ND + 1.
C
C
COMMON /TRPCOC/ DSPMIN,CMSCAL,NPMAX
C
C *** DSPMIN RESOLUTION OF PLOTTING DEVICE IN USE
C (MINIMAL DISTANCE OF TWO POINTS TO BE PLOTTED)
C *** CMSCAL VARIABLE FOR SWITCHING BETWEEN CM AND INCH
C NPMAX MAXIMUM NUMBER OF POINTS PER LINE FOR A TRIANGLE
C
C
C SET INSTALLATION PARAMETERS
NCP= 4
CMSCAL= 1.
C *** SET CMSCAL= 2.54 FOR INCH CALIBRATED PLOTTERS
DSPMIN= 0.02/CMSCAL
NPMAX= 500
C
C SET STARTING ADDRESSES FOR IWK WORK ARRAY
N1= 2
N2= N1 + 6*ND - 15
IF (MODE.GE.1) N2= N1 + 3*IWK(1)
N3= N2 + 6*ND
N4= N3 + 18*ND
C
IF (MODE.EQ.1) GOTO 50
IF (MODE.EQ.2) GOTO 60
IF (MODE.EQ.3) GOTO 70
C
C TRIANGULATION
CALL IDTANG (ND,XD,YD,IWK,IWK(N1),NL,IWK(N2),IWK(N3),IWK(N4),WK)
C IDTANG IS PART OF ACM ALG. 526 BY H. AKIMA
N2= N1+3*IWK(1)
C
C FIND NCP POINTS CLOSEST TO EACH DATA POINT
50 CONTINUE
CALL IDCLDP (ND,XD,YD,NCP,IWK(N2))
C IDCLDP IS PART OF ACM ALG. 526 BY H. AKIMA
C
C COMPUTE PARTIAL DERIVATIVES
60 CONTINUE
CALL IDPDRV (ND,XD,YD,ZD,NCP,IWK(N2),WK)
C IDPDRV IS PART OF ACM ALG. 526 BY H. AKIMA
C
C COMPUTE CONTOURS BY SUCCESSIVE SOLUTION OF
C QUINTIC POLYNOMIAL EQUATIONS
70 CONTINUE
CALL TRP00 (XD,YD,ZD,C,NC,IWK(N1),WK,IWK(1))
C
RETURN
END
SUBROUTINE TRP00 (XD,YD,ZD,CN,NC,IPT,PDD,NT)
C
C T R I ANGLE C ONTOUR P LOTTING
C
C MAIN ROUTINE
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
DIMENSION XD(*),YD(*),ZD(*),CN(NC),IPT(*),PDD(*)
C
C XD,YD,ZD COORDINATES OF DATA POINTS
C CN CONTOUR LEVELS
C NC NUMBER OF CONTOUR LEVELS
C IPT POINT NUMBERS STORED AS 3*I-2, 3*I-1, 3*I TH
C ELEMENT FOR TRIANGLE I
C PDD PARTIAL DERIVATIVES
C (ZX,ZY,ZXX,ZXY,ZYY)
C
C
C FOR PLOTTING THE CALCOMP ROUTINE PLOT(X,Y,IPEN) IS USED
C IPEN= 3 MEANS MOVE TO X,Y WITH PEN IN UPWARD POSITION
C IPEN= 2 MEANS MOVE TO X,Y WITH PEN DOWN
C
C
COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
1, SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
2, AP,BP,CP,DP
C
C X,Y,Z COORDINATES AT THE THREE VERTICES OF A TRIANGLE
C DX,DY DIFFERENCES OF X,Y
C DX2,DY2 (DIFFERENCES OF X,Y)**2
C SL SIDE LENGTHS
C SL2 (SIDE LENGTHS)**2
C ZT(IV,K) FIRST PARTIAL DERIVATIVE WITH RESPECT TO VARIABLE IV
C (U,V,W FOR IV=1,2,3) AT POINT K
C ZTT SECOND PARTIAL DERIVATIVES
C ZUV MIXED PARTIAL DERIVATIVES
C AP,BP,CP,DP CONSTANTS DERIVED FROM DX,DY
C NAMES DUE TO ALG. 526 CACM
C
COMMON /TRPCOC/ DSPMIN,CMSCAL,NPMAX
C FOR EXPLANATION OF VARIABLES IN /TRPCOC/ SEE SUBROUTINE TRICP
C
COMMON /TRPCOF/ KK,KSE,XX4F,YY4F,SIR,COR,CL
C
C /TRPCOF/ CONTAINS VARIABLES WHICH ARE PASSED TO FUNCTION
C TRICPF AS PARAMETERS
C
C KK NUMBER OF FUNCTION TO BE EVALUATED BY TRICPF
C KSE ACTUAL SIDE NUMBER
C XX4F,YY4F COORDINATES FOR POINT P4F (PRELEMINARY POSITION
C OF POINT P4)
C SIR,COR COSINUS OF DIRECTION NORMAL TO CURVE DIRECTION
C CL ACTUAL SCALED CONTOUR LEVEL
C
DIMENSION TS1(6,3),Z1(6,3),ZMAX(3),ZMIN(3),TS2(6),IN(3)
1, TZR(5,3),X0(5,3),Y0(5,3),NZ(3),TPER(3),DT(3),THETAS(3)
2, SI(3),CO(3)
C
C TS1 U,V,W-VALUES AT ENDPOINTS OF INTERVALS
C Z1 Z-VALUES AT ENDPOINTS OF INTERVALS
C ZMAX,ZMIN MAXIMUM AND MINIMUM Z VALUE ALONG A SIDE
C TS2 WORKING ARRAY FOR COMPUTING TS1
C IN NUMBER OF INTERVALS
C TZR U,V,W VALUES FOR ZEROS FOUND ON SIDES
C X0,Y0 X,Y COORDINATES FOR ZEROS FOUND
C NZ NUMBER OF ZEROS FOUND
C TPER PERMITTED POSITION ERROR FOR ZEROS ON SIDES
C DT DIFFERENCE IN U,V,W FOR ESTIMATING STARTING
C DIRECTION OF A CONTOUR LINE
C THETAS ANGLE OF DIRECTION FOR SIDES
C SI,CO COSINUS OF DIRECTION FOR SIDES
C
PI= 4.*ATAN2(1.,1.)
C
C START MAIN LOOP OVER TRIANGLES
DO 5000 IT=1,NT
C
C LOAD COORDINATES AT VERTICES
JI= 3*(IT-1)
DO 50 J=1,3
JI= JI + 1
ID= IPT(JI)
X(J)= XD(ID)
Y(J)= YD(ID)
Z(J)= ZD(ID)
50 CONTINUE
ZS= (Z(1)+Z(2)+Z(3))/3.
C
C SOME BASIC GEOMETRY FOR THE TRIANGLE
DO 60 J=1,3
Z(J)= Z(J) - ZS
NP1= 3
NP2= 2
IF (J.EQ.3) NP1= 1
IF (J.EQ.2) NP2= 1
DX(J)= X(NP2) - X(NP1)
DY(J)= Y(NP2) - Y(NP1)
DX2(J)= DX(J)*DX(J)
DY2(J)= DY(J)*DY(J)
SL2(J)= DX2(J) + DY2(J)
SL(J)= SQRT(SL2(J))
CO(J)= DX(J)/SL(J)
SI(J)= DY(J)/SL(J)
60 CONTINUE
CO(1)= -CO(1)
SI(1)= -SI(1)
SLMAX= AMAX1(SL(1),SL(2),SL(3))
DI1= ABS(DY(3)*CO(1) - DX(3)*SI(1))
DI2= ABS(DY(1)*CO(2) - DX(1)*SI(2))
DI3= ABS(DY(2)*CO(3) - DX(2)*SI(3))
HMIN= AMIN1(DI1,DI2,DI3)
C = SHORTEST DISTANCE BETWEEN A VERTEX AND
C ITS OPPOSITE SIDE (MINIMUM HEIGHT)
C
C *** ISSUE A WARNING MESSAGE HERE, IF HMIN IS LESS THAN,
C SAY 0.01/CMSCAL
C
C DEFINE CONSTANTS DEPENDING ON HMIN
C
RMAX = AMIN1(0.02/CMSCAL,HMIN*0.01)
C = DISTANCE NORMAL TO CURVE DIRECTION WITHIN WHICH
C A ZERO MUST BE FOUND
RMIN = RMAX*0.2
C IF ZERO HAS BEEN FOUND WITHIN A DISTANCE SMALLER
C THAN RMIN, THEN STEPSIZE DS IS MULTIPLIED BY 2.
POSERR= AMIN1(1.E-3/CMSCAL,1.E-3*HMIN*HMIN/SLMAX)
C = PERMITTED POSITION ERROR
DSMAX = HMIN*0.2
C = MAXIMUM STEP SIZE
FSTEP = RMAX*10.
C = STARTING STEP SIZE
EPS = HMIN*0.01
C = DIFFERENCE FOR ESTIMATING STARTING DIRECTION
RMAX2 = RMAX*2.
C = DISTANCE NORMAL TO CURVE DIRECTION. A ZERO MUST BE
C FOUND WITHIN WHEN CROSSING TRIANGLE BORDER.
C
C COMPUTE COEFFICIENTS FOR POLYNOMIALS ALONG SIDES
CALL TRP001 (IT,IPT,PDD)
C
CL= 0.
C
C LOOP OVER SIDES
DO 200 JSE=1,3
KSE= JSE
TPER(KSE)= POSERR/SL(KSE)
C
C COMPUTE ENDPOINTS OF INTERVALS
TS1(1,KSE)= 0.
TS1(2,KSE)= 1.
I= 2
DO 150 K=1,4
KK= K
TS2(1)= 0.
II= 2
TB= 0.
F2= TRICPF(TB)
DO 100 J=2,I
TA= TB
F1= F2
TB= TS1(J,KSE)
F2= TRICPF(TB)
IF (F1*F2.GT.0.) GOTO 100
IF (F1.EQ.0. .AND. F2.EQ.0.) GOTO 100
TS2(II)= TRICPZ(TA,TB,F1,F2,TPER(KSE))
II= II + 1
100 CONTINUE
TS2(II)= 1.
DO 120 J=1,II
120 TS1(J,KSE)= TS2(J)
I= II
150 CONTINUE
C IN(KSE)= NUMBER OF INTERVALS
I= I-1
IN(KSE)= I
C (E.G. IF IN(KSE)=1, THERE IS NO POINT FOR WHICH 1ST DER.=0)
C
C COMPUTE MAXIMA AND MINIMA FOR EACH SIDE
NP1= 3
NP2= 2
IF (KSE.EQ.3) NP1=1
IF (KSE.EQ.2) NP2=1
ZMAX(KSE)= AMAX1(Z(NP1),Z(NP2))
ZMIN(KSE)= AMIN1(Z(NP1),Z(NP2))
Z1(1,KSE)= Z(NP1)
Z1(I+1,KSE)= Z(NP2)
IF (I.EQ.1) GOTO 170
KK= 5
DO 160 J=2,I
Z1(J,KSE)= TRICPF(TS1(J,KSE))
IF (Z1(J,KSE).GT.ZMAX(KSE)) ZMAX(KSE)= Z1(J,KSE)
IF (Z1(J,KSE).LT.ZMIN(KSE)) ZMIN(KSE)= Z1(J,KSE)
160 CONTINUE
170 CONTINUE
200 CONTINUE
C
NPTRI= 0
KA= 1
KE= 1
KC= 0
C
C START LOOP OVER CONTOUR LEVELS
DO 4000 K=1,NC
C
CL= CN(K) - ZS
C
C COMPUTE ZEROS (IF ANY) ON THE THREE SIDES FOR LEVEL CL
KK= 5
DO 1200 JSE=1,3
NZ(JSE)= 0
IF (CL.LT.ZMIN(JSE) .OR. CL.GT.ZMAX(JSE)) GOTO 1200
KSE= JSE
JN=0
NI= IN(KSE)
C CHECK INTERVALS FOR ZEROS
DO 1100 J= 1,NI
F1= Z1(J,KSE) - CL
F2= Z1(J+1,KSE) - CL
F1F2= F1*F2
IF (F1F2.GT.0.) GOTO 1100
IF (F1F2.LT.0.) GOTO 1090
C
C SPECIAL SITUATIONS
IF (NI.EQ.1 .AND. F1.EQ.0. .AND. F2.EQ.0.) GOTO 1070
IF (J.EQ.1 .AND. F1.EQ.0. .OR.
A J.EQ.NI .AND. F2.EQ.0.) GOTO 1080
GOTO 1090
1070 CONTINUE
C CONTOURLINE = SIDE KSE
NP1= 3
NP2= 2
IF (KSE.EQ.3) NP1= 1
IF (KSE.EQ.2) NP2= 1
CALL PLOT (X(NP1),Y(NP1),3)
CALL PLOT (X(NP2),Y(NP2),2)
GOTO 1100
1080 CONTINUE
C LINE PASSES THROUGH DATA POINT
C ONLY ONE ZERO FOR THIS TRIANGLE, SKIP EITHER SIDE
IF (KSE.EQ.3 .OR. (F1.EQ.0. .AND. KSE.EQ.2)) GOTO 1100
JN= JN+1
IF (F2.EQ.0.) TZR(JN,KSE)= 1.-TPER(KSE)*0.5
IF (F1.EQ.0.) TZR(JN,KSE)= TPER(KSE)*0.5
GOTO 1100
C
1090 JN= JN+1
TZR(JN,KSE)= TRICPZ(TS1(J,KSE),TS1(J+1,KSE),F1,F2,TPER(KSE))
1100 CONTINUE
NZ(KSE)= JN
C = NUMBER OF ZEROS FOR SIDE KSE AND CONTOUR LEVEL CL
1200 CONTINUE
C
IF (NZ(1)+NZ(2)+NZ(3).LT.2) GOTO 4000
C IF (.TRUE.) THEN STOP PROCESSING THIS CONTOUR LEVEL
C
C COMPUTE X0,Y0 FOR EACH ZERO (RELATIVE TO X(3),Y(3))
DO 1300 JSE=1,3
JN= NZ(JSE)
IF (JN.EQ.0) GOTO 1300
DO 1280 JZR=1,JN
T= TZR(JZR,JSE)
IF (JSE.EQ.3) GOTO 1230
X0(JZR,JSE)= DX(JSE)*T
Y0(JZR,JSE)= DY(JSE)*T
GOTO 1280
1230 X0(JZR,JSE)= DX(2) + DX(3)*T
Y0(JZR,JSE)= DY(2) + DY(3)*T
1280 CONTINUE
1300 CONTINUE
C
KC= KC+1
IF (KC.NE.1) GOTO 1350
C COMPUTE DIRECTION OF SIDES
THETAS(1)= ATAN2(-DY(1),-DX(1))
THETAS(2)= ATAN2( DY(2), DX(2))
THETAS(3)= ATAN2( DY(3), DX(3))
C COMPUTE DIFFERENCES FOR ESTIMATING START DIRECTION
DT(1)= -EPS/SL(1)
DT(2)= EPS/SL(2)
DT(3)= EPS/SL(3)
C COMPUTE COEFFICIENTS FOR POLYNOMIAL INSIDE TRIANGLE
CALL TRP002
1350 CONTINUE
C
C FOLLOW CONTOURS
C
RMA= RMAX
C LOOP OVER SIDES, KSE= SIDE NUMBER
DO 3000 JSE=1,3
KSE= KA
JN= NZ(KSE)
IF (JN.EQ.0) GOTO 2010
RMA= SIGN (RMAX,-RMA)
C LOOP OVER ZEROS
DO 2000 JZR= 1,JN
IF (TZR(JZR,KSE).GT.5.) GOTO 2000
C IF () THEN THIS ZERO HAS ALREADY BEEN PROCESSED
C
C ESTIMATE STARTING DIRECTION
C
C COMPUTE F1 ON SIDE USING FIRST DERIVATIVE FD1
KK=4
FD1= TRICPF(TZR(JZR,KSE))
F1= FD1 * DT(KSE)
C
C COMPUTE F2 NORMAL TO SIDE
KK= 6
XX4F= X0(JZR,KSE)
YY4F= Y0(JZR,KSE)
SIR= CO(KSE)
COR= - SI(KSE)
F2= TRICPF(EPS)
C
C COMPUTE ANGLE
IF (F2.EQ.0.) GOTO 1470
THETSC= ATAN2 (-F1,F2)
IF (THETSC.LE.0.) THETSC= THETSC + PI
GOTO 1480
1470 CONTINUE
IF (FD1.EQ.0.) GOTO 2000
THETSC= PI*0.5
1480 THETAC= THETAS(KSE) + THETSC
C
C MOVE PEN TO START
CALL PLOT (X0(JZR,KSE)+X(3),Y0(JZR,KSE)+Y(3),3)
C
C INITIALIZE TRACING
DS= FSTEP
DX12= DS*COS(THETAC)
DY12= DS*SIN(THETAC)
XX3= XX4F
YY3= YY4F
XX2= XX3 - DX12
YY2= YY3 - DY12
XX1= XX2 - DX12
YY1= YY2 - DY12
DS23= DS
DS12= DS
C POINTS P1,P2,P3,P4 WITH COORDINATES XX1...XX4, YY1...YY4
C ARE REFERRED TO AS *QUEUE*. DS12...DS34 ARE THE DISTANCES
C BETWEEN POINTS IN THE QUEUE.
C EVERY POINT NORMALLY MUST PASS THROUGH THE QUEUE
C BEFORE IT IS PLOTTED.
C FIRST (OLDEST) POINT IN THE QUEUE IS P1, WHICH IS NORMALLY
C (NPQ = 4) THE FIRST CANDIDATE TO BE PLOTTED.
C P4 IS THE NEXT POINT TO BE COMPUTED.
NPQ= 0
C NPQ= NUMBER OF POINTS IN THE QUEUE COUNTED FROM THE END
C WHICH CAN BE PLOTTED.
C E.G. NPQ=1, P4 IS STILL TO BE PLOTTED
NP=0
C = NUMBER OF POINTS COMPUTED FOR THIS CONTOUR LINE
DSP=0.
C = DISTANCE OF POINT TO BE PLOTTED TO LAST POINT PLOTTED
NSTOP= 0
NFOUND= 0
NOST= 0
C
C COMPUTE NEW POINT FOR CONTOUR LINE
C
C COMPUTE CURVE DIRECTION
1500 CONTINUE
DS13= DS23 + DS12
PL0= DS23/(DS12*DS13)
PL1= -DS13/(DS12*DS23)
PL2= (DS13+DS23)/(DS13*DS23)
DXDS= PL0*XX1 + PL1*XX2 + PL2*XX3
DYDS= PL0*YY1 + PL1*YY2 + PL2*YY3
SQ= SQRT(DXDS*DXDS+DYDS*DYDS)
DXDS= DXDS/SQ
DYDS= DYDS/SQ
COR= -DYDS
SIR= DXDS
C
C SEARCH FOR TWO POINTS WITH OPPOSITE SIGN
RMA= SIGN (RMAX,RMA)
1550 CONTINUE
XX4F= XX3 + DXDS*DS
YY4F= YY3 + DYDS*DS
F1= TRICPF(0.)
F2= TRICPF(RMA)
IF (F1*F2 .LE. 0.) GOTO 1600
IF (ABS(F2) .LT. ABS(F1)) GOTO 1560
RMA= -RMA
F2= TRICPF(RMA)
IF (F1*F2 .LE. 0.) GOTO 1600
C
C DIVIDE STEPSIZE IN CURVE DIRECTION BY 2.
1560 CONTINUE
DS= DS*0.5
IF (DS.GT.POSERR) GOTO 1550
C
C DIVIDE STEPSIZE NORMAL TO CURVE BY 2.
NOST= NOST + 1
DS= DS*2.
RMA= RMA*0.5
1570 F2= TRICPF(RMA)
IF (F1*F2.LE.0.) GOTO 1600
RMA= -RMA
F2= TRICPF(RMA)
IF (F1*F2.LE.0.) GOTO 1600
RMA= RMA*0.5
IF (ABS(RMA).GT.POSERR) GOTO 1570
NSTOP= 1
GOTO 1900
C
C
C FIND ZERO FOR NEW POINT
1600 NPQ= NPQ + 1
NP= NP + 1
IF (NP.GT.NPMAX) GOTO 1990
R= TRICPZ(0.,RMA,F1,F2,POSERR)
DS34= SQRT(DS*DS + R*R)
XX4= XX4F + COR*R
YY4= YY4F + SIR*R
U= AP*XX4 + BP*YY4
V= CP*XX4 + DP*YY4
C
C CHECK IF POINT IS OUTSIDE THE TRIANGLE
KE=1
IF (U.LT.0.) GOTO 1700
KE=2
IF (V.LT.0.) GOTO 1700
KE= 3
IF (V .GT. 1.-U) GOTO 1700
C
C POINT IS INSIDE
C PLOT LAST POINT OF QUEUE AND UPDATE QUEUE
IF (NPQ.NE.4) GOTO 1610
C *** INSERT CALL FOR LABELING ROUTINE HERE
C *** USING XX1...XX4, YY1...YY4, DS12...DS34 AS PARAMETERS
C *** AND REDUCE NPQ BY THE NUMBER OF POINTS WHICH HAVE BEEN
C *** USED FROM THE QUEUE FOR LABELING
NPQ= NPQ -1
DSP= DSP + DS01
IF (DSP.LT.DSPMIN) GOTO 1610
CALL PLOT (XX1+X(3),YY1+Y(3),2)
DSP= 0.
1610 DS01= DS12
DS12= DS23
DS23= DS34
XX1= XX2
YY1= YY2
XX2= XX3
YY2= YY3
XX3= XX4
YY3= YY4
C SET NEW STEP SIZE
IF (ABS(R).LT.RMIN) DS= AMIN1(DSMAX,DS*2.)
GOTO 1500
C
C POINT IS OUTSIDE
C
1700 CONTINUE
C NO SEARCH IF FIRST POINT WAS COMPUTED AND
C START WAS AT DATA POINT
IF (NP.NE.1) GOTO 1710
IF (TZR(JZR,KSE).GT.1.-TPER(KSE) .OR.
A TZR(JZR,KSE).LT.TPER(KSE)) GOTO 2000
1710 CONTINUE
C
C FLAG ZERO WHERE LINE STARTED
TZR(JZR,KSE)= TZR(JZR,KSE) + 10.
C
C SEARCH FOR CORRESPONDING ZERO ON SIDE KE
CHECK= 5.
DIST= 99999.
DO 1790 N=1,2
DO 1770 JESE=1,3
JJ= NZ(KE)
IF (JJ.EQ.0) GOTO 1760
DO 1750 J=1,JJ
IF (TZR(J,KE).GT.CHECK) GOTO 1750
DI1= ABS((Y0(J,KE)-YY3)*DXDS - (X0(J,KE)-XX3)*DYDS)
IF (DI1.GE.DIST) GOTO 1750
J1= J
DIST= DI1
1750 CONTINUE
1760 IF (DIST.LT.RMAX2) GOTO 1800
KE= MOD(KE,3) + 1
1770 CONTINUE
C FOR N=2, DISREGARD FLAGS
CHECK= 15.
1790 CONTINUE
C NO ACCEPTABLE ZERO ON ALL THREE SIDES
NFOUND= 1
GOTO 1900
C
C FLAG ZERO FOUND
1800 TZR(J1,KE)= TZR(J1,KE) + 10.
C
C CLEAR QUEUE AND FINISH CONTOUR LINE
1900 CONTINUE
IF (NPQ.GT.3) GOTO 1910
IF (NPQ.GT.2) GOTO 1920
IF (NPQ.GT.1) GOTO 1930
GOTO 1940
1910 DSP= DSP + DS01
IF (DSP.LT.DSPMIN) GOTO 1920
CALL PLOT(XX1+X(3),YY1+Y(3),2)
DSP= 0.
1920 DSP= DSP + DS12
IF (DSP.LT.DSPMIN) GOTO 1930
CALL PLOT(XX2+X(3),YY2+Y(3),2)
DSP= 0.
1930 DSP= DSP + DS23
IF (DSP.LT.DSPMIN) GOTO 1940
CALL PLOT(XX3+X(3),YY3+Y(3),2)
1940 IF (NFOUND+ NSTOP .EQ.0)
A CALL PLOT (X0(J1,KE)+X(3),Y0(J1,KE)+Y(3),2)
C
1990 CONTINUE
NPTRI= NPTRI + NP
2000 CONTINUE
C END OF LOOP OVER ZEROS
2010 CONTINUE
KA= MOD(KA,3) + 1
3000 CONTINUE
C END OF LOOP OVER SIDES
KA= KE
4000 CONTINUE
C END OF LOOP OVER CONTOUR LEVELS
5000 CONTINUE
C END OF LOOP OVER TRIANGLES
RETURN
END
SUBROUTINE TRP001 (IT0,IPT,PDD)
C
C T R I ANGLE C ONTOUR P LOTTING
C
C COMPUTATION OF COEFFICIENTS FOR POLYNOMIALS ALONG SIDES
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
DIMENSION IPT(*),PDD(*)
C
C IT0 NUMBER OF TRIANGLE IN USE
C IPT POINT NUMBERS STORED AS 3*I-2, 3*I-1, 3*I TH
C ELEMENT FOR TRIANGLE I
C PDD PARTIAL DERIVATIVES AT DATA POINTS
C (ZX,ZY,ZXX,ZXY,ZYY)
C
COMMON /TRPCOP/ P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
1, Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
2, R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
3, P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C P0...P5 COEFFICIENTS OF POLYNOMIALS ALONG THE SIDES
C Q0...Q4 COEFFICIENTS OF 1ST DERIVATIVES ALONG THE SIDES
C R0...R3 COEFFICIENTS OF 2ND DERIVATIVES ALONG THE SIDES
C S0...S2 COEFFICIENTS OF 3RD DERIVATIVES ALONG THE SIDES
C T0...T1 COEFFICIENTS OF 4TH DERIVATIVES ALONG THE SIDES
C P11...P41 COEFFICIENTS FOR BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
1, SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
2, AP,BP,CP,DP
C FOR EXPLANATION OF VARIABLES IN /TRPCOT/ SEE SUBROUTINE TRP00
C
DIMENSION PD(5,3)
C PARTIAL DERIVATIVES IN X-Y DIRECTIONS AT THE THREE VERTICES
C
C LOADS PARTIAL DERIVATIVES AT THE VERTICES
JIPT=3*(IT0-1)
DO 230 K=1,3
JIPT=JIPT+1
IDP=IPT(JIPT)
JPDD= 5*(IDP-1)
DO 220 KPD=1,5
JPDD=JPDD+1
PD(KPD,K)=PDD(JPDD)
220 CONTINUE
230 CONTINUE
C DETERMINES THE COEFFICIENTS FOR THE COORDINATE SYSTEM
C TRANSFORMATION FROM THE X-Y SYSTEM TO THE U-V SYSTEM
AD= DX(2)*DY(1)
BC= DX(1)*DY(2)
DLT= AD-BC
AP= DY(1)/DLT
BP= -DX(1)/DLT
CP= -DY(2)/DLT
DP= DX(2)/DLT
C CONVERTS THE PARTIAL DERIVATIVES AT THE VERTEXES OF THE
C TRIANGLE FOR THE U-V COORDINATE SYSTEM.
AB= DX(2)*DX(1)
ADBC= AD+BC
CD= DY(1)*DY(2)
DXDY1= 2.0*DX(1)*DY(1)
DXDY2= 2.0*DX(2)*DY(2)
DXDY3= 2.0*DX(3)*DY(3)
DO 260 K=1,3
ZT(1,K) = DX(2)*PD(1,K) + DY(2)*PD(2,K)
ZT(2,K) = DX(1)*PD(1,K) + DY(1)*PD(2,K)
ZTT(1,K)= DX2(2)*PD(3,K) + DXDY2*PD(4,K) + DY2(2)*PD(5,K)
ZUV(K) = AB *PD(3,K) + ADBC*PD(4,K) + CD *PD(5,K)
ZTT(2,K)= DX2(1)*PD(3,K) + DXDY1*PD(4,K) + DY2(1)*PD(5,K)
260 CONTINUE
DO 270 K=1,2
ZT(3,K) = DX(3)*PD(1,K) + DY(3)*PD(2,K)
ZTT(3,K)= DX2(3)*PD(3,K) + DXDY3*PD(4,K) + DY2(3)*PD(5,K)
270 CONTINUE
C CALCULATES THE COEFFICIENTS OF THE POLYNOMIALS ALONG
C THE THREE SIDES OF THE TRIANGLE
DO 280 JSE=1,3
NP1= 3
NP2= 2
IF (JSE.EQ.3) NP1= 1
IF (JSE.EQ.2) NP2= 1
IV= NP2
IF (JSE.EQ.3) IV= 3
P0(JSE)= Z(NP1)
P1(JSE)= ZT(IV,NP1)
P2(JSE)= 0.5*ZTT(IV,NP1)
H1= Z(NP2) - P0(JSE) - P1(JSE) - P2(JSE)
H2= ZT(IV,NP2) - P1(JSE) - ZTT(IV,NP1)
H3= ZTT(IV,NP2) - ZTT(IV,NP1)
P3(JSE)= 10.0*H1-4.0*H2+0.5*H3
P4(JSE)=-15.0*H1+7.0*H2 -H3
P5(JSE)= 6.0*H1-3.0*H2+0.5*H3
C CALCULATES COEFFICIENTS FOR DERIVATIVES ALONG SIDES
Q0(JSE)= P1(JSE)
Q1(JSE)= 2. *P2(JSE)
Q2(JSE)= 3. *P3(JSE)
Q3(JSE)= 4. *P4(JSE)
Q4(JSE)= 5. *P5(JSE)
R0(JSE)= Q1(JSE)
R1(JSE)= 2. *Q2(JSE)
R2(JSE)= 3. *Q3(JSE)
R3(JSE)= 4. *Q4(JSE)
S0(JSE)= R1(JSE)
S1(JSE)= 2. *R2(JSE)
S2(JSE)= 3. *R3(JSE)
T0(JSE)= S1(JSE)
T1(JSE)= 2. *S2(JSE)
280 CONTINUE
RETURN
END
SUBROUTINE TRP002
C
C T R I ANGLE C ONTOUR P LOTTING
C
C COMPUTATION OF POLYNOMIAL COEFFICIENTS P11...P41 FOR
C BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
COMMON /TRPCOP/ P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
1, Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
2, R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
3, P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C FOR EXPLANATION OF VARIABLES IN /TRPCOP/ SEE SUBROUTINE TRP001
C
COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
1, SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
2, AP,BP,CP,DP
C FOR EXPLANATION OF VARIABLES IN /TRPCOT/ SEE SUBROUTINE TRP00
C
C
P14= P5(1) * (2.5*(SL2(2)-SL2(3))/SL2(1) + 2.5)
P41= P5(2) * (2.5*(SL2(1)-SL2(3))/SL2(2) + 2.5)
P11= ZUV(3)
H1= ZT(2,1)-P1(1)-P11-P41
H2= ZUV(1)-P11-4.0*P41
P21= 3.0*H1-H2
P31=-2.0*H1+H2
H1= ZT(1,2)-P1(2)-P11-P14
H2= ZUV(2)-P11-4.0*P14
P12= 3.0*H1-H2
P13=-2.0*H1+H2
H1= 0.5*ZTT(2,1)-P2(1)-P12
H2= 0.5*ZTT(1,2)-P2(2)-P21
E1= 2.5*(SL2(2)-SL2(1))/SL2(3) + 2.5
G1= 3. - E1
G2=-2. + E1
G3= E1*(P5(1) - P5(2) + P41 -P14)
1 + P14 - 4.*P41 + 5.*P5(2)
P22= G1*H1 + G2*H2 + G3
P32= H1-P22
P23= H2-P22
RETURN
END
FUNCTION TRICPZ (TA,TB,F1,F2,ER)
C
C T R I ANGLE C ONTOUR P LOTTING
C
C COMPUTE ZERO BETWEEN TA AND TB
C
C F1= FUNCTION VALUE AT TA
C F2= FUNCTION VALUE AT TB
C F1 AND F2 MUST HAVE OPPOSITE SIGN
C THIS MUST BE CHECKED BEFORE ENTRY
C ER= PERMITTED ERROR FOR SOLUTION TRICPZ
C NAME OF FUNCTION = TRICPF
C
C THE METHOD IS A COMBINATION OF THE REGULA FALSI
C AND THE MIDPOINT METHOD
C
C IT IS A MODIFIED VERSION OF THE VIM- (CONTROL DATA
C USER GROUP) ROUTINE WITH CATALOG IDENTIFICATION
C C2BKYZERO
C WRITTEN BY LOREN P. MEISSNER, 1965
C
A=TA
B=TB
FA=F1
FB=F2
C=A
FC=FA
S=C
FS=FC
C
10 CONTINUE
H=0.5*(B+C)
IF(ABS(H-B) .LE.ER) GO TO 110
IF (ABS(FB) .LE. ABS(FC)) GO TO 15
Y=B
FY=FB
G=B
FG=FB
S=C
FS=FC
GO TO 20
15 Y=S
FY=FS
G=C
FG=FC
S=B
FS=FB
20 CONTINUE
IF (FY .NE. FS) GO TO 21
B=H
GO TO 29
21 CONTINUE
E=(S*FY-Y*FS)/(FY-FS)
IF (ABS(E-S) .LE.ER) E=S+SIGN(ER,G-S)
IF ((E-H)*(S-E) .LT. 0.0) GO TO 28
B=E
GO TO 29
28 B=H
C
C *** FUNCTION CALL
29 FB=TRICPF(B)
C
IF (FG*FB .LT. 0.0) GO TO 35
C=S
FC=FS
GO TO 10
35 CONTINUE
C=G
FC=FG
GO TO 10
C
110 TRICPZ= H
C
RETURN
END
FUNCTION TRICPF(T)
C
C T R I ANGLE C ONTOUR P LOTTING
C
C FUNCTION EVALUATION
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
COMMON /TRPCOF/ KK,KSE,XX4F,YY4F,SIR,COR,CL
C
C VARIABELS IN /TRPCOF/ ARE USED AS ARGUMENTS
C FOR AN EXPLANATION SEE SUBROUTINE TRP00
C
C KK NUMBER OF FUNCTION TO BE EVALUATED
C KK=1 4TH DERIVATIVE ALONG SIDE KSE
C KK=2 3RD DERIVATIVE ALONG SIDE KSE
C KK=3 2ND DERIVATIVE ALONG SIDE KSE
C KK=4 1ST DERIVATIVE ALONG SIDE KSE
C KK=5 ORIGINAL POLYNOMIAL ALONG SIDE KSE
C KK=6 BIVARIATE POLYNOMIAL INSIDE TRIANGLE
C
COMMON /TRPCOP/ P0(3),P1(3),P2(3),P3(3),P4(3),P5(3)
1, Q0(3),Q1(3),Q2(3),Q3(3),Q4(3)
2, R0(3),R1(3),R2(3),R3(3),S0(3),S1(3),S2(3),T0(3),T1(3)
3, P11,P12,P13,P14,P21,P22,P23,P31,P32,P41
C FOR EXPLANATION OF VARIABLES IN /TRPCOP/ SEE SUBROUTINE TRP001
C
COMMON /TRPCOT/ X(3),Y(3),Z(3),DX(3),DY(3),DX2(3),DY2(3)
1, SL(3),SL2(3),ZT(3,3),ZTT(3,3),ZUV(3)
2, AP,BP,CP,DP
C FOR AN EXPLANATION OF VARIABLES IN /TRPCOT/ SEE SUBROUTINE TRP00
C
IF (KK.EQ.6) GOTO 60
GOTO (10,20,30,40,50),KK
10 TRICPF= T0(KSE) + T*T1(KSE)
RETURN
20 TRICPF= S0(KSE) + T*(S1(KSE)+T*S2(KSE))
RETURN
30 TRICPF= R0(KSE) + T*(R1(KSE)+T*(R2(KSE)+T*R3(KSE)))
RETURN
40 TRICPF= Q0(KSE)+T*(Q1(KSE)+T*(Q2(KSE)+T*(Q3(KSE)+T*Q4(KSE))))
RETURN
50 TRICPF= P0(KSE)+T*(P1(KSE)+T*(P2(KSE)+T*(P3(KSE)+T*(P4(KSE)+
1 T*P5(KSE))))) - CL
RETURN
60 XX4= XX4F + COR*T
YY4= YY4F + SIR*T
U= AP*XX4 + BP*YY4
V= CP*XX4 + DP*YY4
H0=P0(1)+V*(P1(1)+V*(P2(1)+V*(P3(1)+V*(P4(1)+V*P5(1)))))
H1=P1(2)+V*(P11+V*(P12+V*(P13+V*P14)))
H2=P2(2)+V*(P21+V*(P22+V*P23))
H3=P3(2)+V*(P31+V*P32)
H4=P4(2)+V*P41
TRICPF=H0+U*(H1+U*(H2+U*(H3+U*(H4+U*P5(2))))) - CL
RETURN
END
C PROGRAM TPDEM2
C
C T R I ANGLE C ONTOUR P LOTTING
C
C DEMONSTRATION OF THE TRICP SUBROUTINES
C (MODE=1 AND MODE=3)
C
C COPYRIGHT(C): A. PREUSSER, 1984
C SEN. F. WISS. U. FORSCH.
C PROJEKTGRUPPE PARALLELRECHNEN
C HEILBRONNER STR. 10
C D-1000 BERLIN 31
C
C ONLY (STANDARD)-CALCOMP PLOT COMMANDS ARE USED
C
C ND= NUMBER OF DATA POINTS = 47
C NT= NUMBER OF TRIANGLES = 53
C
C WORK SPACE
DIMENSION PD(5,47),IWK(348)
C IWK(3*NT + 4*ND +1)
C
DIMENSION XD(47),YD(47),ZD(47),C(8)
DIMENSION JP(11),X(3),Y(3)
DATA ND/47/
C
DATA XD
1 /21.95,21.95,19.70,19.70,17.50,17.50,15.30,15.30,13.00,13.00,
2 10.80,10.80, 8.60, 8.60, 7.04, 6.77, 6.48, 6.01, 5.44, 5.59,
3 5.37, 5.51, 4.84, 4.73, 4.77, 4.32, 3.55, 3.56, 3.80, 3.49,
4 3.09, 3.06, 2.57, 2.58, 2.42, 2.44, 2.64, 2.03, 1.24, 1.43,
5 1.96, 1.36, .53, 1.28, .68, .90, .29/
C
DATA YD
1 / 5.10, 7.25, 5.10, 7.25, 5.10, 7.25, 5.10, 7.25, 5.10, 7.25,
2 5.10, 7.25, 5.10, 7.25, 5.10, 6.17, 7.29, 5.10, 5.73, 6.55,
3 7.42, 4.82, 6.66, 7.87, 4.10, 5.45, 7.05, 8.80, 3.45, 4.46,
4 5.03, 6.15, 7.06, 7.84, 8.83, 3.52, 4.76, 7.13, 8.22, 3.83,
5 4.58, 6.96, 7.37, 4.40, 6.78, 4.30, 6.68/
C
DATA ZD
1 /0.1, 0.1, 0.1, 0.1, 1.0, 1.0, 2.1, 2.1, 3.6, 3.6,
2 5.0, 5.0, 6.8, 7.0, 8.5, 0.2, 8.2, 8.5, 1.0, 2.7,
3 6.0, 8.0, 3.0, 1.8, 4.5, 0.9, 2.2, 0.6, 7.0, 4.5,
4 3.0, 2.2, 9.8, 1.8, 1.9, 7.0, 4.8, 9.5, 1.0, 2.2,
5 8.2, 1.8, 0.5,15.0, 1.0, 4.0, 0.0 /
C
C IWK(2...3*NT+1) POINT NUMBERS FOR TRIANGLES
DATA IWK
1 / 53,
2 1, 2, 3, 3, 2, 4, 5, 3, 4,
3 5, 4, 6, 7, 5, 6, 7, 6, 8,
4 9, 7, 8, 9, 8, 10, 11, 9, 10,
5 11, 10, 12, 13, 11, 12, 13, 12, 14,
6 13, 14, 16, 15, 13, 16, 16, 14, 17,
7 18, 15, 16, 18, 16, 19, 19, 16, 20,
8 20, 16, 17, 20, 17, 21, 22, 18, 19,
9 26, 22, 19, 26, 19, 23, 23, 19, 20,
A 23, 20, 21, 23, 21, 24, 26, 23, 27,
B 32, 26, 27, 27, 23, 24, 33, 32, 27,
C 33, 27, 34, 34, 27, 28, 27, 24, 28,
D 38, 33, 34, 34, 28, 35, 38, 34, 39,
E 39, 34, 35, 42, 38, 39, 43, 42, 39,
F 43, 45, 42, 47, 45, 43, 46, 40, 44,
G 40, 41, 44, 40, 36, 41, 36, 37, 41,
H 36, 30, 37, 36, 29, 30, 30, 31, 37,
I 29, 25, 30, 25, 22, 26, 25, 26, 30,
J 30, 26, 31, 31, 26, 32 /
C
C POINT NUMBERS FOR WHICH DERIVATIVES WILL BE SET TO ZERO
DATA JP/7,8,9,10,11,12,13,14,16,19,26/
C
C CONTOUR LEVELS
DATA NC/8/
DATA C / 0., 2., 4., 6., 8., 10., 12., 14./
C
C INITIALIZE PLOTTING
CALL PLOTS
C
C EXAMPLE 3.
C
C PLOT CONTOUR LINES (MODE=1)
CALL TRICP (XD,YD,ZD,ND,C,NC,PD,IWK,1)
C
C PLOT TRIANGLE SIDES
NT= IWK(1)
C = NUMBER OF TRIANGLES
DO 130 I=1,NT
DO 120 J=1,3
N= 3*I+2-J
NH= IWK(N)
X(J)= XD(NH)
Y(J)= YD(NH)
120 CONTINUE
CALL PLOT (X(1),Y(1),3)
CALL PLOT (X(2),Y(2),2)
CALL PLOT (X(3),Y(3),2)
CALL PLOT (X(1),Y(1),2)
130 CONTINUE
C
C PLOT POINT NUMBERS
DO 150 I=1,ND
AI= I
150 CALL NUMBER (XD(I)-0.3, YD(I)+0.1, 0.2, AI, 0., -1)
C
C CLOSE PLOT
CALL PLOT(0.,0.,999)
C
C CHANGE PARTIAL DERIVATIVES (SLOPES AND CURVATURES)
C AT SOME DATA POINTS AND REDRAW THE PLOT WITH MODE= 3
C
C SET THE SLOPES AND CURVATURES AT THE JP-POINTS TO ZERO
DO 200 J=1,11
JPP= JP(J)
DO 200 I=1,5
200 PD(I,JPP)= 0.
C DEFINE SLOPES FOR POINTS 8,10,12,14
DO 220 J=8,14,2
PD(1,J)= -1.
PD(2,J)= 8.5
220 CONTINUE
C DEFINE SLOPES FOR POINTS 7,9,11,13
DO 240 J=7,13,2
PD(1,J)= -1.
PD(2,J)= -8.5
240 CONTINUE
C DEFINE SLOPE IN Y-DIRECTION AT POINT 26
PD(2,26)= 2.
C
C EXAMPLE 4.
C
CALL PLOTS
C
C PLOT CONTOUR LINES (MODE=3)
CALL TRICP (XD,YD,ZD,ND,C,NC,PD,IWK,3)
C
C PLOT TRIANGLE SIDES
DO 330 I=1,NT
DO 320 J=1,3
N= 3*I+2-J
NH= IWK(N)
X(J)= XD(NH)
Y(J)= YD(NH)
320 CONTINUE
CALL PLOT (X(1),Y(1),3)
CALL PLOT (X(2),Y(2),2)
CALL PLOT (X(3),Y(3),2)
CALL PLOT (X(1),Y(1),2)
330 CONTINUE
C
CALL PLOT(0.,0.,999)
C
STOP
END
SUBROUTINE IDTANG(NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
C THIS SUBROUTINE PERFORMS TRIANGULATION. IT DIVIDES THE X-Y
C PLANE INTO A NUMBER OF TRIANGLES ACCORDING TO GIVEN DATA
C POINTS IN THE PLANE, DETERMINES LINE SEGMENTS THAT FORM THE
C BORDER OF DATA AREA, AND DETERMINES THE TRIANGLE NUMBERS
C CORRESPONDING TO THE BORDER LINE SEGMENTS.
C AT COMPLETION, POINT NUMBERS OF THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE. POINT NUMBERS OF THE END POINTS
C OF EACH BORDER LINE SEGMENT ARE LISTED COUNTER-CLOCKWISE,
C LISTING ORDER OF THE LINE SEGMENTS BEING COUNTER-CLOCKWISE.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C THIS SUBROUTINE CALLS THE IDXCHG FUNCTION.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD = ARRAY OF DIMENSION NDP CONTAINING THE
C X COORDINATES OF THE DATA POINTS,
C YD = ARRAY OF DIMENSION NDP CONTAINING THE
C Y COORDINATES OF THE DATA POINTS.
C THE OUTPUT PARAMETERS ARE
C NT = NUMBER OF TRIANGLES,
C IPT = INTEGER ARRAY OF DIMENSION 6*NDP-15, WHERE THE
C POINT NUMBERS OF THE VERTEXES OF THE (IT)TH
C TRIANGLE ARE TO BE STORED AS THE (3*IT-2)ND,
C (3*IT-1)ST, AND (3*IT)TH ELEMENTS,
C IT=1,2,...,NT,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = INTEGER ARRAY OF DIMENSION 6*NDP, WHERE THE
C POINT NUMBERS OF THE END POINTS OF THE (IL)TH
C BORDER LINE SEGMENT AND ITS RESPECTIVE TRIANGLE
C NUMBER ARE TO BE STORED AS THE (3*IL-2)ND,
C (3*IL-1)ST, AND (3*IL)TH ELEMENTS,
C IL=1,2,..., NL.
C THE OTHER PARAMETERS ARE
C IWL = INTEGER ARRAY OF DIMENSION 18*NDP USED
C INTERNALLY AS A WORK AREA,
C IWP = INTEGER ARRAY OF DIMENSION NDP USED
C INTERNALLY AS A WORK AREA,
C WK = ARRAY OF DIMENSION NDP USED INTERNALLY AS A
C WORK AREA.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),IPT(585),IPL(600),
1 IWL(1800),IWP(100),WK(100)
DIMENSION ITF(2)
DATA RATIO/1.0E-6/, NREP/100/, LUN/6/
C STATEMENT FUNCTIONS
DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
SIDE(U1,V1,U2,V2,U3,V3)=(V3-V1)*(U2-U1)-(U3-U1)*(V2-V1)
C PRELIMINARY PROCESSING
10 NDP0=NDP
NDPM1=NDP0-1
IF(NDP0.LT.4) GO TO 90
C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
20 DSQMN=DSQF(XD(1),YD(1),XD(2),YD(2))
IPMN1=1
IPMN2=2
DO 22 IP1=1,NDPM1
X1=XD(IP1)
Y1=YD(IP1)
IP1P1=IP1+1
DO 21 IP2=IP1P1,NDP0
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
IF(DSQI.EQ.0.0) GO TO 91
IF(DSQI.GE.DSQMN) GO TO 21
DSQMN=DSQI
IPMN1=IP1
IPMN2=IP2
21 CONTINUE
22 CONTINUE
DSQ12=DSQMN
XDMP=(XD(IPMN1)+XD(IPMN2))/2.0
YDMP=(YD(IPMN1)+YD(IPMN2))/2.0
C SORTS THE OTHER (NDP-2) DATA POINTS IN ASCENDING ORDER OF
C DISTANCE FROM THE MIDPOINT AND STORES THE SORTED DATA POINT
C NUMBERS IN THE IWP ARRAY.
30 JP1=2
DO 31 IP1=1,NDP0
IF(IP1.EQ.IPMN1.OR.IP1.EQ.IPMN2) GO TO 31
JP1=JP1+1
IWP(JP1)=IP1
WK(JP1)=DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
31 CONTINUE
DO 33 JP1=3,NDPM1
DSQMN=WK(JP1)
JPMN=JP1
DO 32 JP2=JP1,NDP0
IF(WK(JP2).GE.DSQMN) GO TO 32
DSQMN=WK(JP2)
JPMN=JP2
32 CONTINUE
ITS=IWP(JP1)
IWP(JP1)=IWP(JPMN)
IWP(JPMN)=ITS
WK(JPMN)=WK(JP1)
33 CONTINUE
C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
35 AR=DSQ12*RATIO
X1=XD(IPMN1)
Y1=YD(IPMN1)
DX21=XD(IPMN2)-X1
DY21=YD(IPMN2)-Y1
DO 36 JP=3,NDP0
IP=IWP(JP)
IF(ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21).GT.AR)
1 GO TO 37
36 CONTINUE
GO TO 92
37 IF(JP.EQ.3) GO TO 40
JPMX=JP
JP=JPMX+1
DO 38 JPC=4,JPMX
JP=JP-1
IWP(JP)=IWP(JP-1)
38 CONTINUE
IWP(3)=IP
C FORMS THE FIRST TRIANGLE. STORES POINT NUMBERS OF THE VER-
C TEXES OF THE TRIANGLE IN THE IPT ARRAY, AND STORES POINT NUM-
C BERS OF THE BORDER LINE SEGMENTS AND THE TRIANGLE NUMBER IN
C THE IPL ARRAY.
40 IP1=IPMN1
IP2=IPMN2
IP3=IWP(3)
IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
1 .GE.0.0) GO TO 41
IP1=IPMN2
IP2=IPMN1
41 NT0=1
NTT3=3
IPT(1)=IP1
IPT(2)=IP2
IPT(3)=IP3
NL0=3
NLT3=9
IPL(1)=IP1
IPL(2)=IP2
IPL(3)=1
IPL(4)=IP2
IPL(5)=IP3
IPL(6)=1
IPL(7)=IP3
IPL(8)=IP1
IPL(9)=1
C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
50 DO 79 JP1=4,NDP0
IP1=IWP(JP1)
X1=XD(IP1)
Y1=YD(IP1)
C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
IP2=IPL(1)
JPMN=1
DXMN=XD(IP2)-X1
DYMN=YD(IP2)-Y1
DSQMN=DXMN**2+DYMN**2
ARMN=DSQMN*RATIO
JPMX=1
DXMX=DXMN
DYMX=DYMN
DSQMX=DSQMN
ARMX=ARMN
DO 52 JP2=2,NL0
IP2=IPL(3*JP2-2)
DX=XD(IP2)-X1
DY=YD(IP2)-Y1
AR=DY*DXMN-DX*DYMN
IF(AR.GT.ARMN) GO TO 51
DSQI=DX**2+DY**2
IF(AR.GE.(-ARMN).AND.DSQI.GE.DSQMN) GO TO 51
JPMN=JP2
DXMN=DX
DYMN=DY
DSQMN=DSQI
ARMN=DSQMN*RATIO
51 AR=DY*DXMX-DX*DYMX
IF(AR.LT.(-ARMX)) GO TO 52
DSQI=DX**2+DY**2
IF(AR.LE.ARMX.AND.DSQI.GE.DSQMX) GO TO 52
JPMX=JP2
DXMX=DX
DYMX=DY
DSQMX=DSQI
ARMX=DSQMX*RATIO
52 CONTINUE
IF(JPMX.LT.JPMN) JPMX=JPMX+NL0
NSH=JPMN-1
IF(NSH.LE.0) GO TO 60
C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
NSHT3=NSH*3
DO 53 JP2T3=3,NSHT3,3
JP3T3=JP2T3+NLT3
IPL(JP3T3-2)=IPL(JP2T3-2)
IPL(JP3T3-1)=IPL(JP2T3-1)
IPL(JP3T3) =IPL(JP2T3)
53 CONTINUE
DO 54 JP2T3=3,NLT3,3
JP3T3=JP2T3+NSHT3
IPL(JP2T3-2)=IPL(JP3T3-2)
IPL(JP2T3-1)=IPL(JP3T3-1)
IPL(JP2T3) =IPL(JP3T3)
54 CONTINUE
JPMX=JPMX-NSH
C - ADDS TRIANGLES TO THE IPT ARRAY, UPDATES BORDER LINE
C - SEGMENTS IN THE IPL ARRAY, AND SETS FLAGS FOR THE BORDER
C - LINE SEGMENTS TO BE REEXAMINED IN THE IWL ARRAY.
60 JWL=0
DO 64 JP2=JPMX,NL0
JP2T3=JP2*3
IPL1=IPL(JP2T3-2)
IPL2=IPL(JP2T3-1)
IT =IPL(JP2T3)
C - - ADDS A TRIANGLE TO THE IPT ARRAY.
NT0=NT0+1
NTT3=NTT3+3
IPT(NTT3-2)=IPL2
IPT(NTT3-1)=IPL1
IPT(NTT3) =IP1
C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
IF(JP2.NE.JPMX) GO TO 61
IPL(JP2T3-1)=IP1
IPL(JP2T3) =NT0
61 IF(JP2.NE.NL0) GO TO 62
NLN=JPMX+1
NLNT3=NLN*3
IPL(NLNT3-2)=IP1
IPL(NLNT3-1)=IPL(1)
IPL(NLNT3) =NT0
C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
C - - LINE SEGMENTS.
62 ITT3=IT*3
IPTI=IPT(ITT3-2)
IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63
IPTI=IPT(ITT3-1)
IF(IPTI.NE.IPL1.AND.IPTI.NE.IPL2) GO TO 63
IPTI=IPT(ITT3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
63 IF(IDXCHG(XD,YD,IP1,IPTI,IPL1,IPL2).EQ.0) GO TO 64
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
IPT(ITT3-2)=IPTI
IPT(ITT3-1)=IPL1
IPT(ITT3) =IP1
IPT(NTT3-1)=IPTI
IF(JP2.EQ.JPMX) IPL(JP2T3)=IT
IF(JP2.EQ.NL0.AND.IPL(3).EQ.IT) IPL(3)=NT0
C - - SETS FLAGS IN THE IWL ARRAY.
JWL=JWL+4
IWL(JWL-3)=IPL1
IWL(JWL-2)=IPTI
IWL(JWL-1)=IPTI
IWL(JWL) =IPL2
64 CONTINUE
NL0=NLN
NLT3=NLNT3
NLF=JWL/2
IF(NLF.EQ.0) GO TO 79
C - IMPROVES TRIANGULATION.
70 NTT3P3=NTT3+3
DO 78 IREP=1,NREP
DO 76 ILF=1,NLF
ILFT2=ILF*2
IPL1=IWL(ILFT2-1)
IPL2=IWL(ILFT2)
C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
C - - THE FLAGGED LINE SEGMENT.
NTF=0
DO 71 ITT3R=3,NTT3,3
ITT3=NTT3P3-ITT3R
IPT1=IPT(ITT3-2)
IPT2=IPT(ITT3-1)
IPT3=IPT(ITT3)
IF(IPL1.NE.IPT1.AND.IPL1.NE.IPT2.AND.
1 IPL1.NE.IPT3) GO TO 71
IF(IPL2.NE.IPT1.AND.IPL2.NE.IPT2.AND.
1 IPL2.NE.IPT3) GO TO 71
NTF=NTF+1
ITF(NTF)=ITT3/3
IF(NTF.EQ.2) GO TO 72
71 CONTINUE
IF(NTF.LT.2) GO TO 76
C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
C - - ON THE LINE SEGMENT.
72 IT1T3=ITF(1)*3
IPTI1=IPT(IT1T3-2)
IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73
IPTI1=IPT(IT1T3-1)
IF(IPTI1.NE.IPL1.AND.IPTI1.NE.IPL2) GO TO 73
IPTI1=IPT(IT1T3)
73 IT2T3=ITF(2)*3
IPTI2=IPT(IT2T3-2)
IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74
IPTI2=IPT(IT2T3-1)
IF(IPTI2.NE.IPL1.AND.IPTI2.NE.IPL2) GO TO 74
IPTI2=IPT(IT2T3)
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
74 IF(IDXCHG(XD,YD,IPTI1,IPTI2,IPL1,IPL2).EQ.0)
1 GO TO 76
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
IPT(IT1T3-2)=IPTI1
IPT(IT1T3-1)=IPTI2
IPT(IT1T3) =IPL1
IPT(IT2T3-2)=IPTI2
IPT(IT2T3-1)=IPTI1
IPT(IT2T3) =IPL2
C - - SETS NEW FLAGS.
JWL=JWL+8
IWL(JWL-7)=IPL1
IWL(JWL-6)=IPTI1
IWL(JWL-5)=IPTI1
IWL(JWL-4)=IPL2
IWL(JWL-3)=IPL2
IWL(JWL-2)=IPTI2
IWL(JWL-1)=IPTI2
IWL(JWL) =IPL1
DO 75 JLT3=3,NLT3,3
IPLJ1=IPL(JLT3-2)
IPLJ2=IPL(JLT3-1)
IF((IPLJ1.EQ.IPL1.AND.IPLJ2.EQ.IPTI2).OR.
1 (IPLJ2.EQ.IPL1.AND.IPLJ1.EQ.IPTI2))
2 IPL(JLT3)=ITF(1)
IF((IPLJ1.EQ.IPL2.AND.IPLJ2.EQ.IPTI1).OR.
1 (IPLJ2.EQ.IPL2.AND.IPLJ1.EQ.IPTI1))
2 IPL(JLT3)=ITF(2)
75 CONTINUE
76 CONTINUE
NLFC=NLF
NLF=JWL/2
IF(NLF.EQ.NLFC) GO TO 79
C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
JWL=0
JWL1MN=(NLFC+1)*2
NLFT2=NLF*2
DO 77 JWL1=JWL1MN,NLFT2,2
JWL=JWL+2
IWL(JWL-1)=IWL(JWL1-1)
IWL(JWL) =IWL(JWL1)
77 CONTINUE
NLF=JWL/2
78 CONTINUE
79 CONTINUE
C REARRANGES THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE.
80 DO 81 ITT3=3,NTT3,3
IP1=IPT(ITT3-2)
IP2=IPT(ITT3-1)
IP3=IPT(ITT3)
IF(SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3))
1 .GE.0.0) GO TO 81
IPT(ITT3-2)=IP2
IPT(ITT3-1)=IP1
81 CONTINUE
NT=NT0
NL=NL0
RETURN
C ERROR EXIT
90 WRITE (LUN,2090) NDP0
GO TO 93
91 WRITE (LUN,2091) NDP0,IP1,IP2,X1,Y1
GO TO 93
92 WRITE (LUN,2092) NDP0
93 WRITE (LUN,2093)
NT=0
RETURN
C FORMAT STATEMENTS
2090 FORMAT(1X/23H *** NDP LESS THAN 4./8H NDP =,I5)
2091 FORMAT(1X/29H *** IDENTICAL DATA POINTS./
1 8H NDP =,I5,5X,5HIP1 =,I5,5X,5HIP2 =,I5,
2 5X,4HXD =,E12.4,5X,4HYD =,E12.4)
2092 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS./
1 8H NDP =,I5)
2093 FORMAT(35H ERROR DETECTED IN ROUTINE IDTANG/)
END
FUNCTION IDXCHG(X,Y,I1,I2,I3,I4)
C THIS FUNCTION DETERMINES WHETHER OR NOT THE EXCHANGE OF TWO
C TRIANGLES IS NECESSARY ON THE BASIS OF MAX-MIN-ANGLE CRITERION
C BY C. L. LAWSON.
C THE INPUT PARAMETERS ARE
C X,Y = ARRAYS CONTAINING THE COORDINATES OF THE DATA
C POINTS,
C I1,I2,I3,I4 = POINT NUMBERS OF FOUR POINTS P1, P2,
C P3, AND P4 THAT FORM A QUADRILATERAL WITH P3
C AND P4 CONNECTED DIAGONALLY.
C THIS FUNCTION RETURNS AN INTEGER VALUE 1 (ONE) WHEN AN EX-
C CHANGE IS NECESSARY, AND 0 (ZERO) OTHERWISE.
C DECLARATION STATEMENTS
DIMENSION X(100),Y(100)
EQUIVALENCE (C2SQ,C1SQ),(A3SQ,B2SQ),(B3SQ,A1SQ),
1 (A4SQ,B1SQ),(B4SQ,A2SQ),(C4SQ,C3SQ)
C PRELIMINARY PROCESSING
10 X1=X(I1)
Y1=Y(I1)
X2=X(I2)
Y2=Y(I2)
X3=X(I3)
Y3=Y(I3)
X4=X(I4)
Y4=Y(I4)
C CALCULATION
20 IDX=0
U3=(Y2-Y3)*(X1-X3)-(X2-X3)*(Y1-Y3)
U4=(Y1-Y4)*(X2-X4)-(X1-X4)*(Y2-Y4)
IF(U3*U4.LE.0.0) GO TO 30
U1=(Y3-Y1)*(X4-X1)-(X3-X1)*(Y4-Y1)
U2=(Y4-Y2)*(X3-X2)-(X4-X2)*(Y3-Y2)
A1SQ=(X1-X3)**2+(Y1-Y3)**2
B1SQ=(X4-X1)**2+(Y4-Y1)**2
C1SQ=(X3-X4)**2+(Y3-Y4)**2
A2SQ=(X2-X4)**2+(Y2-Y4)**2
B2SQ=(X3-X2)**2+(Y3-Y2)**2
C3SQ=(X2-X1)**2+(Y2-Y1)**2
S1SQ=U1*U1/(C1SQ*AMAX1(A1SQ,B1SQ))
S2SQ=U2*U2/(C2SQ*AMAX1(A2SQ,B2SQ))
S3SQ=U3*U3/(C3SQ*AMAX1(A3SQ,B3SQ))
S4SQ=U4*U4/(C4SQ*AMAX1(A4SQ,B4SQ))
IF(AMIN1(S1SQ,S2SQ).LT.AMIN1(S3SQ,S4SQ)) IDX=1
30 IDXCHG=IDX
RETURN
END
SUBROUTINE IDCLDP(NDP,XD,YD,NCP,IPC)
C THIS SUBROUTINE SELECTS SEVERAL DATA POINTS THAT ARE CLOSEST
C TO EACH OF THE DATA POINT.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD = ARRAYS OF DIMENSION NDP CONTAINING THE X AND Y
C COORDINATES OF THE DATA POINTS,
C NCP = NUMBER OF DATA POINTS CLOSEST TO EACH DATA
C POINTS.
C THE OUTPUT PARAMETER IS
C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP, WHERE THE
C POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C EACH OF THE NDP DATA POINTS ARE TO BE STORED.
C THIS SUBROUTINE ARBITRARILY SETS A RESTRICTION THAT NCP MUST
C NOT EXCEED 25.
C THE LUN CONSTANT IN THE DATA INITIALIZATION STATEMENT IS THE
C LOGICAL UNIT NUMBER OF THE STANDARD OUTPUT UNIT AND IS,
C THEREFORE, SYSTEM DEPENDENT.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),IPC(400)
DIMENSION DSQ0(25),IPC0(25)
DATA NCPMX/25/, LUN/6/
C STATEMENT FUNCTION
DSQF(U1,V1,U2,V2)=(U2-U1)**2+(V2-V1)**2
C PRELIMINARY PROCESSING
10 NDP0=NDP
NCP0=NCP
IF(NDP0.LT.2) GO TO 90
IF(NCP0.LT.1.OR.NCP0.GT.NCPMX.OR.NCP0.GE.NDP0) GO TO 90
C CALCULATION
20 DO 59 IP1=1,NDP0
C - SELECTS NCP POINTS.
X1=XD(IP1)
Y1=YD(IP1)
J1=0
DSQMX=0.0
DO 22 IP2=1,NDP0
IF(IP2.EQ.IP1) GO TO 22
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
J1=J1+1
DSQ0(J1)=DSQI
IPC0(J1)=IP2
IF(DSQI.LE.DSQMX) GO TO 21
DSQMX=DSQI
JMX=J1
21 IF(J1.GE.NCP0) GO TO 23
22 CONTINUE
23 IP2MN=IP2+1
IF(IP2MN.GT.NDP0) GO TO 30
DO 25 IP2=IP2MN,NDP0
IF(IP2.EQ.IP1) GO TO 25
DSQI=DSQF(X1,Y1,XD(IP2),YD(IP2))
IF(DSQI.GE.DSQMX) GO TO 25
DSQ0(JMX)=DSQI
IPC0(JMX)=IP2
DSQMX=0.0
DO 24 J1=1,NCP0
IF(DSQ0(J1).LE.DSQMX) GO TO 24
DSQMX=DSQ0(J1)
JMX=J1
24 CONTINUE
25 CONTINUE
C - CHECKS IF ALL THE NCP+1 POINTS ARE COLLINEAR.
30 IP2=IPC0(1)
DX12=XD(IP2)-X1
DY12=YD(IP2)-Y1
DO 31 J3=2,NCP0
IP3=IPC0(J3)
DX13=XD(IP3)-X1
DY13=YD(IP3)-Y1
IF((DY13*DX12-DX13*DY12)**2/(DSQ0(1)*DSQ0(J3))
A .GT.0.06698) GOTO 50
C ----- 0.06698 CORRESPONDS TO AN ANGLE OF ABOUT 15. DEGREES(=SIN**2)
31 CONTINUE
C - SEARCHES FOR THE CLOSEST NONCOLLINEAR POINT.
40 NCLPT=0
DO 43 IP3=1,NDP0
IF(IP3.EQ.IP1) GO TO 43
DO 41 J4=1,NCP0
IF(IP3.EQ.IPC0(J4)) GO TO 43
41 CONTINUE
DX13=XD(IP3)-X1
DY13=YD(IP3)-Y1
DSQI=DSQF(X1,Y1,XD(IP3),YD(IP3))
IF((DY13*DX12-DX13*DY12)**2/(DSQ0(1)*DSQI)
A .LE.0.06698) GO TO 43
IF(NCLPT.EQ.0) GO TO 42
IF(DSQI.GE.DSQMN) GO TO 43
42 NCLPT=1
DSQMN=DSQI
IP3MN=IP3
43 CONTINUE
IF(NCLPT.EQ.0) GO TO 91
DSQMX=DSQMN
IPC0(JMX)=IP3MN
C - REPLACES THE LOCAL ARRAY FOR THE OUTPUT ARRAY.
50 J1=(IP1-1)*NCP0
DO 51 J2=1,NCP0
J1=J1+1
IPC(J1)=IPC0(J2)
51 CONTINUE
59 CONTINUE
RETURN
C ERROR EXIT
90 WRITE (LUN,2090)
GO TO 92
91 WRITE (LUN,2091)
92 WRITE (LUN,2092) NDP0,NCP0
IPC(1)=0
RETURN
C FORMAT STATEMENTS FOR ERROR MESSAGES
2090 FORMAT(1X/41H *** IMPROPER INPUT PARAMETER VALUE(S).)
2091 FORMAT(1X/33H *** ALL COLLINEAR DATA POINTS.)
2092 FORMAT(8H NDP =,I5,5X,5HNCP =,I5/
1 35H ERROR DETECTED IN ROUTINE IDCLDP/)
END
SUBROUTINE IDPDRV(NDP,XD,YD,ZD,NCP,IPC,PD)
C THIS SUBROUTINE ESTIMATES PARTIAL DERIVATIVES OF THE FIRST AND
C SECOND ORDER AT THE DATA POINTS.
C THE INPUT PARAMETERS ARE
C NDP = NUMBER OF DATA POINTS,
C XD,YD,ZD = ARRAYS OF DIMENSION NDP CONTAINING THE X,
C Y, AND Z COORDINATES OF THE DATA POINTS,
C NCP = NUMBER OF ADDITIONAL DATA POINTS USED FOR ESTI-
C MATING PARTIAL DERIVATIVES AT EACH DATA POINT,
C IPC = INTEGER ARRAY OF DIMENSION NCP*NDP CONTAINING
C THE POINT NUMBERS OF NCP DATA POINTS CLOSEST TO
C EACH OF THE NDP DATA POINTS.
C THE OUTPUT PARAMETER IS
C PD = ARRAY OF DIMENSION 5*NDP, WHERE THE ESTIMATED
C ZX, ZY, ZXX, ZXY, AND ZYY VALUES AT THE DATA
C POINTS ARE TO BE STORED.
C DECLARATION STATEMENTS
DIMENSION XD(100),YD(100),ZD(100),IPC(400),PD(500)
REAL NMX,NMY,NMZ,NMXX,NMXY,NMYX,NMYY
C PRELIMINARY PROCESSING
10 NDP0=NDP
NCP0=NCP
NCPM1=NCP0-1
C ESTIMATION OF ZX AND ZY
20 DO 24 IP0=1,NDP0
X0=XD(IP0)
Y0=YD(IP0)
Z0=ZD(IP0)
NMX=0.0
NMY=0.0
NMZ=0.0
JIPC0=NCP0*(IP0-1)
DO 23 IC1=1,NCPM1
JIPC=JIPC0+IC1
IPI=IPC(JIPC)
DX1=XD(IPI)-X0
DY1=YD(IPI)-Y0
DZ1=ZD(IPI)-Z0
IC2MN=IC1+1
DO 22 IC2=IC2MN,NCP0
JIPC=JIPC0+IC2
IPI=IPC(JIPC)
DX2=XD(IPI)-X0
DY2=YD(IPI)-Y0
DNMZ=DX1*DY2-DY1*DX2
IF(DNMZ.EQ.0.0) GO TO 22
DZ2=ZD(IPI)-Z0
DNMX=DY1*DZ2-DZ1*DY2
DNMY=DZ1*DX2-DX1*DZ2
IF(DNMZ.GE.0.0) GO TO 21
DNMX=-DNMX
DNMY=-DNMY
DNMZ=-DNMZ
21 NMX=NMX+DNMX
NMY=NMY+DNMY
NMZ=NMZ+DNMZ
22 CONTINUE
23 CONTINUE
JPD0=5*IP0
PD(JPD0-4)=-NMX/NMZ
PD(JPD0-3)=-NMY/NMZ
24 CONTINUE
C ESTIMATION OF ZXX, ZXY, AND ZYY
30 DO 34 IP0=1,NDP0
JPD0=JPD0+5
X0=XD(IP0)
JPD0=5*IP0
Y0=YD(IP0)
ZX0=PD(JPD0-4)
ZY0=PD(JPD0-3)
NMXX=0.0
NMXY=0.0
NMYX=0.0
NMYY=0.0
NMZ =0.0
JIPC0=NCP0*(IP0-1)
DO 33 IC1=1,NCPM1
JIPC=JIPC0+IC1
IPI=IPC(JIPC)
DX1=XD(IPI)-X0
DY1=YD(IPI)-Y0
JPD=5*IPI
DZX1=PD(JPD-4)-ZX0
DZY1=PD(JPD-3)-ZY0
IC2MN=IC1+1
DO 32 IC2=IC2MN,NCP0
JIPC=JIPC0+IC2
IPI=IPC(JIPC)
DX2=XD(IPI)-X0
DY2=YD(IPI)-Y0
DNMZ =DX1*DY2 -DY1*DX2
IF(DNMZ.EQ.0.0) GO TO 32
JPD=5*IPI
DZX2=PD(JPD-4)-ZX0
DZY2=PD(JPD-3)-ZY0
DNMXX=DY1*DZX2-DZX1*DY2
DNMXY=DZX1*DX2-DX1*DZX2
DNMYX=DY1*DZY2-DZY1*DY2
DNMYY=DZY1*DX2-DX1*DZY2
IF(DNMZ.GE.0.0) GO TO 31
DNMXX=-DNMXX
DNMXY=-DNMXY
DNMYX=-DNMYX
DNMYY=-DNMYY
DNMZ =-DNMZ
31 NMXX=NMXX+DNMXX
NMXY=NMXY+DNMXY
NMYX=NMYX+DNMYX
NMYY=NMYY+DNMYY
NMZ =NMZ +DNMZ
32 CONTINUE
33 CONTINUE
PD(JPD0-2)=-NMXX/NMZ
PD(JPD0-1)=-(NMXY+NMYX)/(2.0*NMZ)
PD(JPD0) =-NMYY/NMZ
34 CONTINUE
RETURN
END
.