//Utils.cpp, Copyright (c) 2000-2008 R.Lackner
//Collection of utility functions and classes for RLPlot
//
//    This file is part of RLPlot.
//
//    RLPlot is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    RLPlot 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 General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with RLPlot; if not, write to the Free Software
//    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <errno.h>
#include "rlplot.h"

extern GraphObj *CurrGO;			//Selected Graphic Objects
extern Label *CurrLabel;
extern Default defs;
extern UndoObj Undo;

char TmpTxt[TMP_TXT_SIZE];

//----------------------------------------------------------------------------
// Get the rectanpular part of bitmap. Used for screen updates
//----------------------------------------------------------------------------
anyOutput *GetRectBitmap(RECT *rc, anyOutput *src)
{
	RECT cr;
	anyOutput *bm;

	if(!rc || !src) return 0L;
	src->ActualSize(&cr);
	if(rc->left < cr.left) rc->left = cr.left;
	if(rc->right > cr.right) rc->right = cr.right;
	if(rc->top < cr.top) rc->top = cr.top;
	if(rc->bottom > cr.bottom) rc->bottom = cr.bottom;
	if(rc->left == rc->right) return 0L;
	if(rc->top == rc->bottom) return 0L;
	if(!(bm = NewBitmapClass(rc->right - rc->left, rc->bottom - rc->top, 
		src->hres, src->vres)))return 0L;
	bm->CopyBitmap(0, 0, src, rc->left, rc->top, rc->right - rc->left, 
		rc->bottom - rc->top, false);
	return bm;
}

void RestoreRectBitmap(anyOutput **pmo, RECT *mrc, anyOutput *o)
{
	if(pmo && *pmo && o && mrc) {
		o->CopyBitmap(mrc->left, mrc->top, *pmo, 0, 0, mrc->right - mrc->left, 
			mrc->bottom - mrc->top, false);
		DelBitmapClass(*pmo);		*pmo = 0L;
		o->UpdateRect(mrc, false);
		}
}

//----------------------------------------------------------------------------
// Axis utility functions
//----------------------------------------------------------------------------
void NiceAxis(AxisDEF *axis, int nTick)
{
	double diff, logStep, Step, Magn, LoVal, HiVal;
	int i;

	axis->Start = axis->min;
	axis->Step = (axis->max - axis->min)/((double)nTick);
	diff = axis->max - axis->min;
	if(axis->breaks) for(i = 0; i < axis->nBreaks; i++) {
		diff -= fabs(axis->breaks[i].fy - axis->breaks[i].fx);
		}
	if(diff <= 0.0) return;
	logStep = log10(diff/(double)nTick);
	Magn = floor(logStep);			logStep -= Magn;
	if(logStep > 0.8) Step = 10.0;
	else if(logStep > 0.5) Step = 5.0;
	else if(logStep > 0.2) Step = 2.0; 
	else Step = 1.0;
	Step *= pow(10.0, Magn);		HiVal = LoVal = Step * floor(axis->min/Step);
	axis->max += (diff * 0.05);
	while(HiVal < axis->max) HiVal += Step;
	if((axis->flags & AXIS_LOG) == AXIS_LOG) {
		if (LoVal > defs.min4log) axis->min = LoVal;
		if ((LoVal + Step) > defs.min4log && (LoVal + Step) < axis->min) axis->min = LoVal+Step;
		}
	else axis->min = LoVal;
	axis->max = HiVal;		axis->Start = axis->min;		axis->Step = Step;
}

void NiceStep(AxisDEF *axis, int nTick)
{
	double diff, d, logStep, Step, Magn;
	int i;

	diff = axis->max - axis->min;		d = axis->Step != 0.0 ? diff/axis->Step : HUGE_VAL;
	if((d - floor(d)) < 0.1 && axis->Step != 0.0 && diff/axis->Step < 12.0)return;
	if(axis->breaks) for(i = 0; i < axis->nBreaks; i++) {
		diff -= fabs(axis->breaks[i].fy - axis->breaks[i].fx);
		}
	if(diff <= 0.0) return;
	logStep = log10(diff/(double)nTick);
	Magn = floor(logStep);				logStep -= Magn;
	if(logStep > 0.8) Step = 10.0;
	else if(logStep > 0.5) Step = 5.0;
	else if(logStep > 0.2) Step = 2.0; 
	else Step = 1.0;
	Step *= pow(10.0, Magn);			axis->Step = Step;
}

double base4log(AxisDEF *axis, int direc)
{
	double lv, Step = 1.0, Magn;
	int cmd;

	switch (direc) {
	case 0:		cmd = SIZE_BOUNDS_XMIN;		break;
	case 1:		cmd = SIZE_BOUNDS_YMIN;		break;
	case 2:		cmd = SIZE_BOUNDS_ZMIN;		break;
	default:	return 1.0;
		}
	lv = axis->min > defs.min4log ? axis->min : defs.min4log;
	if(lv <= defs.min4log) return defs.min4log;
	lv = log10(lv);
	lv -= (Magn = floor(lv));
	if(lv > 0.301) Step = 2.0; 
	if(lv > 0.699) Step = 5.0;
	Step *= pow(10.0, Magn);
	return Step > defs.min4log ? Step : 1.0;
}

double TransformValue(AxisDEF *axis, double val, bool transform)
{
	int i;
	double f1, f2, RetVal = val;

	if(!axis) return val;
	if(axis->breaks) {
		for (i = 0; i < axis->nBreaks; i++) {
			f1 = axis->breaks[i].fx;	f2 = axis->breaks[i].fy;
			if(val > f2) RetVal -= (f2-f1);
			else if(val > f1 && val <= f2) RetVal -= (val-f1);
			}
		}
	else (axis->nBreaks = 0);
	if(transform) {
		switch(axis->flags & 0x7000L) {
		case AXIS_LINEAR:	break;
		case AXIS_LOG:
			if(axis->flags & AXIS_RADIAL) RetVal = fabs(RetVal);
			RetVal = RetVal > defs.min4log ? log10(RetVal): log10(defs.min4log);
			break;
		case AXIS_RECI:		RetVal = RetVal > defs.min4log || RetVal < - defs.min4log ?
				1.0/RetVal : 0.0;		break;
		case AXIS_SQR:		RetVal = RetVal >= 0.0 ? sqrt(RetVal) : 0.0;					break;
			}
		}
	return RetVal;
}

void SortAxisBreaks(AxisDEF *axis)
{
	int i, j;
	double ftmp;
	bool sorted;

	if(!axis || !axis->nBreaks || !axis->breaks) return;
	//low values first
	for(i = 0; i < axis->nBreaks; i++) {
		if(axis->breaks[i].fy < axis->breaks[i].fx) {
			ftmp = axis->breaks[i].fx;
			axis->breaks[i].fx = axis->breaks[i].fy;
			axis->breaks[i].fy = ftmp;
			}
		}
	//a simple bubble sort should do
	if(axis->nBreaks >1) do {
		sorted = true;
		for(i = 1; i < axis->nBreaks; i++) {
			if(axis->breaks[i-1].fx > axis->breaks[i].fx) {
				ftmp = axis->breaks[i-1].fx;
				axis->breaks[i-1].fx = axis->breaks[i].fx;
				axis->breaks[i].fx = ftmp;
				ftmp = axis->breaks[i-1].fy;
				axis->breaks[i-1].fy = axis->breaks[i].fy;
				axis->breaks[i].fy = ftmp;
				sorted = false;
				}
			}
		}while(!sorted);
	//combine overlapping ranges
	if((j = axis->nBreaks) >1) for(i = j = 1; i < axis->nBreaks; i++) {
		if(axis->breaks[i].fx > axis->breaks[j-1].fy) {
			axis->breaks[j].fx = axis->breaks[i].fx;
			axis->breaks[j].fy = axis->breaks[i].fy;
			j++;
			}
		else {
			j--;
			axis->breaks[j++].fy = axis->breaks[i].fy;
			}
		}
	axis->nBreaks = j;
}

double GetAxisFac(AxisDEF *axis, double delta, int direc)
{
	double da, v1, v2;

	switch(axis->flags & 0x7000L) {
	case AXIS_LOG:		
		axis->max = axis->max > defs.min4log ? log10(axis->max): log10(defs.min4log);
		axis->min = axis->min > defs.min4log ? log10(axis->min): 
			log10(base4log(axis, direc));
		if(axis->max <= axis->min) axis->max = axis->min +1.0;
		break;
	case AXIS_RECI:
		v1 = fabs(axis->min) >defs.min4log ? axis->min :
			base4log(axis, direc);
		if(fabs(v1) > defs.min4log) v1 = 1.0/v1;
		else v1 = 1.0e+34;
		if(fabs(axis->max) >defs.min4log) v2 = 1.0/axis->max;
		else v2 = 0.0;
		if(fabs(v2) < fabs(v1/10.0)) v2 = 0.0;
		axis->min = v2;		axis->max = v1;
		break;
	case AXIS_SQR:
		axis->max = axis->max > defs.min4log ? sqrt(axis->max) : 0.0;
		axis->min = axis->min > defs.min4log ? sqrt(axis->min) : 0.0;
		break;
		}
	v2 = TransformValue(axis, axis->max, false);	//process breaks
	v1 = TransformValue(axis, axis->min, false);
	da = v2 != v1 ? v2 - v1 : 1.0;
	return delta / da;
}

//----------------------------------------------------------------------------
// Text utility functions: internationalization and formats
//----------------------------------------------------------------------------
//remove leading/trailing whitespace
char *str_ltrim(char *str) {
	int i, j;

	if(!str || !str[0]) return str;
	for(i = 0; str[i] && str[i] <= ' '; i++);
	for(j = 0; str[i]; str[j++] = str[i++]);
	str[j++] = '\0';	return str;
	}

char *str_rtrim(char *str) {
	size_t i;

	i = strlen(str);
	while(i > 0 && str[i-1] <= ' ') str[--i] = '\0';
	return str;
}

char *str_trim(char *str) {
	str = str_ltrim(str);			return str_rtrim(str);
}

//remove leading and tailing quotatation
void rmquot(char *str)
{
	size_t i, len;
	char c;

	if(str && str[0] && (*str == '"' || *str == '\'')) {
		len = strlen(str);		c = *str;
		if(str[len-1] == c) {
			str[len-1] = 0;
			for(i = 1; i < len; str[i-1] = str[i++]);
			}
		}
}

int strpos(char *needle, char *haystack)
{
	int i, j;

	if(!needle || !needle[0] || !haystack || !haystack[0]) return -1;
	for(i = j = 0; haystack[i]; i++, j=0) {
		if(haystack[i] == needle[0]) for (j = 1; haystack[i+j]; j++) {
			if(needle[j] != haystack[i+j]) break; 
			}
		if(j && !needle[j]) return i;
		}
	return -1;
}

char *strreplace(char *needle, char *replace, char *haystack)
{
	static char *result = 0L;
	static size_t reslen = 0;
	size_t i, j, k, l;

	if(!needle || !needle[0] || !haystack || !haystack[0]) return result;
	if(!result) result = (char*)malloc(reslen = 100);
	result[0] = 0;		l = strlen(needle);
	for(i = j = k = 0; haystack[i]; i++, j=0) {
		if(haystack[i] == needle[0]) for (j = 1; haystack[i+j]; j++) {
			if(needle[j] != haystack[i+j]) break; 
			}
		if(j && !needle[j]) {
			if(replace && replace[0]) {
				if(reslen < (i + (int)strlen(replace) + 10)) {
					result = (char*)realloc(result, reslen += 100);
					}
				for(j = 0; replace[j]; j++) result[k++] = replace[j];
				}
			i += (l-1);
			}
		else result[k++] = haystack[i];
		}
	result[k++] = 0;
	return result;
}

char *substr(char *text, int pos1, int pos2)
{
	static char *result = 0L;
	static size_t reslen = 0;
	int i, j;
	size_t l;

	if(!text || !text[0]) return 0L;
	l = strlen(text);
	if(pos1 < 0) pos1 = 0;			if(pos2 < pos1) pos2 = (int)(l+1);
	if(!result) result = (char*)malloc(reslen = 100);
	while (reslen < l) result = (char*) realloc(result, reslen += 100);
	for(i = 0; i < pos1 && text[i]; i++);
	for(j = 0; i <= pos2 && text[i]; result[j++] = text[i++]);
	result[j] = 0;
	return result;
}

//copy string in src to dest, returning the stringlength of src
//   seek better solution for long strings
int rlp_strcpy(char*dest, int size, char*src)
{
	int i;

	if(dest && src) {
		for(i = 0; i < size; i++) {
			if(!(dest[i] = src[i])) return i;
			}
		dest[i-1] = 0;
		return i-1;
		}
	return 0;
}

// restyle formula
void ReshapeFormula(char **text)
{
	int i, j, l;

	if(!text || !*text || !**text) return;
	l = (int)strlen(*text);
	for(i = j = 0; i < l; i++) {
		if((*text)[i] == ';') {
			if((*text)[i+1] == ';' || (*text)[i+1] == '\n') i++;
			}
		TmpTxt[j++] = (*text)[i];
		}
	TmpTxt[j] = 0;
	if(j && j <= l && TmpTxt[0]) {
		rlp_strcpy(*text, l+1, TmpTxt);
		}
}

//translate anyResult to output format
void TranslateResult(anyResult *res)
{
	static char tr_text[80];

	switch (res->type) {
	case ET_VALUE:
		if(res->value == HUGE_VAL) rlp_strcpy(tr_text, 80, "inf");
		else if(res->value == -HUGE_VAL) rlp_strcpy(tr_text, 80, "-inf");
#ifdef USE_WIN_SECURE
		else sprintf_s(tr_text, 80, "%g", res->value);
#else
		else sprintf(tr_text, "%g", res->value);
#endif
		res->text = tr_text;				return;
	case ET_BOOL:
		rlp_strcpy(tr_text, 80, ((int)res->value) ? (char*)"true" : (char*)"false"); 
		res->text = tr_text;				return;
	case ET_DATE:
		rlp_strcpy(tr_text, 80, value_date(res->value, defs.fmt_date)); 
		res->text = tr_text;				return;
	case ET_TIME:
		rlp_strcpy(tr_text, 80, value_date(res->value, defs.fmt_time)); 
		res->text = tr_text;				return;
	case ET_DATETIME:
		rlp_strcpy(tr_text, 80, value_date(res->value, defs.fmt_datetime)); 
		res->text = tr_text;				return;
	case ET_TEXT:
		if(res->text && res->text[0])		return;
		}
	if(!(res->text)) res->text="";
}

//remove invalid tag combinations from string
void CleanTags(char *txt, int *i1, int *i2, int *i3)
{
	char *no_tags[] = {"<b></b>", "</b><b>", "<b><b>", "</b></b>",
		"<i></i>", "</i><i>", "<i><i>", "</i></i>", "<u></u>",
		"</u><u>", "<u><u>", "</u></u>", "<sub></sub>",
		"</sub><sub>", "<sup></sup>", "</sup><sup>", 
		0L};
	int i, j, k, l, w;
	bool na;

	for(i = j = 0; txt[i]; i++) {
		if(txt[i] != '<') txt[j++] = txt[i];
		else {
			for(k = 0, na = true; no_tags[k]; k++) {
				for(l=1; no_tags[k][l] && txt[i+l]; l++) 
					if(no_tags[k][l] != txt[i+l]) break;
				if(!no_tags[k][l]){
					na = false;
					i += ((w=(int)strlen(no_tags[k]))-1);
					if(i1 && *i1 > i) *i1 -= w;
					if(i2 && *i2 > i) *i2 -= w;
					if(i3 && *i3 > i) *i3 -= w;
					break;
					}
				}
			if(na) txt[j++] = txt[i];
			}
		}
	txt[j++] = 0;
}

//replace one character in string
void ChangeChar(char *text, char c1, char c2)	//replace one char in string
{
	int i;

	for(i = 0; text[i]; i++) if (text[i] == c1) text[i] = c2;
}

char *Int2Nat(char *text)	//format ASCII number to locale format
{
	int i;

	for(i = 0; text[i]; i++) if (text[i] == '.') text[i] = defs.DecPoint[0];
	return text;
}

char *Nat2Int(char *text)	//format locale number to intranational notation
{
	int i;

	for(i = 0; text[i]; i++) if (text[i] == defs.DecPoint[0]) text[i] = '.';
	return text;
}

void WriteNatFloatToBuff(char *buff, double val)
{
#ifdef USE_WIN_SECURE
	sprintf_s(buff, 20, " %g", val);
#else
	sprintf(buff, " %g", val);
#endif
	Int2Nat(buff);
}

bool Txt2Flt(char *txt, double *val)
{
	char *tmp = 0L;

	if(txt && txt[0] && val) {
		if(!txt[1] && (txt[0] == defs.DecPoint[0] || txt[0] < '0' ||
			txt[0] > '9'))return false;
		if(txt && txt[0] && (tmp = (char*)memdup(txt, (int)strlen(txt)+1, 0))){
			Nat2Int(tmp);
			//the return value of sscanf only roughly identifies a number
#ifdef USE_WIN_SECURE
			if(!sscanf_s(tmp,"%lf", val)){
#else
			if(!sscanf(tmp,"%lf", val)){
#endif
				free(tmp);
				return false;
				}
			free(tmp);
			return true;
			}
		}
	return false;
}

void RmTrail(char *txt)
{
	int i;

	i = (int)strlen(txt);
	while(i >0 && (txt[i-1] == '0' || txt[i-1] < 32)) txt[--i] = 0;
	if(i > 1 && txt[i-1] == '.') txt[i-1] = 0;
}

double NiceValue(double fv)
{
	double sign = fv > 0.0f ? 1.0 : -1.0;
	double fa = fabs(fv);
	double magn = floor(log10(fa));
	int i = iround(fa/pow(10.0, magn-1.0));

	return sign*pow(10.0, magn-1.0) *(double)i;
}

char *NiceTime(double val)
{
	rlp_datetime dt;

	parse_datevalue(&dt, val);
	if(dt.year > 1905) {
		if(dt.hours) return date2text(&dt, defs.fmt_datetime);
		else return date2text(&dt, defs.fmt_date);
		}
	else return date2text(&dt, defs.fmt_time);
}

char *Int2ColLabel(int nr1, bool uc)
{
	static char RetTxt[12];
	int i, j, nr = nr1;
	char base = uc ? 'A' : 'a';

#ifdef USE_WIN_SECURE
	sprintf_s(RetTxt+8, 4, "%c\0", base + (nr %26));
#else
	sprintf(RetTxt+8, "%c\0", base + (nr %26));
#endif
	nr /= 26;
	for (i = 7; nr && i>=0; i--) {
		j = nr %27;
		RetTxt[i] = base + (j ? j-1 : j);
		if (nr == 26) nr = 0;
		else nr = nr/26;
		}
	return RetTxt+i+1;
}

char *mkCellRef(int row, int col)
{
	static char RetTxt[20];

#ifdef USE_WIN_SECURE
	sprintf_s(RetTxt, 20, "%s%d", Int2ColLabel(col, false), row+1);
#else
	sprintf(RetTxt, "%s%d", Int2ColLabel(col, false), row+1);
#endif
	return RetTxt;
}

char *mkRangeRef(int r1, int c1, int r2, int c2)
{
	static char RetTxt[40];
	int cb;

	cb = rlp_strcpy(RetTxt, 20, mkCellRef(r1, c1)); 
	RetTxt[cb++] = ':';
	rlp_strcpy(RetTxt+cb, 40-cb, mkCellRef(r2, c2));
	return RetTxt;
}
//convert strings to XML specifications
//this offers a limited support for special characters
char *str2xml(char *str, bool bGreek)
{
	int i, j;
	wchar_t wc;

	if(!str) return 0L;
	for(i = j = 0; str[i] && j < 4090; i++) {
		switch(str[i]) {
		case '"':
			j += rlp_strcpy(TmpTxt+j, TMP_TXT_SIZE-j, "&quot;");
			break;
		case '&':
			j += rlp_strcpy(TmpTxt+j, TMP_TXT_SIZE-j, "&amp;");
			break;
		case '<':
			j += rlp_strcpy(TmpTxt+j, TMP_TXT_SIZE-j, "&lt;");
			break;
		case '>':
			j += rlp_strcpy(TmpTxt+j, TMP_TXT_SIZE-j, "&gt;");
			break;
		default:
#ifdef USE_WIN_SECURE
			if(bGreek && str[i] >= 'A' && str[i] <= 'Z') j += sprintf_s(TmpTxt+j, 20, "&#%d;", str[i] - 'A' + 0x391);
			else if(bGreek && str[i] >= 'a' && str[i] <= 'z') j += sprintf_s(TmpTxt+j, 20, "&#%d;", str[i] - 'a' + 0x3B1);
#else
			if(bGreek && str[i] >= 'A' && str[i] <= 'Z') j += sprintf(TmpTxt+j, "&#%d;", str[i] - 'A' + 0x391);
			else if(bGreek && str[i] >= 'a' && str[i] <= 'z') j += sprintf(TmpTxt+j, "&#%d;", str[i] - 'a' + 0x3B1);
#endif
			else if((unsigned char)str[i] <= 127) TmpTxt[j++]=str[i];
			else {
				if(mbtowc(&wc, str+i, 1) >0) 
#ifdef USE_WIN_SECURE
					j += sprintf_s(TmpTxt+j, TMP_TXT_SIZE-j, "&#%d;", ((unsigned short)wc));
#else
					j += sprintf(TmpTxt+j, "&#%d;", ((unsigned short)wc));
#endif
				}
			}
		}
	TmpTxt[j] = 0;
	return(TmpTxt);
}

// split string str into array of strings using sep as separator
// return number of lines created in nl
static char  *split_buf = 0L;
static int split_buf_size, split_buf_pos;
char **split(char *str, char sep, int *nl)
{
	int i, j, l, ns;
	char **ptr;

	if(!str || !str[0] || !sep) return 0L;
	split_buf_pos = 0;
	add_to_buff(&split_buf, &split_buf_pos, &split_buf_size, str, 0);
	if(!split_buf || !split_buf_pos) return 0L;
	for(i = ns = 0; i < split_buf_pos; i++) if(split_buf[i] == sep) ns++;
	if(!(ptr = (char**)calloc(ns+2, sizeof(char*)))) return 0L;
	for(i = j = l = 0; i < split_buf_pos; i++) {
		if(split_buf[i] == sep) {
			split_buf[i] = 0;
			ptr[l++] = (char*)memdup(split_buf+j, (int)strlen(split_buf+j)+1, 0);
			j = i+1;
			}
		}
	ptr[l++] = (char*)memdup(split_buf+j, (int)strlen(split_buf+j)+1, 0);
	if(nl) *nl = l;
	return ptr;
}

char *fit_num_rect(anyOutput *o, int max_width, char *num_str)
{
	int i, j, k, w, h, slen;
	char mant[30], expo[30], fmt[20];
	double num;

	o->oGetTextExtent(num_str, slen = (int)strlen(num_str), &w, &h);
	if(w < (max_width-5)) return num_str;
	//first try to remove leading zero from exponent
	for(i = 0; i < slen; i++) if(num_str[i] == 'e') break;
	if(num_str[i] == 'e') while (num_str[i+2] == '0') {
		for(j = i+2; num_str[j]; j++) num_str[j] = num_str[j+1];
		o->oGetTextExtent(num_str, slen = (int)strlen(num_str), &w, &h);
		if(w < (max_width-5)) return num_str;
		}
	//no success: reduce the number of digit by rounding
	for(i = k = 0; i <= slen; i++){
		if(num_str[i] == '.') k = i;
		if((mant[i] = num_str[i]) == 'e' || mant[i] == 0) break;
		}
	if(num_str[i] =='e') mant[i] = 0;						k = i - k-1;
	for(j = 0; num_str[i]; j++, i++) expo[j] = num_str[i];		expo[j] = 0;
#ifdef USE_WIN_SECURE
	sscanf_s(mant, "%lf", &num);
#else
	sscanf(mant, "%lf", &num);	
#endif
	if(k >0) do {
#ifdef USE_WIN_SECURE
		sprintf_s(fmt, 20, "%%.%dlf%s", k, expo);
		slen = sprintf_s(num_str, 60, fmt, num);		//num_str is much longer than 60
#else
		sprintf(fmt, "%%.%dlf%s", k, expo);
		slen = sprintf(num_str, fmt, num);
#endif
		k--;
		o->oGetTextExtent(num_str, slen, &w, &h);
		if(w < (max_width-5)) return num_str;
		}while (k >= 0);
	//cannot fit: return hash marks instead
	for(i = w = 0; w < (max_width-5) && i < 11; i++) {
		rlp_strcpy(num_str, i+2, "##########");
		o->oGetTextExtent(num_str, i+1, &w, &h);
		}
	num_str[i-1] = 0;		return num_str;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// utilities to add a line or number to a text buffer: memory file
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void add_to_buff(char** dest, int *pos, int *csize, char *txt, int len)
{
	if(!len)len = (int)strlen(txt);
	if(!len) return;
	if((*pos+len+1)>= *csize) {
		*csize += 1000;
		while(*csize < (*pos+len+1)) *csize += 1000;
		*dest = (char*)realloc(*dest, *csize);
		}
	if(*dest) {
		*pos += rlp_strcpy(*dest+*pos, len+1, txt);
		}
}

void add_int_to_buff(char** dest, int *pos, int *csize, int value, bool lsp, int ndig)
{
	int len;
	char txt[40];

#ifdef USE_WIN_SECURE
	len = sprintf_s(txt, 40, lsp ? " %d" : "%d", value);
#else
	len = sprintf(txt, lsp ? " %d" : "%d", value);
#endif
	add_to_buff(dest, pos, csize, txt, len);
}

void add_dbl_to_buff(char** dest, int *pos, int *csize, double value, bool lsp)
{
	int len;
	char txt[40];

#ifdef USE_WIN_SECURE
	len = sprintf_s(txt, 40, lsp ? " %g" : "%g", value);
#else
	len = sprintf(txt, lsp ? " %g" : "%g", value);
#endif
	add_to_buff(dest, pos, csize, txt, len);
}

void add_hex_to_buff(char** dest, int *pos, int *csize, DWORD value, bool lsp)
{
	int len;
	char txt[40];

	if(value) {
#ifdef USE_WIN_SECURE
		len = sprintf_s(txt, 40, lsp ? " 0x%08x" : "0x%08x", value);
#else
		len = sprintf(txt, lsp ? " 0x%08x" : "0x%08x", value);
#endif
		add_to_buff(dest, pos, csize, txt, len);
		}
	else if(lsp) add_to_buff(dest, pos, csize, " 0x0", 4);
	else add_to_buff(dest, pos, csize, "0x0", 3);
}

//----------------------------------------------------------------------------
// bounding rectangle utilities
//----------------------------------------------------------------------------
void SetMinMaxRect(RECT *rc, int x1, int y1, int x2, int y2)
{
	if(x1 > x2) {
		rc->left = x2;		rc->right = x1;
		}
	else {
		rc->left = x1;		rc->right = x2;
		}
	if(y1 > y2) {
		rc->top = y2;		rc->bottom = y1;
		}
	else {
		rc->top = y1;		rc->bottom = y2;
		}
}

void UpdateMinMaxRect(RECT *rc, int x, int y)
{
	if(x < rc->left) rc->left = x;
	if(x > rc->right) rc->right = x;
	if(y < rc->top) rc->top = y;
	if(y > rc->bottom) rc->bottom = y;
}

void IncrementMinMaxRect(RECT *rc, int i)
{
	rc->left -= i;			rc->right += i;
	rc->top -= i;			rc->bottom += i;
}

bool IsInRect(RECT *rc, int x, int y)
{
	if(x < rc->left || x > rc->right || 
		y < rc->top || y > rc->bottom) return false;
	return true;
}

//----------------------------------------------------------------------------
// test if point (e.g. mouse position) is close to or iside a shape
//----------------------------------------------------------------------------
bool IsCloseToLine(POINT *p1, POINT *p2, int x, int y)
{
	bool bFound = false;
	int tmp, dx, dy;

	//do not process single point
	if(p1->x == p2->x && p1->y == p2->y) return false;
	//point must be bracketed by p1:p2
	if((x-2) > p1->x && (x-2) > p2->x) return false;
	if((x+2) < p1->x && (x+2) < p2->x) return false;
	if((y-2) > p1->y && (y-2) > p2->y) return false;
	if((y+2) < p1->y && (y+2) < p2->y) return false;
	if(abs(dx = (p2->x - p1->x)) < abs(dy = (p2->y - p1->y))) {	//y dominant
		tmp = (p1->x + ((y - p1->y) * dx)/dy);
		if(x > (tmp-3) && x < (tmp+3)) bFound = true;
		}
	else {							// x dominant
		tmp = (p1->y + ((x - p1->x) * dy)/dx);
		if(y > (tmp-3) && y < (tmp+3)) bFound = true;
		}
	return bFound;
}

bool IsCloseToPL(POINT p, POINT *pts, int cp)
{
	int i;

	for( i = 1; i < cp; i++) if(IsCloseToLine(pts+i-1, pts+i, p.x, p.y)) return true;
	return false;
}

//test if poitn p is within the polygon pts of size cp
//note: the last point of the polypon should be equal to the first point,
//   i.e. the polygon must be closed
int idx_point_on_line;
bool IsInPolygon(POINT *p, POINT *pts, int cp)
{
	int tmp1 = 0, tmp2 = 0, tmp3, i;

	for(i = 1; i < cp; i++) {
		if((pts[i-1].x <= p->x && pts[i].x > p->x) || (pts[i-1].x > p->x && pts[i].x <= p->x)) {
			tmp3 = ((p->x - pts[i-1].x)*(pts[i].y -pts[i-1].y))/(pts[i].x - pts[i-1].x);
			tmp3 += pts[i-1].y;
			if(p->y > tmp3) tmp1++;
			else if(p->y < tmp3) tmp2++;
			else return false;			//ignore points on the line
			}
		}
	return((tmp1 & 0x1) && (tmp2 & 0x1));
}

bool IsInPolygon3D(double x, double y, POINT3D *pts, int cp, int *us, int *ls)
{
	int tmp1 = 0, tmp2 = 0, i, idx1, idx2;
	double tmp3;

	for(i = 1, idx1 = idx2 = -1; i < cp; i++) {
		if((pts[i-1].x <= x && pts[i].x > x) || (pts[i-1].x > x && pts[i].x <= x)) {
			tmp3 = ((x - pts[i-1].x)*(pts[i].y -pts[i-1].y))/(pts[i].x - pts[i-1].x);
			tmp3 += pts[i-1].y;
			if(y > tmp3){
				tmp1++;				idx1 = i;
				}
			else if(y < tmp3) {
				tmp2++;				idx2 = i;
				}
			else {			//points is on the line
				idx1 = idx2 = idx_point_on_line = i;
				if(us) *us = idx1;		if(ls) *ls = idx2;
				return true;
				}
			}
		}
	//return an index to the bracketing segments of the polygon
	if(us) *us = idx1;			if(ls) *ls = idx2;
	return((tmp1 & 0x1) && (tmp2 & 0x1));
}

//return true if two rectangles overlap
bool OverlapRect(RECT *rc1, RECT *rc2)
{
	if((rc1->left < rc2->right && rc1->right > rc2->left) ||
		(rc2->left < rc1->right && rc2->right > rc1->left)) {
		if((rc1->top < rc2->bottom && rc1->bottom > rc2->top) ||
			(rc2->top < rc1->bottom && rc2->bottom > rc1->top))
			return true;
		}
	return false;
}

//----------------------------------------------------------------------------
// collect points to a polygon
// keep number of points low by extrapolation
//----------------------------------------------------------------------------
void AddToPolygon(long *cp, POINT *pts, POINT *np)
{
	int ix, iy;
	long i = *cp;

	if(i && pts[i-1].x == np->x && pts[i-1].y == np->y) return;
	if(i < 2) {					//accept first points of polygon
		pts[i].x = np->x;	pts[i].y = np->y;	*cp = i+1;
		return;
		}
	if(pts[i-1].x == pts[i-2].x && pts[i-1].x == np->x) {
		if(np->y == pts[i-1].y) return;
		if((np->y > pts[i-1].y && pts[i-1].y > pts[i-2].y) ||
			(np->y < pts[i-1].y && pts[i-1].y < pts[i-2].y)) {
			pts[i-1].x = np->x;		pts[i-1].y = np->y;
			return;
			}
		}
	if(pts[i-1].y == pts[i-2].y && pts[i-1].y == np->y) {
		if(np->x == pts[i-1].x) return;
		if((np->x > pts[i-1].x && pts[i-1].x > pts[i-2].x) ||
			(np->x < pts[i-1].x && pts[i-1].x < pts[i-2].x)) {
			pts[i-1].x = np->x;		pts[i-1].y = np->y;
			return;
			}
		}
	//try linear extrapolation
	if((pts[i-1].x > pts[i-2].x && np->x > pts[i-1].x) ||
		(pts[i-1].x < pts[i-2].x && np->x < pts[i-1].x)) {
		ix = (pts[i-1].y != pts[i-2].y) ? (pts[i-2].x + ((np->y - pts[i-2].y) * 
			(pts[i-1].x - pts[i-2].x))/(pts[i-1].y - pts[i-2].y)) : 0;
		iy = (pts[i-1].x != pts[i-2].x) ? (pts[i-2].y + ((np->x - pts[i-2].x) * 
			(pts[i-1].y - pts[i-2].y))/(pts[i-1].x - pts[i-2].x)) : 0;
		if((ix && ix == np->x) && (iy && iy == np->y)) {
			pts[i-1].x = np->x;		pts[i-1].y = np->y;
			return;
			}
		}
	//not explained by extrapolation, accept new point
	pts[i].x = np->x;				pts[i].y = np->y;
	*cp = i+1;
	return;
}

//----------------------------------------------------------------------------
// create a Bezier polygon
#define MIN_SEG 11
#define MAX_DEPTH 5
void DrawBezier(long *cp, POINT *pts, POINT p0, POINT p1, POINT p2, POINT p3, int depth)
{
	int i;
	POINT np, p01, p12, p23, p012, p123, p0123;
	POINT *ap[] = {&p0, &p1, &p2, &p3};

	depth ++;
	if(depth > MAX_DEPTH) {
		for(i= 0; i < 4; i++) {
			np.x = (*ap[i]).x >> 2;		np.y = (*ap[i]).y >> 2;
			AddToPolygon(cp, pts, &np);
			}
		return;
		}
	else if(depth == 1) for(i=0; i < 4; i++) {
		(*ap[i]).x <<= 2;			(*ap[i]).y <<= 2;
		}
	p01.x = (p0.x + p1.x) >> 1;			p01.y = (p0.y + p1.y) >> 1;
	p12.x = (p1.x + p2.x) >> 1;			p12.y = (p1.y + p2.y) >> 1;
	p23.x = (p2.x + p3.x) >> 1;			p23.y = (p2.y + p3.y) >> 1;
	p012.x = (p01.x + p12.x) >> 1;			p012.y = (p01.y + p12.y) >> 1;
	p123.x = (p12.x + p23.x) >> 1;			p123.y = (p12.y + p23.y) >> 1;
	p0123.x = (p012.x + p123.x) >> 1;		p0123.y = (p012.y + p123.y) >> 1;
	if(abs(p0.x - p0123.x)> MIN_SEG || abs(p0.y - p0123.y)> MIN_SEG) {
		DrawBezier(cp, pts, p0, p01, p012, p0123, depth);	//recursion: refine
		}
	else {
		DrawBezier(cp, pts, p0, p01, p012, p0123, MAX_DEPTH);	//recursion: store data
		}
	if(abs(p3.x - p0123.x)> MIN_SEG || abs(p3.y - p0123.y)> MIN_SEG) {
		DrawBezier(cp, pts, p0123, p123, p23, p3, depth);	//recursion: refine
		}
	else {
		DrawBezier(cp, pts, p0123, p123, p23, p3, MAX_DEPTH);	//recursion: store data
		}
}
#undef MAX_DEPTH
#undef MIN_SEG

// create a Bezier polygon clipped to rectangle
static RECT BezClipRec;
#define MIN_SEGCLP 5
#define MAX_DEPTHCLP 6
static void DrawBezierClip(long *cp, POINT *pts, POINT p0, POINT p1, POINT p2, POINT p3, int depth)
{
	int i;
	POINT np, p01, p12, p23, p012, p123, p0123;
	POINT *ap[] = {&p0, &p1, &p2, &p3};

	depth++;
	if(depth > MAX_DEPTHCLP) {
		for(i= 0; i < 4; i++) {
			np.x = (*ap[i]).x >> 2;		np.y = (*ap[i]).y >> 2;
			if(np.x < BezClipRec.left) np.x = BezClipRec.left;
			if(np.x > BezClipRec.right) np.x = BezClipRec.right;
			if(np.y < BezClipRec.top) np.y = BezClipRec.top;
			if(np.y > BezClipRec.bottom) np.y = BezClipRec.bottom;
			AddToPolygon(cp, pts, &np);
			}
		return;
		}
	else if(depth == 1) for(i=0; i < 4; i++) {
		(*ap[i]).x <<= 2;			(*ap[i]).y <<= 2;
		}
	p01.x = (p0.x + p1.x) >> 1;			p01.y = (p0.y + p1.y) >> 1;
	p12.x = (p1.x + p2.x) >> 1;			p12.y = (p1.y + p2.y) >> 1;
	p23.x = (p2.x + p3.x) >> 1;			p23.y = (p2.y + p3.y) >> 1;
	p012.x = (p01.x + p12.x) >> 1;			p012.y = (p01.y + p12.y) >> 1;
	p123.x = (p12.x + p23.x) >> 1;			p123.y = (p12.y + p23.y) >> 1;
	p0123.x = (p012.x + p123.x) >> 1;		p0123.y = (p012.y + p123.y) >> 1;
	if((abs(p0.x - p0123.x)> MIN_SEGCLP || abs(p0.y - p0123.y)> MIN_SEGCLP)) {
		DrawBezierClip(cp, pts, p0, p01, p012, p0123, depth);		//recursion: refine
		}
	else {
		DrawBezierClip(cp, pts, p0, p01, p012, p0123, MAX_DEPTHCLP);	//store data
		}
	if((abs(p3.x - p0123.x)> MIN_SEGCLP || abs(p3.y - p0123.y)> MIN_SEGCLP)) {
		DrawBezierClip(cp, pts, p0123, p123, p23, p3, depth);		//recursion: refine
		}
	else {
		DrawBezierClip(cp, pts, p0123, p123, p23, p3, MAX_DEPTHCLP);	//store data
		}
}

void ClipBezier(long *cp, POINT *pts, POINT p0, POINT p1, POINT p2, POINT p3, POINT *clp1, POINT *clp2)
{
	if(clp1 && clp2){
		SetMinMaxRect(&BezClipRec, clp1->x, clp1->y, clp2->x, clp2->y);
		}
	if(cp && pts) DrawBezierClip(cp, pts, p0, p1, p2, p3, 0);
}
#undef MIN_SEGCLP
#undef MAX_DEPTHCLP

//----------------------------------------------------------------------------
// create a Bezier spline through data points
static void ipol_curve(lfPOINT *p1, lfPOINT *p2, lfPOINT *p3, lfPOINT *cp1, lfPOINT *cp2)
{
	double dx, dy, l;
	lfPOINT tHat;

	tHat.fx = p3->fx - p1->fx;		tHat.fy = p3->fy - p1->fy;
	l = sqrt(tHat.fx * tHat.fx + tHat.fy * tHat.fy);
	tHat.fx /= l;					tHat.fy /= l;
	cp1->fx = cp2->fx = p2->fx;		cp1->fy = cp2->fy = p2->fy;
	dx = p3->fx - p2->fx;			dy = p3->fy - p2->fy;		l = sqrt(dx*dx + dy*dy)/3.0;
	cp2->fx += (tHat.fx * l);		cp2->fy += (tHat.fy *l);
	dx = p2->fx - p1->fx;			dy = p2->fy - p1->fy;		l = sqrt(dx*dx + dy*dy)/3.0;
	cp1->fx -= (tHat.fx * l);		cp1->fy -= (tHat.fy *l);
}

static void mirr_vecvec(lfPOINT *a0, lfPOINT *a1, lfPOINT *v1)
{
	double dx, dy, as, ac, vs, vc, s, c, l;
	lfPOINT sol1, sol2;

	//mirror vector a0 -> v1 by rotation around vector a0 -> a1
	dx = a1->fx - a0->fx;					dy = a1->fy - a0->fy;
	l = sqrt(dx*dx + dy*dy);	as = dy/l;			ac = dx/l;
	dx = v1->fx - a0->fx;					dy = v1->fy - a0->fy;
	l = sqrt(dx*dx + dy*dy);	vs = dy/l;			vc = dx/l;
	//calculate cw and ccw solution
	s = sin((asin(as)-asin(vs))*2.0);		c = cos((acos(ac)-acos(vc))*2.0);
	sol1.fx = dx * c + dy * s + a0->fx;		sol1.fy = dy * c - dx * s + a0->fy;
	s = sin((asin(as)-asin(-vs))*2.0);		c = cos((acos(ac)-acos(-vc))*2.0);
	sol2.fx = dx * c + dy * s + a0->fx;		sol2.fy = dy * c - dx * s + a0->fy;
	//get better solution
	dx = sol1.fx - a1->fx;					dy = sol1.fy - a1->fy;
	l = sqrt(dx*dx + dy*dy);
	dx = sol2.fx - a1->fx;					dy = sol2.fy - a1->fy;
	if( l < sqrt(dx*dx + dy*dy)) {
		v1->fx = sol1.fx;					v1->fy = sol1.fy;
		}
	else {
		v1->fx = sol2.fx;					v1->fy = sol2.fy;
		}
}

int mkCurve(lfPOINT *src, int n1, lfPOINT **dst, bool bClosed)
{
	int i, j, iret;

	if(!src || n1 < 3) return 0;
	if(!(*dst = (lfPOINT*)malloc((n1*3+2) * sizeof(lfPOINT))))return 0;
	for(i = j = iret = 0; i < n1; i++, j += 3) {
		(*dst)[j].fx = src[i].fx;		(*dst)[j].fy = src[i].fy;
		if(i && i < (n1-1)){
			ipol_curve(&src[i-1], &src[i], &src[i+1], *dst+j-1, *dst+j+1);
			}
		else if(i) {
			iret = j+1;					break;
			}
		}
	if(bClosed && iret >2) {
		if((*dst)[j].fx != (*dst)[0].fx || (*dst)[j].fy != (*dst)[0].fy) {
			j += 3;								iret += 3;
			(*dst)[j].fx = (*dst)[0].fx;		(*dst)[j].fy = (*dst)[0].fy;	
			}
		if(j < 6) {
			free(*dst);		*dst = 0L;				return 0;
			}
		ipol_curve(*dst+j-3, *dst+j, *dst+3, *dst+j-1, *dst+1);
		ipol_curve(*dst+j-6, *dst+j-3, *dst+j, *dst+j-4, *dst+j-2);
		}
	else if(!bClosed && iret>3) {
		(*dst)[1].fx = (*dst)[0].fx + (*dst)[3].fx - (*dst)[2].fx;
		(*dst)[1].fy = (*dst)[0].fy + (*dst)[3].fy - (*dst)[2].fy;
		mirr_vecvec(*dst, *dst+3, *dst+1);
		(*dst)[j-1].fx = (*dst)[j].fx + (*dst)[j-3].fx - (*dst)[j-2].fx;
		(*dst)[j-1].fy = (*dst)[j].fy + (*dst)[j-3].fy - (*dst)[j-2].fy;
		mirr_vecvec(*dst+j, *dst+j-3, *dst+j-1);
		}
	else {
		free(*dst);			*dst = 0L;				return 0;
		}
	return iret;
}

//----------------------------------------------------------------------------
// create a circular polygon
//use circular Bresenham's algorithm to draw arcs
//Ref: C. Montani, R. Scopigno (1990) "Speres-To-Voxel Conversion", in:
//   Graphic Gems (A.S. Glassner ed.) Academic Press, Inc.; 
//   ISBN 0-12-288165-5 
//----------------------------------------------------------------------------
POINT *MakeArc(int ix, int iy, int r, int qad, long *npts)
{
	int i, x, y, di, de, lim;
	static POINT *pts, *rpts;
	POINT np;

	if(r < 1) return 0L;
	if(!(pts = (POINT*) malloc((r+1)*8*sizeof(POINT))))return 0L;
	for(i = *npts = 0; i < 4; i++) {
		x = lim = 0;	y = r;	di = 2*(1-r);
		if(qad & (1<<i))while (y >= lim){
			if(di < 0) {
				de = 2*di + 2*y -1;
				if(de > 0) {
					x++;	y--;	di += (2*x -2*y +2);
					}
				else {
					x++;	di += (2*x +1);
					}
				}
			else {
				de = 2*di -2*x -1;
				if(de > 0) {
					y--;	di += (-2*y +1);
					}
				else {
					x++;	y--;	di += (2*x -2*y +2);
					}
				}
			switch(i) {
			case 0:	np.x = ix-y;	np.y = iy+x;	break;
			case 1: np.x = ix+x;	np.y = iy+y;	break;
			case 2: np.x = ix+y;	np.y = iy-x;	break;
			case 3:	np.x = ix-x;	np.y = iy-y;	break;
				}
			AddToPolygon(npts, pts, &np);
			}
		}
	if(*npts < 3) return 0L;
	if(rpts = (POINT*)realloc(pts, sizeof(POINT)*(*npts+4))) return rpts;
	return pts;
}

//----------------------------------------------------------------------------
// display a marked line using complement colors
//----------------------------------------------------------------------------
void InvertLine(POINT *pts, int nPts, LineDEF *Line, RECT *rc, 
	anyOutput *o, bool mark)
{
	int i;
	LineDEF OldLine;

	memcpy(&OldLine, Line, sizeof(LineDEF));
	if(OldLine.width <= 0.0001) OldLine.width = 0.0001;
	for(i = 0; o->un2fiy(OldLine.width) < 1.0 && i < 50; i++) OldLine.width *= 2.0;
	OldLine.color = mark ? (Line->color & 0x00ffffffL) : 0x00ffffffL;
	OldLine.width *= 3.0;				o->SetLine(&OldLine);
	o->oPolyline(pts, nPts);			OldLine.width = Line->width;
	OldLine.color = mark ? (Line->color & 0x00ffffffL) ^ 0x00ffffffL : (Line->color & 0x00ffffff);
	o->SetLine(&OldLine);				o->oPolyline(pts, nPts);
	if(rc) o->UpdateRect(rc, false);
}

//----------------------------------------------------------------------------
// color utilitis
//----------------------------------------------------------------------------
// calculate distance between two colors
unsigned int ColDiff(DWORD col1, DWORD col2)
{
	int ret = 0, d;

	d = (col1 & 0xff) - (col2 & 0xff);
	ret = isqr(d*d);
	d = ((col1>>8) & 0xff) - ((col2>>8) & 0xff);
	ret += isqr(d*d);
	d = ((col1>>16) & 0xff) - ((col2>>16) & 0xff);
	ret += isqr(d*d);
	return ret;	
}

//----------------------------------------------------------------------------
// interpolate between two colors
DWORD IpolCol(DWORD color1, DWORD color2, double fact)
{
	DWORD col1, col2, col3, c1;
	int i;

	col1 = color1;	col2 = color2;	col3 = 0x0L;
	for(i = 0; i < 4; i++) {
		c1 = iround(fabs((((col1 & 0xff000000)>>24) & 0xff) * fact));
		c1 += iround(fabs((((col2 & 0xff000000)>>24) & 0xff) *(1.0-fact)));
		col3 |= c1 < 0xff ? c1 : 0xff;
		col1 <<= 8;	col2 <<= 8; 
		if(i < 3) col3 <<= 8;
		}
	return col3;
}

//----------------------------------------------------------------------------
// Random number generator with low sequential correlations.
// ran2 returns a number betwee 0.0f and 1.0f;
// Ref: W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Veterling
// Numerical Recipe in C, The Art of Scientific Computing
// Cambridge University Press 1988, ISBN 0-521-35465-X
//----------------------------------------------------------------------------
#define M 714025
#define IA 1366
#define IC 150889
double ran2(long *idum)
{
	static long iy, ir[98];
	static int iff = 0;
	int j;

	if(*idum < 0 || iff == 0) {
		iff = 1;
		if((*idum = (IC-(*idum)) % M) < 0) *idum = -(*idum);
		for(j = 1; j <= 97; j++) {
			*idum = (IA*(*idum)+IC) % M;
			ir[j] = (*idum);
			}
		*idum=(IA*(*idum)+IC) % M;
		iy=(*idum);
		}
	j = 1+97 * iy/M;
	if(j > 97 || j< 1) return 0.0f;		//impossible
	iy = ir[j];
	*idum = (IA*(*idum)+IC) % M;
	ir[j] = (*idum);
	return (float) iy/M;
}
#undef IC
#undef IA
#undef M

//----------------------------------------------------------------------------
// integer square root
// calculate the largest number <= sqr(n)
// Christopher J. Musial in Graphics Gems II, page 37ff, page 610
// Academic Press, 1991
// modified
//----------------------------------------------------------------------------

unsigned long int isqr(unsigned long int n)
{
	unsigned long int nextTrial, decrement;

	if (nextTrial = n>>1) {
		for ( ; ; ){
			if(decrement = (nextTrial - n/nextTrial)>>1) nextTrial -= decrement;
			else if(nextTrial * nextTrial > n) return nextTrial-1;
			else return nextTrial;
			}
		}
	return n;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// multiply two rotation matrices
//  see: Graphic Gems, A.S. Glassner ed.; Academic Press Inc.
//  ISBN 0-12-286165-5, p.640
//
bool MatMul(double a[3][3], double b[3][3], double c[3][3])
{
	int i, j, k;
	bool success = true;

	for(i = 0; i < 3; i++) {
		for (j = 0; j < 3; j++) {
			c[i][j] = 0.0;
			for(k = 0; k < 3; k++) c[i][j] += a[i][k] * b[k][j];
			if(c[i][j] < -.99999 || c[i][j] > .99999) success = false;
			}
		}
	return success;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//Return a format string depending on the range
//
char *GetNumFormat(double Magn)
{
	if(Magn < -3.0) return("%0.2le");
	if(Magn == -3.0) return("%0.4lf");
	if(Magn == -2.0) return("%0.3lf");
	if(Magn == -1.0) return("%0.2lf");
	if(Magn == 0.0) return("%0.1lf");
	if(Magn >= 5.0) return("%0.2le");
	return("%0.0lf");
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//Delete a graphic object using the Id member of the class to select the 
//   the proper destructor
void DeleteGO(GraphObj *go)
{
	if(!go) return;
	if(!go->Id) {
		//this will also set the Id member of the class
		go->Command(CMD_SET_DATAOBJ, 0L, 0L);
		}
	if(go == CurrGO) CurrGO = 0L;
	switch(go->Id) {
	case GO_AXIS:			delete((Axis*)go);			break;
	case GO_TICK:			delete((Tick*)go);			break;
	case GO_GRIDLINE:		delete((GridLine*)go);		break;
	case GO_GRIDLINE3D:		delete((GridLine3D*)go);	break;
	case GO_SYMBOL:			delete((Symbol*)go);		break;
	case GO_BUBBLE:			delete((Bubble*)go);		break;
	case GO_BAR:			delete((Bar*)go);			break;
	case GO_ERRBAR:			delete((ErrorBar*)go);		break;
	case GO_ARROW:			delete((Arrow*)go);			break;
	case GO_BOX:			delete((Box*)go);			break;
	case GO_WHISKER:		delete((Whisker*)go);		break;
	case GO_DROPLINE:		delete((DropLine*)go);		break;
	case GO_DATALINE:
		//we call ~DataLine for ~DataPolygon because variables are
		//   initialized in DataLine ??!!??!!
		//   otherwise we would crash with ~DataPolygon.
	case GO_DATAPOLYGON:	delete((DataLine*)go);		break;
	case GO_SPHERE:			delete((Sphere*)go);		break;
	case GO_PLANE:			delete((plane*)go);			break;
	case GO_BRICK:			delete((Brick*)go);			break;
	case GO_LINE3D:			delete((Line3D*)go);		break;
	case GO_LABEL:			delete((Label*)go);			break;
	case GO_MLABEL:			delete((mLabel*)go);		break;
	case GO_SEGMENT:		delete((segment*)go);		break;
	case GO_POLYGON:
	case GO_POLYLINE:		delete((polyline*)go);		break;
	case GO_REGLINE:		delete((RegLine*)go);		break;
	case GO_SDELLIPSE:		delete((SDellipse*)go);		break;
	case GO_ELLIPSE:		
	case GO_ROUNDREC:
	case GO_RECTANGLE:		delete((rectangle*)go);		break;
	case GO_DRAGHANDLE:		delete((dragHandle*)go);	break;
	case GO_DRAGRECT:		delete((dragRect*)go);		break;
	case GO_DRAG3D:			delete((Drag3D*)go);		break;
	case GO_FRAMERECT:		delete((FrmRect*)go);		break;
	case GO_PLOT:			delete((Plot*)go);			break;
	case GO_BARCHART:
	case GO_PLOTSCATT:		delete((PlotScatt*)go);		break;
	case GO_REGRESSION:		delete((Regression*)go);	break;
	case GO_BUBBLEPLOT:		delete((BubblePlot*)go);	break;
	case GO_BOXPLOT:		delete((BoxPlot*)go);		break;
	case GO_DENSDISP:		delete((DensDisp*)go);		break;
	case GO_STACKBAR:
	case GO_WATERFALL:
	case GO_STACKPG:		delete((StackBar*)go);		break;
	case GO_POLARPLOT:		delete((PolarPlot*)go);		break;
	case GO_RINGCHART:
	case GO_PIECHART:		delete((PieChart*)go);		break;
	case GO_GROUP:
	case GO_STARCHART:		delete((GoGroup*)go);		break;
	case GO_SCATT3D:		delete((Scatt3D*)go);		break;
	case GO_PLOT3D:			delete((Plot3D*)go);		break;
	case GO_PAGE:
	case GO_GRAPH:			delete((Graph*)go);			break;
	case GO_SVGOPTIONS:		delete((svgOptions*)go);	break;
	case GO_DROPL3D:		delete((DropLine3D*)go);	break;
	case GO_ARROW3D:		delete((Arrow3D*)go);		break;
	case GO_LIMITS:			delete((Limits*)go);		break;
	case GO_GRIDRADIAL:		delete((GridRadial*)go);	break;
	case GO_DEFRW:			delete((DefsRW*)go);		break;
	case GO_PLANE3D:		delete((Plane3D*)go);		break;
	case GO_RIBBON:			delete((Ribbon*)go);		break;
	case GO_FUNCTION:		delete((Function*)go);		break;
	case GO_FITFUNC:		delete((FitFunc*)go);		break;
	case GO_LEGITEM:		delete((LegItem*)go);		break;
	case GO_LEGEND:			delete((Legend*)go);		break;
	case GO_OBJTREE:		delete((ObjTree*)go);		break;
	case GO_FREQDIST:		delete((FreqDist*)go);		break;
	case GO_GRID3D:			delete((Grid3D*)go);		break;
	case GO_FUNC3D:			delete((Func3D*)go);		break;
	case GO_XYSTAT:			delete((xyStat*)go);		break;
	case GO_FITFUNC3D:		delete((FitFunc3D*)go);		break;
	case GO_BEZIER:			delete((Bezier*)go);		break;
	case GO_TEXTFRAME:		delete((TextFrame*)go);		break;
	case GO_NORMQUANT:		delete((NormQuant*)go);		break;
	case GO_CONTOUR:		delete((ContourPlot*)go);	break;
	default:
#ifdef USE_WIN_SECURE
		sprintf_s(TmpTxt, TMP_TXT_SIZE, "Cannot delete Object\nwith Id %ld", go->Id);
#else
		sprintf(TmpTxt, "Cannot delete Object\nwith Id %ld", go->Id);
#endif
		ErrorBox(TmpTxt);
		//we do not delete the object, probably we recover
		}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//Delete a graphic object from a list 
bool DeleteGOL(GraphObj ***gol, long n, GraphObj *go, anyOutput *o)
{
	long i;
	int c;
	GraphObj **g, *p;

	if(!gol || !(*gol) || !go || !n) return false;
	for (i = 0, c = 0, g = *gol; i < n; i++, g++) {
		if(*g) {
			c++;
			if(*g == go) {
				p = (*g)->parent;
				if(o) o->HideMark();
				Undo.DeleteGO(g, 0L, o);
				if(c == 1) {
					for (g++, i++ ;i < n; i++, g++) {
						if(*g) return true;
						}
					if(p) Undo.DropMemory(p, (void**) gol, UNDO_CONTINUE);
					}
				return true;
				}
			}
		}
	return false;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//backup file before writing a new one
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool BackupFile(char *FileName)
{
	int i;
	FILE *TheFile;
	char TmpFileName[512], Name[512], ext[6];

	//no backup necessary if file does not exist
	if(!(FileExist(FileName))) return true;
	rlp_strcpy(Name, 512,FileName);
	i = (int)strlen(Name)-1;
	if(Name[--i] == '.') Name[i] = 0;
	else if(Name[--i] == '.') Name[i] = 0;
	else if(Name[--i] == '.') Name[i] = 0;
	else return false;
	i = 1;
	do {
#ifdef USE_WIN_SECURE
		sprintf_s(ext, 6, ".%03d", i);
		sprintf_s(TmpFileName, 512, "%s%s", Name, ext);
		if(!(fopen_s(&TheFile, TmpFileName, "r"))) {
			fclose(TheFile);
			}
		else break;
#else
		sprintf(ext, ".%03d", i);
		sprintf(TmpFileName, "%s%s", Name, ext);
		if((TheFile = fopen(TmpFileName, "r"))) fclose(TheFile);
#endif
		i++;
	} while (i < 999 && TheFile);
	if(i >= 999) {								//too many backups exist already
		ErrorBox("Too many backup\files exist already.");
		return false;
		}
	if(-1 == rename(FileName, TmpFileName)) {
		ErrorBox("Error accessing file.");
		return false;
		}
	return true;
}

bool IsRlpFile(char *FileName)
{
	FILE *TheFile;
	char Line[10];
	bool bRet = false;

#ifdef USE_WIN_SECURE
	if(fopen_s(&TheFile, FileName, "r")) return false;
#else
	if(0L ==(TheFile = fopen(FileName, "r"))) return false;
#endif
	fread(Line, 1, 8, TheFile);
	Line[5] = 0;
	if(0 == strcmp(Line, ";RLP "))bRet = true;
	fclose(TheFile);
	return bRet;
}

bool IsXmlFile(char *FileName)
{
	FILE *TheFile;
	char Line[10];
	bool bRet = false;

#ifdef USE_WIN_SECURE
	if(fopen_s(&TheFile, FileName, "r")) return false;
#else
	if(0L ==(TheFile = fopen(FileName, "r"))) return false;
#endif
	fread(Line, 1, 8, TheFile);
	Line[6] = 0;
	if(0 == strcmp(Line, "<?xml "))bRet = true;
	fclose(TheFile);
	return bRet;
}

bool FileExist(char *FileName)
{
	FILE *TheFile;

#ifdef USE_WIN_SECURE
	if(ENOENT == fopen_s(&TheFile, FileName, "r")) return false;
#else
	if(0L ==(TheFile = fopen(FileName, "r"))) {
		if(errno == ENOENT) return false;
		return true;
		}
#endif
	fclose(TheFile);
	return true;
}


//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//check Object for certain properties
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool IsPlot3D(GraphObj *g)
{
	if(g && (g->Id == GO_PLOT3D || g->Id == GO_FUNC3D || g->Id == GO_FITFUNC3D)) return true;
	return false;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//duplicate a block of memory
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void *memdup(void *ptr, int cb_old, int cb_new)
{
	void *p;

	if(cb_new > cb_old) {
		if((p = calloc(cb_new, 1)) && ptr) memcpy(p, ptr, cb_old);
		}
	else {
		if((p = malloc(cb_old)) && ptr) memcpy(p, ptr, cb_old);
		}
	return p;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//calculate angle from sin(angle)
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
double trig2deg(double si, double csi)
{
	double ang;

	ang = acos(csi);
	if(si < 0.0) ang *= -1.0;
	return floor(ang * 57.29577951 +0.5);
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//replace graphic object with new: typically used for undo
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool ReplaceGO(GraphObj **oldobj, GraphObj **newobj)
{
	newobj[1]->parent = newobj[0]->parent;
	newobj[1]->Command(CMD_SET_DATAOBJ, newobj[0]->data, 0L);
	*oldobj = newobj[1];
	newobj[0]->parent = 0L;		//disable messaging
	Undo.InvalidGO(newobj[0]);
	DeleteGO(newobj[0]);
	return true;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//calculate a 'unique' hash value for a string
//Ref: Corman T.H., Leiserson C.E. & Rivest R.L. (1990) Hash Functions.
//   in: Introduction to Algorithms (MIT Press & McGraw-Hill)
//   ISBN 0-262-03141-8 and ISBN 0-07-013143-0, pp. 226ff
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
unsigned int HashValue(unsigned char *str)
{
	unsigned int i = 0, ret = 0;

	if(!str || !str[0]) return 0;
	do {
		if(str[i] > 32) ret = ((str[i]-32) + (ret <<2));
		i++;
		}while(str[i]);
	return ret;
}

unsigned int Hash2(unsigned char * str)
{
	unsigned int i = 0, ret = 0, c;

	if(!str) return 0;
	do {
		c = str[i++];
		ret = ((ret * c)<<2) | c;
		}while(str[i]);
	return ret;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//compare data structures: return true if changed or save undo info
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
bool cmpLineDEF(LineDEF *l1, LineDEF *l2)
{
	if(!l1 || !l2 || l1 == l2) return false;	//oh, oh !
	if(l1->width != l2->width) return true;
	if(l1->patlength != l2->patlength) return true;
	if(l1->color != l2->color) return true;
	if(l1->pattern != l2->pattern) return true;
	return false;
}

bool cmpFillDEF(FillDEF *f1, FillDEF *f2)
{
	if(!f1 || !f2 || f1 == f2) return false;	//oh, oh!
	if(f1->type != f2->type || f1->color != f2->color ||
		f1->scale != f2->scale || f1->color2 != f2->color2) return true;
	//the hatch line is subject to a separate call to cmpLineDEF
	return false;
}

bool cmpAxisDEF(AxisDEF *a1, AxisDEF *a2)
{
	int i;

	if(!a1 || !a2 || a1 == a2) return false;	//oh, oh!
	if(a1->flags != a2->flags || a1->min != a2->min || a1->max != a2->max ||
		a1->loc[0].fx != a2->loc[0].fx || a1->loc[0].fy != a2->loc[0].fy ||
		a1->loc[0].fz != a2->loc[0].fz || a1->loc[1].fx != a2->loc[1].fx ||
		a1->loc[1].fy != a2->loc[1].fy || a1->loc[1].fz != a2->loc[1].fz ||
		a1->Start != a2->Start || a1->Step != a2->Step || a1->Radius != a2->Radius ||
		a1->Center.fx != a2->Center.fx || a1->Center.fy != a2->Center.fy ||
		a1->nBreaks != a2->nBreaks) return true;
	for(i = 0; i < a1->nBreaks; i++) {
		if(a1->breaks[i].fx != a2->breaks[i].fx ||
			a1->breaks[i].fy != a2->breaks[i].fy) return true;
		}
	return false;
}

bool cmpTextDEF(TextDEF *t1, TextDEF *t2)
{
	if(!t1 || !t2) return false;
	if(t1->ColTxt != t2->ColTxt || t1->ColBg != t2->ColBg || t1->fSize != t2->fSize ||
		t1->RotBL != t2->RotBL || t1->RotCHAR != t2->RotCHAR || t1->iSize != t2->iSize ||
		t1->Align != t2->Align || t1->Mode != t2->Mode || t1->Style != t2->Style ||
		t1->Font != t2->Font) return true;
	if((!(t1->text) && (t2->text)) || (!(t2->text) && (t1->text))) return true;
	if(t1->text && t2->text && strcmp(t1->text, t2->text)) return true;
	return false;
}

// Dialog Undo utilitites
DWORD CheckNewFloat(double *loc, double old_v, double new_v, GraphObj *par, DWORD flags)
{
	if(loc && old_v != new_v) {
		Undo.ValFloat(par, loc, flags);
		*loc = new_v;						return (flags | UNDO_CONTINUE);
		}
	return flags;
}

DWORD CheckNewInt(int *loc, int old_v, int new_v, GraphObj *par, DWORD flags)
{
	if(loc && old_v != new_v) {
		Undo.ValInt(par, loc, flags);
		*loc = new_v;						return (flags | UNDO_CONTINUE);
		}
	return flags;
}

DWORD CheckNewDword(DWORD *loc, DWORD old_v, DWORD new_v, GraphObj *par, DWORD flags)
{
	if(loc && old_v != new_v) {
		Undo.ValDword(par, loc, flags);
		*loc = new_v;						return (flags | UNDO_CONTINUE);
		}
	return flags;
}

DWORD CheckNewLFPoint(lfPOINT *loc, lfPOINT *old_v, lfPOINT *new_v, GraphObj *par, DWORD flags)
{
	if(loc && old_v && new_v && (old_v->fx != new_v->fx || old_v->fy != new_v->fy)) {
		Undo.SaveLFP(par, loc, flags);
		loc->fx = new_v->fx;	loc->fy = new_v->fy;	return (flags | UNDO_CONTINUE);
		}
	return flags;
}

DWORD CheckNewString(char **loc, char *s_old, char *s_new, GraphObj *par, DWORD flags)
{
	int ocb, ncb, cb;

	if(s_old && s_new) {
		if(!strcmp(s_old, s_new)) return flags;
		ocb = (int)strlen(s_old);		ncb = (int)strlen(s_new);
		cb = ncb > ocb ? ncb : ocb;
		if(cb > ocb) {
			*loc = (char*)realloc(*loc, cb * sizeof(char)+1);
			}
		Undo.String(par, loc, flags);	flags |= UNDO_CONTINUE;
		if(*loc) rlp_strcpy(*loc, cb+1, s_new); 
		}
	else if(!s_old && s_new && s_new[0]) {
		Undo.String(par, loc, flags);	flags |= UNDO_CONTINUE;
		*loc = (char *)memdup(s_new, (int)strlen(s_new) +1, 0);
		}
	else if(s_old && s_old[0] && !s_new) {
		Undo.String(par, loc, flags);	flags |= UNDO_CONTINUE;
		if(*loc) *loc[0] = 0;
		}
	return flags;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//execute clipping of 3D objects
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
typedef struct {			//structure required by 3D Bresenhan's algorithm
	int d, a, s, *r, l;
	}moveDEF;

long sph_r2;
int sph_x, sph_y, sph_z, nclp_pg, seg_x_seg;
POINT3D *clp_pg, *sphlim1, *sphlim2;
double *vclp_pg = 0L;

int test_sphere(int x, int y, int z)
{
	int d;
	long ld;
	
	ld = (d = x-sph_x) * d;			ld += (d = y-sph_y) * d;
	ld += (d = z-sph_z) * d;
	return (ld > sph_r2) ? 1 : 0;
}

//use a 3D Bresenham alorithm to find z coordinates where x == lx or y == ly
bool line3D_z(POINT3D *p1, POINT3D *p2, int lx, int ly, int *cx, int *cy, int *cz)
{
	moveDEF mx, my, mz, *m1, *m2, *m3;
	int x, y, z, d1, d2;

	mx.d = p2->x - p1->x;		mx.a = mx.d >= 0 ? mx.d : -mx.d;
	mx.s = mx.d >= 0 ? 1 : -1;	mx.r = &x;		mx.l = p2->x;	x = p1->x;
	my.d = p2->y - p1->y;		my.a = my.d >= 0 ? my.d : -my.d;
	my.s = my.d >= 0 ? 1 : -1;	my.r = &y;		my.l = p2->y;	y = p1->y;
	mz.d = p2->z - p1->z;		mz.a = mz.d >= 0 ? mz.d : -mz.d;
	mz.s = mz.d >= 0 ? 1 : -1;	mz.r = &z;		mz.l = p2->z;	z = p1->z;
	if(mx.a > my.a) {
		if(mz.a > mx.a) {
			m1 = &mz;		m2 = &mx;	m3 = &my;
			}
		else if(mz.a > my.a) {
			m1 = &mx;		m2 = &mz;	m3 = &my;
			}
		else {
			m1 = &mx;		m2 = &my;	m3 = &mz;
			}
		}
	else {
		if(mz.a > my.a) {
			m1 = &mz;		m2 = &my;	m3 = &mx;
			}
		else if(mz.a > mx.a) {
			m1 = &my;		m2 = &mz;	m3 = &mx;
			}
		else {
			m1 = &my;		m2 = &mx;	m3 = &mz;
			}
		}
	d1 = m2->a - (m1->a >>1);		d2 = m3->a - (m1->a >>1);
	for(; ;) {
		//process point at (m1.r, m2.r, m3.r);
		if(x == lx || y == ly) {
			*cx = x;	*cy = y;	*cz = z;
			return true;
			}
		if(*(m1->r) == m1->l) return false;
		if(d1 >= 0) {
			*(m2->r) += m2->s;	d1 -= m1->a;
			}
		if(d2 >= 0) {
			*(m3->r) += m3->s;	d2 -= m1->a;
			}
		*(m1->r) += m1->s;	d1 += m2->a;		d2 += m3->a;
		}
}

//test if point is 1) outside, 2) above, 3) on the line, or 0) hidden by a plane 
int test_plane(int x, int y, int z)
{
	int us, ls, us1, ls1, x1, y1, z1, x2, y2, z2;
	static int last = 0;
	POINT3D p1, p2;

	if(IsInPolygon3D(x, y, clp_pg, nclp_pg, &us, &ls)) {
		if(us == ls){
			seg_x_seg = us;
			return last = 3;		//point is on line: visible
			}
		if(vclp_pg) {
			if(iround((x * vclp_pg[0] + y * vclp_pg[1] - vclp_pg[3])/vclp_pg[2]) < z)
				return last = 2;
			else return last = 0;
			}
		if(us < 1) us1 = nclp_pg -1;	else us1 = us -1;
		if(ls < 1) ls1 = nclp_pg -1;	else ls1 = ls -1;
		if(z < clp_pg[us1].z && z < clp_pg[us].z && z < clp_pg[ls1].z &&
			z < clp_pg[ls].z) return last = 0;		//far below plane
		if(z > clp_pg[us1].z && z > clp_pg[us].z && z > clp_pg[ls1].z &&
			z > clp_pg[ls].z) return last = 2;		//far above plane
		if(line3D_z(&clp_pg[us1], &clp_pg[us], x, -1, &x1, &y1, &z1) &&
			line3D_z(&clp_pg[ls1], &clp_pg[ls], x, -1, &x2, &y2, &z2)){
			if(z1 < z && z2 < z) return last = 2;
			if(z1 >= z && z2 >= z) return last = 0;
			p1.x = x1;	p1.y = y1; p1.z = z1;	p2.x = x2;	p2.y = y2; p2.z = z2;
			if(line3D_z(&p1, &p2, -1, y, &x1, &y1, &z1)) {
				if(z > z1) return last = 2;
				else return last = 0;
				}
			}
		return last;
		}
	return last = 1;
}

int test_planeandline(int x, int y, int z)
{
	int ret;

	ret = test_plane(x, y, z);
	if(ret == 3 && clp_pg && nclp_pg && idx_point_on_line < nclp_pg) {
		if(vclp_pg) {
			if(iround((x * vclp_pg[0] + y * vclp_pg[1] - vclp_pg[3])/vclp_pg[2]) < z)
				return 2;
			else return 0;
			}
		else {
			// no equation for plane available
			}
		}
	return ret;
}


void proc_polygon(int vis, POINT3D *pnt, POINT3D * last);	//prototype

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//use a 3D Bresenham algorithm to clip lines
void clip_line_3D(GraphObj *go, POINT3D *p1, POINT3D *p2, int(*proc)(int, int, int))
{
	moveDEF mx, my, mz, *m1, *m2, *m3;
	int x, y, z, d1, d2, vis;
	bool bVisible, bDrawLater = false;
	POINT3D p, lp;

	if(p1->x == p2->x && p1->y == p2->y && p1->z == p2->z) return;
	mx.d = p2->x - p1->x;		mx.a = mx.d >= 0 ? mx.d : -mx.d;
	mx.s = mx.d >= 0 ? 1 : -1;	mx.r = &x;		mx.l = p2->x;	x = p1->x;
	my.d = p2->y - p1->y;		my.a = my.d >= 0 ? my.d : -my.d;
	my.s = my.d >= 0 ? 1 : -1;	my.r = &y;		my.l = p2->y;	y = p1->y;
	mz.d = p2->z - p1->z;		mz.a = mz.d >= 0 ? mz.d : -mz.d;
	mz.s = mz.d >= 0 ? 1 : -1;	mz.r = &z;		mz.l = p2->z;	z = p1->z;
	if(mx.a > my.a) {
		if(mz.a > mx.a) {
			m1 = &mz;		m2 = &mx;	m3 = &my;
			}
		else if(mz.a > my.a) {
			m1 = &mx;		m2 = &mz;	m3 = &my;
			}
		else {
			m1 = &mx;		m2 = &my;	m3 = &mz;
			}
		}
	else {
		if(mz.a > my.a) {
			m1 = &mz;		m2 = &my;	m3 = &mx;
			}
		else if(mz.a > mx.a) {
			m1 = &my;		m2 = &mz;	m3 = &mx;
			}
		else {
			m1 = &my;		m2 = &mx;	m3 = &mz;
			}
		}
	bVisible = (0 != (vis = (*proc)(x, y, z)));
	if((bDrawLater = (vis == 2)) && go) go->Command(CMD_DRAW_LATER, 0L, 0L);
	lp.x = p1->x;	lp.y = p1->y;	lp.z = p1->z;
	if(!go && vis) proc_polygon(vis, p1, p1);
	d1 = m2->a - (m1->a >>1);		d2 = m3->a - (m1->a >>1);
	for(; ;) {
		//process point at (m1.r, m2.r, m3.r);
		vis = (*proc)(x, y, z);
		if(!bDrawLater && vis == 2){
			if(go) go->Command(CMD_DRAW_LATER, 0L, 0L);
			bDrawLater = true;
			}
		if(bVisible) {
			if(!vis) {
				p.x = x;	p.y = y;	p.z = z;
				if(go) go->Command(CMD_ADDTOLINE, &lp, 0L);
				else proc_polygon(vis, &p, &lp);
				bVisible = false;
				}
			}
		else if(vis){
			p.x = x;	p.y = y;	p.z = z;
			if(go) go->Command(CMD_STARTLINE, &p, 0L);
			else proc_polygon(vis, &p, &lp);
			bVisible = true;
			}
		if(*(m1->r) == m1->l){
			if(vis){
				p.x = x;	p.y = y;	p.z = z;
				if(go) go->Command(CMD_ADDTOLINE, &p, 0L);
				else proc_polygon(vis, &p, &lp);
				}
			return;
			}
		lp.x = x;	lp.y = y;	lp.z = z;
		if(d1 >= 0) {
			*(m2->r) += m2->s;	d1 -= m1->a;
			}
		if(d2 >= 0) {
			*(m3->r) += m3->s;	d2 -= m1->a;
			}
		*(m1->r) += m1->s;	d1 += m2->a;		d2 += m3->a;
		}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//use circular Bresenham's algorithm to clip a spherical scanline
//Ref: C. Montani, R. Scopigno (1990) "Speres-To-Voxel Conversion", in:
//   Graphic Gems (A.S. Glassner ed.) Academic Press, Inc.; 
//   ISBN 0-12-288165-5 
void clip_spher_line(GraphObj *go, POINT3D *cent, int r, bool bVert, int(*proc)(int, int, int))
{
	int x, y, q, di, de, lim;
	int vis = 0, o_vis;
	POINT3D cpos;

	if(r < 1) return;
	cpos.x = cent->x;	cpos.y = cent->y;	cpos.z = cent->z;
	if(bVert) cpos.y -= r;
	else cpos.x -= r;
	for(q = 0; q < 2; q++) {
		x = lim = 0;	y = r;	di = 2*(1-r);
		while (y >= lim){
			o_vis = vis;
			if(bVert && (cpos.x < sphlim1->x || cpos.x > sphlim2->x)) vis = 0;
			else if (cpos.y < sphlim1->y || cpos.y > sphlim2->y) vis = 0;
			else if (cpos.x > 0 && cpos.y >0) vis = (*proc)(cpos.x, cpos.y, cpos.z);
			if(vis != o_vis) {
				if(vis) go->Command(CMD_STARTLINE, &cpos, 0L);
				else go->Command(CMD_ADDTOLINE, &cpos, 0L);
				}
			if(di < 0) {
				de = 2*di + 2*y -1;
				if(de > 0) {
					x++;	y--;	di += (2*x -2*y +2);
					}
				else {
					x++;	di += (2*x +1);
					}
				}
			else {
				de = 2*di -2*x -1;
				if(de > 0) {
					y--;	di += (-2*y +1);
					}
				else {
					x++;	y--;	di += (2*x -2*y +2);
					}
				}
			switch(q) {
			case 0:
				if(bVert) cpos.y = cent->y - y;
				else cpos.x = cent->x - y;
				cpos.z = cent->z + x;
				break;
			case 1:
				if(vis && y < lim) go->Command(CMD_ADDTOLINE, &cpos, 0L);
				if(bVert) cpos.y = cent->y + x;
				else cpos.x = cent->x + x;
				cpos.z = cent->z + y;
				break;
				}
			}
		}
	return;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//come here to process polygon clipping events
int ppg_vis, ppg_level, ppg_reason, ppg_firstvis, ppg_nact, ppg_nmask;
bool ppg_nowvis;
POINT *ppg_mask;
POINT3D ppg_first, *ppg_act;
GraphObj *ppg_par;
double *ppg_vec = 0L;

// test if point is on line of 3D polygon
//Ref: P.S. Heckbert (1990) "Digital Line Drawing", in: Graphic Gems
//   (A.S. Glassner, ed.); Academic Press, Inc.,
//   ISBN 0-12-286165-5
bool IsOnLine(POINT *p1, POINT *p2, int x, int y)
{
	int d, ax, ay, sx, sy, dx, dy;
	POINT tr;

	dx = p2->x - p1->x;
	if ( p2->x < p1->x) { 	ax = (-dx)<<1;		sx = -1;		}
	else {					ax = dx <<1;		sx = 1;		}
	dy = p2->y - p1->y;
	if (p2->y < p1->y) {	ay = (-dy)<<1;		sy = -1;		}
	else {					ay = dy<<1;			sy = 1;		}
	tr.x = p1->x;			tr.y = p1->y;
	if (ax > ay) {									// x dominant
		d = ay - (ax >>1);
		for ( ; ; ) {
			if(tr.y == y && tr.x == x) return true;	//match
			if(tr.x == p2->x) return false;
			if (d >= 0) {tr.y += sy;		d -= ax;}
			tr.x += sx;		d += ay;
			}
		}
	else {											// y dominant
		d = ax - (ay >>1);
		for ( ; ; ) {
			if(tr.x == x && tr.y == y) return true;	//match
			if (tr.y == p2->y) return false;
			if (d >= 0) {tr.x += sx;		d -= ay;}
			tr.y += sy;		d += ax;
			}
		}
}

//use Bresenham's algorithm to complete a partially hidden polygon
//   by the shadow of the above surface
//Ref: P.S. Heckbert (1990) "Digital Line Drawing", in: Graphic Gems
//   (A.S. Glassner, ed.); Academic Press, Inc.,
//   ISBN 0-12-286165-5
bool ShadowPolygon(POINT *p1, POINT *p2, POINT *tp, int ntp, POINT *pts, long *cp, POINT *lim)
{
	int d, ax, ay, sx, sy, dx, dy;
	bool bPen, vis;
	POINT tr, mp;
	long d1, ndist, dist = 99999;

	dx = p2->x - p1->x;
	if ( p2->x < p1->x) { 	ax = (-dx)<<1;		sx = -1;		}
	else {					ax = dx <<1;		sx = 1;		}
	dy = p2->y - p1->y;
	if (p2->y < p1->y) {	ay = (-dy)<<1;		sy = -1;		}
	else {					ay = dy<<1;			sy = 1;		}
	tr.x = p1->x;			tr.y = p1->y;
	bPen = IsInPolygon(&tr, tp, ntp) || IsCloseToPL(tr, tp, ntp);
	AddToPolygon(cp, pts, p1);				//first Point is always visible
	if (ax > ay) {							// x dominant
		d = ay - (ax >>1);
		for ( ; ; ) {
			if(abs(tr.y - lim->y) < 3 && abs(tr.x - lim->x) < 3) {
				ndist = (d1=(tr.y - lim->y))*d1;
				ndist += ((d1=(tr.x - lim->x))*d1);
				if(ndist <= dist) {
					mp.x = tr.x;	mp.y = tr.y;	dist = ndist;
					}
				else {
					//mp is the closest point
					AddToPolygon(cp, pts, &mp);
					return false;
					}
				}
			vis = IsInPolygon(&tr, tp, ntp) || IsCloseToPL(tr, tp, ntp);
			if(bPen && !vis) {
				AddToPolygon(cp, pts, &tr);
				return false;				//now leaving the polygon
				}
			if(tr.x == p2->x) return true;	//still inside
			if (d >= 0) {tr.y += sy;		d -= ax;}
			bPen = vis;		tr.x += sx;		d += ay;
			}
		}
	else {									// y dominant
		d = ax - (ay >>1);
		for ( ; ; ) {
			if(abs(tr.x - lim->x) < 3 && abs(tr.y - lim->y) < 3) {
				ndist = (d1=(tr.y - lim->y))*d1;
				ndist += ((d1=(tr.x - lim->x))*d1);
				if(ndist <= dist) {
					mp.x = tr.x;	mp.y = tr.y;	dist = ndist;
					}
				else {
					//mp is the closest point
					AddToPolygon(cp, pts, &mp);
					return false;
					}
				}
			vis = IsInPolygon(&tr, tp, ntp) || IsCloseToPL(tr, tp, ntp);
			if(bPen && !vis) {
				AddToPolygon(cp, pts, &tr);
				return false;				//now leaving the polygon
				}
			if (tr.y == p2->y) return true;	//still inside
			if (d >= 0) {tr.x += sx;		d -= ay;}
			bPen = vis;		tr.y += sy;		d += ax;
			}
		}
}

//find segment which is closest to point
int FindClosestSeg(POINT3D *pg, int npg, int x, int y, int start)
{
	int i, i1, j, tmp, idx = -2, dx, dy, d, dist = 10000;
	POINT p1, p2;

	if(start < 1) start = 1;
	for(j = start + npg, i1 = start; i1 < j; i1++) {
		i = ((i1-1)%npg)+1;
		if( i == npg) {
			p1.x = pg[i-1].x;	p1.y = pg[i-1].y;
			p2.x = pg[0].x;		p2.y = pg[0].y;
			}
		else {
			p1.x = pg[i-1].x;	p1.y = pg[i-1].y;
			p2.x = pg[i].x;		p2.y = pg[i].y;
			}
		if(p1.x != p2.x || p1.y != p2.y) {
			d = dist;
			if(abs(dx = (p2.x - p1.x)) < abs(dy = (p2.y - p1.y))) {	//y dominant
				if(dy && ((p1.y >= y && p2.y < y) || (p2.y > y && p1.y <= y))) {
					tmp = (p1.x + ((y - p1.y) * dx)/dy);	d = abs(x-tmp);
					}
				}
			else {								// x dominant
				if(dx && ((p1.x >= x && p2.x < x) || (p2.x > x && p1.x <= x))) {
					tmp = (p1.y + ((x - p1.x) * dy)/dx);	d = abs(y-tmp);
					}
				}
			if(d < dist) {
				dist = d;	idx = i;
				}
			}
		if(dist < 3) break;
		}
	return idx;
}

//finish a partially visible 3D polygon by its shadow of the above
//   3D polygon
bool AddShadowPolygon(POINT3D *pnt, POINT3D *ep, int cidx) {
	int us, ls, i, j, k, x1, x2, y1, y2, z1, z2, idx_clppg_line;
	long cpg1 = 0, d, d1, d2, ntp=0, ntp1=0, ntp2=0;
	POINT *pg1 = 0L, np, *tp=0L, *tp1, *tp2, lim;	
	POINT3D nep;
	bool bRet = false;

	d = ((d1 = (ep->x - pnt->x))*d1);
	d += ((d1 = (ep->y - pnt->y))*d1);
	if(d < 4) {			//propably too close
		if(d && ppg_par) ppg_par->Command(CMD_ADDTOLINE, ep, 0L);
		return true;	//connect
		}
	lim.x = ep->x;		lim.y = ep->y;
	idx_clppg_line = FindClosestSeg(clp_pg, nclp_pg, pnt->x, pnt->y, cidx);
	//create track from hiding polygon
	//the ppoint 'pnt' is expected to be on the line
	//   clp_pg[idx_clpgg_line] and clp_pg[idx_clpgg_line - 1]
	if(!(pg1 = (POINT*)calloc(nclp_pg +4, sizeof(POINT))))return true;
	np.x = pnt->x;	np.y = pnt->y;	AddToPolygon(&cpg1, pg1, &np);
	j = idx_clppg_line + nclp_pg;
	for(i = idx_clppg_line; i < j; i++) {
		np.x = clp_pg[k = (i%nclp_pg)].x;		np.y = clp_pg[k].y;
		AddToPolygon(&cpg1, pg1, &np);
		}
	//close polygon
	np.x = pnt->x;	np.y = pnt->y;	AddToPolygon(&cpg1, pg1, &np);
	//calculate two possible solutions
	tp1 = (POINT*)calloc(nclp_pg+4, sizeof(POINT));
	tp2 = (POINT*)calloc(nclp_pg+4, sizeof(POINT));
	if(!tp1 || !tp2) {			//memory allocation error
		free(pg1);	free(ppg_mask);	return false;
		}
	ShadowPolygon(&pg1[0], &pg1[1], ppg_mask, (int)ppg_nmask, tp1, &ntp1, &lim);
	if(ntp1 == 1){				//more than one segment
		for(i = 2; i < cpg1; i++) {
			if(!ShadowPolygon(&pg1[i-1], &pg1[i], ppg_mask, (int)ppg_nmask, tp1, &ntp1, &lim))
				break;
			}
		}
	if(i == cpg1) {				//last segment required
		ShadowPolygon(&pg1[i-1], &pg1[0], ppg_mask, (int)ppg_nmask, tp1, &ntp1, &lim);
		}
	ShadowPolygon(&pg1[0], &pg1[cpg1-1], ppg_mask, (int)ppg_nmask, tp2, &ntp2, &lim);
	if(ntp2 == 1){				//more than one segment
		for(i = cpg1-1; i > 1; i--) {
			if(!ShadowPolygon(&pg1[i], &pg1[i-1], ppg_mask, (int)ppg_nmask, tp2, &ntp2, &lim))
				break;
			}
		}
	//find better solution
	d1 = ((d = (ep->x - tp1[ntp1-1].x))*d);		d1 += ((d = (ep->y - tp1[ntp1-1].y))*d);
	d2 = ((d = (ep->x - tp2[ntp2-1].x))*d);		d2 += ((d = (ep->y - tp2[ntp2-1].y))*d);
	if(d1 < d2 && d1 < 5) {			//use solution 1
		tp = tp1;	ntp = ntp1;
		}
	else if(d2 < d1 && d2 < 5) {	//use solution 2
		tp = tp2;	ntp = ntp2;
		}
	else if (d1 == d2 && d1 < 5) {			//ambiguous result: connect stright 
		if(d && ppg_par) ppg_par->Command(CMD_ADDTOLINE, ep, 0L);
		}
	else {							//no valid solution;
		if(cidx >= 0) return AddShadowPolygon(pnt, ep, -2);
		bRet = false;
		}
	if(tp && ntp>1) {				//create shadow line
		bRet = true;
		for(i = 1; i < ntp; i++) {
			if(i == ntp -1) {
				d = ((d1 = tp[i].x - ep->x) * d1);
				d += ((d1 = tp[i].y - ep->y) * d1);
				if(d < 2){						//too close to end point
					nep.x = ep->x = tp[i].x;	nep.y = ep->y = tp[i].y;
					nep.z = ep->z;
					break;
					}
				}
			np.x = nep.x = tp[i].x;	np.y = nep.y = tp[i].y;
			nep.z = pnt->z > ep->z ? ep->z : pnt->z;
			if(IsInPolygon(&np, ppg_mask, ppg_nmask) || 
				IsCloseToPL(np, ppg_mask, ppg_nmask)){
				if(ppg_vec) {					//valid plane eqation
					nep.z = iround((nep.x * ppg_vec[0] + nep.y * ppg_vec[1] - ppg_vec[3])/ppg_vec[2]);
					if(ppg_par) ppg_par->Command(CMD_ADDTOLINE, &nep, 0L);
					}
				else if(IsInPolygon3D(nep.x, nep.y, ppg_act, ppg_nact, &us, &ls)) {
					if(us == ls) if(ls){	//point is on the line
						j = nep.z;
						if(ppg_act[ls].x == ppg_act[ls-1].x && ppg_act[ls].y == ppg_act[ls-1].y){
							nep.z = (ppg_act[ls].z + ppg_act[ls-1].z)>>1; //impropable
							}
						else if(abs(ppg_act[ls].x - ppg_act[ls-1].x) > 
							abs(ppg_act[ls].y - ppg_act[ls-1].y)){	// x dominant
							line3D_z(&ppg_act[ls-1], &ppg_act[ls], nep.x, -1, &k, &k, &j);
							}
						else {										// y dominant
							line3D_z(&ppg_act[ls-1], &ppg_act[ls], -1, nep.y, &k, &k, &j);
							}
						nep.z = j;
						}
					else {
						if(line3D_z(&ppg_act[ls-1], &ppg_act[ls], nep.x, -1, &x1, &y1, &z1) &&
							line3D_z(&ppg_act[us-1], &ppg_act[us], nep.x, -1, &x2, &y2, &z2)){
							nep.z = (z1 + z2)>>1;					//impropable
							}
						}
					if(ppg_par) ppg_par->Command(CMD_ADDTOLINE, &nep, 0L);
					}
				else {
					//point is inside by one algorithm but outside with another
					//try without this point
					}
				}
			}
		if(ppg_par && (nep.x != ep->x || nep.y != ep->y)) ppg_par->Command(CMD_ADDTOLINE, ep, 0L);
		}
	free(pg1);	free(tp1);	free(tp2);
	return bRet;
}

//calculate the clipping line between two planes
bool CuttingEdge(POINT3D* pt, POINT3D* np)
{
	int i, j, us1, ls1, us2, ls2;
	long d, di[2];
	POINT3D res[2];
	double v[3], s[2][3];

	if(!vclp_pg || !ppg_vec) return false;
	v[0] = vclp_pg[1]*ppg_vec[2] - vclp_pg[2]*ppg_vec[1];
	v[1] = vclp_pg[2]*ppg_vec[0] - vclp_pg[0]*ppg_vec[2];
	v[2] = vclp_pg[0]*ppg_vec[1] - vclp_pg[1]*ppg_vec[0];
	if(fabs(v[0]) < 1e-5 || fabs(v[1]) < 1e-5 || fabs(v[2]) < 1e-5) return false;
	v[0] *= (v[2]*2048.0);	v[1] *= (v[2]*2048.0);	v[2] *= (v[2]*2048.0);
	//find two solutions +/- vector
	for(i = 0; i < 2; i++) {
		s[i][0] = (double)(pt->x);	s[i][1] = (double)(pt->y);	s[i][2] = (double)(pt->z);
		for(j = 0; j < 5; j++) {
			do {
				s[i][0] += v[0];			s[i][1] += v[1];			s[i][2] += v[2];
				}while(IsInPolygon3D(s[i][0], s[i][1], ppg_act, ppg_nact, &us1, &ls1) 
						&& IsInPolygon3D(s[i][0], s[i][1], clp_pg, nclp_pg, &us2, &ls2));
			s[i][0] -= v[0];			s[i][1] -= v[1];			s[i][2] -= v[2];
			v[0] /= 4.0;	v[1] /= 4.0;	v[2] /= 4.0;
			}
		s[i][0] += v[0];			s[i][1] += v[1];			s[i][2] += v[2];
		v[0] *= -1024.0;	v[1] *= -1024.0;	v[2] *= -1024.0;
		res[i].x = iround(s[i][0]);	res[i].y = iround(s[i][1]);	res[i].z = iround(s[i][2]);
		di[i] = (d=(res[i].x - pt->x))*d;		di[i] += ((d=(res[i].y - pt->y))*d);
		}
	if(di[0] > 1 && di[0] > di[1]) {
		//first solution has longer projection
		np->x = res[0].x;	np->y = res[0].y;	np->z = res[0].z;
		return true;
		}
	if(di[1] > 1 && di[1] > di[0]) {
		//second solution has longer projection
		np->x = res[1].x;	np->y = res[1].y;	np->z = res[1].z;
		return true;
		}
	return false;
}

//come here from clip_line_3D to process changes in visibility when
//   clipping one polygon with another
void proc_polygon(int vis1, POINT3D *pnt, POINT3D *last)
{
	static POINT3D np, lp;
	long d, d1;
	int vis = vis1;
	bool spg_valid;

	switch(ppg_level){
	case 0:						//searching first visible point of polygon
		if(vis == 3) vis = 1;						//on line is visible
		if(!ppg_vis && vis && !ppg_nowvis){			//found it !
			if(ppg_par && (vis1 == 3 || ppg_vis == 3)) ppg_par->Command(CMD_SIGNAL_POL, pnt, 0L);
			ppg_nowvis = true;			ppg_first.x = pnt->x;
			ppg_first.y = pnt->y;		ppg_first.z = pnt->z;
			ppg_reason = vis;
			}
		else if(!vis && ppg_nowvis) {	//check if too short
			d = (d1 = pnt->x - ppg_first.x) * d1;
			d += (d1 = pnt->y - ppg_first.y) * d1;
			if(d < 3) ppg_nowvis = false;	//cancel first point
			}
		ppg_vis = vis;
		break;
	case 1:
		if(vis == 3) {						//on line: visible
			vis = 1;
			}
		if(ppg_first.x < 0 && ppg_first.y < 0 && vis) {
			memcpy(&ppg_first, pnt, sizeof(POINT3D));
			ppg_firstvis = vis;			lp.x = pnt->x;
			lp.y = pnt->y;				lp.z = pnt->z;
			if(ppg_par) ppg_par->Command(CMD_STARTLINE, pnt, 0L);
			}
		else if (ppg_vis) {			//leaving visibility or continue
			spg_valid = false;
			if(lp.x == pnt->x && lp.y == pnt->y && lp.z == pnt->z) break;
			if(vis1 == 3 && ppg_par) ppg_par->Command(CMD_SIGNAL_POL, pnt, 0L);
			if(vis && ppg_par) ppg_par->Command(CMD_ADDTOLINE, pnt, 0L);
			else if(ppg_par) ppg_par->Command(CMD_ADDTOLINE, last, 0L);
			if(!vis){				//leaving visibility
				ppg_vis = test_plane(last->x, last->y, last->z);
				if(ppg_vis == 3 && ppg_par) ppg_par->Command(CMD_SIGNAL_POL, pnt, 0L);
				if(ppg_vis == 3) ppg_vis = 1;
				if(ppg_firstvis == 1 && ppg_vis == 1) {
					//from below surface to below surface
					if(!(spg_valid = AddShadowPolygon(last, &ppg_first, -2)))
						if(ppg_par) spg_valid = ppg_par->Command(CMD_ADDTOLINE, pnt, 0L);
					}
				else if(ppg_firstvis == 1 && ppg_vis == 2) {
					//from below surface enter inside surface
					if(CuttingEdge(last, &np)){
						if(ppg_par)ppg_par->Command(CMD_ADDTOLINE, &np, 0L);
						spg_valid = AddShadowPolygon(&np, &ppg_first, seg_x_seg);
						}
					else if(!(spg_valid = AddShadowPolygon(last, &ppg_first, seg_x_seg)))
						if(ppg_par) spg_valid = ppg_par->Command(CMD_ADDTOLINE, pnt, 0L);
					}
				else if(ppg_firstvis == 2 && ppg_vis == 1 && ppg_par) {
					//from inside surface to below surface
					if(CuttingEdge(&ppg_first, &np)){
						if(!(spg_valid = AddShadowPolygon(last, &np, seg_x_seg)))
							spg_valid = ppg_par->Command(CMD_REQ_POINT, &ppg_first, 0L);
						ppg_par->Command(CMD_ADDTOLINE, &np, 0L);
						}
					else if(!(spg_valid = AddShadowPolygon(last, &ppg_first, seg_x_seg)))
						spg_valid = ppg_par->Command(CMD_ADDTOLINE, pnt, 0L);
					}
				else if(ppg_firstvis == 2 && ppg_vis == 2) {
					//from inside surface to inside surface
					// nothing to do: connect straight
					spg_valid = true;
					}
				//prepare for new polygon
				if(spg_valid) {
					ppg_first.x = ppg_first.y = ppg_first.z = -1;
					}
				else {
					//we could not find a proper connection between the two points
					//   probably due to high complexity of graph or shape.
					//   We ignore part of the polygon and continue with the
					//   started shape.
					vis = ppg_vis;
					}
				}
			}
		ppg_vis = vis;				lp.x = pnt->x;
		lp.y = pnt->y;				lp.z = pnt->z;
		break;
		}
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//test if line is available in 3D-polygons: determine shared edges
bool LineInPolgon(POINT3D *p1, POINT3D *p2, POINT3D *tpg, int ntpg)
{
	int i;

	for(i = 1; i < ntpg; i++) {
		if(p2->x == tpg[i].x && p2->y == tpg[i].y && p2->z == tpg[i].z) {
			if(p1->x == tpg[i-1].x && p1->y == tpg[i-1].y && p1->z == tpg[i-1].z) return true;
			}
		if(p1->x == tpg[i].x && p1->y == tpg[i].y && p1->z == tpg[i].z) {
			if(p2->x == tpg[i-1].x && p2->y == tpg[i-1].y && p2->z == tpg[i-1].z) return true;
			}
		}
	i = ntpg -1;
	if(p1->x == tpg[i].x && p1->y == tpg[i].y && p1->z == tpg[i].z &&
		p2->x == tpg[0].x && p2->y == tpg[0].y && p2->z == tpg[0].z) return true;
	if(p2->x == tpg[i].x && p2->y == tpg[i].y && p2->z == tpg[i].z &&
		p1->x == tpg[0].x && p1->y == tpg[0].y && p1->z == tpg[0].z) return true;
	return false;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//entry points for clipping requests

void clip_line_sphere(GraphObj *par, POINT3D **li, int r, int cx, int cy, int cz)
{
	sph_r2 = r*r;	sph_x = cx;		sph_y = cy;		sph_z = cz;
	if(test_sphere(li[0]->x, li[0]->y, li[0]->z)) par->Command(CMD_STARTLINE, li[0], 0L);
	clip_line_3D(par, &li[0][0], &li[0][1], test_sphere);
}

void clip_line_plane(GraphObj *par, POINT3D **li, POINT3D *pg, int np, double *vec)
{
	nclp_pg = np;		clp_pg = pg;	vclp_pg = vec;
	if(test_plane(li[0]->x, li[0]->y, li[0]->z)) par->Command(CMD_STARTLINE, li[0], 0L);
	clip_line_3D(par, &li[0][0], &li[0][1], test_plane);
	vclp_pg = 0L;
}

void clip_sphline_sphere(GraphObj *par, POINT3D *lim1, POINT3D *lim2, POINT3D *cent, 
	int r1, int r2, int cx, int cy, int cz)
{
	sphlim1 = lim1;		sphlim2 = lim2;
	sph_r2 = r2*r2;	sph_x = cx;		sph_y = cy;		sph_z = cz;	
	clip_spher_line(par, cent, r1, lim1->x == lim2->x, test_sphere);
}

void clip_sphline_plane(GraphObj *par, POINT3D *lim1, POINT3D *lim2, POINT3D *cent, 
	int r1, POINT3D *pg, int np, double *vec)
{
	sphlim1 = lim1;		sphlim2 = lim2;		nclp_pg = np;
	clp_pg = pg;		vclp_pg = vec;
	clip_spher_line(par, cent, r1, lim1->x == lim2->x, test_planeandline);
	vclp_pg = 0L;
}

void clip_plane_plane(GraphObj *par, POINT3D *pg1, int n1, double *v1, POINT3D *pg2, 
	int n2, double *v2, POINT *mask, int nmask)
{
	int i, j;
	POINT3D sp;

	nclp_pg = n2;	clp_pg = pg2;	ppg_level = 0;	ppg_nowvis = false;
	ppg_reason = 0;	vclp_pg = v2;	ppg_par = par;
    ppg_vis = test_plane(pg1[0].x, pg1[0].y, pg1[0].z);
	ppg_act = pg1;	ppg_nact = n1;	ppg_vec = v1;
	for(i = 1; i < n1 && !ppg_nowvis; i++){		//assume first pnt == last pnt
		if(!LineInPolgon(&pg1[i-1], &pg1[i], pg2, n2))
			clip_line_3D(0L, &pg1[i-1], &pg1[i], test_plane);
		}
	if(!ppg_nowvis && !ppg_vis) {				//complete pg hidden
		ppg_vec = vclp_pg = 0L;		return;
		}
	if(ppg_nowvis) {							//clip this polygon
		ppg_mask = mask;		ppg_nmask = nmask;
		ppg_level = 1;
		sp.x = ppg_first.x;		sp.y = ppg_first.y;
		sp.z = ppg_first.z;		ppg_first.x = ppg_first.y = ppg_first.z = -1;
		ppg_vis = test_plane(sp.x, sp.y, sp.z);
		proc_polygon(ppg_vis, &sp, &sp);
		clip_line_3D(0L, &sp, &pg1[i-1], test_plane);
		for(j = i+n1-1; i < j; i++) {
			seg_x_seg = -2;
			clip_line_3D(0L, &pg1[(i-1)%n1], &pg1[i%n1], test_plane);
			}
		clip_line_3D(0L, &pg1[(i-1)%n1], &sp, test_plane);
		}
	else {										//all visible
		par->Command(CMD_STARTLINE, pg1, 0L);
		for(i = 1; i < n1; i++) {
			par->Command(CMD_ADDTOLINE, &pg1[i], 0L);
			}
		}
	ppg_vec = vclp_pg = 0L;
}