/*
 * ACE/vis - Visualization of Flow and Transport
 *
 * Paul J. Turner and Antonio M. Baptista
 *
 * Copyright 1990-2003 Oregon Health and Science University
 * All Rights Reserved
 *
 */

/*
 *
 * misc utilities
 *
 * Contents:
 *
 * void cxfree() - free and check for NULL pointer
 * void fswap()  - swap doubles
 * void iswap()  - swap ints
 * void lowtoupper() - convert a string to upper case
 * void convertchar() - remove commas and Fortran D format
 * int ilog2() - integer log base 2, for the fft routine
 * double comp_area() - compute the area of a polygon
 * double comp_perimeter() - compute the perimeter
 *
 */

#ifndef lint
static char RCSid[] = "$Id: utils.c,v 1.5 2007/03/12 21:46:28 pturner Exp $";
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef WIN32
#include <pwd.h>
#endif
#include <time.h>
#include <sys/types.h>
#include <sys/stat.h>

#include "defines.h"
#include "globals.h"

static char readbuf[80];

int ibounds(int x, int lower, int upper, char *name)
{
    int test;

    test = ((x >= lower) && (x <= upper));
    if (!test) {
	sprintf(readbuf, " in %s : parameter must be in (%d , %d)", name, lower, upper);
	errwin(readbuf);
    }
    return (test);
}

int fbounds(double x, double lower, double upper, char *name)
{
    int test;

    test = ((x >= lower) && (x <= upper));
    if (!test) {
	sprintf(readbuf, "In %s : parameter must be in [%lf, %lf]", name, lower, upper);
	errwin(readbuf);
    }
    return (test);
}

void minmax2(double *x, int n, double missing, double *xmin, double *xmax, int *imin, int *imax)
{
    int i;

    *xmin = x[0];
    *xmax = x[0];
    *imin = 1;
    *imax = 1;
    for (i = 1; i < n; i++) {
	if (fabs(x[i]) < missing) {
	    if (x[i] < *xmin) {
		*xmin = x[i];
		*imin = i + 1;
	    }
	    if (x[i] > *xmax) {
		*xmax = x[i];
		*imax = i + 1;
	    }
	}
    }
}

void iminmax(int n, int *x, int *xmin, int *xmax)
{
    int i;

    *xmin = x[0];
    *xmax = x[0];
    for (i = 1; i < n; i++) {
	if (x[i] < *xmin) {
	    *xmin = x[i];
	}
	if (x[i] > *xmax) {
	    *xmax = x[i];
	}
    }
}

void minmax(int n, float *x, double *xmin, double *xmax)
{
    int i;

    *xmin = x[0];
    *xmax = x[0];
    for (i = 1; i < n; i++) {
	if (x[i] < *xmin) {
	    *xmin = x[i];
	}
	if (x[i] > *xmax) {
	    *xmax = x[i];
	}
    }
}

void dminmax(int n, double *x, double *xmin, double *xmax)
{
    int i;

    *xmin = x[0];
    *xmax = x[0];
    for (i = 1; i < n; i++) {
	if (x[i] < *xmin) {
	    *xmin = x[i];
	}
	if (x[i] > *xmax) {
	    *xmax = x[i];
	}
    }
}

void minmaxm(int n, float *x, double miss, double *xmin, double *xmax)
{
    int i;

    *xmin = 1.0e9;
    *xmax = -1.0e9;
    for (i = 0; i < n; i++) {
	if (x[i] != miss) {
	    if (x[i] < *xmin) {
		*xmin = x[i];
	    }
	    if (x[i] > *xmax) {
		*xmax = x[i];
	    }
	}
    }
}

/*
int tolower(c)
int c;
{
   if (c >= 'A' && c <= 'Z') {
        return(c - ' ');
    }
    else {
        return c;
    }
}
*/

/*
 * free and check for NULL pointer
 */
void cxfree(void *ptr)
{
    if (ptr != NULL) {
	free(ptr);
    }
}

/*
 * swap doubles and ints
 */
void fswap(double *x, double *y)
{
    double tmp;

    tmp = *x;
    *x = *y;
    *y = tmp;
}

void iswap(int *x, int *y)
{
    int tmp;

    tmp = *x;
    *x = *y;
    *y = tmp;
}

void cswap(char *x, char *y)
{
    char tmp;

    tmp = *x;
    *x = *y;
    *y = tmp;
}

/*
 * convert a string from lower to upper case
void lowtoupper(char *s)
{
    int i;

    for (i = 0; i < strlen(s); i++) {
	if (s[i] >= 'a' && s[i] <= 'z')
	    s[i] -= ' ';
    }
}
 */

/*
 * convert a string from lower to upper case
 * leaving quoted strings alone
 */
void lowtoupper(char *s)
{
    int i, quoteon = 0;

    for (i = 0; i < strlen(s); i++) {
	if (s[i] == '"') {
	    if (!quoteon) {
		quoteon = 1;
	    } else {
		quoteon = 0;
	    }
	}
	if (s[i] >= 'a' && s[i] <= 'z' && !quoteon) {
	    s[i] -= ' ';
	}
    }
}

/*
 * remove all that fortran nastiness
 */
void convertchar(char *s)
{
    while (*s++) {
	if (*s == ',')
	    *s = ' ';
	if (*s == 'D' || *s == 'd')
	    *s = 'e';
    }
}

/*
 * log base 2
 */
int ilog2(int n)
{
    int i = 0;
    int n1 = n;

    while (n1 >>= 1)
	i++;
    if (1 << i != n)
	return -1;
    else
	return i;
}

/*
 * compute the area bounded by the polygon (xi,yi)
 */
double comp_area(int n, double *x, double *y)
{
    int i;
    double sum = 0.0;

    for (i = 0; i < n; i++) {
	sum = sum + x[i] * y[(i + 1) % n] - y[i] * x[(i + 1) % n];
    }
    return sum * 0.5;
}

/*
 * compute the perimeter bounded by the polygon (xi,yi)
 */
double comp_perimeter(int n, double *x, double *y)
{
    int i;
    double sum = 0.0;

    for (i = 0; i < n - 1; i++) {
	sum = sum + hypot(x[i] - x[(i + 1) % n], y[i] - y[(i + 1) % n]);
    }
    return sum;
}

/*
 * interpolate between 2 time steps
 */
int interpolate_vector(int n, double t1, double t2, double t3, double *v1, double *v2, double *v3)
{
    int i;
    double dt, dv, d, v;
    dt = t2 - t1;
    if (dt == 0.0) {
	return 1;
    }
    d = (t3 - t1) / dt;
    for (i = 0; i < n; i++) {
	dv = v2[i] - v1[i];
	v3[i] = v1[i] + dv * d;
    }
    return 0;
}

int fexists(char *to)
{
    struct stat stto;
    char tbuf[256];

    if (stat(to, &stto) == 0) {
	sprintf(tbuf, "Overwrite %s?", to);
	if (!yesno(tbuf, "", "YES", "NO")) {
	    return (1);
	}
	return (0);
    }
    return (0);
}

long fsize(char *to)
{
    struct stat stto;
    char tbuf[256];

    if (stat(to, &stto) == 0) {
	return (long) stto.st_size;
    }
    return 0L;
}

/*
int isdigit(c)
    int c;
{
    return (c >= '0' && c <= '9');
}
*/

int argmatch(char *s1, char *s2, int atleast)
{
    int l1 = strlen(s1);
    int l2 = strlen(s2);

    if (l1 < atleast) {
	return 0;
    }
    if (l1 > l2) {
	return 0;
    }
    return (strncmp(s1, s2, l1) == 0);
}

/*
 * Time and date routines
 */

char *dayofweekstrs[] = { "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat" };
char *dayofweekstrl[] = { "Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday" };
char *months[] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" };
char *monthl[] = { "January", "February", "March", "April", "May", "June",
    "July", "August", "September", "October", "November", "December"
};

static int days1[] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334, 365 };
static int days2[] = { 0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335, 366 };

/*
 * return the Julian day + hms as a real number
 */

/*
** Takes a date, and returns a Julian day. A Julian day is the number of
** days since some base date  (in the very distant past).
** Handy for getting date of x number of days after a given Julian date
** (use jdate to get that from the Gregorian date).
** Author: Robert G. Tantzen, translator: Nat Howard
** Translated from the algol original in Collected Algorithms of CACM
** (This and jdate are algorithm 199).
*/
double julday(int mon, int day, int year, int h, int mi, double se)
{
    long m = mon, d = day, y = year;
    long c, ya, j;
    double seconds = h * 3600.0 + mi * 60 + se;

    if (m > 2)
	m -= 3;
    else {
	m += 9;
	--y;
    }
    c = y / 100L;
    ya = y - (100L * c);
    j = (146097L * c) / 4L + (1461L * ya) / 4L + (153L * m + 2L) / 5L + d + 1721119L;
    if (seconds < 12 * 3600.0) {
	j--;
	seconds += 12.0 * 3600.0;
    } else {
	seconds = seconds - 12.0 * 3600.0;
    }
    return (j + (seconds / 3600.0) / 24.0);
}

/* Julian date converter. Takes a julian date (the number of days since
** some distant epoch or other), and returns an int pointer to static space.
** ip[0] = month;
** ip[1] = day of month;
** ip[2] = year (actual year, like 1977, not 77 unless it was  77 a.d.);
** ip[3] = day of week (0->Sunday to 6->Saturday)
** These are Gregorian.
** Copied from Algorithm 199 in Collected algorithms of the CACM
** Author: Robert G. Tantzen, Translator: Nat Howard
*/
void calcdate(double jd, int *m, int *d, int *y, int *h, int *mi, double *sec)
{
    static int ret[4];

    long j = jd;
    double tmp, frac = jd - j;

    if (frac >= 0.5) {
	frac = frac - 0.5;
	j++;
    } else {
	frac = frac + 0.5;
    }

    ret[3] = (j + 1L) % 7L;
    j -= 1721119L;
    *y = (4L * j - 1L) / 146097L;
    j = 4L * j - 1L - 146097L * *y;
    *d = j / 4L;
    j = (4L * *d + 3L) / 1461L;
    *d = 4L * *d + 3L - 1461L * j;
    *d = (*d + 4L) / 4L;
    *m = (5L * *d - 3L) / 153L;
    *d = 5L * *d - 3 - 153L * *m;
    *d = (*d + 5L) / 5L;
    *y = 100L * *y + j;
    if (*m < 10)
	*m += 3;
    else {
	*m -= 9;
	*y += 1;
    }
    tmp = 3600.0 * (frac * 24.0);
    *h = (int) (tmp / 3600.0);
    tmp = tmp - *h * 3600.0;
    *mi = (int) (tmp / 60.0);
    *sec = tmp - *mi * 60.0;
}

int dayofweek(double j)
{
    j += 0.5;
    return (int) (j + 1) % 7;
}

int leapyear(int year)
{
    if (year % 4 == 0) {
	return (1);
    } else {
	return (0);
    }
}

/*
   get the month and day given the number of days
   from the beginning of the year 'yr'
*/
void getmoday(int days, int yr, int *mo, int *da)
{
    int i;

    if (leapyear(yr)) {
	for (i = 0; i < 13; i++) {
	    if (days <= days2[i]) {
		*mo = i;
		*da = (days - days2[i - 1]);
		goto out1;
	    }
	}
    } else {
	for (i = 0; i < 13; i++) {
	    if (days <= days1[i]) {
		*mo = i;
		*da = (days - days1[i - 1]);
		goto out1;
	    }
	}
    }
  out1:;
}

/*
   return the number of days from the beginning of the year 'yr'
*/
int getndays(double j)
{
    int m, d, y, hh, mm;
    double ss;

    calcdate(j, &m, &d, &y, &hh, &mm, &ss);
    if (leapyear(y)) {
	return days2[m - 1] + d;
    } else {
	return days1[m - 1] + d;
    }
}

/*
   return hms
*/
void gethms(double j, int *h, int *m, int *s)
{
    double rem = j - (long) j;

    *h = (int) (rem * 24);
    *m = (int) ((rem * 24 - *h) * 60);
    *s = (int) (((rem * 24 - *h) * 60 - *m) * 60);
    *h = (*h + 12) % 24;
}

void set_timestamp(void)
{
    struct tm tm;
    long time_value;
    char *str;

    (void) time(&time_value);
    //time_value -=  8 * 3600;
    tm = *gmtime(&time_value);
    str = asctime(&tm);
    strcpy(timestamp.s, str);
    timestamp.s[strlen(timestamp.s) - 1] = 0;
}

void revbytes(unsigned char *ptr, int len)
{
    int i, iadj = len - 1;
    unsigned char ctmp;
    for (i = 0; i < len / 2; i++) {
	ctmp = ptr[i];
	ptr[i] = ptr[iadj - i];
	ptr[iadj - i] = ctmp;
    }
}

int read_double(double *d, int n, FILE * fp)
{
    int err, type, i;
    err = fread(d, sizeof(double), n, fp);
    if (err == n && swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(double));
	}
    }
    return err;
}

int write_double(double *d, int n, FILE * fp)
{
    int err, type, i;
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(double));
	}
    }
    err = fwrite(d, sizeof(double), n, fp);
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(double));
	}
    }
    return err == n ? 0 : err;
}

int read_int(int *d, int n, FILE * fp)
{
    int err, type, i;
    err = fread(d, sizeof(int), n, fp);
    if (err == n && swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(int));
	}
    }
    return err;
}

int write_int(int *d, int n, FILE * fp)
{
    int err, type, i;
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(int));
	}
    }
    err = fwrite(d, sizeof(int), n, fp);
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(int));
	}
    }
    return err == n ? 0 : err;
}

int read_char(char *d, int n, FILE * fp)
{
    int err, type;
    err = fread(d, sizeof(char), n, fp);
    return err;
}

int write_char(char *d, int n, FILE * fp)
{
    int err, type;
    err = fwrite(d, sizeof(char), n, fp);
    return err == n ? 0 : err;
}

int read_short(short *d, int n, FILE * fp)
{
    int err, type, i;
    err = fread(d, sizeof(short), n, fp);
    if (err == n && swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(short));
	}
    }
    return err;
}

int write_short(short *d, int n, FILE * fp)
{
    int err, type, i;
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(short));
	}
    }
    err = fwrite(d, sizeof(short), n, fp);
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(short));
	}
    }
    return err == n ? 0 : err;
}

int read_float(float *d, int n, FILE * fp)
{
    int err, type, i;
    err = fread(d, sizeof(float), n, fp);
    if (err == n && swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(float));
	}
    }
    return err;
}

int write_float(float *d, int n, FILE * fp)
{
    int err, type, i;
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(float));
	}
    }
    err = fwrite(d, sizeof(float), n, fp);
    if (swapBytes) {
	for (i = 0; i < n; i++) {
	    revbytes((unsigned char *) &d[i], sizeof(float));
	}
    }
    return err == n ? 0 : err;
}