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

#define MPF (0.3048)
#define RAD (3.141592654/180.0)
#define NREAD 200
#define NOB 100
#define NCOORD 1000
#define NBENCH 40
#define FLAGCHAR 8

char str[120], lastob[FLAGCHAR]="", tempfl[FLAGCHAR]="";

struct readstruc {
	char flag[FLAGCHAR];
	float azim;
	float vert;
	float slope;
	float rx;
	float ry;
	float rz;
	} ;

struct coordstruc {
	char flag[FLAGCHAR];
	float x;
	float y;
	float z;
	float dx;
	float dy;
	float dz;
	} bench[NBENCH], coord[NCOORD];

struct obstruc {
	char ob[FLAGCHAR];
	float x0;
	float y0;
	float z0;
	float azim0;
	int nread;
	struct readstruc read[NREAD];
	} obs[NOB];

main()
	{
	int iob, tob, iread, ibench, icoord, ncoord, coflg, nbench, nob;
	float azmin, azsec, vertmin, vertsec, horiz, temp;
	float sx, sx2, sy, sy2, sxy, r, m, b;
	float westx, westy;
	
	ncoord = nbench = nob = 0;
	for (iob=0; iob<NOB; iob++)
		{
		obs[iob].nread = 0;
		obs[iob].x0 = 0.0;
		obs[iob].y0 = 0.0;
		obs[iob].z0 = 0.0;
		obs[iob].azim0 = 0.0;
		}

	iob = iread = ibench = 0;
	while (gets(str) != NULL)
		{
		if (strncmp(str+0, "BM", 2) == 0 && nbench < NBENCH)
			{
			sscanf(str, "%s %f %f %f", bench[ibench].flag,
			&bench[ibench].x, &bench[ibench].y, &bench[ibench].z);
			bench[ibench].dx = 0.0;
			bench[ibench].dy = 0.0;
			bench[ibench].dz = 0.0;
			nbench++;
			}
		
		else if (strncmp(str+0, "OB", 2) == 0)
			{
			if (strcmp(lastob, "") == 0)
				{
				sscanf(str+0, "%s", lastob);
				nob = 1;
				}
			sscanf(str+0, "%s", tempfl);
			if (strcmp(tempfl, lastob) != 0)
				{
				if (nob >= NOB)
					{
					fprintf(stderr,
					"surv: max. no. OBs = %d exceeded, ignoring further data\n",
					NOB);
					break;
					}
				sscanf(str+0, "%s", lastob);
				nob++;
				iob++;
				iread = 0;
				}
			else if (obs[iob].nread >= NREAD)
				{
				fprintf(stderr,
				"surv: max. no. flags per OB = %d exceeded for OB = %s, skipping rest of OB\n",
				NREAD, lastob);
				continue;
				}
			
			sscanf(str, "%s %s %f %f %f %f %f %f %f",
			obs[iob].ob, obs[iob].read[iread].flag,
			&obs[iob].read[iread].azim, &azmin, &azsec,
			&obs[iob].read[iread].vert, &vertmin, &vertsec,
			&obs[iob].read[iread].slope);

			if (obs[iob].read[iread].slope != 0.0)
			    {
			    azmin += azsec/60.0;
			    obs[iob].read[iread].azim += azmin/60.0;
			    vertmin += vertsec/60.0;
			    obs[iob].read[iread].vert += vertmin/60.0;
			    obs[iob].read[iread].slope *= MPF;
    
			    horiz = obs[iob].read[iread].slope
			    * sin(RAD*(double)(obs[iob].read[iread].vert));
			    obs[iob].read[iread].rx = horiz
			    * sin(RAD*(double)(obs[iob].read[iread].azim));
			    obs[iob].read[iread].ry = horiz
			    * cos(RAD*(double)(obs[iob].read[iread].azim));
			    obs[iob].read[iread].rz = obs[iob].read[iread].slope
			    * cos(RAD*(double)(obs[iob].read[iread].vert));
			    }
			else
			    {
			    obs[iob].read[iread].azim = 0.0;
			    obs[iob].read[iread].vert = 0.0;
    
			    obs[iob].read[iread].rx = 0.0;
			    obs[iob].read[iread].ry = 0.0;
			    obs[iob].read[iread].rz = 0.0;
			    }

			obs[iob].nread++;
			iread++;
			}
		}
	/*
	for (iob=0; iob<nob; iob++)
		for (iread=0; iread < obs[iob].nread; iread++)
			printf("1 %s %s %f %f %f %f %f %f\n",
			obs[iob].ob, obs[iob].read[iread].flag,
			obs[iob].read[iread].azim, obs[iob].read[iread].vert,
			obs[iob].read[iread].slope/MPF,
			obs[iob].read[iread].rx/MPF,
			obs[iob].read[iread].ry/MPF,
			obs[iob].read[iread].rz/MPF);
	*/
	
	for (iob=0; iob<nob; iob++)
		{
		for (iread=0; iread < obs[iob].nread; iread++)
			{
			for (ibench=0; ibench<nbench; ibench++)
				{
				if (strcmp(obs[iob].read[iread].flag,
				bench[ibench].flag) == 0)
					{
					tob = iob;
					obs[iob].z0 = bench[ibench].z
						- obs[iob].read[iread].rz;
					obs[iob].x0 = bench[ibench].x
						- obs[iob].read[iread].rx;
					obs[iob].y0 = bench[ibench].y
						- obs[iob].read[iread].ry;
					obs[iob].azim0 = 0.0;
					break; break;
					}
				}
			}
	
		for (tob=iob; tob > 0; tob--)
			connect(iob, tob-1);
		for (tob=iob; tob < nob-1; tob++)
			connect(iob, tob+1);
		}
					
	ncoord = 0;
	for (iob=0; iob<nob; iob++)
		for (iread=0; iread < obs[iob].nread; iread++)
			{
			coflg = 0;
			for (icoord=0; icoord<ncoord; icoord++)
				if (strcmp(coord[icoord].flag,
				obs[iob].read[iread].flag) == 0)
					{
					coflg = 1;
					break;
					}
			if (coflg == 0)
				{
				if (ncoord >= NCOORD)
					{
					fprintf(stderr,
					"surv: ncoord > max=%d, ignoring further data\n",
					NCOORD);
					break;
					}
					
				strcpy(coord[ncoord].flag,
					obs[iob].read[iread].flag);
				coord[ncoord].x = 
					obs[iob].read[iread].rx + obs[iob].x0;
				coord[ncoord].y = 
					obs[iob].read[iread].ry + obs[iob].y0;
				coord[ncoord].z = 
					obs[iob].read[iread].rz + obs[iob].z0;
				coord[ncoord].dx = 0.0;
				coord[ncoord].dy = 0.0;
				coord[ncoord].dz = 0.0;
				ncoord++;
				}
			else
				{
				temp = obs[iob].read[iread].rx + obs[iob].x0
					- coord[icoord].x;
				temp = temp<0 ? -temp : temp;
				if (coord[icoord].dx < temp)
					coord[icoord].dx = temp;
				temp = obs[iob].read[iread].ry + obs[iob].y0
					- coord[icoord].y;
				temp = temp<0 ? -temp : temp;
				if (coord[icoord].dy < temp)
					coord[icoord].dy = temp;
				temp = obs[iob].read[iread].rz + obs[iob].z0
					- coord[icoord].z;
				temp = temp<0 ? -temp : temp;
				if (coord[icoord].dz < temp)
					coord[icoord].dz = temp;
				}
			}

	printf("Resolved Station Coordinates, meters\n");
	printf(" UTM-E (x)  UTM-N (y) Elev-z    dx    dy    dz Station\n");
	printf("------------------------------------------------------\n");
	for (iob=0; iob<nob; iob++)
		{
		printf("%10.1f %10.1f %6.1f  ----  ----  ---- %s\n",
		obs[iob].x0, obs[iob].y0, obs[iob].z0, obs[iob].ob);
		}
	for (icoord=0; icoord<ncoord; icoord++)
		{
		if (coord[icoord].dx == 0.0 && coord[icoord].dy == 0.0
		&& coord[icoord].dz == 0.0)
			printf("%10.1f %10.1f %6.1f  ----  ----  ---- %s\n",
			coord[icoord].x, coord[icoord].y, coord[icoord].z,
			coord[icoord].flag);
		else
			printf("%10.1f %10.1f %6.1f %5.2f %5.2f %5.2f %s\n",
			coord[icoord].x, coord[icoord].y, coord[icoord].z,
			coord[icoord].dx, coord[icoord].dy, coord[icoord].dz,
			coord[icoord].flag);
		}
	printf("---------------------------------------------------\n");

	/* Compute and print projection line */
	sx = 0;
	sx2 = 0;
	sy = 0;
	sy2 = 0;
	sxy = 0;
	westx = coord[0].x;
	westy = coord[0].y;
	for (icoord=0; icoord<ncoord; icoord++)
		{
		sx += coord[icoord].x;
		sx2 += coord[icoord].x * coord[icoord].x;
		sy += coord[icoord].y;
		sy2 += coord[icoord].y * coord[icoord].y;
		sxy += coord[icoord].x * coord[icoord].y;
		if (coord[icoord].x < westx)
			{
			westx = coord[icoord].x;
			westy = coord[icoord].y;
			}
		}
	if (sx2 - sx*sx/ncoord == 0)
		{
		printf("S-N Projection Line\n");
		for (icoord=0; icoord<ncoord; icoord++)
			printf("%s %8.1f %7.1f %10.1f %10.1f\n",
			coord[icoord].flag,
			coord[icoord].y - westy, coord[icoord].z,
			coord[icoord].y, coord[icoord].x);
		}
	else
		{
		m = (sxy - sx*sy/ncoord)/(sx2 - sx*sx/ncoord);
		b = sy/ncoord - m*sx/ncoord;
		if (m == 0.0)
			{
			printf("W-E Projection Line\n");
			for (icoord=0; icoord<ncoord; icoord++)
				printf("%s %8.1f %7.1f %10.1f %10.1f\n",
				coord[icoord].flag,
				coord[icoord].x - westx, coord[icoord].z,
				coord[icoord].y, coord[icoord].x);
			}
		else
			{
			temp = (b + westx/m -westy)/(1.0/m - m);
			printf("Projection Origin UTM %10.1f m East %10.1f m North\n",
			temp, m*temp + b);
			printf("Slope %g m N/m E, Direction N %g deg E\n",
			m, (float)(atan2((double)(1.0), (double)(m)))/RAD);
			m = -1.0/m;
			for (icoord=0; icoord<ncoord; icoord++)
				{
				temp = (m * coord[icoord].x -
				coord[icoord].y + westy - m*westx)
				/ -(float)(sqrt((double)(m*m + 1)));
				printf("%s %8.1f %7.1f %10.1f %10.1f\n",
				coord[icoord].flag, temp, coord[icoord].z,
				coord[icoord].y, coord[icoord].x);
				}
			}
		}	
	
	exit(0);
	}

#define NLAP 20
connect(iob0, iob1)
int iob0, iob1;
	{
	int nlap, ir0, ir1;
	struct lapstruc { int i0; int i1;} laps[NLAP];
	
	nlap = 0;
	for (ir0=0; ir0 < obs[iob0].nread && nlap < NLAP; ir0++)
		for (ir1=0; ir1 < obs[iob1].nread; ir1++)
			if (strcmp(obs[iob0].read[ir0].flag, obs[iob1].read[ir1].flag) == 0)
				{
				laps[nlap].i0 = ir0;
				laps[nlap].i1 = ir1;
				nlap++;
				break;
				}
	if (nlap == 0)
		return(nlap);
	/*
	else if (nlap == 1)
		{
	*/
		obs[iob1].z0 = obs[iob0].z0 + obs[iob0].read[laps[0].i0].rz
			- obs[iob1].read[laps[0].i1].rz;
		obs[iob1].x0 = obs[iob0].x0 + obs[iob0].read[laps[0].i0].rx
			- obs[iob1].read[laps[0].i1].rx;
		obs[iob1].y0 = obs[iob0].y0 + obs[iob0].read[laps[0].i0].ry
			- obs[iob1].read[laps[0].i1].ry;
		obs[iob1].azim0 = 0.0;
	/*
		}
	else if (nlap > 1)
		{
		
		}
	*/
	
	return(nlap);
	}