home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 3 / PDCD_3.iso / utilities / utilsp / pgplot / Examples / f77 / PGDemo3 < prev    next >
Encoding:
Text File  |  1994-02-24  |  14.4 KB  |  491 lines
  1.       PROGRAM PGDEM3
  2. C-----------------------------------------------------------------------
  3. C Demonstration program for PGPLOT contouring routines.
  4. C-----------------------------------------------------------------------
  5.       INTEGER PGBEG
  6.       WRITE (*,'(A)') ' Demonstration of PGPLOT contouring routines'
  7. C
  8. C Call PGBEG to initiate PGPLOT and open the output device; PGBEG
  9. C will prompt the user to supply the device name and type.
  10. C
  11.       IF (PGBEG(0,'?',1,1) .NE. 1) STOP
  12. C
  13. C Call the demonstration subroutines.
  14. C
  15.       WRITE (*,'(A)') ' Routine PGCONT'
  16.       CALL PGEX31
  17.       WRITE (*,'(A)') ' Routine PGCONS'
  18.       CALL PGEX32
  19.       WRITE (*,'(A)') ' Routine PGCONB'
  20.       CALL PGEX33
  21.       WRITE (*,'(A)') ' Routine PGCONX'
  22.       CALL PGEX34
  23.       WRITE (*,'(A)') ' Routine PGVECT'
  24.       CALL PGEX35
  25. C
  26. C Finally, call PGEND to terminate things properly.
  27. C
  28.       CALL PGEND
  29. C-----------------------------------------------------------------------
  30.       END
  31.  
  32.       SUBROUTINE PGEX31
  33. C-----------------------------------------------------------------------
  34. C Demonstration of contouring routine PGCONT.
  35. C-----------------------------------------------------------------------
  36.       INTEGER I,J
  37.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6)
  38.       DATA TR/0.,1.,0.,0.,0.,1./
  39. C
  40. C Compute a suitable function.
  41. C
  42.       FMIN = 0.0
  43.       FMAX = 0.0
  44.       DO 20 I=1,40
  45.           DO 10 J=1,40
  46.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  47.      1                     (I-J)/40.0
  48.               FMIN = MIN(F(I,J),FMIN)
  49.               FMAX = MAX(F(I,J),FMAX)
  50.    10     CONTINUE
  51.    20 CONTINUE
  52. C
  53. C Clear the screen. Set up window and viewport.
  54. C
  55.       CALL PGPAGE
  56.       CALL PGSVP(0.05,0.95,0.05,0.95)
  57.       CALL PGSWIN(1.0,40.0,1.0,40.0)
  58.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  59.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONT')
  60. C
  61. C Draw the map.  PGCONT is called once for each contour, using
  62. C different line attributes to distinguish contour levels.
  63. C
  64.       CALL PGBBUF
  65.       DO 30 I=1,21
  66.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  67.           IF (MOD(I,5).EQ.0) THEN
  68.               CALL PGSLW(3)
  69.           ELSE
  70.               CALL PGSLW(1)
  71.           END IF
  72.           IF (I.LT.10) THEN
  73.               CALL PGSCI(2)
  74.               CALL PGSLS(2)
  75.           ELSE
  76.               CALL PGSCI(3)
  77.               CALL PGSLS(1)
  78.           END IF
  79.           CALL PGCONT(F,40,40,1,40,1,40,ALEV,-1,TR)
  80.    30 CONTINUE
  81.       CALL PGSLW(1)
  82.       CALL PGSLS(1)
  83.       CALL PGSCI(1)
  84.       CALL PGEBUF
  85.       END
  86.  
  87.       SUBROUTINE PGEX32
  88. C-----------------------------------------------------------------------
  89. C Demonstration of contouring routine PGCONS.
  90. C-----------------------------------------------------------------------
  91.       INTEGER I,J
  92.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6)
  93.       DATA TR/0.,1.,0.,0.,0.,1./
  94. C
  95. C Compute a suitable function.
  96. C
  97.       FMIN = 0.0
  98.       FMAX = 0.0
  99.       DO 20 I=1,40
  100.           DO 10 J=1,40
  101.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  102.      1                     (I-J)/40.0
  103.               FMIN = MIN(F(I,J),FMIN)
  104.               FMAX = MAX(F(I,J),FMAX)
  105.    10     CONTINUE
  106.    20 CONTINUE
  107. C
  108. C Clear the screen. Set up window and viewport.
  109. C
  110.       CALL PGPAGE
  111.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  112.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONS')
  113. C
  114. C Draw the map.  PGCONS is called once for each contour, using
  115. C different line attributes to distinguish contour levels.
  116. C
  117.       CALL PGBBUF
  118.       DO 40 I=1,21
  119.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  120.           IF (MOD(I,5).EQ.0) THEN
  121.               CALL PGSLW(3)
  122.           ELSE
  123.               CALL PGSLW(1)
  124.           END IF
  125.           IF (I.LT.10) THEN
  126.               CALL PGSCI(2)
  127.               CALL PGSLS(2)
  128.           ELSE
  129.               CALL PGSCI(3)
  130.               CALL PGSLS(1)
  131.           END IF
  132.           CALL PGCONS(F,40,40,1,40,1,40,ALEV,-1,TR)
  133.    40 CONTINUE
  134.       CALL PGSLW(1)
  135.       CALL PGSLS(1)
  136.       CALL PGSCI(1)
  137.       CALL PGEBUF
  138.       END
  139.  
  140.       SUBROUTINE PGEX33
  141. C-----------------------------------------------------------------------
  142. C Demonstration of contouring routine PGCONB.
  143. C-----------------------------------------------------------------------
  144.       REAL BLANK
  145.       PARAMETER (BLANK=-1.2E20)
  146.       INTEGER I,J
  147.       REAL F(40,40),FMIN,FMAX,ALEV,TR(6),X,Y,R
  148.       DATA TR/0.,1.,0.,0.,0.,1./
  149. C
  150. C Compute a suitable function.
  151. C
  152.       FMIN = 0.0
  153.       FMAX = 0.0
  154.       DO 20 I=1,40
  155.           DO 10 J=1,40
  156.               F(I,J) = COS(0.3*SQRT(I*2.)-0.4*J/3.)*COS(0.4*I/3)+
  157.      1                     (I-J)/40.0
  158.               FMIN = MIN(F(I,J),FMIN)
  159.               FMAX = MAX(F(I,J),FMAX)
  160.    10     CONTINUE
  161.    20 CONTINUE
  162. C
  163. C "Blank" the data outside an annulus.
  164. C
  165.       DO 60 I=1,40
  166.           DO 50 J=1,40
  167.               R = SQRT((I-20.5)**2 + (J-20.5)**2)
  168.               IF (R.GT.20.0 .OR. R.LT.3.0) F(I,J) = BLANK
  169.    50     CONTINUE
  170.    60 CONTINUE
  171. C
  172.       CALL PGPAGE
  173.       CALL PGBOX('bcts',0.0,0,'bcts',0.0,0)
  174.       CALL PGMTXT('t',1.0,0.0,0.0,'Contouring using PGCONB')
  175.       CALL PGBBUF
  176.       DO 80 I=1,21
  177.           ALEV = FMIN + (I-1)*(FMAX-FMIN)/20.0
  178.           IF (MOD(I,5).EQ.0) THEN
  179.               CALL PGSLW(3)
  180.           ELSE
  181.               CALL PGSLW(1)
  182.           END IF
  183.           IF (I.LT.10) THEN
  184.               CALL PGSCI(2)
  185.               CALL PGSLS(2)
  186.           ELSE
  187.               CALL PGSCI(3)
  188.               CALL PGSLS(1)
  189.           END IF
  190.           CALL PGCONB(F,40,40,1,40,1,40,ALEV,-1,TR,BLANK)
  191.    80 CONTINUE
  192.       CALL PGEBUF
  193. C
  194. C Mark the blanked points for easy identification.
  195. C
  196.       CALL PGBBUF
  197.       CALL PGSCI(1)
  198.       DO 100 I=1,40
  199.           DO 90 J=1,40
  200.               IF (F(I,J).EQ.BLANK) THEN
  201.                   X = TR(1) + REAL(I)*TR(2) + REAL(J)*TR(3)
  202.                   Y = TR(4) + REAL(I)*TR(5) + REAL(J)*TR(6)
  203.                   CALL PGPT(1, X, Y, -1)
  204.               END IF
  205.    90     CONTINUE
  206.   100 CONTINUE
  207.       CALL PGEBUF
  208.       END
  209.  
  210.       SUBROUTINE PGEX34
  211. C-----------------------------------------------------------------------
  212. C This program is intended to demonstrate the use of the PGPLOT routine
  213. C PGCONX. As an example, we take data defined on a sphere. We want to
  214. C draw a contour map of the data on an equal-area projection of the
  215. C surface of the sphere; we choose the Hammer-Aitoff equal-area
  216. C projection centered on Declination (latitude) 0, Right Ascension
  217. C (longitude) 0. The data are defined at 1-degree intervals in both
  218. C coordinates. We thus need a data array dimensioned 361 by 181; the
  219. C array index runs from -90 to +90 in declination (181 elements) and
  220. C from -180 to +180 in right ascension (361 elements). The data at -180 
  221. C and +180 must be identical, of course, but they need to be duplicated 
  222. C in the array as these two longitudes appear on opposite sides of the
  223. C map. 
  224. C-----------------------------------------------------------------------
  225.       REAL  RPDEG
  226.       PARAMETER (RPDEG=3.1415926/180.0)
  227.       INTEGER I, J
  228.       REAL RA, DEC, B, L, XC(361), YC(361)
  229.       REAL Q(361,181), C(9)
  230.       EXTERNAL PLOT
  231. C
  232. C Call PGENV to create a rectangular window of 4 x 2 units. This is 
  233. C the bounding rectangle of the plot. The JUST argument is 1
  234. C to get equal scales in x and y.  
  235. C
  236.       CALL PGBBUF
  237.       CALL PGENV(-2.0, 2.0, -1.0, 1.0, 1, -2)
  238.       CALL PGLAB('Right Ascension', 'Declination', ' ')
  239.       CALL PGMTXT('t',2.0,0.0,0.0,
  240.      1           'Contouring on a non-Cartesian grid using PGCONX')
  241.       CALL PGSCH(0.6)
  242.       CALL PGMTXT('b',8.0,0.0,0.0,
  243.      1            'Hammer-Aitoff Equal-Area Projection of the Sphere')
  244.       CALL PGSCH(1.0)
  245. C
  246. C Draw 7 lines of constant longitude at longitude 0, 60, 120, ..., 
  247. C 360 degrees. Each line is made up of 180 straight-line segments.
  248. C
  249.       DO 20 J=1,7
  250.           RA = (-180.+(J-1)*60.)*RPDEG
  251.           DO 10 I=1,181
  252.               DEC = (I-91)*RPDEG
  253.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  254.    10     CONTINUE
  255.           CALL PGLINE(181,XC,YC)
  256.    20 CONTINUE
  257. C
  258. C Draw 5 lines of constant latitude at latitudes -60, -30, 0, 30, 
  259. C 60 degrees. Each line is made up of 360 straight-line segments.
  260. C
  261.       DO 40 J=1,5
  262.           DEC = (-60.+(J-1)*30.)*RPDEG
  263.           DO 30 I=1,361
  264.               RA = (I-181)*RPDEG
  265.               CALL AITOFF(DEC,RA,XC(I),YC(I))
  266.    30     CONTINUE
  267.           CALL PGLINE(361,XC,YC)
  268.    40 CONTINUE
  269.       CALL PGEBUF
  270. C
  271. C Compute the data to be contoured. In practice the data might be read
  272. C in from an external file. In this example the data are computed: they
  273. C are the galactic latitudes of the points on the sphere. Thus the 
  274. C contours will be lines of constant galactic latitude.
  275. C
  276.       DO 60 J=1,181
  277.           DEC = (J-91)*RPDEG
  278.           DO 50 I=1,361
  279.               RA = (I-181)*RPDEG
  280.               CALL GALACT(RA, DEC, B,L)
  281.               Q(I,J) = B
  282.    50     CONTINUE
  283.    60 CONTINUE
  284. C
  285. C Draw the contour map using PGCONX. Contours at 0, 20, 40, 60, 80.
  286. C
  287.       DO 70 I=1,9
  288.           C(I) = -100.0 +I*20.0
  289.    70 CONTINUE
  290.       CALL PGBBUF
  291.       CALL PGSCI(2)
  292.       CALL PGCONX(Q, 361, 181, 1, 361, 1, 181, C, 9, PLOT)
  293.       CALL PGSCI(1)
  294.       CALL PGEBUF
  295.       END
  296.  
  297.       SUBROUTINE PLOT(VISBLE, X, Y, Z)
  298.       INTEGER VISBLE
  299.       REAL X,Y,Z
  300. C-----------------------------------------------------------------------
  301. C Plotting subroutine for PGCONX. This routine converts the input
  302. C coordinates (latitude and longitude) into the projected coordinates
  303. C (x and y), and moves or draws as requested by VISBLE.
  304. C-----------------------------------------------------------------------
  305.       REAL  RPDEG
  306.       PARAMETER (RPDEG=3.1415926/180.0)
  307.       REAL B, L, XWORLD, YWORLD
  308.       B = (Y-91.0)*RPDEG
  309.       L = (X-181.0)*RPDEG
  310.       CALL AITOFF(B, L, XWORLD, YWORLD)
  311.       IF (VISBLE.EQ.0) THEN
  312.           CALL PGMOVE(XWORLD, YWORLD)
  313.       ELSE
  314.           CALL PGDRAW(XWORLD, YWORLD)
  315.       END IF
  316.       END
  317.  
  318.       SUBROUTINE AITOFF(B,L,X,Y)
  319. C-----------------------------------------------------------------------
  320. C Hammer-Aitoff projection.
  321. C
  322. C       Input: latitude and longitude (B,L) in radians
  323. C       Output: cartesian (X,Y) in range +/-2, +/-1
  324. C-----------------------------------------------------------------------
  325.       REAL L,B,X,Y,L2,DEN
  326. C
  327.       L2 = L/2.0
  328.       DEN = SQRT(1.0+COS(B)*COS(L2))
  329.       X = 2.0*COS(B)*SIN(L2)/DEN
  330.       Y = SIN(B)/DEN
  331.       END
  332.  
  333.       SUBROUTINE GALACT(RA,DEC,GLAT,GLONG)
  334. C-----------------------------------------------------------------------
  335. C Convert 1950.0 equatorial coordinates (RA, DEC) to galactic
  336. C latitude and longitude (GLAT, GLONG).
  337. C
  338. C Arguments:
  339. C  RA, DEC (input): 1950.0 RA and Dec (radians).
  340. C  GLAT, GLONG (output): galactic latitude and longitude 
  341. C      (degrees).
  342. C
  343. C Reference: e.g., D. R. H. Johnson and D. R. Soderblom, A. J. v93, 
  344. C  p864 (1987).
  345. C-----------------------------------------------------------------------
  346.       REAL RA, RRA, DEC, RDEC, CDEC, R(3,3), E(3), G(3)
  347.       REAL RADDEG, GLAT, GLONG
  348.       INTEGER I, J
  349.       DATA R/-.066988740D0, .492728466D0,-.867600811D0,-.872755766D0,
  350.      $       -.450346958D0,-.188374601D0,-.483538915D0, .744584633D0,
  351.      $        .460199785D0/
  352.       DATA RADDEG/57.29577951D0/
  353. C-----------------------------------------------------------------------
  354.       RRA = RA
  355.       RDEC = DEC
  356.       CDEC = COS(RDEC)
  357.       E(1) = CDEC*COS(RRA)
  358.       E(2) = CDEC*SIN(RRA)
  359.       E(3) = SIN(RDEC)
  360.       DO 20 I=1,3
  361.           G(I) = 0.0
  362.           DO 10 J=1,3
  363.               G(I) = G(I) + E(J)*R(I,J)
  364.    10     CONTINUE
  365.    20 CONTINUE
  366.       GLAT  = ASIN(G(3))*RADDEG
  367.       GLONG = ATAN2(G(2),G(1))*RADDEG
  368.       IF (GLONG.LT.0.0) GLONG = GLONG+360.0
  369.       RETURN
  370. C-----------------------------------------------------------------------
  371.       END
  372.  
  373.       SUBROUTINE PGEX35
  374. C-----------------------------------------------------------------------
  375. C Program to demonstrate the use of PGVECT along with
  376. C PGCONB by illustrating the flow around a cylinder with circulation.
  377. C
  378. C          NX      total # of axial stations
  379. C          NY      total # of grid pts in y (or r) direction
  380. C-----------------------------------------------------------------------
  381.       INTEGER MAXX, MAXY
  382.       PARAMETER (MAXX=101,MAXY=201)
  383.       INTEGER NX, NY, I, J, NC
  384.       REAL PI, A, GAMMA, VINF, XMAX, XMIN, YMAX, YMIN, DX, DY
  385.       REAL CPMIN, R2, BLANK
  386.       REAL CP(MAXX,MAXY),X(MAXX),Y(MAXY),U(MAXX,MAXY),V(MAXX,MAXY),
  387.      1   PSI(MAXX,MAXY)
  388.       REAL TR(6),C(10)
  389.       PARAMETER(PI = 3.14159265)
  390.       DATA BLANK/-1.E10/
  391. C
  392. C compute the flow about a cylinder with circulation
  393. C
  394. C define various quantities
  395. C
  396. C pi
  397. C number of points in the x and y directions
  398.       NX = 31
  399.       NY = 31
  400. C cylinder radius
  401.       A = 1.
  402. C circulation strength
  403.       GAMMA = 2.
  404. C freestream velocity
  405.       VINF = 1.
  406. C max and min x and y
  407.       XMAX = 3.*A
  408.       XMIN = -3.*A
  409.       YMAX = 3.*A
  410.       YMIN = -3.*A
  411. C point spacing
  412.       DX = (XMAX-XMIN)/(NX-1)
  413.       DY = (YMAX-YMIN)/(NY-1)
  414. C compute the stream function, Cp, and u and v velocities
  415.       CPMIN =1.E10
  416.       DO 20 I=1,NX
  417.          X(I) = XMIN+DX*(I-1)
  418.          DO 10 J=1,NY
  419.             Y(J) = YMIN+DY*(J-1)
  420.             R2 = X(I)**2+Y(J)**2
  421.             IF (R2.GT.0.) THEN
  422.                PSI(I,J) = VINF*Y(J)*(1.-A**2/R2)
  423.      1            +GAMMA/(2.*PI)*0.5*ALOG(R2/A)
  424.                U(I,J) = VINF*(1.+A**2/R2-2.*A**2*X(I)**2/R2**2)
  425.      1            +GAMMA/(2.*PI)*Y(J)/R2
  426.                V(I,J) = VINF*X(I)*(-2.*A**2*Y(J)/R2**2)
  427.      1            +GAMMA/(2.*PI)*X(I)/R2
  428.                CP(I,J) = 1.-(U(I,J)**2+V(I,J)**2)/VINF**2
  429.             ELSE
  430.                PSI(I,J) = 0.
  431.                U(I,J) = 0.
  432.                V(I,J) = 0.
  433.                CP(I,J) = 0.
  434.             END IF
  435.             IF (R2.LT.A**2) THEN
  436.                U(I,J) = BLANK
  437.                V(I,J) = BLANK
  438.             ELSE
  439.                CPMIN = MIN(CPMIN,CP(I,J))
  440.             END IF
  441.    10    CONTINUE
  442.    20 CONTINUE
  443. C
  444. C grid to world transformation
  445. C
  446.       TR(1)=X(1)-DX
  447.       TR(2)=DX
  448.       TR(3)=0.0
  449.       TR(4)=Y(1)-DY
  450.       TR(5)=0.0
  451.       TR(6)=DY
  452. C
  453.       CALL PGENV (X(1),X(NX),Y(1),Y(NY),1,0)
  454.       CALL PGIDEN
  455.       CALL PGLAB ('X','Y','Flow About a Cylinder with Circulation')
  456. C
  457. C contour plot of the stream function (streamlines)
  458. C
  459.       NC=5
  460.       C(1)=1.
  461.       C(2)=.5
  462.       C(3)=0.
  463.       C(4)=-.5
  464.       C(5)=-1.
  465.       CALL PGCONT (PSI,MAXX,MAXY,1,NX,1,NY,C,NC,TR)
  466. C
  467. C draw cylinder
  468. C
  469.       CALL PGBBUF
  470.       CALL PGSCI (0)
  471.       CALL PGSFS (1)
  472.       CALL PGCIRC (0.,0.,A*1.1)
  473.       CALL PGSFS (2)
  474.       CALL PGSCI (14)
  475.       CALL PGCIRC (0.0, 0., A)
  476.       CALL PGSCI (1)
  477.       CALL PGEBUF
  478. C
  479. C vector plot
  480. C
  481.       CALL PGSAH (2, 45.0, 0.7)
  482.       CALL PGSCH (0.3)
  483.       CALL PGVECT (U,V,MAXX,MAXY,2,NX-1,2,NY-1,0.0,0,TR,-1.E10)
  484.       CALL PGSCH(1.0)
  485. C
  486. C finished
  487. C
  488.       RETURN
  489. C----------------------------------------------------------------------
  490.       END
  491.