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