pj_somerc.cpp

projapi是一个关于GIS行业投影转换的程序库

#include "stdafx.h"

#ifndef lint
static const char SCCSID[]="@(#)PJ_somerc.c	4.1	95/08/09	GIE	REL";
#endif
#define PROJ_PARMS__ \
	double	K, c, hlf_e, kR, cosp0, sinp0;
#define PJ_LIB__
#include	"projects.h"
#include <math.h>

PROJ_HEAD(somerc, "Swiss. Obl. Mercator") "\n\tCyl, Ell\n\tFor CH1903";
#define EPS	1.e-10
#define NITER 6
FORWARD(e_forward);
	double phip, lamp, phipp, lampp, sp, cp;

	sp = P->e * sin(lp.phi);
	phip = 2.* atan( exp( P->c * (
		log(tan(FORTPI + 0.5 * lp.phi)) - P->hlf_e * log((1. + sp)/(1. - sp)))
		+ P->K)) - HALFPI;
	lamp = P->c * lp.lam;
	cp = cos(phip);
	phipp = aasin(P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
	lampp = aasin(cp * sin(lamp) / cos(phipp));
	xy.x = P->kR * lampp;
	xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp));
	return (xy);
}
INVERSE(e_inverse); /* ellipsoid & spheroid */
	double phip, lamp, phipp, lampp, cp, esp, con, delp;
	int i;

	phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI);
	lampp = xy.x / P->kR;
	cp = cos(phipp);
	phip = aasin(P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
	lamp = aasin(cp * sin(lampp) / cos(phip));
	con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c;
	for (i = NITER; i ; --i) {
		esp = P->e * sin(phip);
		delp = (con + log(tan(FORTPI + 0.5 * phip)) - P->hlf_e *
			log((1. + esp)/(1. - esp)) ) *
			(1. - esp * esp) * cos(phip) * P->rone_es;
		phip -= delp;
		if (fabs(delp) < EPS)
			break;
	}
	if (i) {
		lp.phi = phip;
		lp.lam = lamp / P->c;
	} else
		I_ERROR
	return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(somerc)
	double cp, phip0, sp;

	P->hlf_e = 0.5 * P->e;
	cp = cos(P->phi0);
	cp *= cp;
	P->c = sqrt(1 + P->es * cp * cp * P->rone_es);
	sp = sin(P->phi0);
	P->cosp0 = cos( phip0 = aasin(P->sinp0 = sp / P->c) );
	sp *= P->e;
	P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * (
		log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e *
		log((1. + sp) / (1. - sp)));
	P->kR = P->k0 * sqrt(P->one_es) / (1. - sp * sp);
	P->inv = e_inverse;
	P->fwd = e_forward;
ENDENTRY(P)