Root/nanonote-files/example-files/data/Examples/lua-plplot-examples/x29.lua

1--[[ $Id: x29.lua 10299 2009-08-19 17:54:27Z smekal $
2
3       Sample plots using date / time formatting for axes
4
5  Copyright (C) 2009 Werner Smekal
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 General Library 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
25-- initialise Lua bindings for PLplot examples.
26dofile("plplot_examples.lua")
27
28
29
30-- Plot a model diurnal cycle of temperature
31function plot1()
32  x = {}
33  y = {}
34
35  xerr1 = {}
36  xerr2 = {}
37  yerr1 = {}
38  yerr2 = {}
39  -- Data points every 10 minutes for 1 day
40  npts = 73
41
42  xmin = 0
43  xmax = 60*60*24 -- Number of seconds in a day
44  ymin = 10
45  ymax = 20
46
47  for i = 1, npts do
48    x[i] = xmax*((i-1)/npts)
49    y[i] = 15 - 5*math.cos(2*math.pi*((i-1)/npts))
50    -- Set x error bars to +/- 5 minute
51    xerr1[i] = x[i]-60*5
52    xerr2[i] = x[i]+60*5
53    -- Set y error bars to +/- 0.1 deg C
54    yerr1[i] = y[i]-0.1
55    yerr2[i] = y[i]+0.1
56  end
57  
58  pl.adv(0)
59
60  -- Rescale major ticks marks by 0.5
61  pl.smaj(0, 0.5)
62  -- Rescale minor ticks and error bar marks by 0.5
63  pl.smin(0, 0.5)
64
65  pl.vsta()
66  pl.wind(xmin, xmax, ymin, ymax)
67
68  -- Draw a box with ticks spaced every 3 hour in X and 1 degree C in Y.
69  pl.col0(1)
70  -- Set time format to be hours:minutes
71  pl.timefmt("%H:%M")
72  pl.box("bcnstd", 3*60*60, 3, "bcnstv", 1, 5)
73
74  pl.col0(3)
75  pl.lab("Time (hours:mins)", "Temperature (degC)", "@frPLplot Example 29 - Daily temperature")
76  
77  pl.col0(4)
78
79  pl.line(x, y)
80  pl.col0(2)
81  pl.errx(xerr1, xerr2, y)
82  pl.col0(3)
83  pl.erry(x, yerr1, yerr2)
84
85  -- Rescale major / minor tick marks back to default
86  pl.smin(0, 1)
87  pl.smaj(0, 1)
88end
89
90
91-- Plot the number of hours of daylight as a function of day for a year
92function plot2()
93  x = {}
94  y = {}
95
96  -- Latitude for London
97  lat = 51.5
98
99  npts = 365
100
101  xmin = 0
102  xmax = npts*60*60*24
103  ymin = 0
104  ymax = 24
105  
106  -- Formula for hours of daylight from
107  -- "A Model Comparison for Daylength as a Function of Latitude and
108  -- Day of the Year", 1995, Ecological Modelling, 80, pp 87-95.
109  for j = 1, npts do
110    x[j] = (j-1)*60*60*24
111    p = math.asin(0.39795*math.cos(0.2163108 + 2*math.atan(0.9671396*math.tan(0.00860*(j-1-186)))))
112    d = 24 - (24/math.pi)*
113      math.acos( (math.sin(0.8333*math.pi/180) + math.sin(lat*math.pi/180)*math.sin(p)) /
114        (math.cos(lat*math.pi/180)*math.cos(p)) )
115    y[j] = d
116  end
117
118  pl.col0(1)
119  -- Set time format to be abbreviated month name followed by day of month
120  pl.timefmt("%b %d")
121  pl.prec(1, 1)
122  pl.env(xmin, xmax, ymin, ymax, 0, 40)
123
124
125  pl.col0(3)
126  pl.lab("Date", "Hours of daylight", "@frPLplot Example 29 - Hours of daylight at 51.5N")
127  
128  pl.col0(4)
129
130  pl.line(x, y)
131  
132  pl.prec(0, 0)
133end
134
135
136function plot3()
137  x = {}
138  y = {}
139
140  tstart = 1133395200
141   
142  npts = 62
143
144  xmin = tstart
145  xmax = xmin + npts*60*60*24
146  ymin = 0
147  ymax = 5
148  
149  for i = 1, npts do
150    x[i] = xmin + (i-1)*60*60*24
151    y[i] = 1 + math.sin(2*math.pi*(i-1)/7) + math.exp(math.min(i-1,npts-i+1)/31)
152  end
153  pl.adv(0)
154
155  pl.vsta()
156  pl.wind(xmin, xmax, ymin, ymax)
157
158  pl.col0(1)
159  -- Set time format to be ISO 8601 standard YYYY-MM-DD. Note that this is
160  --equivalent to %f for C99 compliant implementations of strftime.
161  pl.timefmt("%Y-%m-%d")
162  -- Draw a box with ticks spaced every 14 days in X and 1 hour in Y.
163  pl.box("bcnstd", 14*24*60*60,14, "bcnstv", 1, 4)
164
165  pl.col0(3)
166  pl.lab("Date", "Hours of television watched", "@frPLplot Example 29 - Hours of television watched in Dec 2005 / Jan 2006")
167  
168  pl.col0(4)
169
170  -- Rescale symbol size (used by plpoin) by 0.5
171  pl.ssym(0, 0.5)
172  pl.poin(x, y, 2)
173  pl.line(x, y)
174end
175
176
177function plot4()
178  -- TAI-UTC (seconds) as a function of time.
179  -- Use Besselian epochs as the continuous time interval just to prove
180  -- this does not introduce any issues.
181  
182  x = {}
183  y = {}
184
185  -- Use the definition given in http://en.wikipedia.org/wiki/Besselian_epoch
186  -- B = 1900. + (JD -2415020.31352)/365.242198781
187  -- ==> (as calculated with aid of "bc -l" command)
188  -- B = (MJD + 678940.364163900)/365.242198781
189  -- ==>
190  -- MJD = B*365.24219878 - 678940.364163900
191  scale = 365.242198781
192  offset1 = -678940
193  offset2 = -0.3641639
194  pl.configtime(scale, offset1, offset2, 0, 0, 0, 0, 0, 0, 0, 0.)
195
196  for kind = 0, 6 do
197    if kind == 0 then
198      xmin = pl.ctime(1950, 0, 2, 0, 0, 0)
199      xmax = pl.ctime(2020, 0, 2, 0, 0, 0)
200      npts = 70*12 + 1
201      ymin = 0
202      ymax = 36
203      time_format = "%Y%"
204      if_TAI_time_format = 1
205      title_suffix = "from 1950 to 2020"
206      xtitle = "Year"
207      xlabel_step = 10
208    end
209    if kind==1 or kind==2 then
210      xmin = pl.ctime(1961, 7, 1, 0, 0, 1.64757-0.20)
211      xmax = pl.ctime(1961, 7, 1, 0, 0, 1.64757+0.20)
212      npts = 1001
213      ymin = 1.625
214      ymax = 1.725
215      time_format = "%S%2%"
216      title_suffix = "near 1961-08-01 (TAI)"
217      xlabel_step = 0.05/(scale*86400)
218      if kind==1 then
219        if_TAI_time_format = 1
220        xtitle = "Seconds (TAI)"
221      else
222        if_TAI_time_format = 0
223        xtitle = "Seconds (TAI) labelled with corresponding UTC"
224      end
225    end
226    if kind==3 or kind==4 then
227      xmin = pl.ctime(1963, 10, 1, 0, 0, 2.6972788-0.20)
228      xmax = pl.ctime(1963, 10, 1, 0, 0, 2.6972788+0.20)
229      npts = 1001
230      ymin = 2.55
231      ymax = 2.75
232      time_format = "%S%2%"
233      title_suffix = "near 1963-11-01 (TAI)"
234      xlabel_step = 0.05/(scale*86400)
235      if kind==3 then
236                if_TAI_time_format = 1
237                xtitle = "Seconds (TAI)"
238      else
239                if_TAI_time_format = 0
240                xtitle = "Seconds (TAI) labelled with corresponding UTC"
241      end
242    end
243    if kind==5 or kind==6 then
244      xmin = pl.ctime(2009, 0, 1, 0, 0, 34-5)
245      xmax = pl.ctime(2009, 0, 1, 0, 0, 34+5)
246      npts = 1001
247      ymin = 32.5
248      ymax = 34.5
249      time_format = "%S%2%"
250      title_suffix = "near 2009-01-01 (TAI)"
251      xlabel_step = 1/(scale*86400)
252      if kind==5 then
253        if_TAI_time_format = 1
254        xtitle = "Seconds (TAI)"
255      else
256        if_TAI_time_format = 0
257        xtitle = "Seconds (TAI) labelled with corresponding UTC"
258      end
259    end
260
261    for i = 1, npts do
262      x[i] = xmin + (i-1)*(xmax-xmin)/(npts-1)
263      pl.configtime(scale, offset1, offset2, 0, 0, 0, 0, 0, 0, 0, 0)
264      tai = x[i]
265      tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec = pl.btime(tai)
266      pl.configtime(scale, offset1, offset2, 2, 0, 0, 0, 0, 0, 0, 0)
267      utc_year, utc_month, utc_day, utc_hour, utc_min, utc_sec = pl.btime(tai)
268      pl.configtime(scale, offset1, offset2, 0, 0, 0, 0, 0, 0, 0, 0.)
269      utc = pl.ctime(utc_year, utc_month, utc_day, utc_hour, utc_min, utc_sec)
270      y[i]=(tai-utc)*scale*86400.
271    end
272
273    pl.adv(0)
274    pl.vsta()
275    pl.wind(xmin, xmax, ymin, ymax)
276    pl.col0(1)
277    if if_TAI_time_format ~= 0 then
278      pl.configtime(scale, offset1, offset2, 0, 0, 0, 0, 0, 0, 0, 0)
279    else
280      pl.configtime(scale, offset1, offset2, 2, 0, 0, 0, 0, 0, 0, 0)
281    end
282    pl.timefmt(time_format)
283    pl.box("bcnstd", xlabel_step, 0, "bcnstv", 0., 0)
284    pl.col0(3)
285    title = "@frPLplot Example 29 - TAI-UTC " .. title_suffix
286    pl.lab(xtitle, "TAI-UTC (sec)", title)
287    
288    pl.col0(4)
289    
290    pl.line(x, y)
291  end
292end
293
294----------------------------------------------------------------------------
295-- main
296--
297-- Draws several plots which demonstrate the use of date / time formats for
298-- the axis labels.
299-- Time formatting is done using the system strftime routine. See the
300-- documentation of this for full details of the available formats.
301--
302-- 1) Plotting temperature over a day (using hours / minutes)
303-- 2) Plotting
304--
305-- Note: Times are stored as seconds since the epoch (usually 1st Jan 1970).
306--
307----------------------------------------------------------------------------
308
309-- Parse command line arguments
310pl.parseopts(arg, pl.PL_PARSE_FULL)
311
312-- Initialize plplot
313pl.init()
314
315-- Change the escape character to a '@' instead of the default '#'
316pl.sesc('@')
317
318plot1()
319
320plot2()
321
322plot3()
323
324plot4()
325
326-- Don't forget to call plend() to finish off!
327pl.plend()
328

Archive Download this file

Branches:
backfire
master
release_2010-11-17
release_2010-12-14
release_2011-02-23
release_2011-05-28
release_2011-08-27
release_2011-11-13
release_2012-03-18
release_2012-04-09
release_2012-07-11
release_2012-10-11



interactive