PROGRAM PGDEM3 C----------------------------------------------------------------------- C Demonstration program for PGPLOT contouring routines. C----------------------------------------------------------------------- INTEGER PGBEG WRITE (*,'(A)') ' Demonstration of PGPLOT contouring routines' C C Call PGBEG to initiate PGPLOT and open the output device; PGBEG C will prompt the user to supply the device name and type. C IF (PGBEG(0,'?',1,1) .NE. 1) STOP C C Call the demonstration subroutines. C WRITE (*,'(A)') ' Routine PGCONT' CALL PGEX31 WRITE (*,'(A)') ' Routine PGCONS' CALL PGEX32 WRITE (*,'(A)') ' Routine PGCONT with PGCONL labels' CALL PGEX36 WRITE (*,'(A)') ' Routine PGCONB' CALL PGEX33 WRITE (*,'(A)') ' Routine PGCONX with arrow labels' CALL PGEX37 WRITE (*,'(A)') ' Routine PGCONX' CALL PGEX34 WRITE (*,'(A)') ' Routine PGCONF' CALL PGEXX1 C C Finally, call PGEND to terminate things properly. C CALL PGEND C----------------------------------------------------------------------- END SUBROUTINE PGEX31 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONT. C----------------------------------------------------------------------- INTEGER I,J REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6) DATA TR/0.,1.,0.,0.,0.,1./ C C Compute a suitable function. C FMIN = 0.0 FMAX = 0.0 DO 20 I=1,40 DO 10 J=1,40 F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+ 1 (I-J)/40.0 FMIN = MIN(F(I,J),FMIN) FMAX = MAX(F(I,J),FMAX) 10 CONTINUE 20 CONTINUE C C Clear the screen. Set up window and viewport. C CALL PGPAGE CALL PGSVP(0.05,0.95,0.05,0.95) CALL PGSWIN(1.0,40.0,1.0,40.0) CALL PGBOX('bcts',0.0,0,'bcts',0.0,0) CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONT') C C Draw the map. PGCONT is called once for each contour, using C different line attributes to distinguish contour levels. C CALL PGBBUF DO 30 I=1,21 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 IF (MOD(I,5).EQ.0) THEN CALL PGSLW(3) ELSE CALL PGSLW(1) END IF IF (I.LT.10) THEN CALL PGSCI(2) CALL PGSLS(2) ELSE CALL PGSCI(3) CALL PGSLS(1) END IF CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR) 30 CONTINUE CALL PGSLW(1) CALL PGSLS(1) CALL PGSCI(1) CALL PGEBUF END SUBROUTINE PGEX32 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONS. C----------------------------------------------------------------------- INTEGER I,J REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6) DATA TR/0.,1.,0.,0.,0.,1./ C C Compute a suitable function. C FMIN = 0.0 FMAX = 0.0 DO 20 I=1,40 DO 10 J=1,40 F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+ 1 (I-J)/40.0 FMIN = MIN(F(I,J),FMIN) FMAX = MAX(F(I,J),FMAX) 10 CONTINUE 20 CONTINUE C C Clear the screen. Set up window and viewport. C CALL PGPAGE CALL PGBOX('bcts',0.0,0,'bcts',0.0,0) CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONS') C C Draw the map. PGCONS is called once for each contour, using C different line attributes to distinguish contour levels. C CALL PGBBUF DO 40 I=1,21 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 IF (MOD(I,5).EQ.0) THEN CALL PGSLW(3) ELSE CALL PGSLW(1) END IF IF (I.LT.10) THEN CALL PGSCI(2) CALL PGSLS(2) ELSE CALL PGSCI(3) CALL PGSLS(1) END IF CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR) 40 CONTINUE CALL PGSLW(1) CALL PGSLS(1) CALL PGSCI(1) CALL PGEBUF END SUBROUTINE PGEX33 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONB. C----------------------------------------------------------------------- REAL BLANK PARAMETER (BLANK=-1.2E20) INTEGER I,J REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6),X,Y,R DATA TR/0.,1.,0.,0.,0.,1./ C C Compute a suitable function. C FMIN = 0.0 FMAX = 0.0 DO 20 I=1,40 DO 10 J=1,40 F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+ 1 (I-J)/40.0 FMIN = MIN(F(I,J),FMIN) FMAX = MAX(F(I,J),FMAX) 10 CONTINUE 20 CONTINUE C C "Blank" the data outside an annulus. C DO 60 I=1,40 DO 50 J=1,40 R = SQRT((I-20.5)**2 + (J-20.5)**2) IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK 50 CONTINUE 60 CONTINUE C CALL PGPAGE CALL PGBOX('bcts',0.0,0,'bcts',0.0,0) CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONB') CALL PGBBUF DO 80 I=1,21 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 IF (MOD(I,5).EQ.0) THEN CALL PGSLW(3) ELSE CALL PGSLW(1) END IF IF (I.LT.10) THEN CALL PGSCI(2) CALL PGSLS(2) ELSE CALL PGSCI(3) CALL PGSLS(1) END IF CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK) 80 CONTINUE CALL PGEBUF C C Mark the blanked points for easy identification. C CALL PGBBUF CALL PGSCI(1) DO 100 I=1,40 DO 90 J=1,40 IF (F(I,J).EQ.BLANK) THEN X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3) Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6) CALL PGPT1(X, Y, -1) END IF 90 CONTINUE 100 CONTINUE CALL PGEBUF END SUBROUTINE PGEX34 C----------------------------------------------------------------------- C This program is intended to demonstrate the use of the PGPLOT routine C PGCONX. As an example, we take data defined on a sphere. We want to C draw a contour map of the data on an equal-area projection of the C surface of the sphere; we choose the Hammer-Aitoff equal-area C projection centered on Declination (latitude) 0, Right Ascension C (longitude) 0. The data are defined at 2-degree intervals in both C coordinates. We thus need a data array dimensioned 181 by 91; the C array index runs from -90 to +90 in declination (91 elements) and C from -180 to +180 in right ascension (181 elements). The data at -180 C and +180 must be identical, of course, but they need to be duplicated C in the array as these two longitudes appear on opposite sides of the C map. C----------------------------------------------------------------------- REAL PI, RPDEG PARAMETER (PI=3.14159265359) PARAMETER (RPDEG=PI/180.0) INTEGER I, J REAL RA, DEC, B, L, XC(181), YC(181) REAL Q(181,91), C(9) EXTERNAL PLCAIT C C Call PGENV to create a rectangular window of 4 x 2 units. This is C the bounding rectangle of the plot. The JUST argument is 1 C to get equal scales in x and y. C CALL PGBBUF CALL PGENV(-2.0, 2.0, -1.0, 1.0, 1, -2) CALL PGLAB('Right Ascension', 'Declination', ' ') CALL PGMTXT('t',2.0,0.0,0.0, 1 'Contouring on a non-Cartesian grid using PGCONX') CALL PGSCH(0.6) CALL PGMTXT('b',8.0,0.0,0.0, 1 'Hammer-Aitoff Equal-Area Projection of the Sphere') CALL PGSCH(1.0) C C Draw 7 lines of constant longitude at longitude 0, 60, 120, ..., C 360 degrees. Each line is made up of 90 straight-line segments. C DO 20 J=1,7 RA = (-180.+(J-1)*60.)*RPDEG DO 10 I=1,91 DEC = 2*(I-46)*RPDEG CALL AITOFF(DEC,RA,XC(I),YC(I)) 10 CONTINUE CALL PGLINE(91,XC,YC) 20 CONTINUE C C Draw 5 lines of constant latitude at latitudes -60, -30, 0, 30, C 60 degrees. Each line is made up of 360 straight-line segments. C DO 40 J=1,5 DEC = (-60.+(J-1)*30.)*RPDEG DO 30 I=1,181 RA = 2*(I-91)*RPDEG CALL AITOFF(DEC,RA,XC(I),YC(I)) 30 CONTINUE CALL PGLINE(181,XC,YC) 40 CONTINUE CALL PGEBUF C C Compute the data to be contoured. In practice the data might be read C in from an external file. In this example the data are computed: they C are the galactic latitudes of the points on the sphere. Thus the C contours will be lines of constant galactic latitude. C DO 60 J=1,91 DEC = 2*(J-46)*RPDEG DO 50 I=1,181 RA = 2*(I-91)*RPDEG CALL GALACT(RA, DEC, B,L) Q(I,J) = B 50 CONTINUE 60 CONTINUE C C Draw the contour map using PGCONX. Contours at 0, 20, 40, 60, 80. C DO 70 I=1,9 C(I) = -100.0 +I*20.0 70 CONTINUE CALL PGBBUF CALL PGSCI(2) CALL PGCONX(Q, 181, 91, 1, 181, 1, 91, C, 9, PLCAIT) CALL PGSCI(1) CALL PGEBUF END SUBROUTINE PLCAIT(VISBLE, X, Y, Z) INTEGER VISBLE REAL X,Y,Z C----------------------------------------------------------------------- C Plotting subroutine for PGCONX. This routine converts the input C coordinates (latitude and longitude) into the projected coordinates C (x and y), and moves or draws as requested by VISBLE. C----------------------------------------------------------------------- REAL PI, RPDEG PARAMETER (PI=3.14159265359) PARAMETER (RPDEG=PI/180.0) REAL B, L, XWORLD, YWORLD B = 2.0*(Y-46.0)*RPDEG L = 2.0*(X-91.0)*RPDEG CALL AITOFF(B, L, XWORLD, YWORLD) IF (VISBLE.EQ.0) THEN CALL PGMOVE(XWORLD, YWORLD) ELSE CALL PGDRAW(XWORLD, YWORLD) END IF END SUBROUTINE AITOFF(B,L,X,Y) C----------------------------------------------------------------------- C Hammer-Aitoff projection. C C Input: latitude and longitude (B,L) in radians C Output: cartesian (X,Y) in range +/-2, +/-1 C----------------------------------------------------------------------- REAL L,B,X,Y,L2,DEN C L2 = L/2.0 DEN = SQRT(1.0+COS(B)*COS(L2)) X = 2.0*COS(B)*SIN(L2)/DEN Y = SIN(B)/DEN END SUBROUTINE GALACT(RA,DEC,GLAT,GLONG) C----------------------------------------------------------------------- C Convert 1950.0 equatorial coordinates (RA, DEC) to galactic C latitude and longitude (GLAT, GLONG). C C Arguments: C RA, DEC (input): 1950.0 RA and Dec (radians). C GLAT, GLONG (output): galactic latitude and longitude C (degrees). C C Reference: e.g., D. R. H. Johnson and D. R. Soderblom, A. J. v93, C p864 (1987). C----------------------------------------------------------------------- REAL RA, RRA, DEC, RDEC, CDEC, R(3,3), E(3), G(3) REAL RADDEG, GLAT, GLONG INTEGER I, J DATA R/-.066988740D0, .492728466D0,-.867600811D0,-.872755766D0, $ -.450346958D0,-.188374601D0,-.483538915D0, .744584633D0, $ .460199785D0/ DATA RADDEG/57.29577951D0/ C----------------------------------------------------------------------- RRA = RA RDEC = DEC CDEC = COS(RDEC) E(1) = CDEC*COS(RRA) E(2) = CDEC*SIN(RRA) E(3) = SIN(RDEC) DO 20 I=1,3 G(I) = 0.0 DO 10 J=1,3 G(I) = G(I) + E(J)*R(I,J) 10 CONTINUE 20 CONTINUE GLAT = ASIN(G(3))*RADDEG GLONG = ATAN2(G(2),G(1))*RADDEG IF (GLONG.LT.0.0) GLONG = GLONG+360.0 RETURN C----------------------------------------------------------------------- END SUBROUTINE PGEX37 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONX. C----------------------------------------------------------------------- INTEGER I,J REAL F(40,40),FMIN,FMAX,ALEV(1) EXTERNAL PLCARO C C Compute a suitable function. C FMIN = 0.0 FMAX = 0.0 DO 20 I=1,40 DO 10 J=1,40 F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+ 1 (I-J)/40.0 FMIN = MIN(F(I,J),FMIN) FMAX = MAX(F(I,J),FMAX) 10 CONTINUE 20 CONTINUE C C Clear the screen. Set up window and viewport. C CALL PGPAGE CALL PGSVP(0.05,0.95,0.05,0.95) CALL PGSWIN(1.0,40.0,1.0,40.0) CALL PGBOX('bcts',0.0,0,'bcts',0.0,0) CALL PGMTXT('t',1.0,0.0,0.0, : 'Contouring using PGCONX with arrows') C C Draw the map. PGCONX is called once for each contour, using C different line attributes to distinguish contour levels. C CALL PGBBUF DO 30 I=1,21 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 IF (MOD(I,5).EQ.0) THEN CALL PGSLW(3) ELSE CALL PGSLW(1) END IF IF (I.LT.10) THEN CALL PGSCI(2) CALL PGSLS(2) ELSE CALL PGSCI(3) CALL PGSLS(1) END IF CALL PGCONX(F,40,40,1,40,1,40,ALEV,-1,PLCARO) 30 CONTINUE CALL PGSLW(1) CALL PGSLS(1) CALL PGSCI(1) CALL PGEBUF END SUBROUTINE PGEX36 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONT and PGCONL. C----------------------------------------------------------------------- INTEGER I,J REAL F(40,40),FMIN,FMAX,ALEV(1),TR(6) CHARACTER*32 LABEL DATA TR /0.0, 1.0, 0.0, 0.0, 0.0, 1.0/ C C Compute a suitable function. C FMIN = 0.0 FMAX = 0.0 DO 20 I=1,40 DO 10 J=1,40 F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+ 1 (I-J)/40.0 FMIN = MIN(F(I,J),FMIN) FMAX = MAX(F(I,J),FMAX) 10 CONTINUE 20 CONTINUE C C Clear the screen. Set up window and viewport. C CALL PGPAGE CALL PGBOX('bcts',0.0,0,'bcts',0.0,0) CALL PGMTXT('t',1.0,0.0,0.0, 1 'Contouring using PGCONT and PGCONL labels') C C Draw the map. PGCONT is called once for each contour, using C different line attributes to distinguish contour levels. C CALL PGBBUF DO 40 I=1,21 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 IF (MOD(I,5).EQ.0) THEN CALL PGSLW(3) ELSE CALL PGSLW(1) END IF IF (I.LT.10) THEN CALL PGSCI(2) CALL PGSLS(2) ELSE CALL PGSCI(3) CALL PGSLS(1) END IF CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR) 40 CONTINUE CALL PGSLW(1) CALL PGSLS(1) CALL PGEBUF C C Label the contours with PGCONL. Only even-numbered contours C are labelled. C CALL PGBBUF DO 50 I=2,21,2 ALEV(1) = FMIN + (I-1)*(FMAX-FMIN)/20.0 WRITE (LABEL,'(I2)') I C WRITE (LABEL,'(F8.2)') ALEV IF (I.LT.10) THEN CALL PGSCI(2) ELSE CALL PGSCI(3) END IF CALL PGCONL(F,40,40,1,40,1,40,ALEV,TR,LABEL,16,8) 50 CONTINUE CALL PGSCI(1) CALL PGEBUF END SUBROUTINE PLCARO(VISBLE, X, Y, Z) INTEGER VISBLE REAL X,Y,Z C----------------------------------------------------------------------- C Plotting subroutine for PGCONX. This routine labels contours with C arrows. Arrows point clockwise around minima, anticlockwise around C maxima. Arrows are drawn on 1/16 of contour line segments. C----------------------------------------------------------------------- REAL XP, YP INTEGER I SAVE I DATA I /0/ C I = MOD(I+1,16) IF (VISBLE.EQ.0) THEN I = 0 CALL PGMOVE(X, Y) ELSE IF (I.EQ.8) THEN C -- Draw line segment with arrow at midpoint CALL PGQPOS(XP,YP) CALL PGARRO(XP, YP, (X+XP)*0.5, (Y+YP)*0.5) CALL PGDRAW(X, Y) ELSE C -- Draw plain line segment CALL PGDRAW(X, Y) END IF END SUBROUTINE PGEXX1 C----------------------------------------------------------------------- C Demonstration of contouring routine PGCONF. C----------------------------------------------------------------------- INTEGER NX, NY, NC PARAMETER (NX=51, NY=51, NC=9) INTEGER I, J REAL Z(NX,NY),TR(6), R REAL X, Y, XMIN, XMAX, YMIN, YMAX, DX, DY, MU, C(NC) DATA C /3.0, 3.2, 3.5, 3.6, 3.766413, 4.0 ,5.0, 10.0, 100.0/ C C Compute a suitable function. This is the test function used by C W. V. Snyder, Algorithm 531, Contour Plotting, ACM Trans. Math. C Softw. v.4, pp.290-294 (1978). C XMIN = -2.0 XMAX = 2.0 YMIN =-2.0 YMAX = 2.0 MU = 0.3 DX = (XMAX-XMIN)/FLOAT(NX-1) DY = (YMAX-YMIN)/FLOAT(NY-1) TR(1) = XMIN - DX TR(2) = DX TR(3) = 0.0 TR(4) = YMIN - DY TR(5) = 0.0 TR(6) = DY DO 20 I=1,NX X = TR(1) + I*TR(2) DO 10 J=1,NY Y = TR(4) + J*TR(6) Z(I,J) = (1.0-MU)*(2.0/SQRT((X-MU)**2+Y**2)+(X-MU)**2+Y**2) * + MU*(2.0/SQRT((X+1.0-MU)**2+Y**2)+(X+1.0-MU)**2+Y**2) 10 CONTINUE 20 CONTINUE C C Clear the screen. Set up window and viewport. C CALL PGPAGE CALL PGVSTD(0.05,0.95,0.05,0.95) CALL PGWNAD(XMIN, XMAX, YMIN, YMAX) C C Fill contours with PGCONF. C CALL PGSFS(1) DO 30 I=1, NC-1 R = 0.5+0.5*REAL(I-1)/REAL(NC-1) CALL PGSCR(I+10, R, R, R) CALL PGSCI(I+10) CALL PGCONF(Z,NX,NY,1,NX,1,NY,C(I),C(I+1),TR) 30 CONTINUE C C Draw the contour lines with PGCONT. C CALL PGSCI(3) CALL PGCONT(Z,NX,NY,1,NX,1,NY,C,NC,TR) C C Labels and box. C CALL PGSCI(1) CALL PGSCH(0.6) CALL PGBOX('bctsin',1.0,10,'bctsinv',1.0,10) CALL PGSCH(1.0) CALL PGMTXT('t',1.0,0.0,0.0,'Contour filling using PGCONF') C END