-- Sample plots using date / time formatting for axes -- Copyright (C) 2008-2010 Jerry Bauck -- This file is part of PLplot. -- PLplot is free software; you can redistribute it and/or modify -- it under the terms of the GNU Library General Public License as published -- by the Free Software Foundation; either version 2 of the License, or -- (at your option) any later version. -- PLplot is distributed in the hope that it will be useful, -- but WITHOUT ANY WARRANTY; without even the implied warranty of -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -- GNU Library General Public License for more details. -- You should have received a copy of the GNU Library General Public License -- along with PLplot; if not, write to the Free Software -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA with Ada.Strings.Unbounded, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions, PLplot_Auxiliary, PLplot_Standard; use Ada.Strings.Unbounded, Ada.Numerics, Ada.Numerics.Long_Elementary_Functions, PLplot_Auxiliary, PLplot_Standard; -------------------------------------------------------------------------------- -- Draws several plots which demonstrate the use of date / time formats for -- the axis labels. -- Time formatting is done using the strfqsas routine from the qsastime -- library. This is similar to strftime, but works for a broad -- date range even on 32-bit systems. See the -- documentation of strfqsas for full details of the available formats. -- -- 1) Plotting temperature over a day (using hours / minutes) -- 2) Plotting -- -------------------------------------------------------------------------------- -- NOTE: The Ada user is reminded that Ada.Calendar is very capable and complete. -- See especially Time_Of. procedure xstandard29a is -- Plot a model diurnal cycle of temperature procedure plot1 is x, y : Real_Vector(0 .. 72); xerr1, xerr2, yerr1, yerr2 : Real_Vector(0 .. 72); -- Data points every 10 minutes for 1 day xmin, xmax, ymin, ymax : Long_Float; begin xmin := 0.0; xmax := 60.0 * 60.0 * 24.0; -- Number of seconds in a day ymin := 10.0; ymax := 20.0; for i in x'range loop x(i) := xmax * Long_Float(i) / Long_Float(x'length); y(i) := 15.0 - 5.0 * cos( 2.0 * pi * Long_Float(i) / Long_Float(x'length)); -- Set x error bars to +/- 5 minute xerr1(i) := x(i) - Long_Float(60 * 5); xerr2(i) := x(i) + Long_Float(60 * 5); -- Set y error bars to +/- 0.1 deg C */ yerr1(i) := y(i) - 0.1; yerr2(i) := y(i) + 0.1; end loop; Advance_To_Subpage(Next_Subpage); -- Rescale major ticks marks by 0.5 Set_Major_Tick_Length(0.0, 0.5); -- Rescale minor ticks and error bar marks by 0.5 Set_Minor_Tick_Length(0.0, 0.5); Set_Viewport_Standard; Set_Viewport_World(xmin, xmax, ymin, ymax); -- Draw a box with ticks spaced every 3 hours in X and 1 degree C in Y. Set_Pen_Color(Red); -- Set time format to be hours:minutes Set_Date_Time_Label_Format("%H:%M"); Box_Around_Viewport("bcnstd", 3.0 * 60.0 * 60.0, 3, "bcnstv", 1.0, 5); Set_Pen_Color(Green); Write_Labels("Time (hours:mins)", "Temperature (degC)", "@frPLplot Example 29 - Daily temperature"); Set_Pen_Color(Aquamarine); Draw_Curve(x, y); Set_Pen_Color(Yellow); Draw_Error_Bars_X(xerr1, xerr2, y); Set_Pen_Color(Green); Draw_Error_Bars_Y(x, yerr1, yerr2); -- Rescale major / minor tick marks back to default Set_Minor_Tick_Length(0.0, 1.0); Set_Major_Tick_Length(0.0, 1.0); end plot1; -- Plot the number of hours of daylight as a function of day for a year procedure plot2 is xmin, xmax, ymin, ymax : Long_Float; lat, p, d : Long_Float; x, y : Real_Vector(0 .. 364); begin -- Latitude for London lat := 51.5; xmin := 0.0; xmax := Long_Float(x'length) * 60.0 * 60.0 * 24.0; ymin := 0.0; ymax := 24.0; -- Formula for hours of daylight from -- "A Model Comparison for Daylength as a Function of Latitude and -- Day of the Year", 1995, Ecological Modelling, 80, pp 87-95. for j in x'range loop x(j):= Long_Float(j) * 60.0 * 60.0 * 24.0; p := arcsin(0.39795 * cos(0.2163108 + 2.0 * arctan(0.9671396 * tan(0.00860 * Long_Float(j-186))))); d := 24.0 - (24.0 / pi) * arccos((sin(0.8333 * pi / 180.0) + sin(lat * pi / 180.0) * sin(p)) / (cos(lat * pi / 180.0) * cos(p))); y(j):= d; end loop; Set_Pen_Color(Red); -- Set time format to be abbreviated month name followed by day of month Set_Date_Time_Label_Format("%b %d"); Set_Numeric_Label_Precision(True, 1); Set_Environment(xmin, xmax, ymin, ymax, 0, 40); Set_Pen_Color(Green); Write_Labels("Date", "Hours of daylight", "@frPLplot Example 29 - Hours of daylight at 51.5N"); Set_Pen_Color(Aquamarine); Draw_Curve(x, y); Set_Numeric_Label_Precision(False, 0); end plot2; procedure plot3 is xmin, xmax, ymin, ymax, tstart : Long_Float; x, y : Real_Vector(0 .. 61); begin -- Calculate continuous time corresponding to 2005-12-01 UTC. Continuous_From_Broken_Down_Time(2005, 11, 01, 0, 0, 0.0, tstart); xmin := tstart; xmax := xmin + Long_Float(x'length) * 60.0 * 60.0 * 24.0; ymin := 0.0; ymax := 5.0; for i in x'range loop x(i) := xmin + Long_Float(i) * 60.0 * 60.0 * 24.0; y(i) := 1.0 + sin(2.0 * pi * Long_Float(i) / 7.0) + exp((Long_Float(Integer'min(i, x'length - i))) / 31.0); end loop; Advance_To_Subpage(Next_Subpage); Set_Viewport_Standard; Set_Viewport_World(xmin, xmax, ymin, ymax); Set_Pen_Color(Red); -- Set time format to be ISO 8601 standard YYYY-MM-DD. Set_Date_Time_Label_Format("%F"); -- Draw a box with ticks spaced every 14 days in X and 1 hour in Y. Box_Around_Viewport("bcnstd", 14.0 * 24.0 * 60.0 * 60.0, 14, "bcnstv", 1.0, 4); Set_Pen_Color(Green); Write_Labels("Date", "Hours of television watched", "@frPLplot Example 29 - Hours of television watched in Dec 2005 / Jan 2006"); Set_Pen_Color(Aquamarine); -- Rescale symbol size (used by Draw_Points) by 0.5 Set_Symbol_Size(0.0,0.5); Draw_Points(x, y, 2); Draw_Curve(x, y); end plot3; procedure plot4 is -- TAI-UTC (seconds) as a function of time. Scale : Long_Float; xmin, xmax, ymin, ymax, xlabel_step : Long_Float; npts : Integer; if_TAI_time_format : Boolean; time_format : Unbounded_String := To_Unbounded_String(""); title_suffix : Unbounded_String := To_Unbounded_String(""); xtitle : Unbounded_String := To_Unbounded_String(""); title : Unbounded_String := To_Unbounded_String(""); x, y : Real_Vector(0 .. 1000); epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min : Integer; epoch_sec : Long_Float; tai_year, tai_month, tai_day, tai_hour, tai_min : Integer; tai_sec, tai : Long_Float; utc_year, utc_month, utc_day, utc_hour, utc_min : Integer; utc_sec, utc : Long_Float; begin -- Continuous time unit is Besselian years from whatever epoch is -- chosen below. Could change to seconds (or days) from the -- epoch, but then would have to adjust xlabel_step below. scale := 365.242198781; -- MJD epoch (see <https://en.wikipedia.org/wiki/Julian_day>). -- This is only set for illustrative purposes, and is overwritten -- below for the time-representation reasons given in the -- discussion below. epoch_year := 1858; epoch_month := 11; epoch_day := 17; epoch_hour := 0; epoch_min := 0; epoch_sec := 0.0; -- To illustrate the time-representation issues of using the MJD -- epoch, in 1985, MJD was roughly 46000 days which corresponds to -- 4e9 seconds. Thus, for the -DPL_DOUBLE=ON case where PLFLT is -- a double which can represent continuous time to roughly 16 -- decimal digits of precision the time-representation error is -- roughly ~400 nanoseconds. Therefore the MJD epoch would be -- acceptable for the plots below in the -DPL_DOUBLE=ON case. -- However, that epoch is obviously not acceptable for the -- -DPL_DOUBLE=OFF case where PLFLT is a float which can represent -- continuous time to only ~7 decimal digits of precision -- corresponding to a time representation error of 400 seconds (!) -- in 1985. For this reason, we do not use the MJD epoch below -- and instead choose the best epoch for each case to minimize -- time-representation issues. for kind in 0 .. 6 loop if kind = 0 then -- Choose midpoint to maximize time-representation precision. epoch_year := 1985; epoch_month := 0; epoch_day := 2; epoch_hour := 0; epoch_min := 0; epoch_sec := 0.0; Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); Continuous_From_Broken_Down_Time(1950, 0, 2, 0, 0, 0.0, xmin); Continuous_From_Broken_Down_Time(2020, 0, 2, 0, 0, 0.0, xmax); npts := 70 * 12 + 1; ymin := 0.0; ymax := 36.0; time_format := To_Unbounded_String("%Y%"); if_TAI_time_format := True; title_suffix := To_Unbounded_String("from 1950 to 2020"); xtitle := To_Unbounded_String("Year"); xlabel_step := 10.0; elsif kind = 1 or kind = 2 then -- Choose midpoint to maximize time-representation precision. epoch_year := 1961; epoch_month := 7; epoch_day := 1; epoch_hour := 0; epoch_min := 0; epoch_sec := 1.64757; Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); Continuous_From_Broken_Down_Time(1961, 7, 1, 0, 0, 1.64757 - 0.20, xmin); Continuous_From_Broken_Down_Time(1961, 7, 1, 0, 0, 1.64757 + 0.20, xmax); npts := 1001; ymin := 1.625; ymax := 1.725; time_format := To_Unbounded_String("%S%2%"); title_suffix := To_Unbounded_String("near 1961-08-01 (TAI)"); xlabel_step := 0.05 / (scale * 86400.0); if kind = 1 then if_TAI_time_format := True; xtitle := To_Unbounded_String("Seconds (TAI)"); else if_TAI_time_format := False; xtitle := To_Unbounded_String("Seconds (TAI) labelled with corresponding UTC"); end if; elsif kind = 3 or kind = 4 then -- Choose midpoint to maximize time-representation precision. epoch_year := 1963; epoch_month := 10; epoch_day := 1; epoch_hour := 0; epoch_min := 0; epoch_sec := 2.6972788; Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); Continuous_From_Broken_Down_Time(1963, 10, 1, 0, 0, 2.6972788 - 0.20, xmin); Continuous_From_Broken_Down_Time(1963, 10, 1, 0, 0, 2.6972788 + 0.20, xmax); npts := 1001; ymin := 2.55; ymax := 2.75; time_format := To_Unbounded_String("%S%2%"); title_suffix := To_Unbounded_String("near 1963-11-01 (TAI)"); xlabel_step := 0.05 / (scale * 86400.0); if kind = 3 then if_TAI_time_format := True; xtitle := To_Unbounded_String("Seconds (TAI)"); else if_TAI_time_format := False; xtitle := To_Unbounded_String("Seconds (TAI) labelled with corresponding UTC"); end if; elsif kind = 5 or kind = 6 then -- Choose midpoint to maximize time-representation precision. epoch_year := 2009; epoch_month := 0; epoch_day := 1; epoch_hour := 0; epoch_min := 0; epoch_sec := 34.0; Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); Continuous_From_Broken_Down_Time(2009, 0, 1, 0, 0, 34.0 - 5.0, xmin); Continuous_From_Broken_Down_Time(2009, 0, 1, 0, 0, 34.0 + 5.0, xmax); npts := 1001; ymin := 32.5; ymax := 34.5; time_format := To_Unbounded_String("%S%2%"); title_suffix := To_Unbounded_String("near 2009-01-01 (TAI)"); xlabel_step := 1.0 / (scale * 86400.0); if kind = 5 then if_TAI_time_format := True; xtitle := To_Unbounded_String("Seconds (TAI)"); else if_TAI_time_format := False; xtitle := To_Unbounded_String("Seconds (TAI) labelled with corresponding UTC"); end if; end if; for i in 0 .. npts - 1 loop x(i) := xmin + Long_Float(i) * (xmax - xmin) / (Long_Float(npts - 1)); tai := x(i); Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); Broken_Down_From_Continuous_Time(tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec, tai); -- Calculate residual using tai as the epoch to nearly maximize time-representation precision. Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec) ; -- Calculate continuous tai with new epoch. Continuous_From_Broken_Down_Time(tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec, tai); -- Calculate broken-down utc (with leap seconds inserted) from continuous tai with new epoch. Configure_Time_Transformation(scale, 0.0, 0.0, 2, True, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec) ; Broken_Down_From_Continuous_Time(utc_year, utc_month, utc_day, utc_hour, utc_min, utc_sec, tai); -- Calculate continuous utc from broken-down utc using same epoch as for the continuous tai. Configure_Time_Transformation(scale, 0.0, 0.0, 2, True, tai_year, tai_month, tai_day, tai_hour, tai_min, tai_sec) ; Continuous_From_Broken_Down_Time(utc_year, utc_month, utc_day, utc_hour, utc_min, utc_sec, utc); -- Convert residuals to seconds. y(i) := (tai - utc) * scale * 86400.0; end loop; Advance_To_Subpage(Next_Subpage); Set_Viewport_Standard; Set_Viewport_World(xmin, xmax, ymin, ymax); Set_Pen_Color(Red); if if_TAI_time_format then Configure_Time_Transformation(scale, 0.0, 0.0, 0, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); else Configure_Time_Transformation(scale, 0.0, 0.0, 2, True, epoch_year, epoch_month, epoch_day, epoch_hour, epoch_min, epoch_sec); end if; Set_Date_Time_Label_Format(To_String(time_format)); Box_Around_Viewport("bcnstd", xlabel_step, 0, "bcnstv", 0.0, 0); Set_Pen_Color(Green); title := To_Unbounded_String("@frPLplot Example 29 - TAI-UTC "); title := title & title_suffix; Write_Labels(To_String(xtitle), "TAI-UTC (sec)", To_String(title)); Set_Pen_Color(Aquamarine); if kind = 0 then -- Shorter x and y Draw_Curve(x(0 .. 70 * 12), y(0 .. 70 * 12)); else -- Longer x and y Draw_Curve(x, y); end if; end loop; -- kind end plot4; begin -- main -- Parse command line arguments Parse_Command_Line_Arguments(Parse_Full); -- Change the escape character to a '@' instead of the default '#' Set_Escape_Character('@'); -- Initialize plplot Initialize_PLplot; plot1; plot2; plot3; plot4; -- Don't forget to call End_PLplot to finish off! End_PLplot; end xstandard29a;