1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
|
SUBROUTINE CONTNG (NDP,XD,YD,NT,IPT,NL,IPL,IWL,IWP,WK)
C
C +-----------------------------------------------------------------+
C | |
C | Copyright (C) 1986 by UCAR |
C | University Corporation for Atmospheric Research |
C | All Rights Reserved |
C | |
C | NCARGRAPHICS Version 1.00 |
C | |
C +-----------------------------------------------------------------+
C
C
C
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 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 = ARRAY OF DIMENSION 6*NDP-15, WHERE THE POINT
C NUMBERS OF THE VERTEXES OF THE (IT)TH TRIANGLE
C ARE TO BE STORED AS THE (3*IT-2)ND, (3*IT-1)ST,
C AND (3*IT)TH ELEMENTS, IT=1,2,...,NT,
C NL = NUMBER OF BORDER LINE SEGMENTS,
C IPL = ARRAY OF DIMENSION 6*NDP, WHERE THE POINT
C NUMBERS OF THE END POINTS OF THE (IL)TH BORDER
C LINE SEGMENT AND ITS RESPECTIVE TRIANGLE NUMBER
C ARE TO BE STORED AS THE (3*IL-2)ND, (3*IL-1)ST,
C AND (3*IL)TH ELEMENTS, 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
C
SAVE
C
INTEGER CONXCH
COMMON /CONRA3/ IREC
DIMENSION XD(*) ,YD(*) ,IPT(*) ,IPL(*) ,
1 IWL(*) ,IWP(*) ,WK(*)
DIMENSION ITF(2)
CHARACTER*4 IP1C, IP2C
CHARACTER*64 ITEMP
DATA RATIO/1.0E-6/, NREP/100/
C
C STATEMENT FUNCTIONS
C
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
C PRELIMINARY PROCESSING
C
NDPM1 = NDP-1
C
C DETERMINES THE CLOSEST PAIR OF DATA POINTS AND THEIR MIDPOINT.
C
DSQMN = DSQF(XD(1),YD(1),XD(2),YD(2))
IPMN1 = 1
IPMN2 = 2
DO 140 IP1=1,NDPM1
X1 = XD(IP1)
Y1 = YD(IP1)
IP1P1 = IP1+1
DO 130 IP2=IP1P1,NDP
DSQI = DSQF(X1,Y1,XD(IP2),YD(IP2))
IF (DSQI .NE. 0.) GO TO 120
C
C ERROR, IDENTICAL INPUT DATA POINTS
C
ITEMP = ' CONTNG-IDENTICAL INPUT DATA POINTS FOUND
1 AT AND '
C
C + NOAO - FTN internal writes rewritten as calls to encode for IRAF
C
C WRITE(IP1C,'(I4)')IP1
C WRITE(IP2C,'(I4)')IP2
call encode (4, '(I4)', ip1c, ip1)
call encode (4, '(I4)', ip2c, ip2)
C - NOAO
C
CALL SETER (ITEMP,1,1)
ITEMP(46:49) = IP1C
ITEMP(55:58) = IP2C
RETURN
120 IF (DSQI .GE. DSQMN) GO TO 130
DSQMN = DSQI
IPMN1 = IP1
IPMN2 = IP2
130 CONTINUE
140 CONTINUE
DSQ12 = DSQMN
XDMP = (XD(IPMN1)+XD(IPMN2))/2.0
YDMP = (YD(IPMN1)+YD(IPMN2))/2.0
C
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.
C
JP1 = 2
DO 150 IP1=1,NDP
IF (IP1.EQ.IPMN1 .OR. IP1.EQ.IPMN2) GO TO 150
JP1 = JP1+1
IWP(JP1) = IP1
WK(JP1) = DSQF(XDMP,YDMP,XD(IP1),YD(IP1))
150 CONTINUE
DO 170 JP1=3,NDPM1
DSQMN = WK(JP1)
JPMN = JP1
DO 160 JP2=JP1,NDP
IF (WK(JP2) .GE. DSQMN) GO TO 160
DSQMN = WK(JP2)
JPMN = JP2
160 CONTINUE
ITS = IWP(JP1)
IWP(JP1) = IWP(JPMN)
IWP(JPMN) = ITS
WK(JPMN) = WK(JP1)
170 CONTINUE
C
C IF NECESSARY, MODIFIES THE ORDERING IN SUCH A WAY THAT THE
C FIRST THREE DATA POINTS ARE NOT COLLINEAR.
C
AR = DSQ12*RATIO
X1 = XD(IPMN1)
Y1 = YD(IPMN1)
DX21 = XD(IPMN2)-X1
DY21 = YD(IPMN2)-Y1
DO 180 JP=3,NDP
IP = IWP(JP)
IF (ABS((YD(IP)-Y1)*DX21-(XD(IP)-X1)*DY21) .GT. AR) GO TO 190
180 CONTINUE
CALL SETER (' CONTNG - ALL COLLINEAR DATA POINTS',1,1)
190 IF (JP .EQ. 3) GO TO 210
JPMX = JP
JP = JPMX+1
DO 200 JPC=4,JPMX
JP = JP-1
IWP(JP) = IWP(JP-1)
200 CONTINUE
IWP(3) = IP
C
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.
C
210 IP1 = IPMN1
IP2 = IPMN2
IP3 = IWP(3)
IF (SIDE(XD(IP1),YD(IP1),XD(IP2),YD(IP2),XD(IP3),YD(IP3)) .GE.
1 0.0) GO TO 220
IP1 = IPMN2
IP2 = IPMN1
220 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
C ADDS THE REMAINING (NDP-3) DATA POINTS, ONE BY ONE.
C
DO 400 JP1=4,NDP
IP1 = IWP(JP1)
X1 = XD(IP1)
Y1 = YD(IP1)
C
C - DETERMINES THE VISIBLE BORDER LINE SEGMENTS.
C
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 240 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 230
DSQI = DX**2+DY**2
IF (AR.GE.(-ARMN) .AND. DSQI.GE.DSQMN) GO TO 230
JPMN = JP2
DXMN = DX
DYMN = DY
DSQMN = DSQI
ARMN = DSQMN*RATIO
230 AR = DY*DXMX-DX*DYMX
IF (AR .LT. (-ARMX)) GO TO 240
DSQI = DX**2+DY**2
IF (AR.LE.ARMX .AND. DSQI.GE.DSQMX) GO TO 240
JPMX = JP2
DXMX = DX
DYMX = DY
DSQMX = DSQI
ARMX = DSQMX*RATIO
240 CONTINUE
IF (JPMX .LT. JPMN) JPMX = JPMX+NL0
NSH = JPMN-1
IF (NSH .LE. 0) GO TO 270
C
C - SHIFTS (ROTATES) THE IPL ARRAY TO HAVE THE INVISIBLE BORDER
C - LINE SEGMENTS CONTAINED IN THE FIRST PART OF THE IPL ARRAY.
C
NSHT3 = NSH*3
DO 250 JP2T3=3,NSHT3,3
JP3T3 = JP2T3+NLT3
IPL(JP3T3-2) = IPL(JP2T3-2)
IPL(JP3T3-1) = IPL(JP2T3-1)
IPL(JP3T3) = IPL(JP2T3)
250 CONTINUE
DO 260 JP2T3=3,NLT3,3
JP3T3 = JP2T3+NSHT3
IPL(JP2T3-2) = IPL(JP3T3-2)
IPL(JP2T3-1) = IPL(JP3T3-1)
IPL(JP2T3) = IPL(JP3T3)
260 CONTINUE
JPMX = JPMX-NSH
C
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.
C
270 JWL = 0
DO 310 JP2=JPMX,NL0
JP2T3 = JP2*3
IPL1 = IPL(JP2T3-2)
IPL2 = IPL(JP2T3-1)
IT = IPL(JP2T3)
C
C - - ADDS A TRIANGLE TO THE IPT ARRAY.
C
NT0 = NT0+1
NTT3 = NTT3+3
IPT(NTT3-2) = IPL2
IPT(NTT3-1) = IPL1
IPT(NTT3) = IP1
C
C - - UPDATES BORDER LINE SEGMENTS IN THE IPL ARRAY.
C
IF (JP2 .NE. JPMX) GO TO 280
IPL(JP2T3-1) = IP1
IPL(JP2T3) = NT0
280 IF (JP2 .NE. NL0) GO TO 290
NLN = JPMX+1
NLNT3 = NLN*3
IPL(NLNT3-2) = IP1
IPL(NLNT3-1) = IPL(1)
IPL(NLNT3) = NT0
C
C - - DETERMINES THE VERTEX THAT DOES NOT LIE ON THE BORDER
C - - LINE SEGMENTS.
C
290 ITT3 = IT*3
IPTI = IPT(ITT3-2)
IF (IPTI.NE.IPL1 .AND. IPTI.NE.IPL2) GO TO 300
IPTI = IPT(ITT3-1)
IF (IPTI.NE.IPL1 .AND. IPTI.NE.IPL2) GO TO 300
IPTI = IPT(ITT3)
C
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
C
300 IF (CONXCH(XD,YD,IP1,IPTI,IPL1,IPL2) .EQ. 0) GO TO 310
C
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
C
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
C - - SETS FLAGS IN THE IWL ARRAY.
C
JWL = JWL+4
IWL(JWL-3) = IPL1
IWL(JWL-2) = IPTI
IWL(JWL-1) = IPTI
IWL(JWL) = IPL2
310 CONTINUE
NL0 = NLN
NLT3 = NLNT3
NLF = JWL/2
IF (NLF .EQ. 0) GO TO 400
C
C - IMPROVES TRIANGULATION.
C
NTT3P3 = NTT3+3
DO 390 IREP=1,NREP
DO 370 ILF=1,NLF
ILFT2 = ILF*2
IPL1 = IWL(ILFT2-1)
IPL2 = IWL(ILFT2)
C
C - - LOCATES IN THE IPT ARRAY TWO TRIANGLES ON BOTH SIDES OF
C - - THE FLAGGED LINE SEGMENT.
C
NTF = 0
DO 320 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 320
IF (IPL2.NE.IPT1 .AND. IPL2.NE.IPT2 .AND.
1 IPL2.NE.IPT3) GO TO 320
NTF = NTF+1
ITF(NTF) = ITT3/3
IF (NTF .EQ. 2) GO TO 330
320 CONTINUE
IF (NTF .LT. 2) GO TO 370
C
C - - DETERMINES THE VERTEXES OF THE TRIANGLES THAT DO NOT LIE
C - - ON THE LINE SEGMENT.
C
330 IT1T3 = ITF(1)*3
IPTI1 = IPT(IT1T3-2)
IF (IPTI1.NE.IPL1 .AND. IPTI1.NE.IPL2) GO TO 340
IPTI1 = IPT(IT1T3-1)
IF (IPTI1.NE.IPL1 .AND. IPTI1.NE.IPL2) GO TO 340
IPTI1 = IPT(IT1T3)
340 IT2T3 = ITF(2)*3
IPTI2 = IPT(IT2T3-2)
IF (IPTI2.NE.IPL1 .AND. IPTI2.NE.IPL2) GO TO 350
IPTI2 = IPT(IT2T3-1)
IF (IPTI2.NE.IPL1 .AND. IPTI2.NE.IPL2) GO TO 350
IPTI2 = IPT(IT2T3)
C
C - - CHECKS IF THE EXCHANGE IS NECESSARY.
C
350 IF (CONXCH(XD,YD,IPTI1,IPTI2,IPL1,IPL2) .EQ. 0)
1 GO TO 370
C
C - - MODIFIES THE IPT ARRAY WHEN NECESSARY.
C
IPT(IT1T3-2) = IPTI1
IPT(IT1T3-1) = IPTI2
IPT(IT1T3) = IPL1
IPT(IT2T3-2) = IPTI2
IPT(IT2T3-1) = IPTI1
IPT(IT2T3) = IPL2
C
C - - SETS NEW FLAGS.
C
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 360 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)
360 CONTINUE
370 CONTINUE
NLFC = NLF
NLF = JWL/2
IF (NLF .EQ. NLFC) GO TO 400
C
C - - RESETS THE IWL ARRAY FOR THE NEXT ROUND.
C
JWL = 0
JWL1MN = (NLFC+1)*2
NLFT2 = NLF*2
DO 380 JWL1=JWL1MN,NLFT2,2
JWL = JWL+2
IWL(JWL-1) = IWL(JWL1-1)
IWL(JWL) = IWL(JWL1)
380 CONTINUE
NLF = JWL/2
390 CONTINUE
400 CONTINUE
C
C REARRANGE THE IPT ARRAY SO THAT THE VERTEXES OF EACH TRIANGLE
C ARE LISTED COUNTER-CLOCKWISE.
C
DO 410 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)) .GE.
1 0.0) GO TO 410
IPT(ITT3-2) = IP2
IPT(ITT3-1) = IP1
410 CONTINUE
NT = NT0
NL = NL0
RETURN
END
|
