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