1 // Sample plots using date / time formatting for axes
2 //
3 // Copyright (C) 2007 Andrew Ross
4 //
5 // This file is part of PLplot.
6 //
7 // PLplot is free software; you can redistribute it and/or modify
8 // it under the terms of the GNU Library General Public License as published
9 // by the Free Software Foundation; either version 2 of the License, or
10 // (at your option) any later version.
11 //
12 // PLplot is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU Library General Public License for more details.
16 //
17 // You should have received a copy of the GNU Library General Public License
18 // along with PLplot; if not, write to the Free Software
19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20 //
21 //
22
23 #include "plcdemos.h"
24
25 static PLFLT x[365], y[365];
26 static PLFLT xerr1[365], xerr2[365], yerr1[365], yerr2[365];
27
28 // Function prototypes
29
30 void plot1( void );
31 void plot2( void );
32 void plot3( void );
33 void plot4( void );
34
35 //--------------------------------------------------------------------------
36 // main
37 //
38 // Draws several plots which demonstrate the use of date / time formats for
39 // the axis labels.
40 // Time formatting is done using the strfqsas routine from the qsastime
41 // library. This is similar to strftime, but works for a broad
42 // date range even on 32-bit systems. See the
43 // documentation of strfqsas for full details of the available formats.
44 //
45 // 1) Plotting temperature over a day (using hours / minutes)
46 // 2) Plotting
47 //
48 // Note: We currently use the default call for plconfigtime (done in
49 // plinit) which means continuous times are interpreted as seconds since
50 // 1970-01-01, but that may change in future, more extended versions of
51 // this example.
52 //
53 //--------------------------------------------------------------------------
54
55 int
main(int argc,char * argv[])56 main( int argc, char *argv[] )
57 {
58 // Parse command line arguments
59 plparseopts( &argc, argv, PL_PARSE_FULL );
60
61 // Initialize plplot
62 plinit();
63
64 // Change the escape character to a '@' instead of the default '#'
65 plsesc( '@' );
66
67 plot1();
68
69 plot2();
70
71 plot3();
72
73 plot4();
74
75 // Don't forget to call plend() to finish off!
76 plend();
77 exit( 0 );
78 }
79
80 // Plot a model diurnal cycle of temperature
81 void
plot1(void)82 plot1( void )
83 {
84 int i, npts;
85 PLFLT xmin, xmax, ymin, ymax;
86
87 // Data points every 10 minutes for 1 day
88 npts = 73;
89
90 xmin = 0;
91 xmax = 60.0 * 60.0 * 24.0; // Number of seconds in a day
92 ymin = 10.0;
93 ymax = 20.0;
94
95 for ( i = 0; i < npts; i++ )
96 {
97 x[i] = xmax * ( (PLFLT) i / (PLFLT) npts );
98 y[i] = 15.0 - 5.0 * cos( 2 * M_PI * ( (PLFLT) i / (PLFLT) npts ) );
99 // Set x error bars to +/- 5 minute
100 xerr1[i] = x[i] - 60 * 5;
101 xerr2[i] = x[i] + 60 * 5;
102 // Set y error bars to +/- 0.1 deg C
103 yerr1[i] = y[i] - 0.1;
104 yerr2[i] = y[i] + 0.1;
105 }
106
107 pladv( 0 );
108
109 // Rescale major ticks marks by 0.5
110 plsmaj( 0.0, 0.5 );
111 // Rescale minor ticks and error bar marks by 0.5
112 plsmin( 0.0, 0.5 );
113
114 plvsta();
115 plwind( xmin, xmax, ymin, ymax );
116
117 // Draw a box with ticks spaced every 3 hour in X and 1 degree C in Y.
118 plcol0( 1 );
119 // Set time format to be hours:minutes
120 pltimefmt( "%H:%M" );
121 plbox( "bcnstd", 3.0 * 60 * 60, 3, "bcnstv", 1, 5 );
122
123 plcol0( 3 );
124 pllab( "Time (hours:mins)", "Temperature (degC)", "@frPLplot Example 29 - Daily temperature" );
125
126 plcol0( 4 );
127
128 plline( npts, x, y );
129 plcol0( 2 );
130 plerrx( npts, xerr1, xerr2, y );
131 plcol0( 3 );
132 plerry( npts, x, yerr1, yerr2 );
133
134 // Rescale major / minor tick marks back to default
135 plsmin( 0.0, 1.0 );
136 plsmaj( 0.0, 1.0 );
137 }
138
139 // Plot the number of hours of daylight as a function of day for a year
140 void
plot2(void)141 plot2( void )
142 {
143 int j, npts;
144 PLFLT xmin, xmax, ymin, ymax;
145 PLFLT lat, p, d;
146
147 // Latitude for London
148 lat = 51.5;
149
150 npts = 365;
151
152 xmin = 0;
153 xmax = npts * 60.0 * 60.0 * 24.0;
154 ymin = 0;
155 ymax = 24;
156
157 // Formula for hours of daylight from
158 // "A Model Comparison for Daylength as a Function of Latitude and
159 // Day of the Year", 1995, Ecological Modelling, 80, pp 87-95.
160 for ( j = 0; j < npts; j++ )
161 {
162 x[j] = j * 60.0 * 60.0 * 24.0;
163 p = asin( 0.39795 * cos( 0.2163108 + 2 * atan( 0.9671396 * tan( 0.00860 * ( j - 186 ) ) ) ) );
164 d = 24.0 - ( 24.0 / M_PI ) *
165 acos( ( sin( 0.8333 * M_PI / 180.0 ) + sin( lat * M_PI / 180.0 ) * sin( p ) ) /
166 ( cos( lat * M_PI / 180.0 ) * cos( p ) ) );
167 y[j] = d;
168 }
169
170 plcol0( 1 );
171 // Set time format to be abbreviated month name followed by day of month
172 pltimefmt( "%b %d" );
173 plprec( 1, 1 );
174 plenv( xmin, xmax, ymin, ymax, 0, 40 );
175
176
177 plcol0( 3 );
178 pllab( "Date", "Hours of daylight", "@frPLplot Example 29 - Hours of daylight at 51.5N" );
179
180 plcol0( 4 );
181
182 plline( npts, x, y );
183
184 plprec( 0, 0 );
185 }
186
187 void
plot3(void)188 plot3( void )
189 {
190 int i, npts;
191 PLFLT xmin, xmax, ymin, ymax;
192 PLFLT tstart;
193
194 // Calculate continuous time corresponding to 2005-12-01 UTC.
195 plctime( 2005, 11, 01, 0, 0, 0., &tstart );
196 npts = 62;
197
198 xmin = tstart;
199 xmax = xmin + npts * 60.0 * 60.0 * 24.0;
200 ymin = 0.0;
201 ymax = 5.0;
202
203 for ( i = 0; i < npts; i++ )
204 {
205 x[i] = xmin + i * 60.0 * 60.0 * 24.0;
206 y[i] = 1.0 + sin( 2 * M_PI * ( (PLFLT) i ) / 7.0 ) +
207 exp( ( (PLFLT) MIN( i, npts - i ) ) / 31.0 );
208 }
209 pladv( 0 );
210
211 plvsta();
212 plwind( xmin, xmax, ymin, ymax );
213
214 plcol0( 1 );
215 // Set time format to be ISO 8601 standard YYYY-MM-DD.
216 pltimefmt( "%F" );
217 // Draw a box with ticks spaced every 14 days in X and 1 hour in Y.
218 plbox( "bcnstd", 14 * 24.0 * 60.0 * 60.0, 14, "bcnstv", 1, 4 );
219
220 plcol0( 3 );
221 pllab( "Date", "Hours of television watched", "@frPLplot Example 29 - Hours of television watched in Dec 2005 / Jan 2006" );
222
223 plcol0( 4 );
224
225 // Rescale symbol size (used by plpoin) by 0.5
226 plssym( 0.0, 0.5 );
227 plpoin( npts, x, y, 2 );
228 plline( npts, x, y );
229 }
230
231 void
plot4(void)232 plot4( void )
233 {
234 // TAI-UTC (seconds) as a function of time.
235
236 PLFLT scale;
237 PLFLT xmin, xmax, ymin = 0.0, ymax = 0.0, xlabel_step = 0.0;
238 int kind, npts = 0, if_TAI_time_format = 0, i;
239 char time_format[10];
240 char title_suffix[100];
241 char xtitle[100];
242 char title[100];
243 PLFLT xx[1001], yy[1001];
244 PLINT epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min;
245 PLFLT epoch_sec;
246 PLINT tai_year, tai_month, tai_day, tai_hour, tai_min;
247 PLFLT tai_sec, tai;
248 PLINT utc_year, utc_month, utc_day, utc_hour, utc_min;
249 PLFLT utc_sec, utc;
250
251 // Continuous time unit is Besselian years from whatever epoch is
252 // chosen below. Could change to seconds (or days) from the
253 // epoch, but then would have to adjust xlabel_step below.
254 scale = 365.242198781;
255 // MJD epoch (see <https://en.wikipedia.org/wiki/Julian_day>).
256 // This is only set for illustrative purposes, and is overwritten
257 // below for the time-representation reasons given in the
258 // discussion below.
259 epoch_year = 1858;
260 epoch_month = 11;
261 epoch_day = 17;
262 epoch_hour = 0;
263 epoch_min = 0;
264 epoch_sec = 0.;
265 // To illustrate the time-representation issues of using the MJD
266 // epoch, in 1985, MJD was roughly 46000 days which corresponds to
267 // 4e9 seconds. Thus, for the -DPL_DOUBLE=ON case where PLFLT is
268 // a double which can represent continuous time to roughly 16
269 // decimal digits of precision the time-representation error is
270 // roughly ~400 nanoseconds. Therefore the MJD epoch would be
271 // acceptable for the plots below in the -DPL_DOUBLE=ON case.
272 // However, that epoch is obviously not acceptable for the
273 // -DPL_DOUBLE=OFF case where PLFLT is a float which can represent
274 // continuous time to only ~7 decimal digits of precision
275 // corresponding to a time representation error of 400 seconds (!)
276 // in 1985. For this reason, we do not use the MJD epoch below
277 // and instead choose the best epoch for each case to minimize
278 // time-representation issues.
279
280 for ( kind = 0; kind < 7; kind++ )
281 {
282 if ( kind == 0 )
283 {
284 // Choose midpoint to maximize time-representation precision.
285 epoch_year = 1985;
286 epoch_month = 0;
287 epoch_day = 2;
288 epoch_hour = 0;
289 epoch_min = 0;
290 epoch_sec = 0.;
291 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
292 plctime( 1950, 0, 2, 0, 0, 0., &xmin );
293 plctime( 2020, 0, 2, 0, 0, 0., &xmax );
294 npts = 70 * 12 + 1;
295 ymin = 0.0;
296 ymax = 36.0;
297 strncpy( time_format, "%Y%", 10 );
298 if_TAI_time_format = 1;
299 strncpy( title_suffix, "from 1950 to 2020", 100 );
300 strncpy( xtitle, "Year", 100 );
301 xlabel_step = 10.;
302 }
303 else if ( kind == 1 || kind == 2 )
304 {
305 // Choose midpoint to maximize time-representation precision.
306 epoch_year = 1961;
307 epoch_month = 7;
308 epoch_day = 1;
309 epoch_hour = 0;
310 epoch_min = 0;
311 epoch_sec = 1.64757;
312 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
313 plctime( 1961, 7, 1, 0, 0, 1.64757 - .20, &xmin );
314 plctime( 1961, 7, 1, 0, 0, 1.64757 + .20, &xmax );
315 npts = 1001;
316 ymin = 1.625;
317 ymax = 1.725;
318 strncpy( time_format, "%S%2%", 10 );
319 strncpy( title_suffix, "near 1961-08-01 (TAI)", 100 );
320 xlabel_step = 0.05 / ( scale * 86400. );
321 if ( kind == 1 )
322 {
323 if_TAI_time_format = 1;
324 strncpy( xtitle, "Seconds (TAI)", 100 );
325 }
326 else
327 {
328 if_TAI_time_format = 0;
329 strncpy( xtitle, "Seconds (TAI) labelled with corresponding UTC", 100 );
330 }
331 }
332 else if ( kind == 3 || kind == 4 )
333 {
334 // Choose midpoint to maximize time-representation precision.
335 epoch_year = 1963;
336 epoch_month = 10;
337 epoch_day = 1;
338 epoch_hour = 0;
339 epoch_min = 0;
340 epoch_sec = 2.6972788;
341 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
342 plctime( 1963, 10, 1, 0, 0, 2.6972788 - .20, &xmin );
343 plctime( 1963, 10, 1, 0, 0, 2.6972788 + .20, &xmax );
344 npts = 1001;
345 ymin = 2.55;
346 ymax = 2.75;
347 strncpy( time_format, "%S%2%", 10 );
348 strncpy( title_suffix, "near 1963-11-01 (TAI)", 100 );
349 xlabel_step = 0.05 / ( scale * 86400. );
350 if ( kind == 3 )
351 {
352 if_TAI_time_format = 1;
353 strncpy( xtitle, "Seconds (TAI)", 100 );
354 }
355 else
356 {
357 if_TAI_time_format = 0;
358 strncpy( xtitle, "Seconds (TAI) labelled with corresponding UTC", 100 );
359 }
360 }
361 else if ( kind == 5 || kind == 6 )
362 {
363 // Choose midpoint to maximize time-representation precision.
364 epoch_year = 2009;
365 epoch_month = 0;
366 epoch_day = 1;
367 epoch_hour = 0;
368 epoch_min = 0;
369 epoch_sec = 34.;
370 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
371 plctime( 2009, 0, 1, 0, 0, 34. - 5., &xmin );
372 plctime( 2009, 0, 1, 0, 0, 34. + 5., &xmax );
373 npts = 1001;
374 ymin = 32.5;
375 ymax = 34.5;
376 strncpy( time_format, "%S%2%", 10 );
377 strncpy( title_suffix, "near 2009-01-01 (TAI)", 100 );
378 xlabel_step = 1. / ( scale * 86400. );
379 if ( kind == 5 )
380 {
381 if_TAI_time_format = 1;
382 strncpy( xtitle, "Seconds (TAI)", 100 );
383 }
384 else
385 {
386 if_TAI_time_format = 0;
387 strncpy( xtitle, "Seconds (TAI) labelled with corresponding UTC", 100 );
388 }
389 }
390
391 for ( i = 0; i < npts; i++ )
392 {
393 xx[i] = xmin + i * ( xmax - xmin ) / ( (PLFLT) ( npts - 1 ) );
394 tai = xx[i];
395 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
396 plbtime( &tai_year, &tai_month, &tai_day, &tai_hour, &tai_min, &tai_sec, tai );
397 // Calculate residual using tai as the epoch to nearly maximize time-representation precision.
398 plconfigtime( scale, 0., 0., 0x0, 1, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec );
399 // Calculate continuous tai with new epoch.
400 plctime( tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec, &tai );
401 // Calculate broken-down utc (with leap seconds inserted) from continuous tai with new epoch.
402 plconfigtime( scale, 0., 0., 0x2, 1, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec );
403 plbtime( &utc_year, &utc_month, &utc_day, &utc_hour, &utc_min, &utc_sec, tai );
404 // Calculate continuous utc from broken-down utc using same epoch as for the continuous tai.
405 plconfigtime( scale, 0., 0., 0x0, 1, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec );
406 plctime( utc_year, utc_month, utc_day, utc_hour, utc_min, utc_sec, &utc );
407 // Convert residuals to seconds.
408 yy[i] = ( tai - utc ) * scale * 86400.;
409 }
410
411 pladv( 0 );
412 plvsta();
413 plwind( xmin, xmax, ymin, ymax );
414 plcol0( 1 );
415 if ( if_TAI_time_format )
416 plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
417 else
418 plconfigtime( scale, 0., 0., 0x2, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
419 pltimefmt( time_format );
420 plbox( "bcnstd", xlabel_step, 0, "bcnstv", 0., 0 );
421 plcol0( 3 );
422 strncpy( title, "@frPLplot Example 29 - TAI-UTC ", 100 );
423 strncat( title, title_suffix, 100 - strlen( title ) - 1 );
424 pllab( xtitle, "TAI-UTC (sec)", title );
425
426 plcol0( 4 );
427
428 plline( npts, xx, yy );
429 }
430 }
431