datconv.c

这是经典的卫星编程应用程序

			strcpy(fname,ProgramName);
			if((p=strrchr(fname,'\\')) == NULL) {
				strcpy(pathbuffer,"g7datums.txt");
			}
			else {
				strcpy(++p,"g7datums.txt");
				strcpy(pathbuffer,fname);
			}
		}
		strcpy(fname,pathbuffer);
#endif
		if((in=fopen(fname,"r")) == NULL) {  // can't find file g7datums.txt
//
// check internal array if file not there
//
// ====================
checkinternal:;
			done=0;
			i=-1;
			while(done==0) {
				i++;
				if(stricmp("FINISXXY",gDatum[i].name) == 0) {
					done=1;
					found=0;
					continue;
				}
				if(stricmp(datumstr,gDatum[i].name) != 0) continue;
				found=done=1;
				strcpy(df,gDatum[i].name);
				deltaA=gDatum[i].a;
				deltaF=gDatum[i].invf;
				deltaX=(double)gDatum[i].dX;
				deltaY=(double)gDatum[i].dY;
				deltaZ=(double)gDatum[i].dZ;
			}
			if(found==1) goto setval;
// ====================
			fprintf(stderr,"Error: datum %s not recognized. Defaulting to WGS-84\n",datumstr);
			deltaA = deltaF = deltaX = deltaY = deltaZ = 0.0;
			strcpy(df,"WGS-84");
			found = 1;
			goto setval;
		}
		fileopen=1;
	}
//
// read looking for datum
//
	fseek(in,0L,0);
	done=0;
	while(done==0) {
		if(fgets((char*)bin,MAX_LENGTH,in)==NULL) {
		    done=1;
		    break;
		}
		if(bin[0]==';') continue;
		trim((char*)bin);
		p=strtok(bin,",");
		if(stricmp(datumstr,p) != 0) continue;
		found=done=1;
	  	strcpy(df,datumstr);
		deltaA=atof(strtok(NULL,","));
		deltaF=atof(strtok(NULL,","));
		deltaX=atof(strtok(NULL,","));
		deltaY=atof(strtok(NULL,","));
		deltaZ=atof(strtok(NULL,","));
	}
	if(found==0) goto checkinternal;
setval:;
	if(found==0) {
		fprintf(stderr,"Error: %s datum '%s' not found.  Defaulting to WGS-84\n",
			to==1?"'to'":"'from'",datumstr);
		goto notfound;
	}
	if(to==0) {
		strcpy(inputdatum,df);
		from_a=deltaA;
		from_f=deltaF;
		from_x=deltaX;
		from_y=deltaY;
		from_z=deltaZ;
	}
	else {
		strcpy(outputdatum,df);
		to_a=deltaA;
		to_f=deltaF;
		to_x=deltaX;
		to_y=deltaY;
		to_z=deltaZ;
	}
	return;	
} // void get_datum

PROCX void set_datum(char *from, char * to)
{
	get_datum(0,from);
	get_datum(1,to);
//
// find real ?_a values
//
	from_a = WGS84a - from_a;
	to_a   = WGS84a - to_a;
//
// find real ?_f values
//
// ?_f = (1/wgs84f - 1/?_f)/1e4;
// ?_f = (1/298.257223563 - 1/?_f);
//
	from_f = from_f/10000.0;
	from_f = WGS84f-from_f;
	to_f = to_f/10000.0;
	to_f = WGS84f-to_f;
	to_x=-to_x; to_y=-to_y; to_z=-to_z;
	from_x=-from_x; from_y=-from_y; from_z=-from_z;
//
// save for UTM
//
	UTMa_in=from_a;
	UTMf_in=from_f;
	UTMa=UTMa_out=to_a;
	UTMf=UTMf_out=to_f;
//
// calculate second eccentricity squared
//
	from_es=(2.0*from_f-from_f*from_f);
//
// calculate delta's
//
	dA=to_a-from_a; dF=to_f-from_f; dX=to_x-from_x; dY=to_y-from_y; dZ=to_z-from_z;
} // set_datum


#define Pi (double)M_PI
static const double Degree=Pi/180.0;

static const double lat0 = 0.0;	// reference transverse mercator latitude
static const double k0 = 0.9996;
//============================================================================================
//   "WGS 84",                     6378137.0,   298.257223563 },  /* 20   */
//============================================================================================
static void datumParams(double *a, double *es)
{
    double f = UTMf;				// flattening

    *es = 2 * f - f*f;				// eccentricity^2
    *a =  UTMa;						// semimajor axis
}
//============================================================================================

// Equations from "Map Projections -- A Working Manual", USGS Professional Paper 1395

static double M(double phi, double a, double es);

/* -------------------------------------------------------------------------- */

static void toTM(double lat, double lon, double lat0, double lon0, double k0, double *x, double *y)
{
	double m, et2, n, t, c, A, a, m0, es, lambda, phi, lambda0, phi0;
	
	datumParams(&a, &es);
	
	lambda = lon * Degree;
	phi = lat * Degree;

	phi0 = lat0 * Degree;
	lambda0 = lon0 * Degree;
	
	m0 = M(phi0, a, es);
	m = M(phi, a, es);
	
	et2 = es / (1 - es);
	n = a / sqrt(1 - es * pow(sin(phi), 2.0));
	t = pow(tan(phi), 2.0);
	c = et2 * pow(cos(phi), 2.0);
	A = (lambda - lambda0) * cos(phi);
	*x = k0*n*(A + (1.0 - t + c)*A*A*A/6.0 
			+ (5.0 - 18.0*t + t*t + 72.0*c - 58.0*et2)*pow(A, 5.0) / 120.0);
	*y = k0*(m - m0 + n*tan(phi)*(A*A/2.0
			+ (5.0 - t + 9.0*c + 4*c*c)*pow(A, 4.0)/24.0
			+ (61.0 - 58.0*t + t*t + 600.0*c - 330.0*et2)*pow(A, 6.0)/720.0) );

}

/* -------------------------------------------------------------------------- */

static void fromTM(double x, double y, double lat0, double lon0, double k0, double *lat, double *lon)
{
//	extern struct PREFS gFilePrefs;
	double a, m0, es, et2, m, e1, mu, phi1, c1, t1, n1, r1, d, phi0, lambda0;
	
	phi0 = lat0 * Degree;
	lambda0 = lon0 * Degree;
	
	datumParams(&a, &es);
		
	m0 = M(phi0, a, es);

	et2 = es / (1.0 - es);
	m = m0 + y / k0;
	e1 = (1.0 - sqrt(1.0 - es)) / (1.0 + sqrt(1.0 - es));
	mu = m / (a * (1.0 - es/4.0 - 3.0 * es*es/64.0 - 5.0 * es*es*es/256.0));
	phi1 = mu + (3.0 * e1/2.0 - 27.0 * pow(e1, 3.0)/32.0) * sin(2.0 * mu)
			+ (21.0 * e1*e1/16.0 - 55.0 * pow(e1, 4.0)/32.0)
			* sin(4.0 * mu) + 151.0 * pow(e1, 3.0)/96.0 * sin(6.0 * mu)
			+ 1097.0 * pow(e1, 4.0)/512.0 * sin(8.0 * mu);
	c1 = et2 * pow(cos(phi1), 2.0);
	t1 = pow(tan(phi1), 2.0);
	n1 = a / sqrt(1 - es * pow(sin(phi1), 2.0));
	r1 = a * (1.0 - es) / pow(1.0 - es * pow(sin(phi1), 2.0), 1.5);
	d = x / (n1 * k0);
	*lat = (phi1 - n1 * tan(phi1) / r1
			* (d*d / 2.0 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1*c1 - 9.0 * et2)
			* pow(d, 4.0) / 24.0 + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1*t1
			- 252.0 * et2 - 3.0 * c1*c1) * pow(d, 6.0) / 720.0 )) / Degree;
	*lon = (lambda0 + (d - (1.0 + 2.0 * t1 + c1) * pow(d, 3.0)/6.0
			+ (5.0 -2.0 * c1 + 28.0 * t1 - 3.0 * c1*c1 + 8.0 * et2 + 24.0 * t1*t1)
			* pow(d, 5.0)/120.0) / cos(phi1)) / Degree;
}

/* -------------------------------------------------------------------------- */

static double M(double phi, double a, double es)
{
	if (phi == 0.0)
		return 0.0;
	else {
		return a * (
			( 1.0 - es/4.0 - 3.0*es*es/64.0 - 5.0*es*es*es/256.0 ) * phi -
			( 3.0*es/8.0 + 3.0*es*es/32.0 + 45.0*es*es*es/1024.0 ) * sin(2.0 * phi) +
			( 15.0*es*es/256.0 + 45.0*es*es*es/1024.0 ) * sin(4.0 * phi) -
			( 35.0*es*es*es/3072.0 ) * sin(6.0 * phi) );
	}
}


//============================================================================================

static void calcPhi(double *phi, double e, double t);

/* --------------------------------------------------------------------------------- */

static void toUPS(double lat, double lon, double *x, double *y)
{
//	extern struct PREFS gPrefs;
	double a, t, e, es, rho;
	const double k0 = 0.994;

	double lambda = lon * Degree;
	double phi = fabs(lat * Degree);
	
	datumParams(&a, &es);
	e = sqrt(es);
	t = tan(Pi/4.0 - phi/2.0) / pow( (1.0 - e * sin(phi)) / (1.0 + e * sin(phi)), (e/2.0) );
	rho = 2.0 * a * k0 * t / sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e));
	*x = rho * sin(lambda);
	*y = rho * cos(lambda);

	if (lat > 0.0)	// Northern hemisphere
		*y = -(*y);
	*x += 2.0e6;	// Add in false easting and northing
	*y += 2.0e6;
}

/* --------------------------------------------------------------------------------- */

static void fromUPS(short southernHemisphere, double x, double y, double *lat, double *lon)
{
//	extern struct PREFS gFilePrefs;
	double a, es, e, t, rho;
	const double k0 = 0.994;
	
	datumParams(&a, &es);
	e = sqrt(es);
	
	x -= 2.0e6;		// Remove false easting and northing
	y -= 2.0e6;
	
	rho = sqrt(x*x + y*y);
	t = rho * sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e)) / (2.0 * a * k0);
	
	calcPhi(lat, e, t);
	*lat /= Degree;
	
	if (y != 0.0)
		t = atan(fabs(x/y));
	else {
		t = Pi / 2.0;
		if (x < 0.0) t = -t;
	}

	if (!southernHemisphere)
		y = -y;
		
	if (y < 0.0)
		t = Pi - t;
		
	if (x < 0.0)
		t = -t;

	*lon = t / Degree;
}

/* --------------------------------------------------------------------------------- */

static void calcPhi(double *phi, double e, double t)
{
	double old = Pi/2.0 - 2.0 * atan(t);
	short maxIterations = 20;
	
	while ( (fabs((*phi - old) / *phi) > 1.0e-8) && maxIterations-- ) { 
		old = *phi;
		*phi = Pi/ 2.0 - 2.0 * atan( t * pow((1.0 - e * sin(*phi)) / ((1.0 + e * sin(*phi))), (e / 2.0)) );
	}
}
//============================================================================================
PROCX void DegToUTM(double lat, double lon, char *zone, double *x, double *y)
{
	char nz;
	double lon0;
	
	if ((lat >= -80.0) && (lat <= 84.0)) {
		nz = 'C'+((short)(lat + 80.0)) / 8;
		if (nz > 'H') ++nz;					// skip 'I' and 'O'
		if (nz > 'N') ++nz;
		lon0 = 6.0 * floor(lon / 6.0) + 3.0;	
		sprintf(zone, "%02d%c", ((short)lon0 +183) / 6, nz);
		toTM(lat, lon, lat0, lon0, k0, x, y);
		*x += 500000.0;				// false easting
		if (lat < 0.0)				// false northing for southern hemisphere
			*y = 10000000.0 + *y;
	}
	else {
		strcpy(zone, "00x");		
		if (lat > 0.0)
			if (lon < 0.0) zone[2] = 'Y';
					  else zone[2] = 'Z';
		else
			if (lon < 0.0) zone[2] = 'A';
					  else zone[2] = 'B';
		toUPS(lat, lon, x, y);
	}
} // DegToUTM

/* --------------------------------------------------------------------------------- */

PROCX void UTMtoDeg(short zone, short southernHemisphere, double x, double y, double *lat, double *lon)
{
	double lon0;
	
	if (zone != 0) {
		lon0 = (double)((-183 + 6 * zone));
		if (southernHemisphere)
			y = 1.0e7 - y;			// remove false northing for southern hemisphere
		x -= 500000.0;				//   and false easting
		fromTM(x, y, lat0, lon0, k0, lat, lon);
	}
	else
		fromUPS(southernHemisphere, x, y, lat, lon);
	*lat=fabs(*lat);
	if(southernHemisphere) *lat=-(*lat);
} // UTMtoDEG