home *** CD-ROM | disk | FTP | other *** search
/ ARM Club 3 / TheARMClub_PDCD3.iso / hensa / maths / plplot / plplot_2 / examples / f77 / x16f.f < prev   
Encoding:
Text File  |  1994-05-26  |  7.6 KB  |  300 lines
  1. ! $Id: x16f.f,v 1.2 1994/05/26 19:34:28 mjl Exp $
  2. ! $Log: x16f.f,v $
  3. ! Revision 1.2  1994/05/26  19:34:28  mjl
  4. ! Inserted missing CVS Id and Log fields for all Fortran demos.  Comment
  5. ! character changed to "!" everywhere, to work well with font-lock in Lucid
  6. ! emacs (requires a small change to fortran-mode.el).
  7. !
  8. !
  9. !     program example16
  10. !     =================
  11. !
  12. ! Demonstration of plshade plotting      
  13.  
  14. ! Reduce colors in cmap 0 so that cmap 1 is useful on a 16-color display
  15.  
  16.       call plscmap0n(3)
  17.  
  18. ! Initialize plplot
  19.  
  20.       call plinit()
  21.  
  22. ! Rectangular coordinate plot
  23.  
  24.       call rect()
  25.  
  26. ! Polar coordinate plot
  27.  
  28.       call polar()
  29.  
  30.       call plend
  31.       end
  32.  
  33. !----------------------------------------------------------------------------!
  34. ! Plot function using the identity transform
  35.  
  36.       subroutine rect()
  37.  
  38.       implicit character(a-z)
  39.       integer    NX, NY, NCONTR
  40.       parameter (NX = 35, NY = 46, NCONTR = 14)
  41.  
  42.       real    z(NX, NY), w(NX, NY), clevel(NCONTR)
  43.       real    xmin, xmax, ymin, ymax, zmin, zmax, x, y
  44.       real    shade_min, shade_max, sh_color
  45.       integer    i, j, sh_cmap, sh_width
  46.       integer    min_color, min_width, max_color, max_width
  47.  
  48.       xmin = -1.
  49.       ymin = -1.
  50.       xmax =  1.
  51.       ymax =  1.
  52.  
  53. ! Set up for plshade call
  54.  
  55.       sh_cmap = 1
  56.       min_color = 1
  57.       min_width = 0
  58.       max_color = 0
  59.       max_width = 0
  60.  
  61. ! Set up data arrays
  62.  
  63.       do 12 i = 1, NX
  64.          x = (i - 1 - (NX/2)) / real(NX/2)
  65.          do 10 j = 1, NY
  66.             y = (j - 1 - (NY/2)) / real(NY/2) - 1.0
  67.             z(i,j) = x*x - y*y + (x - y) / (x*x + y*y + 0.1)
  68.             w(i,j) = 2*x*y
  69.  10      continue
  70.  12   continue
  71.  
  72.       call a2mnmx(z, NX, NY, zmin, zmax)
  73.       do 20 i = 1, NCONTR
  74.          clevel(i) = zmin + (zmax - zmin) * (i + 0.5) / real(NCONTR)
  75.  20   continue
  76.  
  77. ! Plot using identity transform 
  78.  
  79.       call pladv(0)
  80.       call plvpor(0.1, 0.9, 0.1, 0.9)
  81.       call plwind(-1.0, 1.0, -1.0, 1.0)
  82.  
  83.       do 100 i = 1, NCONTR
  84.          shade_min = zmin + (zmax - zmin) * real(i - 1) / real(NCONTR)
  85.          shade_max = zmin + (zmax - zmin) * real(i)     / real(NCONTR)
  86.          sh_color = real(i - 1) / real(NCONTR - 1)
  87.          sh_width = 2
  88.          call plpsty(0)
  89.  
  90.          call plshade0(z, NX, NY, ' ', 
  91.      &        -1., 1., -1., 1., 
  92.      &        shade_min, shade_max, 
  93.      &        sh_cmap, sh_color, sh_width, 
  94.      &        min_color, min_width, max_color, max_width)
  95.  
  96.  100  continue
  97.  
  98.       call plcol(1)
  99.       call plbox('bcnst', 0.0, 0, 'bcnstv', 0.0, 0)
  100.       call plcol(2)
  101.       call pllab('distance', 'altitude', 'Bogon flux')
  102.  
  103.       return
  104.       end
  105.  
  106. !----------------------------------------------------------------------------!
  107. ! Routine for demonstrating use of transformation arrays in contour plots.
  108.  
  109.       subroutine polar()
  110.  
  111.       implicit character(a-z)
  112.       integer    NX, NY, NCONTR, NBDRY
  113.       real    TWOPI
  114.       parameter (NX = 40, NY = 64, NCONTR = 14, NBDRY=200)
  115.       parameter (TWOPI=6.2831853071795864768)
  116.  
  117.       real    z(NX, NY), ztmp(NX, NY+1)
  118.       real    xg(NX, NY+1), yg(NX, NY+1), xtm(NBDRY), ytm(NBDRY)
  119.       real    clevel(NCONTR)
  120.       real    xmin, xmax, ymin, ymax, zmin, zmax
  121.       real    xpmin, xpmax, ypmin, ypmax
  122.       real    x, y, r, theta, rmax, x0, y0
  123.       real    eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i
  124.       real    div1, div1i, div2, div2i
  125.  
  126.       real    shade_min, shade_max, sh_color
  127.       real    xtick, ytick
  128.       integer    nxsub, nysub
  129.       integer    ncolbox, ncollab
  130.       integer    i, j, kx, lx, ky, ly
  131.       integer    sh_cmap, sh_width
  132.       integer    min_color, min_width, max_color, max_width
  133.       character*8 xopt, yopt
  134.  
  135. ! Set up for plshade call
  136.  
  137.       sh_cmap = 1
  138.       min_color = 1
  139.       min_width = 0
  140.       max_color = 0
  141.       max_width = 0
  142.  
  143.       kx = 1
  144.       lx = NX
  145.       ky = 1
  146.       ly = NY
  147.  
  148. ! Set up r-theta grids
  149. ! Tack on extra cell in theta to handle periodicity.
  150.  
  151.       do 12 i = 1, NX
  152.          r = i - 0.5
  153.          do 10 j = 1, NY
  154.             theta = TWOPI/float(NY) * (j-0.5)
  155.             xg(i,j) = r * cos(theta)
  156.             yg(i,j) = r * sin(theta)
  157.  10      continue
  158.          xg(i, NY+1) = xg(i, 1)
  159.          yg(i, NY+1) = yg(i, 1)
  160.  12   continue
  161.       call a2mnmx(xg, NX, NY, xmin, xmax)
  162.       call a2mnmx(yg, NX, NY, ymin, ymax)
  163.  
  164.       rmax = r
  165.       x0 = (xmin + xmax)/2.
  166.       y0 = (ymin + ymax)/2.
  167.  
  168. ! Potential inside a conducting cylinder (or sphere) by method of images.
  169. ! Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
  170. ! Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
  171. ! Also put in smoothing term at small distances.
  172.  
  173.       eps = 2.
  174.  
  175.       q1 = 1.
  176.       d1 = r/4.
  177.  
  178.       q1i = - q1*r/d1
  179.       d1i = r**2/d1
  180.  
  181.       q2 = -1.
  182.       d2 = r/4.
  183.  
  184.       q2i = - q2*r/d2
  185.       d2i = r**2/d2
  186.  
  187.       do 22 i = 1, NX
  188.          do 20 j = 1, NY
  189.             div1 = sqrt((xg(i,j)-d1)**2 + (yg(i,j)-d1)**2 + eps**2)
  190.             div1i = sqrt((xg(i,j)-d1i)**2 + (yg(i,j)-d1i)**2 + eps**2)
  191.  
  192.             div2 = sqrt((xg(i,j)-d2)**2 + (yg(i,j)+d2)**2 + eps**2)
  193.             div2i = sqrt((xg(i,j)-d2i)**2 + (yg(i,j)+d2i)**2 + eps**2)
  194.  
  195.             z(i,j) = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
  196.  20      continue
  197.  22   continue
  198.  
  199. ! Tack on extra cell in theta to handle periodicity.
  200.  
  201.       do 32 i = 1, NX
  202.          do 30 j = 1, NY
  203.             ztmp(i,j) = z(i,j)
  204.  30      continue
  205.          ztmp(i, NY+1) = z(i, 1)
  206.  32   continue
  207.       call a2mnmx(z, NX, NY, zmin, zmax)
  208.  
  209. ! Set up contour levels.
  210.  
  211.       do 40 i = 1, NCONTR
  212.          clevel(i) = zmin + (i-0.5)*abs(zmax - zmin)/float(NCONTR)
  213.  40   continue
  214.  
  215. ! Advance graphics frame and get ready to plot.
  216.  
  217.       ncolbox = 1
  218.       ncollab = 2
  219.  
  220.       call pladv(0)
  221.       call plcol(ncolbox)
  222.  
  223. ! Scale window to user coordinates.
  224. ! Make a bit larger so the boundary does not get clipped.
  225.  
  226.       eps = 0.05
  227.       xpmin = xmin - abs(xmin)*eps
  228.       xpmax = xmax + abs(xmax)*eps
  229.       ypmin = ymin - abs(ymin)*eps
  230.       ypmax = ymax + abs(ymax)*eps
  231.  
  232.       call plvpas(0.1, 0.9, 0.1, 0.9, 1.0)
  233.       call plwind(xpmin, xpmax, ypmin, ypmax)
  234.  
  235.       xopt = ' '
  236.       yopt = ' '
  237.       xtick = 0.
  238.       nxsub = 0
  239.       ytick = 0.
  240.       nysub = 0
  241.  
  242.       call plbox(xopt, xtick, nxsub, yopt, ytick, nysub)
  243.  
  244. ! Call plotter once for z < 0 (dashed), once for z > 0 (solid lines).
  245.  
  246.       do 100 i = 1, NCONTR
  247.          shade_min = zmin + (zmax - zmin) * real(i - 1) / real(NCONTR)
  248.          shade_max = zmin + (zmax - zmin) * real(i)     / real(NCONTR)
  249.          sh_color = real(i - 1) / real(NCONTR - 1)
  250.          sh_width = 2
  251.          call plpsty(0)
  252.  
  253.          call plshade2(z, NX, NY, ' ', 
  254.      &        -1., 1., -1., 1., 
  255.      &        shade_min, shade_max, 
  256.      &        sh_cmap, sh_color, sh_width, 
  257.      &        min_color, min_width, max_color, max_width, xg, yg)
  258.  
  259.  100  continue
  260.  
  261. ! Draw boundary.
  262.  
  263.       do 110 i = 1, NBDRY
  264.          theta = (TWOPI)/(NBDRY-1) * float(i-1)
  265.          xtm(i) = x0 + rmax * cos(theta)
  266.          ytm(i) = y0 + rmax * sin(theta)
  267.  110  continue
  268.       call plcol(ncolbox)
  269.       call plline(NBDRY, xtm, ytm)
  270.  
  271.       call plcol(ncollab)
  272.       call pllab(' ', ' ', 
  273.      &'Shielded potential of charges in a conducting sphere')
  274.  
  275.       return
  276.       end
  277.  
  278. !----------------------------------------------------------------------------!
  279. ! Subroutine a2mnmx
  280. !----------------------------------------------------------------------------!
  281. ! Minimum and the maximum elements of a 2-d array.
  282. !----------------------------------------------------------------------------!
  283.  
  284.       subroutine a2mnmx(f, nx, ny, fmin, fmax)
  285.  
  286.       integer    nx, ny
  287.       real    f(nx, ny), fmin, fmax
  288.  
  289.       fmax = f(1, 1)
  290.       fmin = fmax
  291.       do 12 j = 1, ny
  292.          do 10 i = 1, nx
  293.             fmax = max(fmax, f(i, j))
  294.             fmin = min(fmin, f(i, j))
  295.  10      continue
  296.  12   continue
  297.  
  298.       return 
  299.       end
  300.