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