/* $Id: pldtik.c 8082 2007-12-12 15:44:17Z andrewross $
Modified 30 July 2008 by Steve Schwartz (s.schwartz@imperial.ac.uk) to
replace C mktime by utc alternative.
Determines tick spacing and mode (fixed or floating) of
numeric axis labels.
Copyright (C) 2004 Alan W. Irwin
This file is part of PLplot.
PLplot is free software; you can redistribute it and/or modify
it under the terms of the GNU General Library 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
*/
#include "plplotP.h"
/* Patch submitted by Steve Schwartz s.schwartz@imperial.ac.uk Aug 2008 includes
* following routines to replace C's mktime. The routine mktimeutc is used in pldtfac
* in place of mktime, and then tdiff isn't needed because it is zero by construction.
* Prototypes placed in plplotP.h, but mktimeutc might be useful for the general
* user, to construct times. */
/*
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 pldtik()
*
* Determine tick spacing: works out a "nice" interval (if tick == 0) such
* that there are between 3 and 7.5 major tick intervals in the input
* range vmin to vmax. The recommended number of subticks is returned in
* "nsubt" unless the routine is entered with a non-zero value of "nsubt".
* n.b. big change: now returns only positive values of tick and nsubt
\*----------------------------------------------------------------------*/
void
pldtik(PLFLT vmin, PLFLT vmax, PLFLT *tick, PLINT *nsubt, PLBOOL ld)
{
PLFLT t1, t2, tick_reasonable;
PLINT np, ns;
PLFLT diff, factor;
if (ld) {
/* Check suitable units for tick spacing */
pldtfac(vmin, vmax, &factor, NULL);
*tick = *tick / factor;
vmin = vmin / factor;
vmax = vmax / factor;
}
/* Magnitude of min/max difference to get tick spacing */
t1 = (PLFLT) log10(ABS(vmax - vmin));
np = (PLINT) floor(t1);
t1 = t1 - np;
/* Get tick spacing. */
if (t1 > 0.7781512503) {
t2 = 2.0;
ns = 4;
}
else if (t1 > 0.4771212549) {
t2 = 1.0;
ns = 5;
}
else if (t1 > 0.1760912591) {
t2 = 5.0;
ns = 5;
np = np - 1;
}
else {
t2 = 2.0;
ns = 4;
np = np - 1;
}
/* Now compute reasonable tick spacing */
tick_reasonable = t2 * pow(10.0, (double) np);
if (*tick == 0) {
*tick = t2 * pow(10.0, (double) np);
}
else {
*tick = ABS(*tick);
if(*tick < 1.e-4*tick_reasonable) {
plexit("pldtik: magnitude of specified tick spacing is much too small");
return;
}
}
if (*nsubt == 0)
*nsubt = ns;
*nsubt = ABS(*nsubt);
if (ld) {
*tick = *tick*factor;
}
}
/*----------------------------------------------------------------------*\
* PLFLT pldtfac()
*
* Calculate factor to convert a date/time interval in seconds
* into a more natural units (minutes, hours, days, week, years).
* Also optionally calculate the sensible start time for counting ticks
* from (e.g. beginning of day, beginning of year).
* Used to calculate sensible tick and label spacings.
\*----------------------------------------------------------------------*/
void
pldtfac(PLFLT vmin, PLFLT vmax, PLFLT *factor, PLFLT *start) {
PLFLT diff, tdiff;
time_t t, t2;
struct tm tm;
diff = vmax - vmin;
if (start != NULL) {
t = (time_t) vmin;
tm = *gmtime(&t);
/* with mktimeutc routine, mktimeutc(&tm) will always = t, so no need to
calculate tdiff. */
/* Here was the old way: ******************
t2 = mktime(&tm); */
/* Arg! This is because mktime is in local time and we need to
correct for the offset. C time handling really is broken... */
/* tdiff = difftime(t,t2); *******************/
}
if (diff < 3.0*60.0) {
/* Seconds */
*factor = 1.0;
if (start != NULL) {
tm.tm_sec = 0;
/* Again the old way *****************
t = mktime(&tm);
*start = (PLFLT)t+tdiff; sim edits done in rest of else ifs*/
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
else if (diff < 3.0*60.0*60.0) {
/* Minutes */
*factor = 60.0;
if (start != NULL) {
tm.tm_sec = 0;
tm.tm_min = 0;
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
else if (diff < 3.0*60.0*60.0*24.0) {
/* Hours */
*factor = 60.0*60.0;
if (start != NULL) {
tm.tm_sec = 0;
tm.tm_min = 0;
tm.tm_hour = 0;
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
else if (diff < 3.0*60.0*60.0*24.0*7.0) {
/* Days */
*factor = 60.0*60.0*24.0;
if (start != NULL) {
tm.tm_sec = 0;
tm.tm_min = 0;
tm.tm_hour = 0;
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
else if (diff < 3.0*60.0*60.0*24.0*365) {
/* Weeks */
*factor = 60.0*60.0*24.0*7.0;
if (start != NULL) {
tm.tm_sec = 0;
tm.tm_min = 0;
tm.tm_hour = 0;
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
else {
/* Years */
*factor = 60.0*60.0*24.0*365.25;
if (start != NULL) {
tm.tm_sec = 0;
tm.tm_min = 0;
tm.tm_hour = 0;
tm.tm_mday = 1;
tm.tm_mon = 0;
t = mktimeutc(&tm);
*start = (PLFLT)t;
}
}
}
/*----------------------------------------------------------------------*\
* void pldprec()
*
* Determine precision: the output variable "mode" is set to 0 if labels
* are to be written in floating-point format, or to 1 if they are to be
* written in scientific format. For mode = 1, the exponent will be
* placed at:
*
* top left for vertical axis on left
* top right for vertical axis on right
* bottom right for horizontal axis
*
* The digmax flag can be set by the user, and represents the maximum
* number of digits a label may occupy including sign and decimal point.
* digmin, calculated internally, is the maximum number of digits
* labels at vmin and vmax would occupy if floating point.
* If digmax<0, it is disregarded,
* and if digmax=0 the default value is used. For digmax>0, mode=1 is
* chosen if there is insufficient room for the label within the specified
* # of digits (digmin > digfix, where digfix is determined from digmax with
* fuzz factors).
*
* In the case of mode=0, the actual # of digits will become too large
* when the magnitude of the labels become too large. The mode=1 case
* offers the greatest precision for the smallest field length.
*
* The determination of maximum length for fixed point quantities is
* complicated by the fact that very long fixed point representations look
* much worse than the same sized floating point representation. Further,
* a fixed point number with a large negative exponent will actually gain
* in precision when written as floating point. Thus we use certain fuzz
* factors to get 'digfix' from 'digmax', however it will always be true
* that digfix<=digmax.
*
* Finally, if 'digmax' is set, 'prec' is reduced in size if necessary so
* that the labels fit the requested field length, where prec is the number of
* places after the decimal place.
\*----------------------------------------------------------------------*/
#define MIN_FLTDIG 3 /* disregarded if fractional part is 0 */
#define MAX_FIXDIG_POS 6
#define MAX_FIXDIG_NEG 4
#define DIGMAX_DEF 5
void
pldprec(PLFLT vmin, PLFLT vmax, PLFLT tick, PLINT lf,
PLINT *mode, PLINT *prec, PLINT digmax, PLINT *scale)
{
PLFLT chosen, notchosen, vmod, t0;
PLINT msd, notmsd, np, digmin, digfix;
*mode = 0;
*scale = 0;
if (digmax == 0)
digmax = DIGMAX_DEF;
/* Choose vmin or vmax depending on magnitudes of vmin and vmax. */
chosen = (ABS(vmax) >= ABS(vmin))? vmax: vmin;
notchosen = (ABS(vmax) >= ABS(vmin))? vmin: vmax;
/* Magnitute of chosen to get number of significant digits */
if(ABS(chosen) > 0.) {
vmod = ABS(chosen);
t0 = (PLFLT) log10(vmod);
msd = (PLINT) floor(t0);
}
else {
/* this branch occurs only when 0. --- 0. range put in */
vmod = 1.;
t0 = (PLFLT) log10(vmod);
msd = (PLINT) floor(t0);
}
if(ABS(notchosen) > 0.)
notmsd = (PLINT) floor( (PLFLT) log10(ABS(notchosen)));
else
notmsd = msd;
/* Autoselect the mode flag */
/* 'digmin' is the minimum number of places taken up by the label */
if (msd >= 0) {
/* n.b. no decimal point in the minimal case */
digmin = msd + 1;
digfix = MAX_FIXDIG_POS;
if (digmax > 0)
digfix = MIN(digmax, MAX_FIXDIG_POS);
}
else {
/* adjust digmin to account for leading 0 and decimal point */
digmin = -msd + 2;
digfix = MAX_FIXDIG_NEG;
if (digmax > 0)
digfix = MIN(digmax, MAX_FIXDIG_NEG);
}
/* adjust digmin to account for sign on the chosen end of axis or sign on the
* nonchosen end of axis if notmsd = msd or (msd <= 0 and notmsd < 0)
* For the latter case the notchosen label starts with "-0."
* For checking for the latter case, the notmsd < 0 condition is redundant
* since notmsd <= msd always and the equal part is selected by the first
* condition.
*/
if(chosen < 0.||(notchosen < 0. && (notmsd == msd || msd <= 0)))
digmin = digmin + 1;
if (digmin > digfix && !lf) {
*mode = 1;
*scale = msd;
}
/* Establish precision. */
/* It must be fine enough to resolve the tick spacing */
np = (PLINT) floor(log10(ABS(tick)));
if (*mode != 0)
*prec = msd - np;
else
*prec = MAX(-np, 0);
/* One last hack required: if exponent < 0, i.e. number has leading '0.',
* it's better to change to floating point form if the number of digits
* is insufficient to represent the tick spacing.
*/
if (*mode == 0 && digmax > 0 && !lf) {
if (t0 < 0.0) {
if (digmax - 2 - *prec < 0) {
*mode = 1;
*scale = msd;
}
}
else
*prec = MAX(MIN(*prec, digmax - msd - 1), 0);
}
if (*mode != 0) {
*prec = msd - np;
*prec = MAX(MIN(*prec, MAX(digmax-1, MIN_FLTDIG)), 0);
}
}
/* 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 s.schwartz@imperial.ac.uk 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 at 0 calendar time */
tm0 = *gmtime(&t); /* and find its UTC */
/* convert both tm structs to seconds since 0AD; could error check here, as
* mktimeEpochSecs returns -1 if inputs are out of range, but no sensible way
* for calling routine to test error in returned time_t value ?*/
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);
}
/******************************************************************************
* mktimeEpochSecs.
* 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 mktimeEpochSecs( int year, int month, int day,
int hour, int minute, int second )
{
double epoch;
long lyear, lmonth, lday;
/* Julian day needs to be long; inputs typically from struct tm, and are ints
* so convert for safe portability */
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 {
lyear = (long) year;
lmonth = (long) month;
lday = (long) day;
epoch = (double) JulianDay (year, month, day) - 1721060;
epoch = epoch * 24.0 + hour;
epoch = epoch * 60.0 + minute;
epoch = epoch * 60.0 + second;
}
return epoch;
}