1      PROGRAM PGDEM3
2C-----------------------------------------------------------------------
3C Demonstration program for PGPLOT contouring routines.
4C-----------------------------------------------------------------------
5      INTEGER PGBEG
6      WRITE (*,'(A)') ' Demonstration of PGPLOT contouring routines'
7C
8C Call PGBEG to initiate PGPLOT and open the output device; PGBEG
9C will prompt the user to supply the device name and type.
10C
11      IF (PGBEG(0,'?',1,1) .NE. 1) STOP
12C
13C Call the demonstration subroutines.
14C
15      WRITE (*,'(A)') ' Routine PGCONT'
16      CALL PGEX31
17      WRITE (*,'(A)') ' Routine PGCONS'
18      CALL PGEX32
19      WRITE (*,'(A)') ' Routine PGCONT with PGCONL labels'
20      CALL PGEX36
21      WRITE (*,'(A)') ' Routine PGCONB'
22      CALL PGEX33
23      WRITE (*,'(A)') ' Routine PGCONX with arrow labels'
24      CALL PGEX37
25      WRITE (*,'(A)') ' Routine PGCONX'
26      CALL PGEX34
27      WRITE (*,'(A)') ' Routine PGCONF'
28      CALL PGEXX1
29C
30C Finally, call PGEND to terminate things properly.
31C
32      CALL PGEND
33C-----------------------------------------------------------------------
34      END
35
36      SUBROUTINE PGEX31
37C-----------------------------------------------------------------------
38C Demonstration of contouring routine PGCONT.
39C-----------------------------------------------------------------------
40      INTEGER I,J
41      REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
42      DATA TR/0.,1.,0.,0.,0.,1./
43C
44C Compute a suitable function.
45C
46      FMIN = 0.0
47      FMAX = 0.0
48      DO 20 I=1,40
49          DO 10 J=1,40
50              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
51     1                     (I-J)/40.0
52              FMIN = MIN(F(I,J),FMIN)
53              FMAX = MAX(F(I,J),FMAX)
54   10     CONTINUE
55   20 CONTINUE
56C
57C Clear the screen. Set up window and viewport.
58C
59      CALL PGPAGE
60      CALL PGSVP(0.05,0.95,0.05,0.95)
61      CALL PGSWIN(1.0,40.0,1.0,40.0)
62      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
63      CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONT')
64C
65C Draw the map.  PGCONT is called once for each contour, using
66C different line attributes to distinguish contour levels.
67C
68      CALL PGBBUF
69      DO 30 I=1,21
70          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
71          IF (MOD(I,5).EQ.0) THEN
72              CALL PGSLW(3)
73          ELSE
74              CALL PGSLW(1)
75          END IF
76          IF (I.LT.10) THEN
77              CALL PGSCI(2)
78              CALL PGSLS(2)
79          ELSE
80              CALL PGSCI(3)
81              CALL PGSLS(1)
82          END IF
83          CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
84   30 CONTINUE
85      CALL PGSLW(1)
86      CALL PGSLS(1)
87      CALL PGSCI(1)
88      CALL PGEBUF
89      END
90
91      SUBROUTINE PGEX32
92C-----------------------------------------------------------------------
93C Demonstration of contouring routine PGCONS.
94C-----------------------------------------------------------------------
95      INTEGER I,J
96      REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
97      DATA TR/0.,1.,0.,0.,0.,1./
98C
99C Compute a suitable function.
100C
101      FMIN = 0.0
102      FMAX = 0.0
103      DO 20 I=1,40
104          DO 10 J=1,40
105              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
106     1                     (I-J)/40.0
107              FMIN = MIN(F(I,J),FMIN)
108              FMAX = MAX(F(I,J),FMAX)
109   10     CONTINUE
110   20 CONTINUE
111C
112C Clear the screen. Set up window and viewport.
113C
114      CALL PGPAGE
115      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
116      CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONS')
117C
118C Draw the map.  PGCONS is called once for each contour, using
119C different line attributes to distinguish contour levels.
120C
121      CALL PGBBUF
122      DO 40 I=1,21
123          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
124          IF (MOD(I,5).EQ.0) THEN
125              CALL PGSLW(3)
126          ELSE
127              CALL PGSLW(1)
128          END IF
129          IF (I.LT.10) THEN
130              CALL PGSCI(2)
131              CALL PGSLS(2)
132          ELSE
133              CALL PGSCI(3)
134              CALL PGSLS(1)
135          END IF
136          CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR)
137   40 CONTINUE
138      CALL PGSLW(1)
139      CALL PGSLS(1)
140      CALL PGSCI(1)
141      CALL PGEBUF
142      END
143
144      SUBROUTINE PGEX33
145C-----------------------------------------------------------------------
146C Demonstration of contouring routine PGCONB.
147C-----------------------------------------------------------------------
148      REAL BLANK
149      PARAMETER (BLANK=-1.2E20)
150      INTEGER I,J
151      REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6),X,Y,R
152      DATA TR/0.,1.,0.,0.,0.,1./
153C
154C Compute a suitable function.
155C
156      FMIN = 0.0
157      FMAX = 0.0
158      DO 20 I=1,40
159          DO 10 J=1,40
160              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
161     1                     (I-J)/40.0
162              FMIN = MIN(F(I,J),FMIN)
163              FMAX = MAX(F(I,J),FMAX)
164   10     CONTINUE
165   20 CONTINUE
166C
167C "Blank" the data outside an annulus.
168C
169      DO 60 I=1,40
170          DO 50 J=1,40
171              R = SQRT((I-20.5)**2 + (J-20.5)**2)
172              IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK
173   50     CONTINUE
174   60 CONTINUE
175C
176      CALL PGPAGE
177      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
178      CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONB')
179      CALL PGBBUF
180      DO 80 I=1,21
181          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
182          IF (MOD(I,5).EQ.0) THEN
183              CALL PGSLW(3)
184          ELSE
185              CALL PGSLW(1)
186          END IF
187          IF (I.LT.10) THEN
188              CALL PGSCI(2)
189              CALL PGSLS(2)
190          ELSE
191              CALL PGSCI(3)
192              CALL PGSLS(1)
193          END IF
194          CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK)
195   80 CONTINUE
196      CALL PGEBUF
197C
198C Mark the blanked points for easy identification.
199C
200      CALL PGBBUF
201      CALL PGSCI(1)
202      DO 100 I=1,40
203          DO 90 J=1,40
204              IF (F(I,J).EQ.BLANK) THEN
205                  X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3)
206                  Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6)
207                  CALL PGPT1(X, Y, -1)
208              END IF
209   90     CONTINUE
210  100 CONTINUE
211      CALL PGEBUF
212      END
213
214      SUBROUTINE PGEX34
215C-----------------------------------------------------------------------
216C This program is intended to demonstrate the use of the PGPLOT routine
217C PGCONX. As an example, we take data defined on a sphere. We want to
218C draw a contour map of the data on an equal-area projection of the
219C surface of the sphere; we choose the Hammer-Aitoff equal-area
220C projection centered on Declination (latitude) 0, Right Ascension
221C (longitude) 0. The data are defined at 2-degree intervals in both
222C coordinates. We thus need a data array dimensioned 181 by 91; the
223C array index runs from -90 to +90 in declination (91 elements) and
224C from -180 to +180 in right ascension (181 elements). The data at -180
225C and +180 must be identical, of course, but they need to be duplicated
226C in the array as these two longitudes appear on opposite sides of the
227C map.
228C-----------------------------------------------------------------------
229      REAL  PI, RPDEG
230      PARAMETER (PI=3.14159265359)
231      PARAMETER (RPDEG=PI/180.0)
232      INTEGER I, J
233      REAL RA, DEC, B, L, XC(181), YC(181)
234      REAL Q(181,91), C(9)
235      EXTERNAL PLCAIT
236C
237C Call PGENV to create a rectangular window of 4 x 2 units. This is
238C the bounding rectangle of the plot. The JUST argument is 1
239C to get equal scales in x and y.
240C
241      CALL PGBBUF
242      CALL PGENV(-2.0, 2.0, -1.0, 1.0, 1, -2)
243      CALL PGLAB('Right Ascension', 'Declination', ' ')
244      CALL PGMTXT('t',2.0,0.0,0.0,
245     1           'Contouring on a non-Cartesian grid using PGCONX')
246      CALL PGSCH(0.6)
247      CALL PGMTXT('b',8.0,0.0,0.0,
248     1            'Hammer-Aitoff Equal-Area Projection of the Sphere')
249      CALL PGSCH(1.0)
250C
251C Draw 7 lines of constant longitude at longitude 0, 60, 120, ...,
252C 360 degrees. Each line is made up of 90 straight-line segments.
253C
254      DO 20 J=1,7
255          RA = (-180.+(J-1)*60.)*RPDEG
256          DO 10 I=1,91
257              DEC = 2*(I-46)*RPDEG
258              CALL AITOFF(DEC,RA,XC(I),YC(I))
259   10     CONTINUE
260          CALL PGLINE(91,XC,YC)
261   20 CONTINUE
262C
263C Draw 5 lines of constant latitude at latitudes -60, -30, 0, 30,
264C 60 degrees. Each line is made up of 360 straight-line segments.
265C
266      DO 40 J=1,5
267          DEC = (-60.+(J-1)*30.)*RPDEG
268          DO 30 I=1,181
269              RA = 2*(I-91)*RPDEG
270              CALL AITOFF(DEC,RA,XC(I),YC(I))
271   30     CONTINUE
272          CALL PGLINE(181,XC,YC)
273   40 CONTINUE
274      CALL PGEBUF
275C
276C Compute the data to be contoured. In practice the data might be read
277C in from an external file. In this example the data are computed: they
278C are the galactic latitudes of the points on the sphere. Thus the
279C contours will be lines of constant galactic latitude.
280C
281      DO 60 J=1,91
282          DEC = 2*(J-46)*RPDEG
283          DO 50 I=1,181
284              RA = 2*(I-91)*RPDEG
285              CALL GALACT(RA, DEC, B,L)
286              Q(I,J) = B
287   50     CONTINUE
288   60 CONTINUE
289C
290C Draw the contour map using PGCONX. Contours at 0, 20, 40, 60, 80.
291C
292      DO 70 I=1,9
293          C(I) = -100.0 +I*20.0
294   70 CONTINUE
295      CALL PGBBUF
296      CALL PGSCI(2)
297      CALL PGCONX(Q, 181, 91, 1, 181, 1, 91, C, 9, PLCAIT)
298      CALL PGSCI(1)
299      CALL PGEBUF
300      END
301
302      SUBROUTINE PLCAIT(VISBLE, X, Y, Z)
303      INTEGER VISBLE
304      REAL X,Y,Z
305C-----------------------------------------------------------------------
306C Plotting subroutine for PGCONX. This routine converts the input
307C coordinates (latitude and longitude) into the projected coordinates
308C (x and y), and moves or draws as requested by VISBLE.
309C-----------------------------------------------------------------------
310      REAL  PI, RPDEG
311      PARAMETER (PI=3.14159265359)
312      PARAMETER (RPDEG=PI/180.0)
313      REAL B, L, XWORLD, YWORLD
314      B = 2.0*(Y-46.0)*RPDEG
315      L = 2.0*(X-91.0)*RPDEG
316      CALL AITOFF(B, L, XWORLD, YWORLD)
317      IF (VISBLE.EQ.0) THEN
318          CALL PGMOVE(XWORLD, YWORLD)
319      ELSE
320          CALL PGDRAW(XWORLD, YWORLD)
321      END IF
322      END
323
324      SUBROUTINE AITOFF(B,L,X,Y)
325C-----------------------------------------------------------------------
326C Hammer-Aitoff projection.
327C
328C       Input: latitude and longitude (B,L) in radians
329C       Output: cartesian (X,Y) in range +/-2, +/-1
330C-----------------------------------------------------------------------
331      REAL L,B,X,Y,L2,DEN
332C
333      L2 = L/2.0
334      DEN = SQRT(1.0+COS(B)*COS(L2))
335      X = 2.0*COS(B)*SIN(L2)/DEN
336      Y = SIN(B)/DEN
337      END
338
339      SUBROUTINE GALACT(RA,DEC,GLAT,GLONG)
340C-----------------------------------------------------------------------
341C Convert 1950.0 equatorial coordinates (RA, DEC) to galactic
342C latitude and longitude (GLAT, GLONG).
343C
344C Arguments:
345C  RA, DEC (input): 1950.0 RA and Dec (radians).
346C  GLAT, GLONG (output): galactic latitude and longitude
347C      (degrees).
348C
349C Reference: e.g., D. R. H. Johnson and D. R. Soderblom, A. J. v93,
350C  p864 (1987).
351C-----------------------------------------------------------------------
352      REAL RA, RRA, DEC, RDEC, CDEC, R(3,3), E(3), G(3)
353      REAL RADDEG, GLAT, GLONG
354      INTEGER I, J
355      DATA R/-.066988740D0, .492728466D0,-.867600811D0,-.872755766D0,
356     $       -.450346958D0,-.188374601D0,-.483538915D0, .744584633D0,
357     $        .460199785D0/
358      DATA RADDEG/57.29577951D0/
359C-----------------------------------------------------------------------
360      RRA = RA
361      RDEC = DEC
362      CDEC = COS(RDEC)
363      E(1) = CDEC*COS(RRA)
364      E(2) = CDEC*SIN(RRA)
365      E(3) = SIN(RDEC)
366      DO 20 I=1,3
367          G(I) = 0.0
368          DO 10 J=1,3
369              G(I) = G(I) + E(J)*R(I,J)
370   10     CONTINUE
371   20 CONTINUE
372      GLAT  = ASIN(G(3))*RADDEG
373      GLONG = ATAN2(G(2),G(1))*RADDEG
374      IF (GLONG.LT.0.0) GLONG = GLONG+360.0
375      RETURN
376C-----------------------------------------------------------------------
377      END
378
379      SUBROUTINE PGEX37
380C-----------------------------------------------------------------------
381C Demonstration of contouring routine PGCONX.
382C-----------------------------------------------------------------------
383      INTEGER I,J
384      REAL F(40,40),FMIN,FMAX,ALEV(1)
385      EXTERNAL PLCARO
386C
387C Compute a suitable function.
388C
389      FMIN = 0.0
390      FMAX = 0.0
391      DO 20 I=1,40
392          DO 10 J=1,40
393              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
394     1                     (I-J)/40.0
395              FMIN = MIN(F(I,J),FMIN)
396              FMAX = MAX(F(I,J),FMAX)
397   10     CONTINUE
398   20 CONTINUE
399C
400C Clear the screen. Set up window and viewport.
401C
402      CALL PGPAGE
403      CALL PGSVP(0.05,0.95,0.05,0.95)
404      CALL PGSWIN(1.0,40.0,1.0,40.0)
405      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
406      CALL PGMTXT('t',1.0,0.0,0.0,
407     :            'Contouring using PGCONX with arrows')
408C
409C Draw the map.  PGCONX is called once for each contour, using
410C different line attributes to distinguish contour levels.
411C
412      CALL PGBBUF
413      DO 30 I=1,21
414          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
415          IF (MOD(I,5).EQ.0) THEN
416              CALL PGSLW(3)
417          ELSE
418              CALL PGSLW(1)
419          END IF
420          IF (I.LT.10) THEN
421              CALL PGSCI(2)
422              CALL PGSLS(2)
423          ELSE
424              CALL PGSCI(3)
425              CALL PGSLS(1)
426          END IF
427          CALL PGCONX(F,40,40,1,40,1,40,ALEV,-1,PLCARO)
428   30 CONTINUE
429      CALL PGSLW(1)
430      CALL PGSLS(1)
431      CALL PGSCI(1)
432      CALL PGEBUF
433      END
434
435      SUBROUTINE PGEX36
436C-----------------------------------------------------------------------
437C Demonstration of contouring routine PGCONT and PGCONL.
438C-----------------------------------------------------------------------
439      INTEGER I,J
440      REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6)
441      CHARACTER*32 LABEL
442      DATA TR /0.0, 1.0, 0.0, 0.0, 0.0, 1.0/
443C
444C Compute a suitable function.
445C
446      FMIN = 0.0
447      FMAX = 0.0
448      DO 20 I=1,40
449          DO 10 J=1,40
450              F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
451     1                     (I-J)/40.0
452              FMIN = MIN(F(I,J),FMIN)
453              FMAX = MAX(F(I,J),FMAX)
454   10     CONTINUE
455   20 CONTINUE
456C
457C Clear the screen. Set up window and viewport.
458C
459      CALL PGPAGE
460      CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
461      CALL PGMTXT('t',1.0,0.0,0.0,
462     1            'Contouring using PGCONT and PGCONL labels')
463C
464C Draw the map.  PGCONT is called once for each contour, using
465C different line attributes to distinguish contour levels.
466C
467      CALL PGBBUF
468      DO 40 I=1,21
469          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
470          IF (MOD(I,5).EQ.0) THEN
471              CALL PGSLW(3)
472          ELSE
473              CALL PGSLW(1)
474          END IF
475          IF (I.LT.10) THEN
476              CALL PGSCI(2)
477              CALL PGSLS(2)
478          ELSE
479              CALL PGSCI(3)
480              CALL PGSLS(1)
481          END IF
482          CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
483   40 CONTINUE
484      CALL PGSLW(1)
485      CALL PGSLS(1)
486      CALL PGEBUF
487C
488C Label the contours with PGCONL. Only even-numbered contours
489C are labelled.
490C
491      CALL PGBBUF
492      DO 50 I=2,21,2
493          ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0
494          WRITE (LABEL,'(I2)') I
495C         WRITE (LABEL,'(F8.2)') ALEV
496          IF (I.LT.10) THEN
497              CALL PGSCI(2)
498          ELSE
499              CALL PGSCI(3)
500          END IF
501          CALL PGCONL(F,40,40,1,40,1,40,ALEV,TR,LABEL,16,8)
502 50   CONTINUE
503      CALL PGSCI(1)
504      CALL PGEBUF
505      END
506
507      SUBROUTINE PLCARO(VISBLE, X, Y, Z)
508      INTEGER VISBLE
509      REAL X,Y,Z
510C-----------------------------------------------------------------------
511C Plotting subroutine for PGCONX. This routine labels contours with
512C arrows. Arrows point clockwise around minima, anticlockwise around
513C maxima. Arrows are drawn on 1/16 of contour line segments.
514C-----------------------------------------------------------------------
515      REAL XP, YP
516      INTEGER I
517      SAVE I
518      DATA I /0/
519C
520      I = MOD(I+1,16)
521      IF (VISBLE.EQ.0) THEN
522          I = 0
523          CALL PGMOVE(X, Y)
524      ELSE IF (I.EQ.8) THEN
525C         -- Draw line segment with arrow at midpoint
526          CALL PGQPOS(XP,YP)
527          CALL PGARRO(XP, YP, (X+XP)*0.5, (Y+YP)*0.5)
528          CALL PGDRAW(X, Y)
529      ELSE
530C         -- Draw plain line segment
531          CALL PGDRAW(X, Y)
532      END IF
533      END
534
535      SUBROUTINE PGEXX1
536C-----------------------------------------------------------------------
537C Demonstration of contouring routine PGCONF.
538C-----------------------------------------------------------------------
539      INTEGER NX, NY, NC
540      PARAMETER (NX=51, NY=51, NC=9)
541      INTEGER I, J
542      REAL Z(NX,NY),TR(6), R
543      REAL X, Y, XMIN, XMAX, YMIN, YMAX, DX, DY, MU, C(NC)
544      DATA C /3.0, 3.2, 3.5, 3.6, 3.766413, 4.0 ,5.0, 10.0, 100.0/
545C
546C Compute a suitable function. This is the test function used by
547C W. V. Snyder, Algorithm 531, Contour Plotting, ACM Trans. Math.
548C Softw. v.4, pp.290-294 (1978).
549C
550      XMIN = -2.0
551      XMAX = 2.0
552      YMIN =-2.0
553      YMAX = 2.0
554      MU = 0.3
555      DX = (XMAX-XMIN)/FLOAT(NX-1)
556      DY = (YMAX-YMIN)/FLOAT(NY-1)
557      TR(1) = XMIN - DX
558      TR(2) = DX
559      TR(3) = 0.0
560      TR(4) = YMIN - DY
561      TR(5) = 0.0
562      TR(6) = DY
563      DO 20 I=1,NX
564         X = TR(1) + I*TR(2)
565         DO 10 J=1,NY
566            Y = TR(4) + J*TR(6)
567            Z(I,J) = (1.0-MU)*(2.0/SQRT((X-MU)**2+Y**2)+(X-MU)**2+Y**2)
568     *           + MU*(2.0/SQRT((X+1.0-MU)**2+Y**2)+(X+1.0-MU)**2+Y**2)
569 10      CONTINUE
570   20 CONTINUE
571C
572C Clear the screen. Set up window and viewport.
573C
574      CALL PGPAGE
575      CALL PGVSTD(0.05,0.95,0.05,0.95)
576      CALL PGWNAD(XMIN, XMAX, YMIN, YMAX)
577C
578C Fill contours with PGCONF.
579C
580      CALL PGSFS(1)
581      DO 30 I=1, NC-1
582         R = 0.5+0.5*REAL(I-1)/REAL(NC-1)
583         CALL PGSCR(I+10, R, R, R)
584         CALL PGSCI(I+10)
585         CALL PGCONF(Z,NX,NY,1,NX,1,NY,C(I),C(I+1),TR)
586 30   CONTINUE
587C
588C Draw the contour lines with PGCONT.
589C
590      CALL PGSCI(3)
591      CALL PGCONT(Z,NX,NY,1,NX,1,NY,C,NC,TR)
592C
593C Labels and box.
594C
595      CALL PGSCI(1)
596      CALL PGSCH(0.6)
597      CALL PGBOX('bctsin',1.0,10,'bctsinv',1.0,10)
598      CALL PGSCH(1.0)
599      CALL PGMTXT('t',1.0,0.0,0.0,'Contour filling using PGCONF')
600C
601      END
602