/*----------------------------------------------------------------------*\
 | Program to read Digital Line Graph data in optional format (.lgo)	|
 |																		|
 | Peter N. Schweitzer (USGS-GD-OCG) 									|
 | MS 955, National Center												|
 | U.S. Geological Survey												|
 | Reston, VA 22092														|
 | Tel: (703) 648-6533													|
 | FAX: (703) 648-6647													|
 | email: pschweitzer@usgs.gov										|
\*----------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "lgo.h"
#include "mif.h"
#include "bna.h"
#include "attr.h"

struct state_line {
	int id;
	int index;
	int bb,be,eb,ee;
	int links;
	};

struct state {
	int FIPS_code;
	int line_cnt;
	struct state_line *line;
	int *boundary;
	};

#define TRUE	1
#define FALSE	0

/*----------------------------------------------------------------------*\
 | External procedures
\*----------------------------------------------------------------------*/

extern char *state_name (int FIPS_code);
extern char *state_abbr (int FIPS_code);
extern char *county_name (int state_code, int county_code);

/*----------------------------------------------------------------------*\
 | External data
\*----------------------------------------------------------------------*/
/*----------------------------------------------------------------------*\
 | Public procedures
\*----------------------------------------------------------------------*/

struct node *find_node (struct category *c, int id);
struct area *find_area (struct category *c, int id);
struct line *find_line (struct category *c, int id);

/*----------------------------------------------------------------------*\
 | Public data
\*----------------------------------------------------------------------*/

struct attr search;
int verbose = FALSE;
char file_name[MAXLEN];
int target = NoTarget;

/*----------------------------------------------------------------------*\
 | Static procedures
\*----------------------------------------------------------------------*/

static void read_lgo_header (FILE *in, struct lgo_header *h);
static void read_ctrl_pts (FILE *in, struct ctrl_pt *p[], int n);
static void read_categories (FILE *in, struct category *c[], int n);
static void read_nodes (FILE *in, struct node *p[], int n);
static void read_areas (FILE *in, struct area *p[], int n);
static void read_lines (FILE *in, struct line *p[], int n);
static void find_state_boundaries (struct category *c);
static void validate_lgo (struct category *c);
static int set_const (double lat0, double lat1, double lat2, double lon0);
static int axy_ll_e (double x, double y, double *lat, double *lon);
static int all_xy_e (double *x, double *y, double lat, double lon);

/*----------------------------------------------------------------------*\
 | Static data
\*----------------------------------------------------------------------*/

/*----------------------------------------------------------------------*\
 | strncpy does not null-terminate the destination string; this does.	|
\*----------------------------------------------------------------------*/

static char *strncpy0 (char *dst, const char *src, size_t n) {
	int m;
	m = strlen (src);
	if (m < n) n = m;
	memcpy (dst,src,n);
	dst[n] = 0;				/* the real strncpy doesn't do this */
	}

/*----------------------------------------------------------------------*\
 | C functions cannot read double-precision Fortran exponents, so change|
 | the 'D' or 'd' to 'e'.												|
\*----------------------------------------------------------------------*/

static void d_to_e (char *s) {
	while (*s) {
		if (toupper(*s) == 'D') *s = 'e';
		s++;
		}
	}

/*----------------------------------------------------------------------*\
 | Given an attribute, return a pointer to its name, from the table.	|
\*----------------------------------------------------------------------*/

static char *description (struct attr a) {
	char *s = "";
	struct attr_table *t = NULL;

	if (t = attribute (a)) s = t->name;

	return (s);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/

main (int argc, char *argv[]) {
	int i,j,k,m,n;
	FILE *in, *f_geo, *f_attr;
	char input_file[MAXLEN], base_name[MAXLEN];
	char geo_file[MAXLEN],attr_file[MAXLEN];
	char string[81],*s;
	char number[21];
	double *px,*py;
	double lat0,lat1,lat2,lon0;
	struct lgo_header *header;
	struct category **category;
	struct output output;

	if (argc > 1) {

		for (i=0; i < argc; i++) {
			if (argv[i][0] == '-') {
				if (toupper(argv[i][1]) == 'M') target = MapInfo;
				if (toupper(argv[i][1]) == 'A') target = AtlasPro;
				if (toupper(argv[i][1]) == 'T') target = Tactician;
				if (argv[i][1] == 'v') verbose = TRUE;
				}
			else strcpy (input_file,argv[i]);
			}

		if (in = fopen (input_file,"r")) {

			/*----------------------------------------------------------*\
			 | Determine DLG file name.									|
			\*----------------------------------------------------------*/

			if (s = strrchr (input_file,'/')) s++;
			else s = input_file;
			strcpy (file_name,s);
			if (s = strstr (file_name,".lgo")) *s = 0;

			/*----------------------------------------------------------*\
			 | Compose output file names from the input file name.		|
			\*----------------------------------------------------------*/

			strcpy (base_name,input_file);
			if (s = strstr (base_name,".lgo")) *s = 0;
			else {
				printf ("Error: I expected the input file (%s) to have extension \".lgo\"\n",input_file);
				exit (1);
				}
			strcpy (geo_file,base_name);
			strcpy (attr_file,base_name);
			switch (target) {
				case MapInfo:
					strcat (geo_file,".mif");
					strcat (attr_file,".mid");
					break;
				case AtlasPro:
				case Tactician:
					strcat (geo_file,".bna");
					strcat (attr_file,".attr");
					break;
				default:
					strcat (geo_file,".geo");
					strcat (attr_file,".attr");
					break;
				}

			/*----------------------------------------------------------*\
			 | Read file header											|
			\*----------------------------------------------------------*/

			header = (struct lgo_header *) malloc (sizeof (struct lgo_header));
			if (header) read_lgo_header (in,header);
			else {
				puts ("Error: could not allocate memory for file header.");
				exit (1);
				}

			#ifdef DEBUG
				printf ("header:\n");
				printf ("  banner: \"%s\"\n",header->banner);
				printf ("  mapname: \"%s\"\n",header->mapname);
				printf ("  source_year: %d\n",header->source_year);
				printf ("  source_scale: %lu\n",header->source_scale);
				printf ("  DLG_level: %d\n",header->DLG_level);
				printf ("  ref_system: %d\n",header->ref_system);
				printf ("  zone: %d\n",header->zone);
				printf ("  units: %d\n",header->units);
				printf ("  resolution: %lf\n",header->resolution);
				printf ("  trans_param_cnt: %d\n",header->trans_param_cnt);
				printf ("  acc_rec_cnt: %d\n",header->acc_rec_cnt);
				printf ("  ctrl_pt_cnt: %d\n",header->ctrl_pt_cnt);
				printf ("  category_cnt: %d\n",header->category_cnt);
				printf ("  a    = %lf\n",header->a);
				printf ("  e2   = %lf\n",header->e2);
				printf ("  lat1 = %lf\n",header->lat1);
				printf ("  lat2 = %lf\n",header->lat2);
				printf ("  lon0 = %lf\n",header->lon0);
				printf ("  lat0 = %lf\n",header->lat0);
				printf ("  easting  = %lf\n",header->easting);
				printf ("  northing = %lf\n",header->northing);
				printf ("  a1   = %lf\n",header->a1);
				printf ("  a2   = %lf\n",header->a2);
				printf ("  a3   = %lf\n",header->a3);
				printf ("  a4   = %lf\n",header->a4);
			#endif

			/*----------------------------------------------------------*\
			 | Read control-point identification records				|
			\*----------------------------------------------------------*/

			header->ctrl_pt = (struct ctrl_pt **) malloc (header->ctrl_pt_cnt * sizeof (struct ctrl_pt *));
			if (header->ctrl_pt) read_ctrl_pts (in,header->ctrl_pt,header->ctrl_pt_cnt);
			else {
				puts ("Error: could not allocate control point array");
				exit (1);
				}

			#ifdef DEBUG
				printf ("Control points:\n");
				for (i=0; i < header->ctrl_pt_cnt; i++) {
					printf ("%s %12.6lf %12.6lf %12.6lf %12.6lf\n",
						header->ctrl_pt[i]->corner,
						header->ctrl_pt[i]->latitude,
						header->ctrl_pt[i]->longitude,
						header->ctrl_pt[i]->x,
						header->ctrl_pt[i]->y);
					}
			#endif

			lat0 = header->lat0;
			lat1 = header->lat1;
			lon0 = header->lon0;
			lat2 = header->lat2;
			set_const (lat0,lat1,lat2,lon0);

			#ifdef DEBUG
				for (i=0; i < header->ctrl_pt_cnt; i++) {
					px = &header->ctrl_pt[i]->x;
					py = &header->ctrl_pt[i]->y;
					printf ("Before: %12.6lf %12.6lf ...",*px,*py);
					axy_ll_e (*px,*py,py,px);
					printf ("After: %12.6lf %12.6lf\n",*px,*py);
					}
			#endif

			/*----------------------------------------------------------*\
			 | Read data category identification records				|
			\*----------------------------------------------------------*/

			category = (struct category **) malloc (header->category_cnt * sizeof (struct category *));
			if (category) read_categories (in,category,header->category_cnt);
			else {
				puts ("Error: could not allocate category array");
				exit (1);
				}

			/*----------------------------------------------------------*\
			 | Read data for each category								|
			\*----------------------------------------------------------*/

			for (i=0; i < header->category_cnt; i++) {

				/*------------------------------------------------------*\
				 | Read nodes											|
				\*------------------------------------------------------*/

				category[i]->node = (struct node **) malloc (category[i]->node_cnt * sizeof (struct node *));
				if (category[i]->node) read_nodes (in,category[i]->node,category[i]->node_cnt);
				else {
					printf ("Error: could not allocate node array for category %s\n",category[i]->name);
					exit (1);
					}

				/*------------------------------------------------------*\
				 | Read areas											|
				\*------------------------------------------------------*/

				category[i]->area = (struct area **) malloc (category[i]->area_cnt * sizeof (struct area *));
				if (category[i]->area) read_areas (in,category[i]->area,category[i]->area_cnt);
				else {
					printf ("Error: could not allocate area array for category %s\n",category[i]->name);
					exit (1);
					}

				/*------------------------------------------------------*\
				 | Read lines											|
				\*------------------------------------------------------*/

				category[i]->line = (struct line **) malloc (category[i]->line_cnt * sizeof (struct line *));
				if (category[i]->line) read_lines (in,category[i]->line,category[i]->line_cnt);
				else {
					printf ("Error: could not allocate line array for category %s\n",category[i]->name);
					exit (1);
					}

				/*------------------------------------------------------*\
				 | Convert to latitude and longitude					|
				\*------------------------------------------------------*/

				for (j=0; j < category[i]->node_cnt; j++) {
					px = &category[i]->node[j]->x;
					py = &category[i]->node[j]->y;
					axy_ll_e (*px,*py,py,px);
					}

				for (j=0; j < category[i]->area_cnt; j++) {
					px = &category[i]->area[j]->x;
					py = &category[i]->area[j]->y;
					axy_ll_e (*px,*py,py,px);
					}

				for (j=0; j < category[i]->line_cnt; j++) {
					for (k=0; k < category[i]->line[j]->coord_cnt; k++) {
						px = &category[i]->line[j]->coord[k].x;
						py = &category[i]->line[j]->coord[k].y;
						axy_ll_e (*px,*py,py,px);
						}
					}

				/*------------------------------------------------------*\
				 | Find state boundaries								|
				\*------------------------------------------------------*/

				find_state_boundaries (category[i]);

				/*------------------------------------------------------*\
				 | Show data if desired									|
				\*------------------------------------------------------*/

				#ifdef DEBUG
					for (j=0; j < category[i]->node_cnt; j++) {
						m = category[i]->node[j]->attr_cnt;
						if (m > 0) {
							for (k=0; k < m; k++)
								printf ("Node %3d: %s\n",category[i]->node[j]->id,
									description (category[i]->node[j]->attr[k]));
							}
						}

					for (j=0; j < category[i]->area_cnt; j++) {
						m = category[i]->area[j]->attr_cnt;
						if (m > 0) {
							for (k=0; k < m; k++)
								printf ("Area %3d: %s\n",category[i]->area[j]->id,
									description (category[i]->area[j]->attr[k]));
							}
						}

					for (j=0; j < category[i]->line_cnt; j++) {
						m = category[i]->line[j]->attr_cnt;
						if (m > 0) {
							for (k=0; k < m; k++)
								printf ("Line %3d: %s\n",category[i]->line[j]->id,
									description (category[i]->line[j]->attr[k]));
							}
						}
				#endif

				}

			/*----------------------------------------------------------*\
			 | Close input file											|
			\*----------------------------------------------------------*/

			fclose (in);

			/*----------------------------------------------------------*\
			 | Summarize data and determine what should be output		|
			\*----------------------------------------------------------*/

			printf ("File %s\n",input_file);
			printf ("%s\n",header->banner);
			printf ("%d categor%s\n",header->category_cnt,((header->category_cnt > 1) ? "ies" : "y"));
			for (i=0; i < header->category_cnt; i++) {
				printf ("Category: %s\n",category[i]->name);
				validate_lgo (category[i]);
				for (j=0, k=0; j < category[i]->node_cnt; j++) k += (category[i]->node[j]->attr_cnt > 0);
				printf ("  %4d nodes, %4d with attributes\n",category[i]->node_cnt,k);
				for (j=0, k=0; j < category[i]->area_cnt; j++) k += (category[i]->area[j]->attr_cnt > 0);
				printf ("  %4d areas, %4d with attributes\n",category[i]->area_cnt,k);
				for (j=0, k=0; j < category[i]->line_cnt; j++) k += (category[i]->line[j]->attr_cnt > 0);
				printf ("  %4d lines, %4d with attributes\n",category[i]->line_cnt,k);

				/*----------------------------------------------------------*\
				 | Find out what the user wants to write out.				|
				\*----------------------------------------------------------*/

				do {
					printf ("  Write out (A) all nodes,\n");
					printf ("            (S) only nodes matching search,\n");
					printf ("            (O) only nodes with attributes,\n");
					printf ("         or (N) no nodes?: ");
					gets (string);
					for (s=string; *s && isspace (*s); s++);
					switch (toupper (*s)) {
						case 'A':
							output.nodes = ALL;
							break;
						case 'S':
							output.nodes = SEARCH;
							printf ("Enter search criterion:\n");
							search.major = search.minor = 0;
							while (!search.major) {
								printf ("Major attribute code: ");
								gets (string);
								search.major = (int) strtol (string,0,0);
								}
							while (!search.minor) {
								printf ("Minor attribute code: ");
								gets (string);
								search.minor = (int) strtol (string,0,0);
								}
							break;
						case 'O':
							output.nodes = ATTR_ONLY;
							break;
						case 'N':
							output.nodes = NONE;
							break;
						default:
							output.nodes = 0;
							break;
						}
					} while (output.nodes == 0);

				do {
					printf ("  Write out (A) all areas,\n");
					printf ("            (S) only areas matching search,\n");
					printf ("            (O) only areas with attributes,\n");
					printf ("         or (N) no areas?: ");
					gets (string);
					for (s=string; *s && isspace (*s); s++);
					switch (toupper (*s)) {
						case 'A':
							output.areas = ALL;
							break;
						case 'S':
							output.areas = SEARCH;
							printf ("Enter search criterion:\n");
							search.major = search.minor = 0;
							while (!search.major) {
								printf ("Major attribute code: ");
								gets (string);
								search.major = (int) strtol (string,0,0);
								}
							while (!search.minor) {
								printf ("Minor attribute code: ");
								gets (string);
								search.minor = (int) strtol (string,0,0);
								}
							break;
						case 'O':
							output.areas = ATTR_ONLY;
							break;
						case 'N':
							output.areas = NONE;
							break;
						default:
							output.areas = 0;
							break;
						}
					} while (output.areas == 0);

				do {
					printf ("  Write out (A) all lines,\n");
					printf ("            (S) only lines matching search,\n");
					printf ("            (O) only lines with attributes,\n");
					printf ("         or (N) no lines?: ");
					gets (string);
					for (s=string; *s && isspace (*s); s++);
					switch (toupper (*s)) {
						case 'A':
							output.lines = ALL;
							break;
						case 'S':
							output.lines = SEARCH;
							printf ("Enter search criterion:\n");
							search.major = search.minor = 0;
							while (!search.major) {
								printf ("Major attribute code: ");
								gets (string);
								search.major = (int) strtol (string,0,0);
								}
							while (!search.minor) {
								printf ("Minor attribute code: ");
								gets (string);
								search.minor = (int) strtol (string,0,0);
								}
							break;
						case 'O':
							output.lines = ATTR_ONLY;
							break;
						case 'N':
							output.lines = NONE;
							break;
						default:
							output.lines = 0;
							break;
						}
					} while (output.lines == 0);

				/*------------------------------------------------------*\
				 | Write desired output									|
				\*------------------------------------------------------*/

				if (output.areas != NONE ||
					output.lines != NONE ||
					output.nodes != NONE) {
					if (!(f_geo = fopen (geo_file,"w"))) {
						printf ("Error: could not open output file %s\n",geo_file);
						exit (1);
						}
					if (!(f_attr = fopen (attr_file,"w"))) {
						printf ("Error: could not open output file %s\n",attr_file);
						exit (1);
						}
					printf ("Output files: %s %s\n",geo_file,attr_file);
					}
				else {
					f_geo = NULL;
					f_attr = NULL;
					}

				switch (target) {
					case MapInfo:
						mif_write_header (f_geo);
						mif_write_nodes (f_geo,f_attr,category[i],output);
						mif_write_areas (f_geo,f_attr,category[i],output);
						mif_write_lines (f_geo,f_attr,category[i],output);
						break;
					case AtlasPro:
					case Tactician:
						bna_write_header (f_attr);
						bna_write_nodes (f_geo,f_attr,category[i],output);
						bna_write_areas (f_geo,f_attr,category[i],output);
						bna_write_lines (f_geo,f_attr,category[i],output);
						break;
					default:
						break;
					}
				}

			/*----------------------------------------------------------*\
			 | Close output files										|
			\*----------------------------------------------------------*/

			fclose (f_geo);
			fclose (f_attr);
			}
		else {
			printf ("Error: could not open input file %s\n",input_file);
			exit (1);
			}
		}
	else {
		printf ("Usage: %s [-{a|m|T}] [-v] infile\n",argv[0]);
		exit (1);
		}
	}

/*----------------------------------------------------------------------*\
 | Read the header of the DLG optional-format file						|
\*----------------------------------------------------------------------*/

static void read_lgo_header (FILE *in, struct lgo_header *h) {
	int i;
	char string[MAXLEN];
	char number[25];
	char *s;
	int deg,min,sec;
	long d;

	/*-- Line 1 --------------------------------------------------------*/

	fgets (string,MAXLEN,in);
	if (s = strrchr (string,'\n')) *s = 0;
	strcpy (h->banner,string);

	/*-- Line 2 --------------------------------------------------------*/

	fgets (string,MAXLEN,in);
	*h->mapname = 0;
	strncpy0 (h->mapname,string,40);

	strncpy0 (number,string+41,10);
	h->source_year = (int) strtol (number,0,0);

	strncpy0 (number,string+52,8);
	h->source_scale = strtoul (number,0,0);

	/*-- Line 3 --------------------------------------------------------*/

	fgets (string,MAXLEN,in);

	/*-- Line 4 --------------------------------------------------------*/

	fgets (string,MAXLEN,in);

	strncpy0 (number,string,6);
	h->DLG_level = (int) strtol (number,0,0);

	strncpy0 (number,string+6,6);
	h->ref_system = (int) strtol (number,0,0);

	strncpy0 (number,string+12,6);
	h->zone = (int) strtol (number,0,0);

	strncpy0 (number,string+18,6);
	h->units = (int) strtol (number,0,0);

	strncpy0 (number,string+24,18);
	d_to_e (number);
	h->resolution = strtod (number,0);

	strncpy0 (number,string+42,6);
	h->trans_param_cnt = (int) strtol (number,0,0);

	strncpy0 (number,string+48,6);
	h->acc_rec_cnt = (int) strtol (number,0,0);

	strncpy0 (number,string+54,6);
	h->ctrl_pt_cnt = (int) strtol (number,0,0);

	strncpy0 (number,string+60,6);
	h->category_cnt = (int) strtol (number,0,0);

	/*-- Lines 5 .. 9 --------------------------------------------------*/

	fgets (string,MAXLEN,in);

	strncpy0 (number,string,   24);
	d_to_e (number);
	h->a = strtod (number,0);

	strncpy0 (number,string+24,24);
	d_to_e (number);
	h->e2 = strtod (number,0);

	strncpy0 (number,string+48,24);
	d_to_e (number);
	d = (long) strtod (number,0);
	sec = d % 1000;
	min = (d/1000) % 1000;
	deg = d/1000000;
	h->lat1 = ((double)deg) + ((double)min)/60 + ((double)sec)/3600;

	fgets (string,MAXLEN,in);

	strncpy0 (number,string   ,24);
	d_to_e (number);
	d = (long) strtod (number,0);
	sec = d % 1000;
	min = (d/1000) % 1000;
	deg = d/1000000;
	h->lat2 = ((double)deg) + ((double)min)/60 + ((double)sec)/3600;

	strncpy0 (number,string+24,24);
	d_to_e (number);
	d = (long) strtod (number,0);
	sec = d % 1000;
	min = (d/1000) % 1000;
	deg = d/1000000;
	h->lon0 = ((double)deg) + ((double)min)/60 + ((double)sec)/3600;

	strncpy0 (number,string+48,24);
	d_to_e (number);
	d = (long) strtod (number,0);
	sec = d % 1000;
	min = (d/1000) % 1000;
	deg = d/1000000;
	h->lat0 = ((double)deg) + ((double)min)/60 + ((double)sec)/3600;

	fgets (string,MAXLEN,in);

	strncpy0 (number,string   ,24);
	d_to_e (number);
	h->easting = strtod (number,0);

	strncpy0 (number,string+24,24);
	d_to_e (number);
	h->northing = strtod (number,0);

	fgets (string,MAXLEN,in);
	fgets (string,MAXLEN,in);

	/*-- Line 10 -------------------------------------------------------*/

	fgets (string,MAXLEN,in);

	strncpy0 (number,string,   18);
	d_to_e (number);
	h->a1 = strtod (number,0);

	strncpy0 (number,string+18,18);
	d_to_e (number);
	h->a2 = strtod (number,0);

	strncpy0 (number,string+36,18);
	d_to_e (number);
	h->a3 = strtod (number,0);

	strncpy0 (number,string+54,18);
	d_to_e (number);
	h->a4 = strtod (number,0);
	}

/*----------------------------------------------------------------------*\
 | Read the control points of the DLG optional-format file				|
\*----------------------------------------------------------------------*/

static void read_ctrl_pts (FILE *in, struct ctrl_pt *p[], int n) {
	int i;
	char string[MAXLEN];
	char number[25];

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

		p[i] = (struct ctrl_pt *) malloc (sizeof (struct ctrl_pt));
		if (!p[i]) {
			printf ("Error: Could not allocate ctrl_pt structure\n");
			exit (1);
			}

		fgets (string,MAXLEN,in);
		*p[i]->corner = 0;
		strncpy0 (p[i]->corner,string,2);

		strncpy0 (number,string+6, 12);
		d_to_e (number);
		p[i]->latitude = strtod (number,0);

		strncpy0 (number,string+18,12);
		d_to_e (number);
		p[i]->longitude = strtod (number,0);

		strncpy0 (number,string+36,12);
		d_to_e (number);
		p[i]->x = strtod (number,0);

		strncpy0 (number,string+48,12);
		d_to_e (number);
		p[i]->y = strtod (number,0);
		}
	}

/*----------------------------------------------------------------------*\
 | Read the data category identification records						|
\*----------------------------------------------------------------------*/

static void read_categories (FILE *in, struct category *c[], int n) {
	int i;
	char string[MAXLEN];
	char number[25];

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

		c[i] = (struct category *) malloc (sizeof (struct category));
		if (!c[i]) {
			printf ("Error: Could not allocate category structure\n");
			exit (1);
			}

		fgets (string,MAXLEN,in);
		*c[i]->name = 0;
		strncpy0 (c[i]->name,string,20);

		strncpy0 (number,string+20,4);
		c[i]->attr_format = (int) strtol (number,0,0);

		strncpy0 (number,string+24,6);
		c[i]->node_ref_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+30,6);
		c[i]->node_cnt = (int) strtol (number,0,0);

		c[i]->link_node_area = ((string[37] == '1') ? 1 : 0);
		c[i]->link_node_line = ((string[38] == '1') ? 1 : 0);

		strncpy0 (number,string+40,6);
		c[i]->area_ref_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+46,6);
		c[i]->area_cnt = (int) strtol (number,0,0);

		c[i]->link_area_node = ((string[53] == '1') ? 1 : 0);
		c[i]->link_area_line = ((string[54] == '1') ? 1 : 0);
		c[i]->area_coord_lists = ((string[55] == '1') ? 1 : 0);

		strncpy0 (number,string+56,6);
		c[i]->line_ref_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+62,6);
		c[i]->line_cnt = (int) strtol (number,0,0);

		c[i]->line_coord_lists = ((string[71] == '1') ? 1 : 0);
		}
	}

/*----------------------------------------------------------------------*\
 | Read the node records												|
\*----------------------------------------------------------------------*/

static void read_nodes (FILE *in, struct node *p[], int n) {
	int i,j,k,L,m,mm;
	char string[MAXLEN];
	char number[25];
	char *s;

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

		p[i] = (struct node *) malloc (sizeof (struct node));
		if (!p[i]) {
			printf ("Error: Could not allocate node structure\n");
			exit (1);
			}

		fgets (string,MAXLEN,in);
		if (*string != 'N') {
			if (s = strrchr (string,'\n')) *s = 0;
			printf ("Error: Expected node record, got \"%s\"\n",string);
			exit (1);
			}

		strncpy0 (number,string+1,5);
		p[i]->id = (int) strtol (number,0,0);
		if (p[i]->id != i+1) printf ("Warning: Node %d has id %d\n",i,p[i]->id);

		strncpy0 (number,string+6,12);
		d_to_e (number);
		p[i]->x = strtod (number,0);

		strncpy0 (number,string+18,12);
		d_to_e (number);
		p[i]->y = strtod (number,0);

		strncpy0 (number,string+30,6);
		p[i]->area_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+36,6);
		p[i]->line_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+42,6);
		p[i]->area_coord_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+48,6);
		p[i]->attr_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+54,6);
		p[i]->text_cnt = (int) strtol (number,0,0);

		#ifdef DEBUG
			printf ("node %3d (%3d) x: %15.4lf y:%15.4lf lines: %3d attrs: %2d\n",
				i,
				p[i]->id,
				p[i]->x,
				p[i]->y,
				p[i]->line_list_cnt,
				p[i]->attr_cnt);
		#endif

		m = p[i]->line_list_cnt;
		if (m > 0) {
			p[i]->line = (int *) malloc (m * sizeof (int));
			if (p[i]->line) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/12; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,12);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*6,6);
						p[i]->line[k++] = (int) strtol (number,0,0);
						}
					}
				}
			else {
				printf ("Error: could not allocate line list for node %d\n",i);
				exit (1);
				}
			}

		m = p[i]->attr_cnt;
		if (m > 0) {
			p[i]->attr = (struct attr *) malloc (m * sizeof (struct attr));
			if (p[i]->attr) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/6; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,6);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*12,6);
						p[i]->attr[k].major = (int) strtol (number,0,0);
						strncpy0 (number,string+L*12+6,6);
						p[i]->attr[k].minor = (int) strtol (number,0,0);
						k++;
						}
					}
				}
			else {
				printf ("Error: could not allocate line list for node %d\n",i);
				exit (1);
				}
			}
		}
	}

/*----------------------------------------------------------------------*\
 | Read the area records												|
\*----------------------------------------------------------------------*/

static void read_areas (FILE *in, struct area *p[], int n) {
	int i,j,k,L,m,mm;
	char string[MAXLEN];
	char number[25];
	char *s;

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

		p[i] = (struct area *) malloc (sizeof (struct area));
		if (!p[i]) {
			printf ("Error: Could not allocate area structure %d\n",i);
			exit (1);
			}

		fgets (string,MAXLEN,in);
		if (*string != 'A') {
			if (s = strrchr (string,'\n')) *s = 0;
			printf ("Error: Expected area record, got \"%s\"\n",string);
			exit (1);
			}

		strncpy0 (number,string+1,5);
		p[i]->id = (int) strtol (number,0,0);
		if (p[i]->id != i+1) printf ("Warning: Area %d has id %d\n",i,p[i]->id);

		strncpy0 (number,string+6,12);
		d_to_e (number);
		p[i]->x = strtod (number,0);

		strncpy0 (number,string+18,12);
		d_to_e (number);
		p[i]->y = strtod (number,0);

		strncpy0 (number,string+30,6);
		p[i]->node_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+36,6);
		p[i]->line_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+42,6);
		p[i]->area_coord_list_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+48,6);
		p[i]->attr_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+54,6);
		p[i]->text_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+60,6);
		p[i]->island_cnt = (int) strtol (number,0,0);

		#ifdef DEBUG
			printf ("area %3d (%3d) x: %15.4lf y:%15.4lf lines: %3d attrs: %2d\n",
				i,
				p[i]->id,
				p[i]->x,
				p[i]->y,
				p[i]->line_list_cnt,
				p[i]->attr_cnt);
		#endif

		m = p[i]->line_list_cnt;
		if (m > 0) {
			p[i]->line = (int *) malloc (m * sizeof (int));
			if (p[i]->line) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/12; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,12);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*6,6);
						p[i]->line[k++] = (int) strtol (number,0,0);
						}
					}
				}
			else {
				printf ("Error: could not allocate line list for node %d\n",i);
				exit (1);
				}
			}

		m = p[i]->attr_cnt;
		if (m > 0) {
			p[i]->attr = (struct attr *) malloc (m * sizeof (struct attr));
			if (p[i]->attr) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/6; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,6);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*12,6);
						p[i]->attr[k].major = (int) strtol (number,0,0);
						strncpy0 (number,string+L*12+6,6);
						p[i]->attr[k].minor = (int) strtol (number,0,0);
						k++;
						}
					}
				}
			else {
				printf ("Error: could not allocate attr list for node %d\n",i);
				exit (1);
				}
			}
		}
	}

/*----------------------------------------------------------------------*\
 | Read the line records												|
\*----------------------------------------------------------------------*/

static void read_lines (FILE *in, struct line *p[], int n) {
	int i,j,k,L,m,mm;
	char string[MAXLEN];
	char number[25];
	char *s;

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

		p[i] = (struct line *) malloc (sizeof (struct line));
		if (!p[i]) {
			printf ("Error: Could not allocate line structure\n");
			exit (1);
			}

		fgets (string,MAXLEN,in);
		if (*string != 'L') {
			if (s = strrchr (string,'\n')) *s = 0;
			printf ("Error: Expected line record, got \"%s\"\n",string);
			exit (1);
			}

		strncpy0 (number,string+1, 5);
		p[i]->id = (int) strtol (number,0,0);
		if (p[i]->id != i+1) printf ("Warning: Line %d has id %d\n",i,p[i]->id);

		strncpy0 (number,string+6, 6);
		p[i]->beg_node = (int) strtol (number,0,0);

		strncpy0 (number,string+12,6);
		p[i]->end_node = (int) strtol (number,0,0);

		strncpy0 (number,string+18,6);
		p[i]->left_area = (int) strtol (number,0,0);

		strncpy0 (number,string+24,6);
		p[i]->right_area = (int) strtol (number,0,0);

		strncpy0 (number,string+42,6);
		p[i]->coord_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+48,6);
		p[i]->attr_cnt = (int) strtol (number,0,0);

		strncpy0 (number,string+54,6);
		p[i]->text_cnt = (int) strtol (number,0,0);

		#ifdef DEBUG
			printf ("line %3d (%3d) beg: %3d end:%3d coords: %3d attrs: %2d\n",
				i,
				p[i]->id,
				p[i]->beg_node,
				p[i]->end_node,
				p[i]->coord_cnt,
				p[i]->attr_cnt);
		#endif

		m = p[i]->coord_cnt;
		if (m > 0) {
			p[i]->coord = (struct point *) malloc (m * sizeof (struct point));
			if (p[i]->coord) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/3; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,3);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*24,12);
						d_to_e (number);
						p[i]->coord[k].x = strtod (number,0);
						strncpy0 (number,string+L*24+12,12);
						d_to_e (number);
						p[i]->coord[k].y = strtod (number,0);
						k++;
						}
					}
				}
			else {
				printf ("Error: could not allocate line list for node %d\n",i);
				exit (1);
				}
			}

		m = p[i]->attr_cnt;
		if (m > 0) {
			p[i]->attr = (struct attr *) malloc (m * sizeof (struct attr));
			if (p[i]->attr) {
				k = 0;
				for (j=0; j < 1 + (m - 1)/6; j++) {
					fgets (string,MAXLEN,in);
					mm = min (m-k,6);
					for (L=0; L < mm; L++) {
						strncpy0 (number,string+L*12,6);
						p[i]->attr[k].major = (int) strtol (number,0,0);
						strncpy0 (number,string+L*12+6,6);
						p[i]->attr[k].minor = (int) strtol (number,0,0);
						k++;
						}
					}
				}
			else {
				printf ("Error: could not allocate line list for node %d\n",i);
				exit (1);
				}
			}
		}
	}

/*----------------------------------------------------------------------*\
 | Procedures to find a given node, area, or line by its id code.		|
\*----------------------------------------------------------------------*/

struct node *find_node (struct category *c, int id) {
	int i;
	if (id - 1 < c->node_cnt)
		if (c->node[id-1]->id == id) return (c->node[id-1]);
	for (i=0; i < c->node_cnt; i++)
		if (c->node[i]->id == id) return (c->node[i]);
	return (NULL);
	}

struct area *find_area (struct category *c, int id) {
	int i;
	if (id - 1 < c->area_cnt)
		if (c->area[id-1]->id == id) return (c->area[id-1]);
	for (i=0; i < c->area_cnt; i++)
		if (c->area[i]->id == id) return (c->area[i]);
	return (NULL);
	}

struct line *find_line (struct category *c, int id) {
	int i;
	if (id - 1 < c->line_cnt)
		if (c->line[id-1]->id == id) return (c->line[id-1]);
	for (i=0; i < c->line_cnt; i++)
		if (c->line[i]->id == id) return (c->line[i]);
	return (NULL);
	}

/*----------------------------------------------------------------------*\
 | Find state boundaries												|
\*----------------------------------------------------------------------*/

static void find_state_boundaries (struct category *c) {
	int i,j,k,m,n,nc;
	int jj,kk;
	int state_cnt;
	struct state **state;
	struct line *L;
	struct area *a;
	struct area **aa;
	int left_state,right_state;
	int beg,end;
	struct state_line *s,*t;
	char output_file[MAXLEN];
	FILE *geo, *attr;
	int forward, error;
	int implicit_state_line,explicit_state_line;

	int code,section;
	char *name, *abbr, *q;
	char geo_file[MAXLEN],attr_file[MAXLEN];
	char primary_id[MAXLEN];

	state_cnt = 0;
	state = (struct state **) malloc (16 * sizeof (struct state *));
	if (!state) {
		printf ("Error: could not allocate state boundary array\n");
		exit (1);
		}

	for (i=0; i < c->line_cnt; i++) {

		/*--------------------------------------------------------------*\
		 | For each line in this category								|
		\*--------------------------------------------------------------*/

		L = c->line[i];

		/*--------------------------------------------------------------*\
		 | Find the state code of the area to the left of the line.		|
		\*--------------------------------------------------------------*/

		left_state = 0;

		if (a = find_area (c,L->left_area)) {
			if (a->id != L->left_area)
				for (a=NULL, k=0; k < c->area_cnt; k++)
					if (c->area[k]->id == L->left_area) {
						a = c->area[k];
						break;
						}

			/*----------------------------------------------------------*\
			 | If the area could be found, scan its attributes looking	|
			 | for state codes.											|
			\*----------------------------------------------------------*/

			if (a->attr_cnt) {
				for (k=0; k < a->attr_cnt; k++) {
					if (a->attr[k].major == 91) {
						left_state = a->attr[k].minor;
						break;
						}
					}
				}
			}

		/*--------------------------------------------------------------*\
		 | Find the state code of the area to the right of the line.	|
		\*--------------------------------------------------------------*/

		right_state = 0;

		if (a = find_area(c,L->right_area)) {
			if (a->id != L->right_area)
				for (a=NULL, k=0; k < c->area_cnt; k++)
					if (c->area[k]->id == L->right_area) {
						a = c->area[k];
						break;
						}

			/*----------------------------------------------------------*\
			 | If the area could be found, scan its attributes looking	|
			 | for state codes.											|
			\*----------------------------------------------------------*/

			if (a->attr_cnt) {
				for (k=0; k < a->attr_cnt; k++) {
					if (a->attr[k].major == 91) {
						right_state = a->attr[k].minor;
						break;
						}
					}
				}
			}

		/*--------------------------------------------------------------*\
		 | Compare left area with right area; if they are not the same,	|
		 | add this line to the state line list of each state.			|
		\*--------------------------------------------------------------*/

		if (left_state != right_state) {

			/*----------------------------------------------------------*\
			 | This line is a state boundary. Check its attributes.		|
			 |															|
			 | Issue a warning if the line is marked neither as a state	|
			 | boundary nor as a national boundary.  This indicates a	|
			 | probable DLG attribute error in any states other than	|
			 | Texas, California, Montana, and Alaska, where it might	|
			 | instead be part of the section boundary.					|
			 |															|
			 | If verbose is true, issue a warning if the line is marked|
			 | as a national boundary but is not also marked as a state	|
			 | boundary.												|
			\*----------------------------------------------------------*/

			explicit_state_line = 0;
			implicit_state_line = 0;
			if (L->attr_cnt)
				for (k=0; k < L->attr_cnt; k++)
					if (L->attr[k].major == 290) {
						if (L->attr[k].minor == 6005 ||
							L->attr[k].minor == 6006) explicit_state_line = 1;
						if (L->attr[k].minor >= 6000 &&
							L->attr[k].minor <= 6002) implicit_state_line = 1;
						}
			if (!explicit_state_line && !implicit_state_line)
				printf ("Warning: line %d is neither explicitly nor implicitly marked as a state line\n",L->id);
			else
				if (!explicit_state_line)
					if (verbose) printf ("Warning: line %d is not explicitly marked as a state line\n",L->id);

			if (left_state) {

				/*------------------------------------------------------*\
				 | Search through the state array to find the state		|
				 | whose FIPS_code matches left_state.					|
				\*------------------------------------------------------*/

				for (k=0; k < state_cnt; k++)
					if (left_state == state[k]->FIPS_code) break;

				if (k == state_cnt) {

					/*--------------------------------------------------*\
					 | If you fell through the loop above, this state is|
					 | not in the list, and should be added.			|
					\*--------------------------------------------------*/

					state[k] = (struct state *) malloc (sizeof (struct state));
					if (!state[k]) {
						printf ("Error: could not allocate state %d\n",k);
						exit (1);
						}
					state[k]->FIPS_code = left_state;
					state[k]->line_cnt = 1;
					state[k]->line = (struct state_line *) malloc (c->line_cnt * sizeof(struct state_line));
					state_cnt++;
					state[k]->line[0].index = i;
					state[k]->line[0].id = c->line[i]->id;
					}
				else {

					/*--------------------------------------------------*\
					 | This state is in the list, and the line should be|
					 | added to it.										|
					\*--------------------------------------------------*/

					state[k]->line[state[k]->line_cnt].index = i;
					state[k]->line[state[k]->line_cnt].id = c->line[i]->id;
					state[k]->line_cnt++;
					}
				}
			if (right_state) {

				/*------------------------------------------------------*\
				 | Search through the state array to find the state		|
				 | whose FIPS_code matches right_state.					|
				\*------------------------------------------------------*/

				for (k=0; k < state_cnt; k++)
					if (right_state == state[k]->FIPS_code) break;

				if (k == state_cnt) {

					/*--------------------------------------------------*\
					 | If you fell through the loop above, this state	|
					 | is not in the list, and should be added.			|
					\*--------------------------------------------------*/

					state[k] = (struct state *) malloc (sizeof (struct state));
					if (!state[k]) {
						printf ("Error: could not allocate state %d\n",k);
						exit (1);
						}
					state[k]->FIPS_code = right_state;
					state[k]->line_cnt = 1;
					state[k]->line = (struct state_line *) malloc (c->line_cnt * sizeof(struct state_line));
					state_cnt++;
					state[k]->line[0].index = i;
					state[k]->line[0].id = c->line[i]->id;
					}
				else {

					/*--------------------------------------------------*\
					 | This state is in the list, and the line should be|
					 | added to it.										|
					\*--------------------------------------------------*/

					state[k]->line[state[k]->line_cnt].index = i;
					state[k]->line[state[k]->line_cnt].id = c->line[i]->id;
					state[k]->line_cnt++;
					}
				}
			}
		}

	/*------------------------------------------------------------------*\
	 | Reallocate the area pointer array, to add the states.			|
	\*------------------------------------------------------------------*/

	if (!(aa = (struct area **) realloc (c->area,(c->area_cnt + state_cnt) * sizeof (struct area *)))) {
		printf ("Error: could not reallocate area pointer array to hold states.\n");
		printf ("  Consequently, the states cannot be output.\n");
		return;
		}
	else
		c->area = aa;

	/*------------------------------------------------------------------*\
	 | Determine connections among lines								|
	\*------------------------------------------------------------------*/

	for (i=0; i < state_cnt; i++) {
		for (j=0; j < state[i]->line_cnt; j++) {
			s = &state[i]->line[j];
			s->bb = s->be = s->eb = s->ee = -1;
			s->links = 0;
			L = c->line[s->index];
			beg = L->beg_node;
			end = L->end_node;
			for (k=0; k < state[i]->line_cnt; k++) {
				if (k != j) {
					L = c->line [state[i]->line[k].index];
					if (L->beg_node == beg) {
						s->bb = k;
						s->links++;
						}
					if (L->end_node == beg) {
						s->be = k;
						s->links++;
						}
					if (L->beg_node == end) {
						s->eb = k;
						s->links++;
						}
					if (L->end_node == end) {
						s->ee = k;
						s->links++;
						}
					}
				}
			}

		/*--------------------------------------------------------------*\
		 | Check the links; all lines should be 2-connected. If this is	|
		 | not true, dump a complete line list to stdout and try to go	|
		 | on to the other states in the file without writing output.	|
		\*--------------------------------------------------------------*/

		error = 0;
		for (j=0; j < state[i]->line_cnt; j++) {
			s = &state[i]->line[j];
			if (s->links != 2) {
				error = 1;
				break;
				}
			}

		if (error) {
			printf ("Warning: State %d has a serious topological problem, and cannot be assembled.\n",state[i]->FIPS_code);
			printf ("State %2d has FIPS_code %2d, and %3d lines:\n",i,state[i]->FIPS_code,state[i]->line_cnt);
			for (j=0; j < state[i]->line_cnt; j++) {
				s = &state[i]->line[j];
				if (s->links != 2) printf ("->");
				else printf ("  ");
				printf ("line %3d (index %3d) has bb %3d, be %3d, eb %3d, ee %3d\n",j,s->index,s->bb,s->be,s->eb,s->ee);
				}
			continue;
			}

		/*--------------------------------------------------------------*\
		 | Allocate space for the boundary array						|
		\*--------------------------------------------------------------*/

		state[i]->boundary = (int *) malloc (state[i]->line_cnt * sizeof (int));
		if (!state[i]->boundary) {
			printf ("Error: could not allocate boundary array\n");
			exit (1);
			}

		/*--------------------------------------------------------------*\
		 | Rearrange the list of lines so that connected lines are		|
		 | consecutive.													|
		 |																|
		 | forward = 1													|
		 |		B--E:b--e	add s->eb positive, keep forward true		|
		 |		B--E:e--b	add s->ee negative, make forward false		|
		 |																|
		 | forward = 0													|
		 |		E--B:b--e	add s->bb positive, make forward true		|
		 |		E--B:e--b	add s->be negative, keep forward false		|
		 |																|
		\*--------------------------------------------------------------*/

		k = 0;
		s = &state[i]->line[k];
		state[i]->boundary[k] = k;
		k++;

		forward = 1;
		while (k < state[i]->line_cnt)
			if (forward)
				if (s->ee < 0) {
					state[i]->boundary[k++] = s->eb;
					s = &state[i]->line[s->eb];
					}
				else {
					state[i]->boundary[k++] = -s->ee;
					s = &state[i]->line[s->ee];
					forward = 0;
					}
			else
				if (s->bb < 0) {
					state[i]->boundary[k++] = -s->be;
					s = &state[i]->line[s->be];
					}
				else {
					state[i]->boundary[k++] = s->bb;
					s = &state[i]->line[s->bb];
					forward = 1;
					}

		/*--------------------------------------------------------------*\
		 | Replace the boundary array elements with the corresponding	|
		 | line id numbers.												|
		\*--------------------------------------------------------------*/

		for (j=0; j < state[i]->line_cnt; j++) {
			k = state[i]->boundary[j];
			if (k < 0) state[i]->boundary[j] = -state[i]->line[-k].id;
			else       state[i]->boundary[j] =  state[i]->line[ k].id;
			}

		/*--------------------------------------------------------------*\
		 | Check for duplicate line entries.							|
		\*--------------------------------------------------------------*/

		for (j=1; j < state[i]->line_cnt; j++)
			if (state[i]->boundary[j] == state[i]->boundary[0]) {
				printf ("Warning: state %d has improper topology.\n",state[i]->FIPS_code);
				for (k=0; k < state[i]->line_cnt; k++) {
					m = state[i]->line[k].id;
					for (n=0; n < j; n++)
						if (m == abs (state[i]->boundary[n])) break;
					if (n == j)
						if (L = find_line (c,m))
							printf ("  - Line %d has left area %d, right area %d\n",m,L->left_area,L->right_area);
					}
				printf ("  Attempting to generate a plausible boundary anyway ...\n");
				state[i]->line_cnt = j;
				break;
				}

		/*--------------------------------------------------------------*\
		 | Determine whether the boundary is clockwise. If it is not,	|
		 | reverse the sequence of the lines.							|
		\*--------------------------------------------------------------*/

		L = find_line (c,state[i]->boundary[0]);;
		a = find_area (c,L->left_area);
		left_state = 0;
		if (a->attr_cnt) {
			for (k=0; k < a->attr_cnt; k++) {
				if (a->attr[k].major == 91) {
					left_state = a->attr[k].minor;
					break;
					}
				}
			}
		if (left_state == state[i]->FIPS_code) {
			/* reverse the boundary, then negate the id's */
			for (j=0; j < state[i]->line_cnt/2; j++) {
				k = state[i]->line_cnt - j - 1;
				m                     = state[i]->boundary[j];
				state[i]->boundary[j] = state[i]->boundary[k];
				state[i]->boundary[k] = m;
				}
			for (j=0; j < state[i]->line_cnt; j++) state[i]->boundary[j] *= -1;
			}

		/*--------------------------------------------------------------*\
		 | Get state name and abbreviation, to be used in output.		|
		\*--------------------------------------------------------------*/

		code = state[i]->FIPS_code;
		name = state_name (code);
		abbr = state_abbr (code);

		/*--------------------------------------------------------------*\
		 | Add the state boundary to the area list						|
		\*--------------------------------------------------------------*/

		if (!(c->area[c->area_cnt] = (struct area *) malloc (sizeof (struct area)))) {
			printf ("Error: could not allocate space for area %d (state %d)\n",c->area_cnt+i+1,state[i]->FIPS_code);
			return;
			}

		a = c->area[c->area_cnt];
		a->id = c->area_cnt+1;
		a->x = 0.0;
		a->y = 0.0;
		a->node_list_cnt = 0;
		a->line_list_cnt = state[i]->line_cnt;
		a->area_coord_list_cnt = 0;
		a->attr_cnt = 1;
		a->text_cnt = 0;
		a->island_cnt = 0;
		a->line = state[i]->boundary;
		a->attr = (struct attr *) malloc (sizeof (struct attr));
		if (a->attr) a->attr->major = 91;
		if (a->attr) a->attr->minor = code;
		c->area_cnt++;

		/*--------------------------------------------------------------*\
		 | Optionally show the line list in the format it would have	|
		 | in the DLG file.												|
		\*--------------------------------------------------------------*/

		if (verbose) {
			printf ("DLG file entry for %s:\n",name);
			printf ("A%5d",a->id);
			printf ("%12.2lf%12.2lf",a->x,a->y);
			printf ("%6d",a->node_list_cnt);
			printf ("%6d",a->line_list_cnt);
			printf ("%6d",a->area_coord_list_cnt);
			printf ("%6d",a->attr_cnt);
			printf ("%6d",a->text_cnt);
			printf ("%6d",a->island_cnt);
			printf ("\n");
			for (j=0; j < a->line_list_cnt; j++) {
				printf ("%6d",a->line[j]);
				if (((j+1) % 12) == 0) printf ("\n");
				}
			if (j % 12) printf ("\n");
			printf ("%6d%6d\n",a->attr[0].major,a->attr[0].minor);
			}
		}
	}

/*----------------------------------------------------------------------*\
 | Validation of DLG data												|
\*----------------------------------------------------------------------*/

static void validate_lgo (struct category *c) {
	int i,j,k,m;
	struct node *n;
	struct area *a;
	struct line *L;
	struct line *L0;
	struct line *L1;
	int j0;
	double x0,y0,x1,y1;

	/*------------------------------------------------------------------*\
	 | Node entry validation											|
	\*------------------------------------------------------------------*/

	printf ("Checking nodes\n");

	for (i=0; i < c->node_cnt; i++) {
		n = c->node[i];

		/*--------------------------------------------------------------*\
		 | Check to see that the lines listed for this node do in fact	|
		 | begin or end here.											|
		\*--------------------------------------------------------------*/

		for (j=0; j < n->line_list_cnt; j++) {

			m = abs(n->line[j]);
			if (L = find_line (c,m)) {
				if (n->line[j] > 0) {
					if (L->coord[0].x != n->x || L->coord[0].y != n->y)
						printf ("Warning: node %d lists starting line %d, which does not start here\n",n->id,L->id);
					}
				else
					if (L->coord[L->coord_cnt-1].x != n->x || L->coord[L->coord_cnt-1].y != n->y)
						printf ("Warning: node %d lists ending line %d, which does not end here\n",n->id,L->id);
				}
			else
				printf ("Warning: node %d refers to non-existent line %d\n",n->id,m);
			}
		}

	/*------------------------------------------------------------------*\
	 | Area entry validation											|
	\*------------------------------------------------------------------*/

	printf ("Checking areas\n");

	for (i=0; i < c->area_cnt; i++) {
		a = c->area[i];

		/*--------------------------------------------------------------*\
		 | Check to see that the lines listed for this area exist and	|
		 | are contiguous.												|
		\*--------------------------------------------------------------*/

		L0 = L1 = NULL;
		j0 = 0;
		for (j=0; j < a->line_list_cnt; j++) {
			if (m = abs(a->line[j])) {
				if (L = find_line (c,m)) {
					if (L1) {
						if (a->line[j-1] > 0) {
							/* previous line forward; look at its last point */
							x0 = L1->coord[L1->coord_cnt-1].x;
							y0 = L1->coord[L1->coord_cnt-1].y;
							}
						else {
							/* previous line reverse; look at its first point */
							x0 = L1->coord[0].x;
							y0 = L1->coord[0].y;
							}
						if (a->line[j] > 0) {
							/* current line forward; look at its first point */
							x1 = L->coord[0].x;
							y1 = L->coord[0].y;
							}
						else {
							/* current line reverse; look at its last point */
							x1 = L->coord[L->coord_cnt-1].x;
							y1 = L->coord[L->coord_cnt-1].y;
							}
						if (x0 != x1 || y0 != y1)
							printf ("Warning: area %d contains non-contiguous lines %d and %d\n",a->id,a->line[j-1],a->line[j]);
						L1 = L;
						}
					else {
						L0 = L;
						j0 = j;
						L1 = L;
						}
					}
				else
					printf ("Warning: area %d refers to non-existent line %d\n",a->id,m);
				}
			else {
				if (a->line[j-1] > 0) {
					/* previous line forward; look at its last point */
					x0 = L1->coord[L1->coord_cnt-1].x;
					y0 = L1->coord[L1->coord_cnt-1].y;
					}
				else {
					/* previous line reverse; look at its first point */
					x0 = L1->coord[0].x;
					y0 = L1->coord[0].y;
					}
				if (a->line[j0] > 0) {
					/* current line forward; look at its first point */
					x1 = L0->coord[0].x;
					y1 = L0->coord[0].y;
					}
				else {
					/* current line reverse; look at its last point */
					x1 = L0->coord[L0->coord_cnt-1].x;
					y1 = L0->coord[L0->coord_cnt-1].y;
					}
				if (x0 != x1 || y0 != y1)
					printf ("Warning: area %d contains non-contiguous lines %d and %d\n",a->id,a->line[j-1],a->line[j0]);
				L0 = L1 = NULL;
				}
			}
		}

	/*------------------------------------------------------------------*\
	 | Line entry validation											|
	\*------------------------------------------------------------------*/

	printf ("Checking lines\n");

	for (i=0; i < c->line_cnt; i++) {
		L = c->line[i];

		/*--------------------------------------------------------------*\
		 | Check to see that beginning and ending nodes coincide with	|
		 | the endpoints of the line.									|
		\*--------------------------------------------------------------*/

		if (n = find_node (c,L->beg_node)) {
			if (n->x != L->coord[0].x || n->y != L->coord[0].y)
				printf ("Warning: line %d lists beg_node %d (%lf,%lf), but starts at (%lf,%lf).\n",
					L->id,n->id,n->x,n->y,L->coord[0].x,L->coord[0].y);
			}
		else
			printf ("Warning: line %d lists non-existent beg_node %d\n",L->id,L->beg_node);

		if (n = find_node (c,L->end_node)) {
			k = L->coord_cnt - 1;
			if (n->x != L->coord[k].x || n->y != L->coord[k].y)
				printf ("Warning: line %d lists beg_node %d (%lf,%lf), but starts at (%lf,%lf).\n",
					L->id,n->id,n->x,n->y,L->coord[k].x,L->coord[k].y);
			}
		else
			printf ("Warning: line %d lists non-existent end_node %d\n",L->id,L->end_node);

		/*--------------------------------------------------------------*\
		 | Check to see that left and right areas contain this line.	|
		\*--------------------------------------------------------------*/

		if (a = find_area (c,L->left_area)) {
			for (j=0; j < a->line_list_cnt; j++) if (L->id == abs (a->line[j])) break;
			if (j == a->line_list_cnt)
				printf ("Warning: line %d lists left area %d, but is not in the area line list.\n",
					L->id,L->left_area);
			}
		else
			printf ("Warning: line %d lists non-existent left area %d\n",L->id,L->left_area);

		if (a = find_area (c,L->right_area)) {
			for (j=0; j < a->line_list_cnt; j++) if (L->id == abs (a->line[j])) break;
			if (j == a->line_list_cnt)
				printf ("Warning: line %d lists right area %d, but is not in the area line list.\n",
					L->id,L->right_area);
			}
		else
			printf ("Warning: line %d lists non-existent right area %d\n",L->id,L->left_area);
		}
	}

/*----------------------------------------------------------------------*\
 | Convert x,y to lat,lon												|
\*----------------------------------------------------------------------*/

#define DEG_RAD 0.01745329251

static double rho0,C,n,a,e,e2,q_max,Lon0;

static int set_const (double lat0, double lat1, double lat2, double lon0) {
	double m1,m2,q0,q1,q2;
	double b;

	a = 6378206.4;
	b = 6356583.8;

	lat0 *= DEG_RAD;
	lat1 *= DEG_RAD;
	lat2 *= DEG_RAD;
	lon0 *= DEG_RAD;
	Lon0 = lon0;

	e2 = 1.0 - (b*b)/(a*a);
	e = sqrt (e2);
	q0 = (1.0 - e*e)*(sin(lat0)/(1.0 - e*e*sin(lat0)*sin(lat0)) - 1.0/(2.0*e)*
		(log((1.0 - e*sin(lat0))/(1.0 + e*sin(lat0)))));
	q1 = (1.0 - e*e)*(sin(lat1)/(1.0 - e*e*sin(lat1)*sin(lat1)) - 1.0/(2.0*e)*
		(log((1.0 - e*sin(lat1))/(1.0 + e*sin(lat1)))));
	q2 = (1.0 - e*e)*(sin(lat2)/(1.0 - e*e*sin(lat2)*sin(lat2)) - 1.0/(2.0*e)*
		(log((1.0 - e*sin(lat2))/(1.0 + e*sin(lat2)))));
	m1 = cos(lat1)/sqrt(1.0 - e*e*sin(lat1)*sin(lat1));
	m2 = cos(lat2)/sqrt(1.0 - e*e*sin(lat2)*sin(lat2));
	if (lat1 == lat2) n = sin(lat1);
	else n = (m1*m1 - m2*m2)/(q2 - q1);

	C = m1*m1 + n*q1;
	rho0 = a*sqrt(C - n*q0)/n;
	q_max = 1.0 - ((1.0 - e2)/(2.0*e))*log((1.0 - e)/(1.0 + e));
	}

static int axy_ll_e (double x, double y, double *lat, double *lon) {
	double q,rho,theta;
	double ry=rho0-y,err,sign=1.0,tlat,sin_tlat;
	int time=0,maxtime=100;

	if (n < 0.0) sign=-1.0;
	theta = atan2 (x*sign,ry*sign);
	rho = sqrt (x*x+ry*ry);
	q = (C - rho*rho*n*n/(a*a))/n;

	if(q > 2.0 || q < -2.0) return (-1);
	tlat = asin (q/2.0);

	if (q < -q_max || q > q_max) return (-1);
	else
		if (q == q_max) *lat=90.0*DEG_RAD;
		else
			if (q == -q_max) *lat=-90.0*DEG_RAD;
			else
				do {
					sin_tlat = sin (tlat);
					*lat = tlat + (1.0 - e2*sin_tlat*sin_tlat)*(1.0 - e2*sin_tlat*sin_tlat)/
						(2.0*cos(tlat))*
						(q/(1.0-e2)-sin(tlat)/(1.0-e2*sin_tlat*sin_tlat)+1/(2.0*e)*
						log((1.0-e*sin_tlat)/(1.0+e*sin_tlat)));
					err = (*lat-tlat)*(*lat-tlat);
					tlat = *lat;
					} while (err > 0.000000001 && time++ < maxtime);

	*lon = Lon0 + theta/n;
	*lon /= DEG_RAD;
	*lat /= DEG_RAD;

	return (1);
	}

static int all_xy_e (double *x, double *y, double lat, double lon) {
	double q,m,rho,theta;

	if (lat > 90.0 || lat < -90.0) return (-1);

	lat *= DEG_RAD;
	lon *= DEG_RAD;

	q= (1.0 - e*e)*(sin(lat)/(1.0 - e*e*sin(lat)*sin(lat)) - 1.0/(2.0*e)*
		(log((1.0 - e*sin(lat))/(1.0 + e*sin(lat)))));
	m = cos (lat)/sqrt (1.0 - e*e*sin(lat)*sin(lat));
	rho = a * sqrt (C - n*q)/n;
	theta = n * (lon - Lon0);
	*x = rho * sin (theta);
	*y = rho0 - rho * cos (theta);

	return (1);
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/