1!      Demonstration of plshade plotting
2!      Reduce colors in cmap 0 so that cmap 1 is useful on a 16-color display
3!
4!      Copyright (C) 2004-2016  Alan W. Irwin
5!
6!      This file is part of PLplot.
7!
8!      PLplot is free software; you can redistribute it and/or modify
9!      it under the terms of the GNU Library General Public License as
10!      published by the Free Software Foundation; either version 2 of the
11!      License, or (at your option) any later version.
12!
13!      PLplot is distributed in the hope that it will be useful,
14!      but WITHOUT ANY WARRANTY; without even the implied warranty of
15!      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16!      GNU Library General Public License for more details.
17!
18!      You should have received a copy of the GNU Library General Public
19!      License along with PLplot; if not, write to the Free Software
20!      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21
22!     N.B. the pl_test_flt parameter used in this code is only
23!     provided by the plplot module to allow convenient developer
24!     testing of either kind(1.0) or kind(1.0d0) floating-point
25!     precision regardless of the floating-point precision of the
26!     PLplot C libraries.  We do not guarantee the value of this test
27!     parameter so it should not be used by users, and instead user
28!     code should replace the pl_test_flt parameter by whatever
29!     kind(1.0) or kind(1.0d0) precision is most convenient for them.
30!     For further details on floating-point precision issues please
31!     consult README_precision in this directory.
32!
33program x16af
34    use plplot
35    implicit none
36    integer :: plparseopts_rc
37
38    !      Process command-line arguments
39    plparseopts_rc = plparseopts(PL_PARSE_FULL)
40    if(plparseopts_rc .ne. 0) stop "plparseopts error"
41
42    call plscmap0n(3)
43
44    !      Initialize plplot
45
46    call plinit()
47
48    !      Rectangular coordinate plot
49
50    call rect()
51
52    !      Polar coordinate plot
53
54    call polar()
55
56    call plend
57
58contains
59
60    !      Plot function using the identity transform
61
62    subroutine rect()
63
64        use plplot
65        implicit none
66        integer   xdim, ydim, NX, NY, NCONTR
67        !      xdim and ydim are the static dimensions of the 2D arrays while
68        !      NX and NY are the defined area.
69        parameter (xdim = 99, NX = 35, ydim = 100, NY = 46, NCONTR = 14)
70
71        real(kind=pl_test_flt)    z(xdim, ydim)
72        real(kind=pl_test_flt)    zmin, zmax, x, y
73        real(kind=pl_test_flt)    shade_min, shade_max, sh_color
74        integer   i, j, sh_cmap
75        integer   min_color, max_color
76        real(kind=pl_test_flt)    sh_width, min_width, max_width
77
78        !      Set up for plshade call
79
80        sh_cmap   = 1
81        min_color = 1
82        min_width = 0
83        max_color = 0
84        max_width = 0
85
86        !      Set up data arrays
87
88        do i = 1, NX
89            x = (i - 1 - (NX/2)) / real(NX/2,kind=pl_test_flt)
90            do j = 1, NY
91                y = (j - 1 - (NY/2)) / real(NY/2,kind=pl_test_flt) - 1.0_pl_test_flt
92                z(i,j) = x*x - y*y + (x - y) / (x*x + y*y + 0.1_pl_test_flt)
93            enddo
94        enddo
95
96        call a2mnmx(z, NX, NY, zmin, zmax, xdim)
97
98        !      Plot using identity transform
99
100        call pladv(0)
101        call plvpor(0.1_pl_test_flt, 0.9_pl_test_flt, 0.1_pl_test_flt, 0.9_pl_test_flt)
102        call plwind(-1.0_pl_test_flt, 1.0_pl_test_flt, -1.0_pl_test_flt, 1.0_pl_test_flt)
103
104        do i = 1, NCONTR
105            shade_min = zmin + (zmax - zmin) * real(i - 1,kind=pl_test_flt) / &
106                   real(NCONTR,kind=pl_test_flt)
107            shade_max = zmin + (zmax - zmin) * real(i,kind=pl_test_flt)     / &
108                   real(NCONTR,kind=pl_test_flt)
109            sh_color = real(i - 1,kind=pl_test_flt) / real(NCONTR - 1,kind=pl_test_flt)
110            sh_width = 2
111            call plpsty(0)
112            call plshade(z(:NX,:NY), &
113                   -1._pl_test_flt, 1.0_pl_test_flt, -1.0_pl_test_flt, 1.0_pl_test_flt, &
114                   shade_min, shade_max, &
115                   sh_cmap, sh_color, sh_width, &
116                   min_color, min_width, max_color, max_width, .true. )
117        enddo
118
119        call plcol0(1)
120        call plbox('bcnst', 0.0_pl_test_flt, 0, 'bcnstv', 0.0_pl_test_flt, 0)
121        call plcol0(2)
122        call pllab('distance', 'altitude', 'Bogon flux')
123
124    end subroutine rect
125
126    !      Routine for demonstrating use_ of transformation arrays in contour plots.
127
128    subroutine polar()
129
130        use plplot, double_TWOPI => PL_TWOPI
131        implicit none
132        real(kind=pl_test_flt), parameter :: TWOPI = double_TWOPI
133        integer   xdim, ydim, NX, NY, NCONTR, NBDRY
134        !      xdim and ydim are the static dimensions of the 2D arrays while
135        !      NX and NY are the defined area.
136        parameter (xdim = 99, NX = 40, ydim = 100, NY = 64)
137        parameter (NCONTR = 14, NBDRY=200)
138
139        real(kind=pl_test_flt)    z(xdim, ydim)
140        real(kind=pl_test_flt)    xg(xdim, ydim+1), yg(xdim, ydim+1), &
141               xtm(NBDRY), ytm(NBDRY)
142        real(kind=pl_test_flt)    xmin, xmax, ymin, ymax, zmin, zmax
143        real(kind=pl_test_flt)    xpmin, xpmax, ypmin, ypmax
144        real(kind=pl_test_flt)    r, theta, rmax, x0, y0
145        real(kind=pl_test_flt)    eps, q1, d1, q1i, d1i, q2, d2, q2i, d2i
146        real(kind=pl_test_flt)    div1, div1i, div2, div2i
147
148        real(kind=pl_test_flt)    shade_min, shade_max, sh_color
149        real(kind=pl_test_flt)    xtick, ytick
150        integer   nxsub, nysub
151        integer   ncolbox, ncollab
152        integer   i, j
153        integer   sh_cmap
154        integer   min_color, max_color
155        real(kind=pl_test_flt) sh_width, min_width, max_width
156        character(len=8) xopt, yopt
157
158        !      Set up for plshade call
159
160        sh_cmap = 1
161        min_color = 1
162        min_width = 0
163        max_color = 0
164        max_width = 0
165
166        !      Set up r-theta grids
167        !      Tack on extra cell in theta to handle periodicity.
168
169        do i = 1, NX
170            r = i - 0.5_pl_test_flt
171            do j = 1, NY
172                theta = TWOPI/real(NY,kind=pl_test_flt) * (j-0.5_pl_test_flt)
173                xg(i,j) = r * cos(theta)
174                yg(i,j) = r * sin(theta)
175            enddo
176            xg(i, NY+1) = xg(i, 1)
177            yg(i, NY+1) = yg(i, 1)
178        enddo
179        call a2mnmx(xg, NX, NY, xmin, xmax, xdim)
180        call a2mnmx(yg, NX, NY, ymin, ymax, xdim)
181
182        rmax = r
183        x0 = (xmin + xmax)/2._pl_test_flt
184        y0 = (ymin + ymax)/2._pl_test_flt
185
186        !      Potential inside a conducting cylinder (or sphere) by method of images.
187        !      Charge 1 is placed at (d1, d1), with image charge at (d2, d2).
188        !      Charge 2 is placed at (d1, -d1), with image charge at (d2, -d2).
189        !      Also put in smoothing term at small distances.
190
191        eps = 2._pl_test_flt
192
193        q1 = 1._pl_test_flt
194        d1 = r/4._pl_test_flt
195
196        q1i = - q1*r/d1
197        d1i = r**2/d1
198
199        q2 = -1._pl_test_flt
200        d2 = r/4._pl_test_flt
201
202        q2i = - q2*r/d2
203        d2i = r**2/d2
204
205        do i = 1, NX
206            do j = 1, NY
207                div1 = sqrt((xg(i,j)-d1)**2 + (yg(i,j)-d1)**2 + eps**2)
208                div1i = sqrt((xg(i,j)-d1i)**2 + (yg(i,j)-d1i)**2 + eps**2)
209
210                div2 = sqrt((xg(i,j)-d2)**2 + (yg(i,j)+d2)**2 + eps**2)
211                div2i = sqrt((xg(i,j)-d2i)**2 + (yg(i,j)+d2i)**2 + eps**2)
212
213                z(i,j) = q1/div1 + q1i/div1i + q2/div2 + q2i/div2i
214            enddo
215        enddo
216
217        call a2mnmx(z, NX, NY, zmin, zmax, xdim)
218
219        !      Advance graphics frame and get ready to plot.
220
221        ncolbox = 1
222        ncollab = 2
223
224        call pladv(0)
225        call plcol0(ncolbox)
226
227        !      Scale window to user coordinates.
228        !      Make a bit larger so the boundary does not get clipped.
229
230        eps = 0.05_pl_test_flt
231        xpmin = xmin - abs(xmin)*eps
232        xpmax = xmax + abs(xmax)*eps
233        ypmin = ymin - abs(ymin)*eps
234        ypmax = ymax + abs(ymax)*eps
235
236        call plvpas(0.1_pl_test_flt, 0.9_pl_test_flt, 0.1_pl_test_flt, 0.9_pl_test_flt, 1.0_pl_test_flt)
237        call plwind(xpmin, xpmax, ypmin, ypmax)
238
239        xopt = ' '
240        yopt = ' '
241        xtick = 0._pl_test_flt
242        nxsub = 0
243        ytick = 0._pl_test_flt
244        nysub = 0
245
246        call plbox(xopt, xtick, nxsub, yopt, ytick, nysub)
247
248        !      Call plotter once for z < 0 (dashed), once for z > 0 (solid lines).
249
250        do i = 1, NCONTR
251            shade_min = zmin + (zmax - zmin) * real(i - 1,kind=pl_test_flt) / &
252                   real(NCONTR,kind=pl_test_flt)
253            shade_max = zmin + (zmax - zmin) * real(i,kind=pl_test_flt)     / &
254                   real(NCONTR,kind=pl_test_flt)
255            sh_color = real(i - 1,kind=pl_test_flt) / real(NCONTR - 1,kind=pl_test_flt)
256            sh_width = 2
257            call plpsty(0)
258
259            call plshade(z(:NX,:NY), &
260                   -1.0_pl_test_flt, 1.0_pl_test_flt, -1.0_pl_test_flt, 1.0_pl_test_flt, &
261                   shade_min, shade_max, &
262                   sh_cmap, sh_color, sh_width, &
263                   min_color, min_width, max_color, max_width, &
264                   .false., xg(:NX,:NY), yg(:NX,:NY) )
265        enddo
266
267        !      Draw boundary.
268
269        do i = 1, NBDRY
270            theta = (TWOPI)/(NBDRY-1) * real(i-1,kind=pl_test_flt)
271            xtm(i) = x0 + rmax * cos(theta)
272            ytm(i) = y0 + rmax * sin(theta)
273        enddo
274        call plcol0(ncolbox)
275        call plline(xtm, ytm)
276
277        call plcol0(ncollab)
278        call pllab(' ', ' ', &
279               'Shielded potential of charges in a conducting sphere')
280
281    end subroutine polar
282
283    !----------------------------------------------------------------------------
284    !      Subroutine a2mnmx
285    !      Minimum and the maximum elements of a 2-d array.
286
287    subroutine a2mnmx(f, nx, ny, fmin, fmax, xdim)
288        use plplot
289        implicit none
290
291        integer   i, j, nx, ny, xdim
292        real(kind=pl_test_flt)    f(xdim, ny), fmin, fmax
293
294        fmax = f(1, 1)
295        fmin = fmax
296        do j = 1, ny
297            do  i = 1, nx
298                fmax = max(fmax, f(i, j))
299                fmin = min(fmin, f(i, j))
300            enddo
301        enddo
302    end subroutine a2mnmx
303end program x16af
304