← Back to branch summary

~uhh-ssd/+junk/humidity_readout

« back to all changes in this revision

Viewing changes to plplot/plplot-5.9.9/examples/d/x21d.d

  • Committer: Joachim Erfle
  • Date: 2013-07-24 13:53:41 UTC
  • Revision ID: joachim.erfle@desy.de-20130724135341-1qojpp701zsn009p
initial commit

Show diffs side-by-side

added added

removed removed

Lines of Context:
plplot/plplot-5.9.9/examples/d/x21d.d
 
1
// $Id: x21d.d 11684 2011-03-31 04:15:32Z airwin $
 
2
//      Grid data demo
 
3
//
 
4
// Copyright (C) 2009  Werner Smekal
 
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 published
 
10
// by the Free Software Foundation; either version 2 of the License, or
 
11
// (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 License
 
19
// 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
//
 
23
 
 
24
import std.math;
 
25
 
 
26
import plplot;
 
27
 
 
28
// Options data structure definition.
 
29
PLINT pts       = 500;
 
30
PLINT xp        = 25;
 
31
PLINT yp        = 20;
 
32
PLINT nl        = 16;
 
33
int   knn_order = 20;
 
34
PLFLT threshold = 1.001;
 
35
PLFLT wmin      = -1e3;
 
36
int   randn     = 0;
 
37
int   rosen     = 0;
 
38
 
 
39
 
 
40
PLFLT xm, xM, ym, yM;
 
41
 
 
42
int main( char[][] args )
 
43
{
 
44
    string[] title = [ "Cubic Spline Approximation",
 
45
                       "Delaunay Linear Interpolation",
 
46
                       "Natural Neighbors Interpolation",
 
47
                       "KNN Inv. Distance Weighted",
 
48
                       "3NN Linear Interpolation",
 
49
                       "4NN Around Inv. Dist. Weighted" ];
 
50
 
 
51
    xm = ym = -0.2;
 
52
    xM = yM = 0.6;
 
53
 
 
54
    // plMergeOpts(options, "x21c options", NULL);
 
55
    plparseopts( args, PL_PARSE_FULL );
 
56
 
 
57
    PLFLT[] opt = [ 0.0, 0.0, wmin, knn_order, threshold, 0.0 ];
 
58
 
 
59
    // Initialize plplot
 
60
    plinit();
 
61
 
 
62
    // Initialise random number generator
 
63
    plseed( 5489 );
 
64
 
 
65
    PLFLT[] x, y, z;
 
66
    x.length = y.length = z.length = pts;
 
67
    create_data( x, y, z ); // the sampled data
 
68
    PLFLT zmin = z[0];
 
69
    PLFLT zmax = z[0];
 
70
    for ( int i = 1; i < pts; i++ )
 
71
    {
 
72
        if ( z[i] > zmax )
 
73
            zmax = z[i];
 
74
        if ( z[i] < zmin )
 
75
            zmin = z[i];
 
76
    }
 
77
 
 
78
    PLFLT[] xg, yg;
 
79
    xg.length = xp;
 
80
    yg.length = yp;
 
81
    create_grid( xg, yg ); // grid the data at
 
82
 
 
83
    PLFLT[][] zg = new PLFLT[][xp];
 
84
    for ( int i = 0; i < xp; i++ )
 
85
        zg[i] = new PLFLT[yp];
 
86
 
 
87
    PLFLT[] clev = new PLFLT[nl];
 
88
 
 
89
    plcol0( 1 );
 
90
    plenv( xm, xM, ym, yM, 2, 0 );
 
91
    plcol0( 15 );
 
92
    pllab( "X", "Y", "The original data sampling" );
 
93
    plcol0( 2 );
 
94
    plpoin( x, y, 5 );
 
95
    pladv( 0 );
 
96
 
 
97
    plssub( 3, 2 );
 
98
 
 
99
    for ( int k = 0; k < 2; k++ )
 
100
    {
 
101
        pladv( 0 );
 
102
        for ( int alg = 1; alg < 7; alg++ )
 
103
        {
 
104
            plgriddata( x, y, z, xg, yg, zg, alg, opt[alg - 1] );
 
105
 
 
106
            // - CSA can generate NaNs (only interpolates?!).
 
107
            // - DTLI and NNI can generate NaNs for points outside the convex hull
 
108
            //      of the data points.
 
109
            // - NNLI can generate NaNs if a sufficiently thick triangle is not found
 
110
            //
 
111
            // PLplot should be NaN/Inf aware, but changing it now is quite a job...
 
112
            // so, instead of not plotting the NaN regions, a weighted average over
 
113
            // the neighbors is done.
 
114
            //
 
115
 
 
116
            if ( alg == GRID_CSA || alg == GRID_DTLI || alg == GRID_NNLI || alg == GRID_NNI )
 
117
            {
 
118
                PLFLT dist, d;
 
119
 
 
120
                for ( int i = 0; i < xp; i++ )
 
121
                {
 
122
                    for ( int j = 0; j < yp; j++ )
 
123
                    {
 
124
                        if ( isnan( zg[i][j] ) )                       // average (IDW) over the 8 neighbors
 
125
                        {
 
126
                            zg[i][j] = 0.0;
 
127
                            dist     = 0.0;
 
128
 
 
129
                            for ( int ii = i - 1; ii <= i + 1 && ii < xp; ii++ )
 
130
                            {
 
131
                                for ( int jj = j - 1; jj <= j + 1 && jj < yp; jj++ )
 
132
                                {
 
133
                                    if ( ii >= 0 && jj >= 0 && !isnan( zg[ii][jj] ) )
 
134
                                    {
 
135
                                        d         = ( abs( ii - i ) + abs( jj - j ) ) == 1 ? 1.0 : 1.4142;
 
136
                                        zg[i][j] += zg[ii][jj] / ( d * d );
 
137
                                        dist     += d;
 
138
                                    }
 
139
                                }
 
140
                            }
 
141
                            if ( dist != 0.0 )
 
142
                                zg[i][j] /= dist;
 
143
                            else
 
144
                                zg[i][j] = zmin;
 
145
                        }
 
146
                    }
 
147
                }
 
148
            }
 
149
 
 
150
            PLFLT lzM, lzm;
 
151
            plMinMax2dGrid( zg, lzM, lzm );
 
152
 
 
153
            lzm = fmin( lzm, zmin );
 
154
            lzM = fmax( lzM, zmax );
 
155
 
 
156
            // Increase limits slightly to prevent spurious contours
 
157
            // due to rounding errors
 
158
            lzm = lzm - 0.01;
 
159
            lzM = lzM + 0.01;
 
160
 
 
161
            plcol0( 1 );
 
162
 
 
163
            pladv( alg );
 
164
 
 
165
            if ( k == 0 )
 
166
            {
 
167
                for ( int i = 0; i < nl; i++ )
 
168
                    clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i;
 
169
 
 
170
                plenv0( xm, xM, ym, yM, 2, 0 );
 
171
                plcol0( 15 );
 
172
                pllab( "X", "Y", title[alg - 1] );
 
173
                plshades( zg, null, xm, xM, ym, yM, clev, 1, 0, 1, 1 );
 
174
                plcol0( 2 );
 
175
            }
 
176
            else
 
177
            {
 
178
                for ( int i = 0; i < nl; i++ )
 
179
                    clev[i] = lzm + ( lzM - lzm ) / ( nl - 1 ) * i;
 
180
 
 
181
                cmap1_init();
 
182
                plvpor( 0.0, 1.0, 0.0, 0.9 );
 
183
                plwind( -1.1, 0.75, -0.65, 1.20 );
 
184
                //
 
185
                // For the comparison to be fair, all plots should have the
 
186
                // same z values, but to get the max/min of the data generated
 
187
                // by all algorithms would imply two passes. Keep it simple.
 
188
                //
 
189
                // plw3d(1., 1., 1., xm, xM, ym, yM, zmin, zmax, 30, -60);
 
190
                //
 
191
 
 
192
                plw3d( 1., 1., 1., xm, xM, ym, yM, lzm, lzM, 30, -40 );
 
193
                plbox3( "bntu", "X", 0., 0,
 
194
                    "bntu", "Y", 0., 0,
 
195
                    "bcdfntu", "Z", 0.5, 0 );
 
196
                plcol0( 15 );
 
197
                pllab( "", "", title[alg - 1] );
 
198
                plot3dc( xg, yg, zg, DRAW_LINEXY | MAG_COLOR | BASE_CONT, clev );
 
199
            }
 
200
        }
 
201
    }
 
202
 
 
203
    plend();
 
204
 
 
205
    return 0;
 
206
}
 
207
 
 
208
 
 
209
void create_grid( PLFLT[] x, PLFLT[] y )
 
210
{
 
211
    int px = x.length;
 
212
    int py = y.length;
 
213
 
 
214
    for ( int i = 0; i < px; i++ )
 
215
        x[i] = xm + ( xM - xm ) * i / ( px - 1.0 );
 
216
 
 
217
    for ( int i = 0; i < py; i++ )
 
218
        y[i] = ym + ( yM - ym ) * i / ( py - 1.0 );
 
219
}
 
220
 
 
221
 
 
222
void create_data( PLFLT[] x, PLFLT[] y, PLFLT[] z )
 
223
{
 
224
    int pts = x.length;
 
225
    assert( pts == y.length, "create_data(): Arrays must be of same length" );
 
226
    assert( pts == z.length, "create_data(): Arrays must be of same length" );
 
227
 
 
228
    PLFLT xt, yt, r;
 
229
    for ( int i = 0; i < pts; i++ )
 
230
    {
 
231
        xt = ( xM - xm ) * plrandd();
 
232
        yt = ( yM - ym ) * plrandd();
 
233
        if ( !randn )
 
234
        {
 
235
            x[i] = xt + xm;
 
236
            y[i] = yt + ym;
 
237
        }
 
238
        else // std=1, meaning that many points are outside the plot range
 
239
        {
 
240
            x[i] = sqrt( -2.0 * log( xt ) ) * cos( 2. * PI * yt ) + xm;
 
241
            y[i] = sqrt( -2.0 * log( xt ) ) * sin( 2. * PI * yt ) + ym;
 
242
        }
 
243
        if ( !rosen )
 
244
        {
 
245
            r    = sqrt( x[i] * x[i] + y[i] * y[i] );
 
246
            z[i] = exp( -r * r ) * cos( 2.0 * PI * r );
 
247
        }
 
248
        else
 
249
            z[i] = log( pow( 1. - x[i], 2.9 ) + 100.0 * pow( y[i] - pow( x[i], 2.0 ), 2.0 ) );
 
250
    }
 
251
}
 
252
 
 
253
 
 
254
void cmap1_init()
 
255
{
 
256
    PLFLT[] i = [ 0.0, 1.0 ];           // boundaries
 
257
 
 
258
    PLFLT[] h = [ 240.0, 0.0 ];         // blue -> green -> yellow -> red
 
259
    PLFLT[] l = [ 0.6, 0.6 ];
 
260
    PLFLT[] s = [ 0.8, 0.8 ];
 
261
 
 
262
    plscmap1n( 256 );
 
263
    plscmap1l( 0, i, h, l, s );
 
264
}