This file is indexed.

/usr/share/doc/libplplot11/examples/f95/x21f.f90 is in libplplot-dev 5.9.9-2ubuntu2.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
!      $Id: x21f.f90 11680 2011-03-27 17:57:51Z airwin $
!      Grid data demo.
!
!      Copyright (C) 2004  Joao Cardoso
!      Copyright (C) 2008  Andrew Ross
!
!      This file is part of PLplot.
!
!      PLplot is free software; you can redistribute it and/or modify
!      it under the terms of the GNU Library General Public License as
!      published by the Free Software Foundation; either version 2 of the
!      License, or (at your option) any later version.
!
!      PLplot is distributed in the hope that it will be useful,
!      but WITHOUT ANY WARRANTY; without even the implied warranty of
!      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!      GNU Library General Public License for more details.
!
!      You should have received a copy of the GNU Library General Public
!      License along with PLplot; if not, write to the Free Software
!      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA

      use plplot, PI => PL_PI

      implicit none

      external myisnan
      logical myisnan

      integer pts, xp, yp, nl, knn_order, randn, rosen
      real(kind=plflt) threshold, wmin
      parameter (pts = 500)
      parameter (xp = 25)
      parameter (yp = 20)
      parameter (nl = 16)
      parameter (knn_order = 20)
      parameter (threshold = 1.001_plflt)
      parameter (wmin = -1e3_plflt)
      parameter (randn = 0)
      parameter (rosen = 0)

      real(kind=plflt) xmin, xmax, ymin, ymax
      
      real(kind=plflt) x(pts), y(pts), z(pts), clev(nl)
      real(kind=plflt) xg(xp), yg(yp), zg(xp,yp)
      real(kind=plflt) zmin, zmax, lzmin, lzmax
      integer i, j, k
      integer alg
      character*80 title(6)
      data title /'Cubic Spline Approximation', &
                 'Delaunay Linear Interpolation', &
     		  'Natural Neighbors Interpolation', &
     		  'KNN Inv. Distance Weighted', &
     		  '3NN Linear Interpolation', &
                 '4NN Around Inv. Dist. Weighted'/

      real(kind=plflt) opt(6)
      data opt /0._plflt, 0._plflt, 0._plflt, 0._plflt, 0._plflt, 0._plflt/

      real(kind=plflt) xt, yt

      real(kind=plflt) r
      integer ii, jj
      real(kind=plflt) dist, d

      character*1 defined

      xmin = -0.2_plflt
      ymin = -0.2_plflt
      xmax = 0.6_plflt
      ymax = 0.6_plflt

!      call plMergeOpts(options, "x21c options", NULL);
      call plparseopts(PL_PARSE_FULL)

      opt(3) = wmin
      opt(4) = dble(knn_order)
      opt(5) = threshold

! Initialize plplot

      call plinit

      call plseed(5489)

      do i=1,pts
         xt = (xmax-xmin)*plrandd()
         yt = (ymax-ymin)*plrandd()
         if (randn.eq.0) then
            x(i) = xt + xmin
            y(i) = yt + ymin
         else
            x(i) = sqrt(-2._plflt*log(xt)) * cos(2._plflt*PI*yt) + xmin
            y(i) = sqrt(-2._plflt*log(xt)) * sin(2._plflt*PI*yt) + ymin
         endif
         if (rosen.eq.0) then
            r = sqrt(x(i)*x(i) + y(i)*y(i))
            z(i) = exp(-r*r)*cos(2._plflt*PI*r)
         else
            z(i) = log((1._plflt-x(i))**2 + 100._plflt*(y(i)-x(i)**2)**2)
         endif
      enddo
      
      zmin = z(1)
      zmax = z(1)
      do i=2,pts
         zmax = max(zmax,z(i))
         zmin = min(zmin,z(i))
      enddo

      do i=1,xp
         xg(i) = xmin + (xmax-xmin)*(i-1._plflt)/(xp-1._plflt)
      enddo
      do i=1,yp
         yg(i) = ymin + (ymax-ymin)*(i-1._plflt)/(yp-1._plflt)
      enddo
      
      call plcol0(1)
      call plenv(xmin, xmax, ymin, ymax, 2, 0)
      call plcol0(15)
      call pllab("X", "Y", "The original data sampling")
      call plcol0(2)
      call plpoin(x, y, 5)
      call pladv(0)

      call plssub(3,2)

      do k=1,2
         call pladv(0)
         do alg=1,6

            call plgriddata(x, y, z, xg, yg, zg, alg, opt(alg))

!     - CSA can generate NaNs (only interpolates? !).
!     - DTLI and NNI can generate NaNs for points outside the convex hull
!     of the data points.
!     - NNLI can generate NaNs if a sufficiently thick triangle is not found
!
!     PLplot should be NaN/Inf aware, but changing it now is quite a job...
!     so, instead of not plotting the NaN regions, a weighted average over
!     the neighbors is done.
!

            if ((alg.eq.GRID_CSA).or.(alg.eq.GRID_DTLI).or. &
                 (alg.eq.GRID_NNLI).or.(alg.eq.GRID_NNI)) then
               
               do i=1,xp
                  do j=1,yp
                     if (myisnan(zg(i,j))) then
!     average (IDW) over the 8 neighbors
                        
                        zg(i,j) = 0._plflt
                        dist = 0._plflt
                        
                        ii=i-1
                        do while ((ii.le.i+1).and.(ii.le.xp))
                           jj = j-1
                           do while ((jj.le.j+1).and.(jj.le.yp))
                              if ((ii.ge.1) .and. (jj.ge.1) .and. &
                                   (.not.myisnan(zg(ii,jj))) ) then
                                 if (abs(ii-i) + abs(jj-j) .eq. 1) then
                                    d = 1._plflt
                                 else
                                    d = 1.4142_plflt
                                 endif
                                 zg(i,j) = zg(i,j) + zg(ii,jj)/(d*d)
                                 dist = dist + d
                              endif
                              jj = jj+1
                           enddo
                           ii = ii+1
                        enddo
                        if (dist.ne.0._plflt) then
                           zg(i,j) = zg(i,j) / dist
                        else
                           zg(i,j) = zmin
                        endif
                     endif
                  enddo
               enddo
            endif
            
            call a2mnmx(zg, xp, yp, lzmin, lzmax, xp)

            lzmin = min(lzmin, zmin)
            lzmax = max(lzmax, zmax)

            lzmin = lzmin - 0.01_plflt
            lzmax = lzmax + 0.01_plflt

            call plcol0(1)
            call pladv(alg)
            
            if (k.eq.1) then
               
               do i=1,nl
                  clev(i) = lzmin + (lzmax-lzmin)/(nl-1._plflt)*(i-1._plflt)
               enddo
               call plenv0(xmin, xmax, ymin, ymax, 2, 0)
               call plcol0(15)
               call pllab("X", "Y", title(alg))
               call plshades(zg, defined, xmin, xmax, ymin, &
                    ymax, clev, 1, 0, 1)
               call plcol0(2)
            else
               
               do i = 1,nl
                  clev(i) = lzmin + (lzmax-lzmin)/(nl-1._plflt)*(i-1._plflt)
               enddo
               call cmap1_init()
               call plvpor(0._plflt, 1._plflt, 0._plflt, 0.9_plflt)
               call plwind(-1.1_plflt, 0.75_plflt, -0.65_plflt, 1.20_plflt)
!     
!     For the comparison to be fair, all plots should have the
!     same z values, but to get the max/min of the data generated
!     by all algorithms would imply two passes. Keep it simple.
!
!     plw3d(1., 1., 1., xmin, xmax, ymin, ymax, zmin, zmax, 30, -60);
!
               
               call plw3d(1._plflt, 1._plflt, 1._plflt, xmin, xmax, ymin, ymax,  &
                    lzmin, lzmax, 30._plflt, -40._plflt)
               call plbox3("bntu", "X", 0._plflt, 0, &
                   "bntu", "Y", 0._plflt, 0, &
                   "bcdfntu", "Z", 0.5_plflt, 0)
               call plcol0(15)
               call pllab("", "", title(alg))
               call plot3dc(xg, yg, zg, ior(ior(DRAW_LINEXY, &
                   MAG_COLOR), BASE_CONT), clev)
            endif
         enddo
      enddo

      call plend
            
      end

      subroutine cmap1_init
        use plplot
        implicit none
        real(kind=plflt) i(2), h(2), l(2), s(2)
        
        i(1) = 0._plflt
        i(2) = 1._plflt
        
        h(1) = 240._plflt
        h(2) = 0._plflt
        
        l(1) = 0.6_plflt
        l(2) = 0.6_plflt
        
        s(1) = 0.8_plflt
        s(2) = 0.8_plflt
        
        call plscmap1n(256)
        call plscmap1l(.false., i, h, l, s)
      end subroutine cmap1_init


!----------------------------------------------------------------------------
!      Subroutine a2mnmx
!      Minimum and the maximum elements of a 2-d array.

      subroutine a2mnmx(f, nx, ny, fmin, fmax, xdim)
        use plplot
        implicit none

        integer   i, j, nx, ny, xdim
        real(kind=plflt) f(xdim, ny), fmin, fmax
        
        fmax = f(1, 1)
        fmin = fmax
        do j = 1, ny
           do  i = 1, nx
              fmax = max(fmax, f(i, j))
              fmin = min(fmin, f(i, j))
           enddo
        enddo
      end subroutine a2mnmx

      include 'plf95demos.inc'