external/NonLinLoc/src/GridGraphLib.c

/*
 * Copyright (C) 1999 Anthony Lomax <lomax@faille.unice.fr>
 *
 * This program 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.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 */


/* Grid graphics functions

	07OCT1997	AJL	modified from wls graphics functions

*/

#define EXTERN_MODE 1
#define INPUT_DEV 1
#include "GridLib.h"
#include "GridGraphLib.h"
#include "vector/vector.h"


// define globals
int gr_num_sta;		/* number of stations read */
struct GrStation sta_array[MAX_NUM_STA];
int gr_num_txt;		/* number of text strings read */
struct GrTxt txt_array[MAX_NUM_TXT];
int gr_num_line;		/* number of lines read */
struct GrLine line_array[MAX_NUM_LINE];
int gr_num_map_files;	/* number of map files */
struct GrMapLines mapfile[MAX_NUM_MAP];



/*** function to set constants */

void SetGraphicConstants(void)
{
	gr_num_sta = 0;			/* number of stations to draw */
	gr_num_txt = 0;			/* number of text strings to draw */
	gr_num_line = 0;		/* number of lines to draw */
	gr_num_map_files = 0;		/* number of map lines to draw */
}


/* function to read input file */

int ReadGraphicsInput(FILE* fp_input, char* param, char* line, int istat)
{


	istat = -1;

	/* read parameter line */


	/* skip comment line or white space */

	if (strncmp(param, "#", 1) == 0 || isspace(param[0]))
		istat = 0;


#if 0
	/* read station params */
	if (strncmp(param, "STATION", 7) == 0)
		if ((istat = get_stations(strchr(line, ' ') + 1, fp_input)) < 0)
			fprintf(stderr, "Error reading station parameters.\n");

	/* read graph text params */
	if (strncmp(param, "GTEXT", 5) == 0)
		if ((istat = get_gtext(strchr(line, ' ') + 1)) < 0)
			fprintf(stderr, "Error reading graphics text parameters.\n");

	/* read graph line params */
	if (strncmp(param, "GLINE", 5) == 0)
		if ((istat = get_gline(strchr(line, ' ') + 1, fp_input)) < 0)
			fprintf(stderr, "Error reading graphics line parameters.\n");

#endif

	/* read mapline params */
	if (strncmp(param, "MAPLINE", 7) == 0)
		if ((istat = get_maplines(strchr(line, ' ') + 1)) < 0)
			fprintf(stderr, "Error reading mapline parameters.\n");

	return(istat);


}




/*** function to read mapline parameters from input line ***/

int get_maplines(char* input_line)
{

		/* add map file to mapfile array */

	if (gr_num_map_files == MAX_NUM_MAP) {
		fprintf(stderr, "Error - Too many map line files.\n");
		return (-1);
	}

		/* read mapline input line */

	sscanf(input_line, "%s %s %lf %lf %lf %s",
		mapfile[gr_num_map_files].format,
		mapfile[gr_num_map_files].name,
		&(mapfile[gr_num_map_files].rgb.r),
		&(mapfile[gr_num_map_files].rgb.g),
		&(mapfile[gr_num_map_files].rgb.b),
		mapfile[gr_num_map_files].line_style);

	sprintf(MsgStr, "MAPLINE  fmt %s  `%s' ",
      		mapfile[gr_num_map_files].format,
		mapfile[gr_num_map_files].name);
 	nll_putmsg(2, MsgStr);
	sprintf(MsgStr, "  col r %f g %f b %f   line_style %s",
       		mapfile[gr_num_map_files].rgb.r,
		mapfile[gr_num_map_files].rgb.g,
		mapfile[gr_num_map_files].rgb.b,
		mapfile[gr_num_map_files].line_style);
 	nll_putmsg(2, MsgStr);

	gr_num_map_files++;

    return (0);
}




/*** function to draw map lines */

#if 0 /*!!!!!!!!!!*/
#define MAXLINEPTS 10000
int PlotMapLines()
{
	int  nmapfile, npts;
	char chardummy;
	double maplong, maplat, xtemp, ytemp, xloc, yloc;
	double c111, angle, cosang, sinang;
	double xpt[MAXLINEPTS], ypt[MAXLINEPTS];
	FILE *fp_map;

	if (!prog_mode_3d)
		return;

	c111 = 10000.0 / 90.0;

	      /* draw map lines */

 	set_drawing_window(norm_trans_map);

	for (nmapfile = 0; nmapfile < gr_num_map_files; nmapfile++) {

    		set_line_color(mapfile[nmapfile].rgb, );
		set_line_style(mapfile[nmapfile].line_style);

			/* open input file */

		if ((fp_map = fopen(mapfile[nmapfile].name, "r")) == NULL) {
			fprintf(stderr, "Error - Cannot open map file `%s'.\n",
				mapfile[nmapfile].name);
    	}

		angle = -cRPD * mapfile[nmapfile].map_rot;
		cosang = cos(angle);
		sinang = sin(angle);

			/* draw each segment */

    	while (fgets(line, MAXLINE, fp_map) != NULL) {

			npts = 0;
			for ( ; ; ) {
				fscanf(fp_map, "%lf %lf \n", &maplong, &maplat);
				if (maplong > 200.0 || maplong < -200.0) {
					if (npts > 1)
						drawpolyline(npts, xpt, ypt);
					break;
				}
				if (npts < MAXLINEPTS) {
					xtemp = (maplong - mapfile[nmapfile].map_orig_long)
						* c111 * cos(cRPD * maplat);
					ytemp = (maplat - mapfile[nmapfile].map_orig_lat) * c111;
					xpt[npts] = xtemp * cosang - ytemp * sinang;
					ypt[npts++] = ytemp * cosang + xtemp * sinang;
				}
			}
		}
	}

	set_drawing_window(norm_trans_all);
    setline(SOLID);

}
#endif  /*!!!!!!!!!!*/



/*** function to draw stations on model */

#if 0  /*!!!!!!!!!!*/
int Stations2GMT()
{
	int  nsta;
	int slen;
	static double xpos[8] = {-1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 1.0, 0.0};
	static double ypos[8] = {-1.0, 1.0, 1.0, -1.0, 0.0, 1.5, 0.0, -1.5};
	static double len[8] = {-1.0, -1.0, 0.0, 0.0, -1.0, -0.5, 0.0, -0.5};
	double side;

	/* draw and label each station */

	side = (xmodmax - xmodmin) / 80.0;

	for (nsta = 0; nsta < gr_num_sta; nsta++) {
		set_text_color(sta_array[nsta].col);
		reset_char_hgt(sta_array[nsta].txt_size);
		drawtriangle(sta_array[nsta].xloc, sta_array[nsta].zloc + 1.0,
			side, 0, get_col(sta_array[nsta].col), 1,
			get_col(sta_array[nsta].col));
		if (prog_mode_3d) {
			set_drawing_window(norm_trans_map);
			drawtriangle(sta_array[nsta].xloc, sta_array[nsta].yloc,
				side, 0, get_col(sta_array[nsta].col), 1,
				get_col(sta_array[nsta].col));
				gtextw(sta_array[nsta].xloc - chrwthw / 4.0 +
					xpos[sta_array[nsta].txt_pos] * chrwthw +
					len[sta_array[nsta].txt_pos] * chrwthw *
				(double) (slen = strlen(sta_array[nsta].name)),
					sta_array[nsta].yloc - chrhtw / 2.0 +
					ypos[sta_array[nsta].txt_pos] * chrhtw,
					sta_array[nsta].name);
			set_drawing_window(norm_trans_all);
		}
		reset_char_hgt(1.0);
	}
}
#endif  /*!!!!!!!!!!*/


/*** function to draw user text on model */

#if 0  /*!!!!!!!!!!*/
int Text2GMT()
{
	int  ntxt;

	for (ntxt = 0; ntxt < gr_num_txt; ntxt++) {
		set_text_color(txt_array[ntxt].col);
		reset_char_hgt(txt_array[ntxt].size);
		gcharup(-sin(cRPD * txt_array[ntxt].rot),
			cos(cRPD * txt_array[ntxt].rot));
		gchar_rot(txt_array[ntxt].rot);
		if (prog_mode_3d && strncmp(txt_array[ntxt].plot, "MAP", 3) == 0) {
			set_drawing_window(norm_trans_map);
			gtextw(txt_array[ntxt].xloc, txt_array[ntxt].yloc,
				txt_array[ntxt].text);
			set_drawing_window(norm_trans_all);
		} else if (strncmp(txt_array[ntxt].plot, "SECT", 4) == 0) {
			gtextw(txt_array[ntxt].xloc, -txt_array[ntxt].zloc,
				txt_array[ntxt].text);
		} else if (strncmp(txt_array[ntxt].plot, "TIME", 4) == 0) {
			set_clip(0);
			gtextw(txt_array[ntxt].xloc,
				ttowind(0.0, (double) txt_array[ntxt].yloc),
				txt_array[ntxt].text);
			set_clip(1);
		}
		gcharup(0.0, 1.0);
		gchar_rot(0.0);
		reset_char_hgt(1.0);
	}

}
#endif  /*!!!!!!!!!!*/


/*** function to draw user lines on model */

#if 0  /*!!!!!!!!!!*/
int Line2GMT()
{
	int  nline;

	for (nline = 0; nline < gr_num_line; nline++) {
		set_line_color(line_array[nline].col);
		set_line_style(line_array[nline].style);
		if (prog_mode_3d && strncmp(line_array[nline].plot, "MAP", 3) == 0) {
			set_drawing_window(norm_trans_map);
   			 drawline(line_array[nline].x1, line_array[nline].y1,
				line_array[nline].x2, line_array[nline].y2);
			set_drawing_window(norm_trans_all);
		} else if (strncmp(line_array[nline].plot, "SECT", 4) == 0) {
   			 drawline(line_array[nline].x1, line_array[nline].z1,
				line_array[nline].x2, line_array[nline].z2);
		}
	}

	setline(SOLID);
}

#endif  /*!!!!!!!!!!*/



/*** function to get current date and time in character format ***/

char *get_timestamp()
{
	static char timestr[21];
	time_t tp;
	struct tm *ltime;

	(void) time(&tp);

	ltime = localtime(&tp);

	strftime(timestr, (size_t) 20, "%d %b %y  %H:%M", ltime);

	return(timestr);

}



/** function to write GMT psxy data to draw an error cross */

int Err2GMT(FILE* fp_io, double xcent, double ycent,
	double xerrlen, double xbarlen, double yerrlen, double ybarlen)
{
	double xloc, yloc;


	fprintf(fp_io, ">\n");

	/* cross */

	xloc = xcent;
	yloc = ycent - yerrlen / 2.0;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	yloc = ycent + yerrlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	xloc = xcent - xerrlen / 2.0;
	yloc = ycent;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	xloc = xcent + xerrlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	/* bars */

	xloc = xcent - xerrlen / 2.0;
	yloc = ycent - ybarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	yloc = ycent + ybarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	xloc = xcent + xerrlen / 2.0;
	yloc = ycent - ybarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	yloc = ycent + ybarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	xloc = xcent - xbarlen / 2.0;
	yloc = ycent + yerrlen / 2.0;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	xloc = xcent + xbarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	xloc = xcent - xbarlen / 2.0;
	yloc = ycent - yerrlen / 2.0;
	fprintf(fp_io, "%lf %lf\n", xloc, yloc);
	xloc = xcent + xbarlen / 2.0;
	fprintf(fp_io, "%lf %lf\n>\n", xloc, yloc);

	return(0);

}



/** function to convert covariance matrix to an array of points defining  an
		error ellipse (covariance ellipse) */

/* from "Bivariate Normal Distribution Web page
	http://www.cern.ch/Physics/DataAnalysis/BriefBook/AN10pp/node17.html */

Vect2D* Cov2Ellipse(double covxx,  double covyy, double covxy,
		double xcent, double ycent, Vect2D* xyarray, int npts)
{
	int n;

	double sigma1_2, sigma1, sigma2_2, sigma2, rho, rho_2;
	double angle, cosa, cosa_2, sina, sina_2, p1, p2;
	double theta, dtheta, xloc, yloc;


	/* sigma**2, sigma */
	sigma1_2 = covxx;
	sigma1 = sqrt(sigma1_2);
	sigma2_2 = covyy;
	sigma2 = sqrt(sigma2_2);

	/* corellation coefficient */
	rho = covxy / (sigma1 * sigma2);
	rho_2 = rho * rho;

	/* angle between x-axis and p1 */
	angle = 0.5 * atan2(2.0 * rho * sigma1 * sigma2, sigma1_2 - sigma2_2);
	cosa = cos(angle);
	cosa_2 = cosa * cosa;
	sina = sin(angle);
	sina_2 = sina * sina;

	/* semi-diameters */
	p1 = sqrt( sigma1_2 * sigma2_2 * (1.0 - rho_2) /
		(sigma2_2 * cosa_2 - 2.0 * rho * sigma1 * sigma1 * sina * cosa
			+ sigma1_2 * sina_2) );
	p2 = sqrt( sigma1_2 * sigma2_2 * (1.0 - rho_2) /
		(sigma2_2 * sina_2 + 2.0 * rho * sigma1 * sigma1 * sina * cosa
			+ sigma1_2 * cosa_2) );

	/* prob inside ellipse is ~39%, multiply by 2 to get ~??% */

	p1 *= 2.0;
	p2 *= 2.0;

	/* generate points defining ellipse */

	dtheta = 2.0 * cPI / (double) (npts - 1);
	for (n = 0, theta = 0.0; n < npts; n++, theta += dtheta)
	{
		/* unrotated ellipse */
		xloc = p1 * cos(theta);
		yloc = p2 * sin(theta);
		/* rotate by angle */
		xyarray[n].x = xcent + xloc * cosa - yloc * sina;
		xyarray[n].y = ycent + yloc * cosa + xloc * sina;
	}

	return(xyarray);

}



/** function to write array of Vect2D to GMT psxy data */

int Vect2DArray2GMT(FILE* fp_io, Vect2D* array, int npts)
{
	int n;


	for (n = 0; n < npts; n++)
		fprintf(fp_io, "%lf %lf\n", array[n].x, array[n].y);

	fprintf(fp_io, ">\n");

	return(0);

}



/** method to convert 2 semi-axis in 3D to an array of points
		defining  an error ellipse (covariance ellipse) */

/* converted from method toEllipsoid3D in Java class Ellipsoid3D (04DEC1998) */

/* !!! remember to free ellArray */

Vect3D* toEllipsoid3D(Vect3D ax1, Vect3D ax2, Vect3D center, int npts)
{
	Vect3D *ellArray;
	double d_angle;
	double angle;
	int n;
	double cosang, sinang;


	ellArray = malloc((size_t) npts * sizeof(Vect3D));

	d_angle = 2.0 * cPI / (double) (npts - 1);
	angle = 0.0;

	for (n = 0; n < npts; n++) {

		cosang = cos(angle);
		sinang = sin(angle);
		ellArray[n].x = center.x
			+ ax1.x * cosang + ax2.x * sinang;
		ellArray[n].y = center.y
			+ ax1.y * cosang + ax2.y * sinang;
		ellArray[n].z = center.z
			+ ax1.z * cosang + ax2.z * sinang;
			angle += d_angle;
	}

	return(ellArray);

}



/** function to project 3D array of points to 2D */

Vect2D* Vect3D2To2D(Vect3D *array3d, Vect2D* array2d, int npts, int iComps)
{
	int n;

	if (iComps == 12) {
		for (n = 0; n < npts; n++)
		{
			array2d[n].x = array3d[n].x;
			array2d[n].y = array3d[n].y;
		}
	} else if (iComps == 13) {
		for (n = 0; n < npts; n++)
		{
			array2d[n].x = array3d[n].x;
			array2d[n].y = array3d[n].z;
		}
	} else if (iComps == 32) {
		for (n = 0; n < npts; n++)
		{
			array2d[n].x = array3d[n].z;
			array2d[n].y = array3d[n].y;
		}
	}

	return(array2d);

}


/*------------------------------------------------------------/ */
/* Anthony Lomax           | email: lomax@faille.unice.fr     / */
/* UMR Geosciences Azur    | web: www-geoazur.unice.fr/~lomax / */
/* 250 Rue Albert Einstein | tel: 33 (0) 4 93 95 43 25        / */
/* 06560 Valbonne, FRANCE  | fax: 33 (0) 4 93 65 27 17        / */
/*------------------------------------------------------------/ */