-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathCONTOR.SRC
480 lines (480 loc) · 12.3 KB
/
CONTOR.SRC
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
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
PROGRAM CONTOR
C
C CONTOR CONTOURS A GRID OF DATA GENERATED BY GRID.
C
CHARACTER*70 LINE
CHARACTER*40 WGRD
CHARACTER*4 FGRD /'.grd'/
REAL*8 GRD,A,CR,SCAL
REAL*8 XRAN,YRAN,XX,YY,ZZ
C
DIMENSION GRD(400,400),CTR(19)
C
COMMON /GRID/ GRID(400,400)
COMMON /TRANS/ CR(2),SCAL
C
DATA CTR / 1.0D-3, 2.0D-3, 4.0D-3, 8.0D-3,
+ 2.0D-2, 4.0D-2, 8.0D-2,
+ 2.0D-1, 4.0D-1, 8.0D-1,
+ 2.0D+0, 4.0D+0, 8.0D+0,
+ 2.0D+1, 4.0D+1, 8.0D+1,
+ 2.0D+2, 4.0D+2, 8.0D+2/
DATA NCNTR /19/
DATA XMIN / 0.0/, XMAX /10.0/, YMIN / 0.0/, YMAX /10.0/
C
CALL MAKNAME(1,WGRD,ILEN,FGRD)
IF (ILEN .EQ. 0) STOP 'usage: contor grdfile'
C
OPEN (30,FILE=WGRD)
C
C INITIALIZE PLOTTER AND LINE TYPES
C
CALL PLOTS (53,0,6)
CALL DASHDF (0,0.,0.,0.)
C
READ (30,100) LINE
READ (30,100) LINE
READ (30,*) NX,NY,XRAN,YRAN
READ (30,*) IFUN,NAT
READ (30,*) XX,YY,ZZ
CR(1) = XRAN/2.0D0
CR(2) = XRAN/2.0D0
SCAL = 10.0D0/XRAN
DO 30 I = 1,NAT
READ (30,*) XX,YY,ZZ
CALL NUCLEI(XX,YY,ZZ)
30 CONTINUE
C
DO 10 I = 1,NX
READ (30,*) (GRD(J,I),J=1,NY)
10 CONTINUE
C
DO 20 I = 1,NX
DO 20 J = 1,NY
GRID(I,J) = SNGL(GRD(I,J))
20 CONTINUE
C
CALL CTPLOT (NX,NY,1,NX,1,NY,XMIN,XMAX,YMIN,YMAX,
+ NCNTR,CTR,IFUN,1)
C
CALL FLUSH (6)
CALL PLOT (0.,0.,999)
C
C FORMATS
C
100 FORMAT(A70)
END
SUBROUTINE CONTUR ( IV, JV, KX, NX, KY, NY, CLP)
C
INTEGER UX1, UY1, COM, REC
C
COMMON /GRID/ A(400,400)
COMMON /SCLDAT/ XMIN, YMIN, XINCR, YINCR, KC
COMMON /CRPTPM/ CLA,KKX,MX,KKY,MY,IZ,NR,NP,IXO,IYO,ISO
COMMON /TEMP/ X(1000), Y(1000), REC(1000)
C
C THE 3 ROUTINES, CONTUR, CURVE AND INTERP FORM A
C GENERAL CONTOUR MAPPING PACKAGE THAT SIMPLY REQUIRES
C A PLOTTING ROUTINE COUPLR(N) TO PLOT N POINTS
C WITH COORDINATES IN THE ARRAYS X AND Y.
C IF CONTOURS WITH MORE THAN 400 POINTS ARE TO BE
C DRAWN (THIS IS VERY LARGE), THE DIMENSION STATEMENTS
C MUST BE ADJUSTED.
C
C SEE IMPORTANT NOTES IN INTERP AND CURVE
C
C THIS PROGRAM IS BASED ON M.O. DAYHOFF'S PAPER
C 'A CONTOUR MAP PROGRAM FOR X-RAY CRYSTALLOGRAPHY'
C OCT. 1963 COMMUNICATIONS OF THE ACM. FIRST
C PROGRAMMED BY BRUCE LANGDON, PLASMA PHYSICS
C LAB, PRINCETON UNIVERSITY, NOV. 1966.
C
IZ = MAX(IV,JV)
KKX = KX
KKY = KY
MX = NX
MY = NY
CLA = CLP
NR = 0
C
C SCAN RIGHT EDGE, BOTTOM TO TOP.
C
DO 100 J = KY+1, MY
IF (A(MX,J-1) .LT. CLA .AND. A(MX,J) .GE. CLA)
1 CALL CURVE( IV, JV, MX, J, 7, IRET)
IF (IRET .EQ. 1) RETURN
100 CONTINUE
C
C SCAN TOP EDGE, RIGHT TO LEFT
C
DO 110 I = MX-1, KX, -1
IF (A(I+1,MY) .LT. CLA .AND. A(I,MY) .GE. CLA)
1 CALL CURVE( IV, JV, I, MY, 5, IRET)
IF (IRET .EQ. 1) RETURN
110 CONTINUE
C
C SCAN LEFT EDGE TOP TO BOTTOM
C
DO 120 J = MY-1, KY, -1
IF (A(KX,J+1) .LT. CLA .AND. A(KX,J) .GE. CLA)
1 CALL CURVE( IV, JV, KX, J, 3, IRET)
IF (IRET .EQ. 1) RETURN
120 CONTINUE
C
C SCAN BOTTOM EDGE AND INTERIOR POINTS
C
DO 130 J = KY, MY-1
DO 130 I = KX+1, MX
IF (A(I-1,J) .GE. CLA .OR. A(I,J) .LT. CLA) GOTO 130
IF (NR .EQ. 0) GOTO 140
COM = IZ * I + J
DO 150 ID = 1, NR
IF (REC(ID) .EQ. COM) GOTO 130
150 CONTINUE
140 CALL CURVE( IV, JV, I, J, 1, IRET)
IF (IRET .EQ. 1) RETURN
130 CONTINUE
RETURN
END
SUBROUTINE COUPLR(NN)
C
INTEGER REC
C
REAL*8 CR,SCAL
COMMON /TEMP/ XX(1000), YY(1000), REC(1000)
COMMON /SCLDAT/ XMIN, YMIN, XINCR, YINCR, KK
COMMON /TRANS/ CR(2),SCAL
DIMENSION AXX(1000),AYY(1000)
C
AXX(1)=XX(1)-CR(1)*SCAL+CR(1)
AYY(1)=YY(1)-CR(2)*SCAL+CR(2)
C
CALL PLOT( AXX(1), AYY(1), 3)
C
C CHOOSE LINE TYPE
C
IF (KK .EQ. 1) IPEN = 2
IF (KK .EQ. 2) IPEN = 4
IF (KK .EQ. 3) IPEN = 5
IF (KK .EQ. 4) IPEN = 6
C
DO 10 I = 2, NN
AXX(I)=XX(I)-CR(1)*SCAL+CR(1)
AYY(I)=YY(I)-CR(2)*SCAL+CR(2)
CALL PLOT( AXX(I), AYY(I), IPEN)
10 CONTINUE
C
RETURN
END
SUBROUTINE CTPLOT
1 ( IV, ! ROW DIMENSION OF MATRIX
1 JV, ! COLUMN DIMENSION OF MATRIX
1 KX, ! LOWEST ROW CONSIDERED
1 NX, ! HIGHEST ROW CONSIDERED
1 KY, ! LOWER COLUMN CONSIDERED
1 NY, ! HIGHEST COLUMN CONSIDERED
1 XXMIN,
1 XXMAX,
1 YYMIN,
1 YYMAX,
1 NCTR, ! NUMBER OF CONTOURS
1 CTR, ! CONTOUR VALUES
1 IFUN, ! FUNCTION
1 ICONN )
C
C THIS PROGRAM PRODUCES CONTOUR PLOTS FROM THE MATRIX A WHICH HAS
C DIMENSIONS IV AND JV USING THE CONTOURS GIVEN IN THE ARRAY CTR.
C
C KX, NX, KY, AND NY DELINEATE THE PORTION OF THE ARRAY TO BE MAPPED.
C
C XXMIN, XXMAX, YYMIN, AND YYMAX INDICATES THE RANGE OF X AND Y VALUES FOR
C THIS REGION OF THE ARRAY.
C
C CTR IS AN ARRAY HOLDING THE VALUES OF THE CONTOURS TO BE DRAWN.
C NCTR THE NUMBER OF CONTOURS AND THE ARRAY
C ICONN 1 IF NEGATIVE CONTOURS ARE DESIRED THAT ARE THE SAME AS THE POSITIVE
C CONTOURS (ONLY THE POSITIVE ONES NEED BE SPECIFIED IN CTR AND NCTR).
C ======================================================================-======
C
COMMON /SCLDAT/ XMIN, YMIN, XINCR, YINCR, KC
DIMENSION CTR(NCTR)
C
XMAX = XXMAX
YMAX = YYMAX
XMIN = XXMIN
YMIN = YYMIN
XINCR = (XMAX - XMIN)/ (NX - KX + 1)
YINCR = (YMAX - YMIN)/ (NY - KY + 1)
C
C CALL CONTUR FOR EACH CONTOUR.
C
DO 100 NCL = 1, NCTR
CL = CTR(NCL)
KC = 1
IF (IFUN .EQ. 2 ) KC = 2
IF (IFUN .EQ. 4 ) KC = 4
CALL CONTUR( IV, JV, KX, NX, KY, NY, CL)
C
C FOR ICONN = 1, DO AGAIN WITH -VE CONTOURS.
C
IF (ICONN .EQ. 1) THEN
C DO 110 NCL = 1, NCTR
CL = -CTR(NCL)
KC = 1
IF (IFUN .EQ. 3) KC = 2
CALL CONTUR( IV, JV, KX, NX, KY, NY, CL)
C110 CONTINUE
C
END IF
100 CONTINUE
RETURN
END
SUBROUTINE CURVE( IV, JV, IXP, IYP, ISP, IRT)
C
INTEGER REC, DX, DY
LOGICAL PC, LOPC, LPC
DIMENSION INX(8), INY(8)
C
COMMON /GRID/ A(400,400)
COMMON /CRPTPM/ CL, KKX,MX,KKY,MY, IZ,NR,NP, IXO,IYO,ISO
COMMON /LOGIC/ PC, LOPC
COMMON /TEMP/ X(1000), Y(1000), REC(1000)
C
DATA INX /-1, -1, 0, 3*1, 0, -1/
DATA INY /0, 3*1, 0, 3*-1/
C
C RECMAX SHOULD BE EQUAL TO THE DIMENSION OF X, Y, REC IN
C COMMON /TEMP/
C
RECMAX=1000
C
IS = ISP
IX = IXP
IY = IYP
IXO = IX
IYO = IY
ISO = IS
NP = 1
LOPC = .FALSE.
CALL INTERP( IX, IY, IS, IRET, IV, JV)
IF (IRET .EQ. 1 .OR. IRET .EQ. 2) CALL
1 EXIT1(' BAD RETURN FROM INTERP')
IF (IRET .EQ. 3) GO TO 3
C
C ROTATE
C
1 IS = MOD( IS, 8) + 1
DX = INX(IS)
DY = INY(IS)
C
C FIND PLOT POINT
C
2 IRET = 1
CALL INTERP( IX, IY, IS, IRET, IV, JV)
IF (IRET .EQ. 1) THEN
C
C DIAG FAIL
C
IS = MOD( IS, 8) + 1
LOPC = .TRUE.
IRET = 1
CALL INTERP( IX+INX(IS), IY+INY(IS), IS-3, IRET, IV, JV)
IF (IRET .EQ. 1 .OR. IRET .EQ. 2)
1 CALL EXIT1(' BAD RETURN FROM INTERP')
IF (IRET .EQ. 3) GO TO 10
IX = IX + DX
IY = IY + DY
IS = MOD( IS+3, 8) + 1
GO TO 2
ELSE IF (IRET .EQ. 2) THEN
C
C NON-DIAG FAIL
C
IX = IX + DX
IY = IY + DY
IS = MOD( IS+4, 8) + 1
LPC = PC
CALL INTERP( IX, IY, IS, IRET, IV, JV)
IF (IRET .EQ. 2) CALL EXIT1(' BAD RETURN FROM INTERP')
IF (IRET .EQ. 3) GO TO 10
IF (.NOT. (LPC .AND. PC)) GO TO 1
C
C PLOT POINT SWITCH
C
IF( NP .LE. 2) GO TO 1
TEM = X(NP-2)
X(NP-2) = X(NP-1)
X(NP-1) = TEM
TEM = Y(NP-2)
Y(NP-2) = Y(NP-1)
Y(NP-1) = TEM
GO TO 1
ENDIF
IF (IRET .EQ. 3) GO TO 10
IF (IS .NE. 1) GO TO 1
C
C RECORD A 'CENTER'
C
3 IF (NR .GE. RECMAX) THEN
C
C REC ARRAY OVERFLOW
C
WRITE(2,4)
4 FORMAT(' TOO MANY POINTS IN CONTOUR. PROGRAM MUST BE
+RE-DIMENSIONED')
CALL COUPLR(NP-1)
IRT = 1
IRET = 1
RETURN
ELSE
NR = NR + 1
REC(NR) = IZ * IX + IY
GO TO 1
ENDIF
10 CALL COUPLR(NP-1)
IRET = 0
RETURN
END
SUBROUTINE EXIT1( AMESS)
C
C TERMINATES PLOTTING IN CASE OF AN ABORT.
C
CHARACTER*(*) AMESS
WRITE(6,*) AMESS
STOP ' CONTOR STOPS '
END
SUBROUTINE INTERP( IX, IY, ISP, IRET, IV, JV)
C
INTEGER REC, DX, DY
LOGICAL PC, LOPC
DIMENSION INX(8), INY(8)
C
COMMON /GRID/ A(400,400)
COMMON /CRPTPM/ CL, KKX,MX,KKY,MY, IZ,NR,NP, IXO,IYO,ISO
COMMON /LOGIC/ PC, LOPC
COMMON /SCLDAT/ XMIN, YMIN, XINCR, YINCR, KK
COMMON /TEMP/ X(1000), Y(1000), REC(1000)
C
DATA INX /-1, -1, 0, 3*1, 0, -1/
DATA INY /0, 3*1, 0, 3*-1/
C
IRET = 0
IS = ISP
IF (IS .LT. 1) IS = IS + 8
DX = INX(IS)
DY = INY(IS)
FDX = FLOAT(DX) * XINCR
FDY = FLOAT(DY) * YINCR
FIX = FLOAT( IX - KKX) * XINCR + XMIN
FIY = FLOAT( IY - KKY) * YINCR + YMIN
IX1 = IX + DX
IY1 = IY + DY
IF (IX1 .LT. KKX .OR. IX1 .GT. MX .OR.
1 IY1 .LT. KKY .OR. IY1 .GT. MY) THEN
IRET = 3
RETURN
ENDIF
IF ( DX * DY .EQ. 0) THEN
C
C NON-DIAGONAL CASE
C CHECK FOR FAIL
C
IF (A(IX1,IY1) .GE. CL) THEN
IRET = 2
RETURN
ENDIF
IF (DX .EQ. 0) THEN
X(NP) = FIX
Y(NP) = (A(IX,IY) - CL) / (A(IX,IY) - A(IX,IY1)) * FDY + FIY
ELSE
Y(NP) = FIY
X(NP) = (A(IX,IY) - CL) / (A(IX,IY) - A(IX1,IY)) * FDX + FIX
ENDIF
PC = .FALSE.
ELSE
C
C DIAGONAL CASE
C
CP = (A(IX,IY) + A(IX1,IY) + A(IX,IY1) + A(IX1,IY1)) / 4.
IF (CP .GE. CL .OR. LOPC) THEN
C
C CONTOUR PASSES ON FAR SIDE OF CENTER POINT
C
IF (A(IX1,IY1) .GE. CL .AND. CP .GE. CL) THEN
IRET = 1
RETURN
ENDIF
V = (A(IX1,IY1) - CL) / (A(IX1,IY1) - CP) / 2.
X(NP) = (1. - V) * FDX + FIX
Y(NP) = (1. - V) * FDY + FIY
PC = .TRUE.
LOPC = .FALSE.
ELSE
C
C CONTOUR PASSES ON NEAR SIDE OF CENTER POINT
C
V = (A(IX,IY) - CL) / (A(IX,IY) - CP) / 2.
X(NP) = V * FDX + FIX
Y(NP) = V * FDY + FIY
PC = .FALSE.
ENDIF
ENDIF
NP = NP + 1
IF (IX .EQ. IXO .AND. IY .EQ. IYO .AND. IS .EQ. ISO) THEN
IRET = 3
RETURN
ENDIF
IF (NP .LE. 400) RETURN
C
C PLOT PART OF CURVE AND CONTINUE
C
CALL COUPLR(NP-1)
X(1) = X(NP-1)
Y(1) = Y(NP-1)
NP = 2
RETURN
END
SUBROUTINE MAKNAME(I,STRING,L,EXT)
CHARACTER*(*) STRING,EXT
INTEGER I,J
CALL GETARG(I,STRING)
J = LEN(STRING)
DO 10 N = 1,J
IF(STRING(N:N) .EQ. ' ') THEN
L = N - 1
STRING = STRING(1:L)//EXT
RETURN
ENDIF
10 CONTINUE
STOP ' FAILED TO MAKE A FILE NAME '
RETURN
END
SUBROUTINE NUCLEI(XN,YN,ZN)
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
REAL*8 CR,SCAL
COMMON /TRANS/ CR(2),SCAL
C
DATA SN1,SN2 /.06,.04/
C
XN = (XN*SCAL + CR(1))
YN = (YN*SCAL + CR(2))
C
CALL PLOT (SNGL(XN-SCAL*SN1),SNGL(YN),3)
C
IF (DABS(ZN) .GE. 1.0D-01) THEN
CALL PLOT (SNGL(XN-SCAL*SN2),SNGL(YN),2)
CALL PLOT (SNGL(XN+SCAL*SN2),SNGL(YN),3)
CALL PLOT (SNGL(XN+SCAL*SN1),SNGL(YN),2)
CALL PLOT (SNGL(XN),SNGL(YN-SCAL*SN1),3)
CALL PLOT (SNGL(XN),SNGL(YN-SCAL*SN2),2)
CALL PLOT (SNGL(XN),SNGL(YN+SCAL*SN2),3)
CALL PLOT (SNGL(XN),SNGL(YN+SCAL*SN1),2)
ELSE
CALL PLOT (SNGL(XN+SCAL*SN1),SNGL(YN),2)
CALL PLOT (SNGL(XN),SNGL(YN-SCAL*SN1),3)
CALL PLOT (SNGL(XN),SNGL(YN+SCAL*SN1),2)
END IF
RETURN
END