--- Begin Message ---
Hi Alan,
Many thanks for the quick email. Apart from other reasons, one of my
motivations for posting the patch was to see if there was a receptive
plplot community out there. As I indicated, we are in the process of
moving a large project from pgplot to plplot and I wanted to reassure
myself that this was a sound investment of effort. So you passed that
test! You will probably see more profound offerings and queries from the
programmers who are responsible for the clean, professional bits of
code; I'm the hacker who wants to use the software, and who fronts the
funding proposals.
I haven't looked at your svn trunk for x29c.c - is that something at the
front of an elephant? I hadn't spent a lot of time playing with the C
time routines before a week or two ago when I started this, because
we've used our own based on the NASA NSSDC CDF ones for all our internal
work. Pgplot as you know I guess is written in Fortran. It's actually
quite clever at handling time labelling, but it works in seconds of the
day; this avoids most of the issues.
I can see the logic and ease of relying on the C routines to do most of
the work, indcluding strftime to deal with the labels. I thought I had
the tricks to handle the well-documented quirks of mktime, but then hit
them again by accident. They behaved differently on different platforms
(windows vs linux, and even different flavours of linux), and different
depending on the time zone, summer/winter dst months, etc.
Given how central time is to our work, I just lost all confidence in
mktime and its portability. Maybe you could demonstrate some suitable
trick, from your svn trunk, that appears to be robust, but my
frustration over the past two weeks suggests to me that if it calls
mktime then I'm reluctant to have 100% confidence in it.
When you do get around to looking at my code, you'll realise that
there's almost nothing clever in it at all apart from the calculation of
the number of days since 0 AD, for which I take zero credit. Given that
gmtime and strftime will be used, it's easy to hold enough accuracy in a
double to count seconds; the CDF routines actually work in msecs, which
is about as much accuracy as can be squeezed. They have new ways of
dealing with finer time resolution that, like yours, divide the
representation across two pieces to get down below nanosecond accuracy I
believe. A purist would do all this in integer arithmetic (you'll see
that my code counts floating point seconds and then truncates), but
double arithmetic is very good these days and, as I noted above, the
accuracy needed here is a long way from representation errors.
There are other things I'd think about "fixing" in this area, but they
are harder. For example, the automatic tick spacing merely calculates an
appropriate scaling factor (to translate to minutes, hours, days, weeks,
etc.) and then uses the standard (quite simple, elegant, and otherwise
robust) plplot algorithms to find "nice" starting points and (sub)tick
divisions. This results, for example, in a time axis being divided into
intervals of 50 seconds (with 5 subticks), and hence with labels that
don't fall on natural minute boundaries. Similar things happen with 50
minutes, and dividing a week into 5 sub-units is quite un-nerving. For
our own purposes we may need to address some of these issues, but it
might well be easier to do so in a layer that sits on top of plplot
(i.e., calculate specific non-zero entries for pltbox) rather than
re-works the internals of pldtik. If you want to play, I attach a little
test programme that will let you enter intervals of time and watch the
labelling - which I've massaged a bit to avoid the places the default
behaviour seemed to be the worst [but this is hard to do in a generic
way, e.g., for small plots, large font sizes, etc.]
As I indicated, you're likely to see more stuff from my programmer,
Alban Rochel, who has been digging deep into the X-window and wingcc
drivers. There are issues here in terms of handling multiple open plot
windows, and being able to trap a user clicking a close box (which in
the present versions result in a broken pipe, abort, or leaving the
stream unreachable. We have an additional hurdle in that our already
large application, which uses Qt widgets, relies on gnu autoconfigure
tools, so we're working on integrating at least a sub-set of plplot
within that environment rather than force our end-users to invest in
coaxing cmake to work for them. The less we have to do to the official
plplot source distribution the less maintenance this solution will be
for us, so we have a vested interest in working with you to incorporate
features we need as part of the core distribution. I do hope that will
also result in benefits to the plplot project, either in terms of
enhanced features or enhanced robustness against platform and software
environment flavours.
Watch this space...
Again my thanks for your reply.
Yours
Steve
NOTE: The attached trials aren't the prettiest code, and the helper
routines aren't quite the final ones I submitted in my patch, but for
the purposes of watching the labelling algorithms they suffice. It's
based on x29c.c, and I guess you know how to build it on your system so
there's not much point sending my hand-crafted (=hacked) makefile. I've
run this for time ranges of 0.1 secs (for which I break out of the
time-handling and treat the x-axis as a normal variable) to 10^9 secs.
NOTE 2: None of this addresses the pre-1970 date issues that arises on
some platforms - like windows/mingw for example - that don'l like
negative time_t's. I think basically, since ther is no ANSI standard for
the time_t t=0 epoch, this is virually impossible. If anyone thinks this
is really an issue, the only sensible solution is to abandon the C-time
handling and adopt another. We are quite versed in the CDF routines
(less so in the newer variants that use two separate variables to get
down to nanosecs, but still knowledgeable). Since there are virtually no
space missions pre-1970, this is not a priority for us :-)
On Sun, 2008-08-03 at 10:27 -0700, Alan W. Irwin wrote:
> On 2008-08-03 15:56+0100 Steve Schwartz wrote:
>
> > Dear Plplot-devel,
> >
> > I attach a patch that immunizes plplot from the peculiarities of C's
> > mktime function, that arise because mktime requires input in localtime.
> > The various recipes to calculate offsets from local time to utc, such as
> > those in example x29, are not foolproof or portable. I have not explored
> > the equivalent timegm function that I've seen in one of the plplot-devel
> > threads, but some of them still rely on mktime and then attempt to
> > adjust the result.
>
> Hi Steve:
>
> Thank you for your interest in improving PLplot's treatment of time. I was
> involved in time research in the late 90's (calculation of time ephemerides
> for the Earth where I was aiming for sub ns numerical accuracy over
> millenia) so I have an interest in making sure PLplot's treatment of time is
> the best possible. I fully agree with your comments about the peculiarities
> of local time so it should be avoided completely.
>
> Have you looked at our recent svn trunk revision 8530 of x29c.c which
> temporarily sets TZ to the null string to avoid those peculiarities? One
> concern with that approach is what to do for Windows platforms, but do you
> think there would be any other concerns? Note, we haven't yet tried that TZ
> approach for the calls to mktime in src/pldtik.c, but that is one possible
> approach to avoid the peculiarities of local time there. Thus, your comments
> on applying the current approach for example 29 to our core library would be
> appreciated. We may end up just adopting your patch (for which many thanks,
> by the way), but I want to make sure we are not missing a simpler solution
> to the problem.
>
> Ephemeris calculations have to contend with the issue of the best way to
> represent time over millenia while minimizing rounding errors. IIRC, the
> JPL ephemeris (upon which my time ephemeris research was based) solved the
> problem by representing time (in Julian days) as two 64-bit floating point
> numbers. The first one was exact since it represented the integer portion
> of the Julian date, and the second represented the fractional part to
> something like 1 part in 10^17 which translates to an error of roughly 0.001
> ns at worst.
>
> I haven't had a chance to look at your patch yet, but I am curious about how
> you solved this same "time representation" problem.
>
> Alan
> __________________________
> Alan W. Irwin
>
> Astronomical research affiliation with Department of Physics and Astronomy,
> University of Victoria (astrowww.phys.uvic.ca).
>
> Programming affiliations with the FreeEOS equation-of-state implementation
> for stellar interiors (freeeos.sf.net); PLplot scientific plotting software
> package (plplot.org); the libLASi project (unifont.org/lasi); the Loads of
> Linux Links project (loll.sf.net); and the Linux Brochure Project
> (lbproject.sf.net).
> __________________________
>
> Linux-powered Science
> __________________________
--
+-------------------------------------------------------------------+
Professor Steven J Schwartz Phone: +44-(0)20-7594-7660
Space and Atmospheric Physics Fax: +44-(0)20-7594-7772
The Blackett Laboratory E-mail: [EMAIL PROTECTED]
Imperial College London Office: Huxley 6M70
London SW7 2BW, U.K. Web: http://www.sp.ph.ic.ac.uk/~sjs
+-------------------------------------------------------------------+
/* To DO:
calculate offset for t in case is_t is false; ok down to dt = few msecs
limit to date ~ 2050?
*/
#include "plcdemos.h"
#include "qsas_plplot_time_helpers.c"
/* Function prototypes */
void plot1();
void plot2();
void plot3();
void plot4();
void mkt(PLFLT tmin, PLFLT tmax, int npts, PLFLT *t);
double toffset(struct tm *tm);
/*--------------------------------------------------------------------------*\
* main
*
* Draws several plots which demonstrate the use of date / time formats for
* the axis labels.
* Time formatting is done using the system strftime routine. See the
* documentation of this for full details of the available formats.
*
* 1) Plotting temperature over a day (using hours / minutes)
* 2) Plotting
*
* Note: Times are stored as seconds since the epoch (usually 1st Jan 1970).
*
\*--------------------------------------------------------------------------*/
#include <stdio.h>
int
main(int argc, const char *argv[])
{
/* Parse command line arguments */
plparseopts(&argc, argv, PL_PARSE_FULL);
/* Initialize plplot */
/* strm 0 is null, so can't kill application by close button */
plsdev("null");
plinit(); /* initialise - once only ?*/
plsstrm(1);
plspause(0); /* no terminal input between pages */
plsdev("tk");
plscolbg(255,255,255);
plinit(); /* see example x14, which does it this way */
plot2();
plflush();
printf("hit any key");
getchar();
plend();
exit(0);
}
/* time handling trials */
void mkt(PLFLT tmin, PLFLT tmax, int npts, PLFLT *t){
PLFLT dt;
int i;
dt = (tmax - tmin)/(PLFLT)(npts-1);
for (i=0; i<npts; i++) t[i] = tmin + i*dt;
}
void
plot2()
{
int i, npts=100;
PLFLT tmin, tmax, ymin, ymax, t[100],y[100];
PLINT ntsub;
time_t tstart, tend;
/* time_t are secs since epoch, e.g., 1970
tm structs are broken down time since 1900! */
double toff, dt, absdt;
struct tm tm1;
char timestr[50], datefmt[20],labfmt[20];
int is_t;
ymin = 0.0;
ymax = 5.0;
for (i = 0; i<npts; i++) {
y[i] = 2.5 + sin(3.14159*i/25.);
}
tm1.tm_year = 108; /* Years since 1900 */
tm1.tm_mon = 6; /* 0 == January, 11 = December */
tm1.tm_mday = 2; /* 1 = 1st of month */
tm1.tm_hour = 3;
tm1.tm_min = 4;
tm1.tm_sec = 5;
tmin = mktimeutc(&tm1);
strftime(timestr,49,"%Y-%m-%d %H:%M:%S", &tm1);
printf("tm1 = %s\n",timestr);
/* ok now let's try some different time intervals */
plcol0(1);
plssub(1,5);
while (1<2){
printf("Input dt in secs ");
scanf("%lg",&dt);
tmax = tmin + dt;
mkt(tmin, tmax, npts, t); /* fill t array */
pladv(0);
plvsta();
plwind(t[0], t[npts-1], ymin, ymax);
/* use tmin and tmax to work out format for labels */
timefmts(tmin, tmax, labfmt, datefmt, &is_t);
pltimefmt(labfmt);
if(is_t) {
/* Attempt to specify sensible sub-tick ranges. This isn't perfect as it takes
no account of the major ticks, but is better than default for some ranges,
notably where plplot counts by weeks. Since plplot simply scales the time
to be mins, hours, etc. and then uses its auto-tick routines it will happily
plot in major tick intervals of 50 mins, which is pretty terrible. The pgplot
algorithms are better, but much more complicated.
*/
absdt = abs(tmax-tmin);
if (absdt < 10800) {ntsub = 0;} /* secs & mins ~ ok except 50 mins yuk */
else if (absdt < 259200) {ntsub = 0;} /* hrs ok except 10, 20 hrs yuk */
else if (absdt < 1200000) {ntsub = 4;} /* days split 4 */
else if (absdt < 1814400) {ntsub = 0;} /* several days */
else if (absdt < 94608000){ntsub = 7;} /* plplot counting in weeks */
else {ntsub = 0;}
plbox("bcnstd", 0,ntsub, "bcnstv", 0,0);
} else { /* range too small; convert to seconds */
if (abs(floor(tmax)-floor(tmin)) < 0.1){ /* this is = 0 test */
toff = floor(tmin);
} else { /* use secs in minute unless sub-second range */
toff = 60.*floor(tmin/60);
}
for(i=0;i<npts; i++) t[i] -= toff;
plwind(t[0],t[npts-1],ymin, ymax);
plbox("bcnst", 0,0, "bcnstv", 0,0);
}
strftime(timestr,49,"%Y-%m-%d %H:%M:%S", &tm1);
pllab("UT", "Y variable", timestr);
tstart = tmin; tend = tmax; /* convert to time_t */
strftime(timestr,49,datefmt, gmtime(&tstart));
plmtex("b",3.0,0.0,0.0,timestr);
if(tend != tstart){
/* this avoids overlap with exponent in rare cases where sub-second
range used at beginning of a minute */
strftime(timestr,49,datefmt, gmtime(&tend));
plmtex("b",3.0,1.0,1.0,timestr);
}
plline(npts, t, y);
plflush();
}
}
/* plplot time stuff
* sjs 31 July 2008 */
/* Included in this file are routines for mktimeutc, which is a fixed
alternative to the C mktime and uses utc. Note that time_t has a Year2038
"bug" like the Y2k one, as it hits the limit of what can be held in a
32-bit int at that time. See the wikipedia page.
Also included is a routine to set up the time format strings that will
get used to format time labels, and also to insert the date stamp at
either end of the time axis. These are based on some experimentation but
might not be appropriate in all cases. There are some oddities in plplot
time handling in that it doesn't know, for example, that hours are better
divided by 4 or 6 and not 5. It also counts in weeks for some ranges,
which is less than ideal.
It would be easy to allow the qsas user to specify their own time label
format, since it is just a string. They might like this.
*/
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/*** Prototypes *********/
time_t mktimeutc(struct tm *tm);
static long JulianDay( long y, long m, long d );
double mktimeEpochSecs( int year, int month, int day,
int hour, int minute, int second );
void timefmts(PLFLT tmin, PLFLT tmax, char *labfmt, char *datefmt, int *is_t);
/* Steve Schwartz additions to replace C mktime by mktimeutc. This avoids
* problems due to mktime expecting its input in local time, which can also
* depend on daylight saving etc. The original recipe for determining the
* offset from time_t epoch was found to be not always effective. */
time_t mktimeutc(struct tm *tm){
/* convert tm struct with entries in UTC to the calendar time secs from a
* fixed system epoch. Note that different implementations of C can use different
* epochs, so can't hardwire the value. Instead we calculate it. This is tedious
* but robust and portable. Relies on algorithms borrowed from NASA NSSDC cdf
* project to assemble a time since 0AD.
*
* The returned time_t will be such that gmtime(returned value) gives tm
*
* It needs to be robust against local time zones, daylight saving time, etc.
* and is necessary because there is a known issue with C-time handling. See
* http://www.cl.cam.ac.uk/~mgk25/time/c/ for example
*
* History: Steve Schwartz [EMAIL PROTECTED] 31 July 2008
*
* Bugs/features:
* tm.tm_isdst isn't checked; this is deliberate since tm is assumed UTC
* tm members must be within valid ranges; the C mktime will adjust its
* members, but this one requires:
* (year+1900) 0->9999
* (mon+1) 1->12 note tm struct uses 0-11
* mday 1-31 no check against month
* hour 0->23
* min 0-59
* sec 0-59 no leap seconds
*
* all other members (wday, yday, isdst) ignored. If you need them, you
* can call this routine to find the time_t, then pass that to gmtime to return
* tm.
*
*/
time_t t;
struct tm tm0;
double tsecs, tsecs0;
t=0; /* start somewhere */
tm0 = *gmtime(&t); /* and find its UTC */
/* convert both tm structs to seconds since 0AD */
tsecs = mktimeEpochSecs(tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday, tm->tm_hour, tm->tm_min, tm->tm_sec);
tsecs0 = mktimeEpochSecs(tm0.tm_year+1900, tm0.tm_mon+1, tm0.tm_mday, tm0.tm_hour, tm0.tm_min, tm0.tm_sec);
return (time_t)(tsecs - tsecs0);
}
/* *NSSDC cdf epoch algorithms - in modified form - follow */
/* ****************************************************************************
* JulianDay.
* The year, month, and day are assumed to have already been validated. This
* is the day since 0 AD (julian day may not be the proper term).
******************************************************************************/
static long JulianDay( long y, long m, long d )
{
long jd;
jd = (long)
(367*y-7*(y+(m+9)/12)/4-3*((y+(m-9)/7)/100+1)/4+275*m/9+d+1721029);
return (jd);
}
/******************************************************************************
* computeEPOCH.
* Computes (and returns) an EPOCH in **SECS** since 0AD
* based on its component parts. -1.0
* is returned if an illegal component part is detected.
double computeEPOCH (year, month, day, hour, minute, second, msec)
long year, month, day, hour, minute, second; cdf msec removed
******************************************************************************/
double mktimeEpochSecs( int year, int month, int day,
int hour, int minute, int second )
{
double epoch;
/* printf("QEpochSecs got %d %d %d %d %d %d %d\n", year, month, day, hour,
minute, second, msec);
*/
if ((day < 1 || day > 31) ||
(month < 1 || month > 12) ||
(year < 0 || year > 9999) ||
(hour < 0 || hour > 23) ||
(minute < 0 || minute > 59) ||
(second < 0 || second > 59) )
epoch = -1.0;
else {
epoch = (double) JulianDay (year, month, day) - 1721060;
epoch = epoch * 24.0 + hour;
epoch = epoch * 60.0 + minute;
epoch = epoch * 60.0 + second;
}
return epoch;
}
void timefmts(PLFLT tmin, PLFLT tmax, char *labfmt, char *datefmt, int *is_t){
/* set up strings for time series plots.
* Input:
* tmin and tmax time limits in seconds
* don't need to be compatible with gmtime() to UTC
* Outputs:
*
* labfmt format for axis tick labels (in strftime format)
* datefmt format for start and end to give full date
* is_t if true (!=0), use plplot time handling, otherwise
* treat as x with time in seconds
* Note: fmts are at most 16 chars long + terminating null. Calling routine must
* pass pointer to char at least 17 long.
*/
PLFLT dt;
*is_t = 1; /* default */
dt = abs(tmax-tmin); /* can plot backward time */
if (dt < 3.0){ /* treat as x variable in secs */
*is_t = 0;
strncpy(datefmt, "%Y-%m-%d %H:%M:%S",17); datefmt[17]='\0';
return;
} else if (dt < 180.){ /* 3 minutes */
strncpy(datefmt, "%Y-%m-%d %H:",12); datefmt[12]='\0';
strncpy(labfmt, ":%M:%S",6); labfmt[6]='\0';
return;
} else if (dt < 10800.){ /* 3 hours */
strncpy(datefmt, "%Y-%m-%d",8); datefmt[8]='\0';
strncpy(labfmt, "%H:%M",5); labfmt[5]='\0';
return;
} else if (dt < 259200.){ /* 3 days */
strncpy(datefmt, "%Y-%m-",6); datefmt[6]='\0';
strncpy(labfmt, "%d %H:",6); labfmt[6]='\0';
return;
} else if (dt <7776000.){ /* 3 months */
strncpy(datefmt, "%Y",2); datefmt[2]='\0';
strncpy(labfmt, "%m-%d",5); labfmt[5]='\0';
return;
} else if (dt <94608000.){ /* 3 years */
strncpy(datefmt, "%Y-%m-%d",8); datefmt[8]='\0';
strncpy(labfmt, "%Y-%j",5); labfmt[5]='\0';
return;
} else { /* > 3 years */
strncpy(datefmt, "%Y-%m-%d",8); datefmt[8]='\0';
strncpy(labfmt, "%Y-%j",5); labfmt[5]='\0';
return;
}
/* shouldn't get here */
}
--- End Message ---
