1## Sample plots using date / time formatting for axes
2##
3## Copyright (C) 2007, 2008 Andrew Ross
4##
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
231;
24
25##
26
27function ix29c
28
29
30    ## Initialize plplot
31    plinit();
32
33    plsesc('@');
34
35    plot1();
36
37    plot2();
38
39    plot3();
40
41    plot4();
42
43    plend1();
44
45endfunction
46
47## Plot a model diurnal cycle of temperature
48function plot1
49
50## Data points every 10 minutes for 1 day
51    npts = 73;
52
53    xmin = 0;
54    xmax = 60.0*60.0*24.0;    ## Number of seconds in a day
55    ymin = 10.0;
56    ymax = 20.0;
57
58    i = 0:npts-1;
59    x = i*xmax/npts;
60    y = 15.0 - 5.0*cos( 2*pi*i/npts);
61    xerr1 = x-60.0*5.0;
62    xerr2 = x+60.0*5.0;
63    yerr1 = y-0.1;
64    yerr2 = y+0.1;
65
66    pladv(0);
67
68    ## Rescale major ticks marks by 0.5
69    plsmaj(0.0,0.5);
70    ## Rescale minor ticks and error bar marks by 0.5
71    plsmin(0.0,0.5);
72
73    plvsta();
74    plwind(xmin, xmax, ymin, ymax);
75
76    ## Draw a box with ticks spaced every 3 hour in X and 1 degree C in Y.
77    plcol0(1);
78    ## Set time format to be hours:minutes
79    pltimefmt("%H:%M");
80    plbox("bcnstd", 3.0*60*60, 3, "bcnstv", 1, 5);
81
82    plcol0(3);
83    pllab("Time (hours:mins)", "Temperature (degC)", "@frPLplot Example 29 - Daily temperature");
84
85    plcol0(4);
86
87    plline(x', y');
88    plcol0(2);
89    plerrx(xerr1', xerr2', y');
90    plcol0(3);
91    plerry(x', yerr1', yerr2');
92
93    ## Rescale major / minor tick marks back to default
94    plsmin(0.0,1.0);
95    plsmaj(0.0,1.0);
96
97endfunction
98
99## Plot the number of hours of daylight as a function of day for a year
100function plot2
101
102## Latitude for London
103    lat = 51.5;
104
105    npts = 365;
106
107    xmin = 0;
108    xmax = npts*60.0*60.0*24.0;
109    ymin = 0;
110    ymax = 24;
111
112    ## Formula for hours of daylight from
113    ## "A Model Comparison for Daylength as a Function of Latitude and
114    ## Day of the Year", 1995, Ecological Modelling, 80, pp 87-95.
115    i = 0:npts-1;
116    x = i*60.0*60.0*24.0;
117    p = asin(0.39795*cos(0.2163108 + 2*atan(0.9671396*tan(0.00860*(i-186)))));
118    d = 24.0 - (24.0/pi)*acos( (sin(0.8333*pi/180.0) + ...
119                                sin(lat*pi/180.0)*sin(p)) ...
120                               ./(cos(lat*pi/180.0)*cos(p)) );
121    y = d;
122
123    plcol0(1);
124    ## Set time format to be abbreviated month name followed by day of month
125    pltimefmt("%b %d");
126    plprec(1,1);
127    plenv(xmin, xmax, ymin, ymax, 0, 40);
128
129
130    plcol0(3);
131    pllab("Date", "Hours of daylight", "@frPLplot Example 29 - Hours of daylight at 51.5N");
132
133    plcol0(4);
134
135    plline(x', y');
136
137    plprec(0,0);
138
139endfunction
140
141function plot3
142
143    ## Calculate seconds since the Unix epoch for 2005-12-01 UTC.
144    tm.year = 105; ## Years since 1900
145    tm.mon = 11;   ## 0 == January, 11 = December
146    tm.mday = 1;   ## 1 = 1st of month
147    tm.hour = 0;
148    tm.min = 0;
149    tm.sec = 0;
150
151    ## NB - no need to call tzset in octave - it doesn't exist.
152    tz = getenv("TZ");
153    putenv("TZ","");
154
155    ## tstart is a time_t value (cast to PLFLT) which represents the number
156    ## of seconds elapsed since 00:00:00 on January 1, 1970, Coordinated
157    ## Universal Time (UTC).
158    tstart = mktime(tm);
159
160    ## Note currently octave appears to have no way to unset a env
161    ## variable.
162    putenv("TZ",tz);
163
164    npts = 62;
165
166    xmin = tstart;
167    xmax = xmin + npts*60.0*60.0*24.0;
168    ymin = 0.0;
169    ymax = 5.0;
170
171    i = 0:npts-1;
172    x = xmin + i*60.0*60.0*24.0;
173    y = 1.0 + sin( 2*pi*i / 7.0 ) + exp( min(i,npts-i) / 31.0);
174
175    pladv(0);
176
177    plvsta();
178    plwind(xmin, xmax, ymin, ymax);
179    plcol0(1);
180    ## Set time format to be ISO 8601 standard YYYY-MM-DD. Note that this is
181    ## equivalent to %f for C99 compliant implementations of strftime.
182    pltimefmt("%Y-%m-%d");
183    ## Draw a box with ticks spaced every 14 days in X and 1 hour in Y.
184    plbox("bcnstd", 14*24.0*60.0*60.0,14, "bcnstv", 1, 4);
185
186    plcol0(3);
187    pllab("Date", "Hours of television watched", "@frPLplot Example 29 - Hours of television watched in Dec 2005 / Jan 2006");
188
189    plcol0(4);
190
191    ## Rescale symbol size (used by plpoin) by 0.5
192    plssym(0.0,0.5);
193    plpoin(x', y', 2);
194    plline(x', y');
195
196endfunction
197
198function plot4
199
200    ## TAI-UTC (seconds) as a function of time.
201
202    ## Continuous time unit is Besselian years from whatever epoch is
203    ## chosen below.  Could change to seconds (or days) from the
204    ## epoch, but then would have to adjust xlabel_step below.
205    scale = 365.242198781;
206    ## MJD epoch (see <https://en.wikipedia.org/wiki/Julian_day>).
207    ## This is only set for illustrative purposes, and is overwritten
208    ## below for the time-representation reasons given in the
209    ## discussion below.
210    epoch_year  = 1858;
211    epoch_month = 11;
212    epoch_day   = 17;
213    epoch_hour  = 0;
214    epoch_min   = 0;
215    epoch_sec   = 0.;
216    ## To illustrate the time-representation issues of using the MJD
217    ## epoch, in 1985, MJD was roughly 46000 days which corresponds to
218    ## 4e9 seconds.  Thus, for the -DPL_DOUBLE=ON case where PLFLT is
219    ## a double which can represent continuous time to roughly 16
220    ## decimal digits of precision the time-representation error is
221    ## roughly ~400 nanoseconds.  Therefore the MJD epoch would be
222    ## acceptable for the plots below in the -DPL_DOUBLE=ON case.
223    ## However, that epoch is obviously not acceptable for the
224    ## -DPL_DOUBLE=OFF case where PLFLT is a float which can represent
225    ## continuous time to only ~7 decimal digits of precision
226    ## corresponding to a time representation error of 400 seconds (!)
227    ## in 1985.  For this reason, we do not use the MJD epoch below
228    ## and instead choose the best epoch for each case to minimize
229    ## time-representation issues.
230
231    for kind=0:6
232        if (kind == 0)
233            ## Choose midpoint to maximize time-representation precision.
234            epoch_year  = 1985;
235            epoch_month = 0;
236            epoch_day   = 2;
237            epoch_hour  = 0;
238            epoch_min   = 0;
239            epoch_sec   = 0.;
240            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
241            xmin = plctime(1950,0,2,0,0,0.);
242            xmax = plctime(2020,0,2,0,0,0.);
243            npts = 70*12 + 1;
244            ymin = 0.0;
245            ymax = 36.0;
246            time_format = "%Y%";
247            if_TAI_time_format = 1;
248            title_suffix = "from 1950 to 2020";
249            xtitle =  "Year";
250            xlabel_step = 10.;
251        elseif (kind == 1 || kind ==2)
252            ## Choose midpoint to maximize time-representation precision.
253            epoch_year  = 1961;
254            epoch_month = 7;
255            epoch_day   = 1;
256            epoch_hour  = 0;
257            epoch_min   = 0;
258            epoch_sec   = 1.64757;
259            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
260            xmin = plctime(1961,7,1,0,0,1.64757-.20);
261            xmax = plctime(1961,7,1,0,0,1.64757+.20);
262            npts = 1001;
263            ymin = 1.625;
264            ymax = 1.725;
265            time_format = "%S%2%";
266            title_suffix = "near 1961-08-01 (TAI)";
267            xlabel_step = 0.05/(scale*86400.);
268            if (kind == 1)
269                if_TAI_time_format = 1;
270                xtitle = "Seconds (TAI)";
271            else
272                if_TAI_time_format = 0;
273                xtitle = "Seconds (TAI) labelled with corresponding UTC";
274            endif
275        elseif (kind == 3 || kind ==4)
276            ## Choose midpoint to maximize time-representation precision.
277            epoch_year  = 1963;
278            epoch_month = 10;
279            epoch_day   = 1;
280            epoch_hour  = 0;
281            epoch_min   = 0;
282            epoch_sec   = 2.6972788;
283            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
284            xmin = plctime(1963,10,1,0,0,2.6972788-.20);
285            xmax = plctime(1963,10,1,0,0,2.6972788+.20);
286            npts = 1001;
287            ymin = 2.55;
288            ymax = 2.75;
289            time_format = "%S%2%";
290            title_suffix = "near 1963-11-01 (TAI)";
291            xlabel_step = 0.05/(scale*86400.);
292            if (kind == 3)
293                if_TAI_time_format = 1;
294                xtitle = "Seconds (TAI)";
295            else
296                if_TAI_time_format = 0;
297                xtitle = "Seconds (TAI) labelled with corresponding UTC";
298            endif
299        elseif (kind == 5 || kind == 6)
300            ## Choose midpoint to maximize time-representation precision.
301            epoch_year  = 2009;
302            epoch_month = 0;
303            epoch_day   = 1;
304            epoch_hour  = 0;
305            epoch_min   = 0;
306            epoch_sec   = 34.;
307            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
308            xmin = plctime(2009,0,1,0,0,34.-5.);
309            xmax = plctime(2009,0,1,0,0,34.+5.);
310            npts = 1001;
311            ymin = 32.5;
312            ymax = 34.5;
313            time_format = "%S%2%";
314            title_suffix = "near 2009-01-01 (TAI)";
315            xlabel_step = 1./(scale*86400.);
316            if (kind == 5)
317                if_TAI_time_format = 1;
318                xtitle = "Seconds (TAI)";
319            else
320                if_TAI_time_format = 0;
321                xtitle = "Seconds (TAI) labelled with corresponding UTC";
322            endif
323        endif
324
325        x = xmin + (0:npts-1)*(xmax-xmin)/(npts-1);
326        tai = x;
327        for i=1:npts
328            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
329            [tai_year(i) , tai_month(i), tai_day(i), tai_hour(i), tai_min(i), tai_sec(i)] = plbtime(tai(i));
330            ## Calculate residual using tai as the epoch to nearly maximize time-representation precision.
331            plconfigtime( scale, 0., 0., 0x0, 1, tai_year(i) , tai_month(i), tai_day(i), tai_hour(i), tai_min(i), tai_sec(i));
332            ## Calculate continuous tai with new epoch.
333            tai(i) = plctime(tai_year(i), tai_month(i), tai_day(i), tai_hour(i), tai_min(i), tai_sec(i));
334            ## Calculate broken-down utc (with leap seconds inserted) from continuous tai with new epoch.
335            plconfigtime( scale, 0., 0., 0x2, 1, tai_year(i) , tai_month(i), tai_day(i), tai_hour(i), tai_min(i), tai_sec(i));
336            [utc_year(i), utc_month(i), utc_day(i), utc_hour(i), utc_min(i), utc_sec(i)] = plbtime(tai(i));
337            ## Calculate continuous utc from broken-down utc using same epoch as for the continuous tai.
338            plconfigtime( scale, 0., 0., 0x0, 1, tai_year(i) , tai_month(i), tai_day(i), tai_hour(i), tai_min(i), tai_sec(i));
339            utc(i) = plctime(utc_year(i), utc_month(i), utc_day(i), utc_hour(i), utc_min(i), utc_sec(i));
340        endfor
341        ## Convert residuals to seconds.
342        y = (tai-utc)*scale*86400.0;
343
344        pladv(0);
345        plvsta();
346        plwind(xmin, xmax, ymin, ymax);
347        plcol0(1);
348        if (if_TAI_time_format)
349            plconfigtime( scale, 0., 0., 0x0, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
350        else
351            plconfigtime( scale, 0., 0., 0x2, 1, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec );
352        endif
353        pltimefmt(time_format);
354        plbox("bcnstd", xlabel_step, 0, "bcnstv", 0., 0);
355        plcol0(3);
356        title = ["@frPLplot Example 29 - TAI-UTC ", title_suffix];
357        pllab(xtitle, "TAI-UTC (sec)", title);
358
359        plcol0(4);
360
361        plline(x', y');
362    endfor
363endfunction
364
365ix29c
366