home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / pgplot_1 / Examples / f77 / PGDEMO15 < prev    next >
Text File  |  1997-06-10  |  4KB  |  134 lines
  1.       PROGRAM PGDEM3
  2. C-----------------------------------------------------------------------
  3. C Demonstration program for PGPLOT vector field plot.
  4. C-----------------------------------------------------------------------
  5.       INTEGER PGOPEN
  6.       WRITE (*,'(A)') ' Demonstration of routine PGVECT'
  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 (PGOPEN('?') .LE. 0) STOP
  12.       CALL PGEX35
  13.       CALL PGEND
  14. C-----------------------------------------------------------------------
  15.       END
  16.  
  17.       SUBROUTINE PGEX35
  18. C-----------------------------------------------------------------------
  19. C Program to demonstrate the use of PGVECT along with
  20. C PGCONB by illustrating the flow around a cylinder with circulation.
  21. C
  22. C          NX      total # of axial stations
  23. C          NY      total # of grid pts in y (or r) direction
  24. C-----------------------------------------------------------------------
  25.       INTEGER MAXX, MAXY
  26.       PARAMETER (MAXX=101,MAXY=201)
  27.       INTEGER NX, NY, I, J, NC
  28.       REAL PI, A, GAMMA, VINF, XMAX, XMIN, YMAX, YMIN, DX, DY
  29.       REAL CPMIN, R2, BLANK
  30.       REAL CP(MAXX,MAXY),X(MAXX),Y(MAXY),U(MAXX,MAXY),V(MAXX,MAXY),
  31.      1   PSI(MAXX,MAXY)
  32.       REAL TR(6),C(10)
  33.       PARAMETER (PI=3.14159265359)
  34.       DATA BLANK/-1.E10/
  35. C
  36. C compute the flow about a cylinder with circulation
  37. C
  38. C define various quantities
  39. C
  40. C number of points in the x and y directions
  41.       NX = 31
  42.       NY = 31
  43. C cylinder radius
  44.       A = 1.
  45. C circulation strength
  46.       GAMMA = 2.
  47. C freestream velocity
  48.       VINF = 1.
  49. C max and min x and y
  50.       XMAX = 3.*A
  51.       XMIN = -3.*A
  52.       YMAX = 3.*A
  53.       YMIN = -3.*A
  54. C point spacing
  55.       DX = (XMAX-XMIN)/(NX-1)
  56.       DY = (YMAX-YMIN)/(NY-1)
  57. C compute the stream function, Cp, and u and v velocities
  58.       CPMIN =1.E10
  59.       DO 20 I=1,NX
  60.          X(I) = XMIN+DX*(I-1)
  61.          DO 10 J=1,NY
  62.             Y(J) = YMIN+DY*(J-1)
  63.             R2 = X(I)**2+Y(J)**2
  64.             IF (R2.GT.0.) THEN
  65.                PSI(I,J) = VINF*Y(J)*(1.-A**2/R2)
  66.      1            +GAMMA/(2.*PI)*0.5*ALOG(R2/A)
  67.                U(I,J) = VINF*(1.+A**2/R2-2.*A**2*X(I)**2/R2**2)
  68.      1            +GAMMA/(2.*PI)*Y(J)/R2
  69.                V(I,J) = VINF*X(I)*(-2.*A**2*Y(J)/R2**2)
  70.      1            +GAMMA/(2.*PI)*X(I)/R2
  71.                CP(I,J) = 1.-(U(I,J)**2+V(I,J)**2)/VINF**2
  72.             ELSE
  73.                PSI(I,J) = 0.
  74.                U(I,J) = 0.
  75.                V(I,J) = 0.
  76.                CP(I,J) = 0.
  77.             END IF
  78.             IF (R2.LT.A**2) THEN
  79.                U(I,J) = BLANK
  80.                V(I,J) = BLANK
  81.             ELSE
  82.                CPMIN = MIN(CPMIN,CP(I,J))
  83.             END IF
  84.    10    CONTINUE
  85.    20 CONTINUE
  86. C
  87. C grid to world transformation
  88. C
  89.       TR(1)=X(1)-DX
  90.       TR(2)=DX
  91.       TR(3)=0.0
  92.       TR(4)=Y(1)-DY
  93.       TR(5)=0.0
  94.       TR(6)=DY
  95. C
  96.       CALL PGENV (X(1),X(NX),Y(1),Y(NY),1,0)
  97.       CALL PGIDEN
  98.       CALL PGLAB ('X','Y','Flow About a Cylinder with Circulation')
  99. C
  100. C contour plot of the stream function (streamlines)
  101. C
  102.       NC=5
  103.       C(1)=1.
  104.       C(2)=.5
  105.       C(3)=0.
  106.       C(4)=-.5
  107.       C(5)=-1.
  108.       CALL PGCONT (PSI,MAXX,MAXY,1,NX,1,NY,C,NC,TR)
  109. C
  110. C draw cylinder
  111. C
  112.       CALL PGBBUF
  113.       CALL PGSCI (0)
  114.       CALL PGSFS (1)
  115.       CALL PGCIRC (0.,0.,A*1.1)
  116.       CALL PGSFS (2)
  117.       CALL PGSCI (14)
  118.       CALL PGCIRC (0.0, 0., A)
  119.       CALL PGSCI (1)
  120.       CALL PGEBUF
  121. C
  122. C vector plot
  123. C
  124.       CALL PGSAH (2, 45.0, 0.7)
  125.       CALL PGSCH (0.3)
  126.       CALL PGVECT (U,V,MAXX,MAXY,2,NX-1,2,NY-1,0.0,0,TR,-1.E10)
  127.       CALL PGSCH(1.0)
  128. C
  129. C finished
  130. C
  131.       RETURN
  132. C----------------------------------------------------------------------
  133.       END
  134.