/*====================== MiscPlanetCoordConverter.m =========================*/
/* MiscPlanetCoordConverter class supports the calculations required for
   converting between World and Universal Transverse Mercator coordinates.

   There is only one instance ever, so unless changes are made, this class
   is NON REENTRANT.

   For information on the underlying mathematics, refer to:

	 1- Ordinance Survey Information, "Transverse Mercator Projection,
	    Constants, Formula and Methods", March 1983

	 2- Ordinance Survey, "Tables for the Transverse Mercator Projection
	    of Ireland", 1953, reprinted 1971

IMPORTANT: These equations are accurate to within 1 millimeter. Calculations
	   requiring greater accuracy must use formulae in:

		Redfern, JCB, "Transverse Mercator Formulae", 1948,
		Empire Survey Review, 9(69) pg318-322

	    Extra decimal places are stored only for the purpose of 
	    slowing error propagation that affects the numbers at the
	    millimeter scale, not because they are meaningful in and of 
	    themselves.

   DMA Release 0.8, Copyright @1993 by Genesis Project, Ltd. All Rights
   Reserved. For further information on terms and conditions see:
		Documentation/GISKit/Agreements-Legal-README

HISTORY
12-Mar-93  Dale Amon at GPL
	   Split UTMGrid into constants and converter parts.
22-Feb-93  Dale Amon at GPL
	   Created.
*/

#import <math.h>
#import "_MiscError.h"
#import <MiscFoundation/MiscPlanetCoordConverter.h>
#import <MiscFoundation/MiscWorldCoord.h>
#import <MiscFoundation/MiscUTMCoord.h>
#import <MiscFoundation/MiscUKUTMCoord.h>
#import <MiscFoundation/MiscZoneUTMCoord.h>
#import <MiscFoundation/MiscUTMConstants.h>

@implementation MiscPlanetCoordConverter

/*=========================================================================*/
/* Internal methods */
/*=========================================================================*/
/* calculations required for both directions of conversion. Caches values for
   later use.

   Uses UTMConstants:
	eSqrd = eccentricity squared
	aF0   = major semiaxis

*/

- (void)blackboardCalc: (double) phi
  {	double	tmp;

	sinPhi  = sin(phi);
	sin2Phi = sinPhi * sinPhi;

	tmp    = 1.0 - xlate->eSqrd * sin2Phi;
	nu     =               xlate->aF0 / sqrt(tmp);
	rho    = (nu * (1.0 - xlate->eSqrd)) / tmp;
	etaSqrd= (nu/rho) - 1.0;
}


/*---------------------------------------------------------------------------*/
/* Calculate the Developed Arc of a Meridian from phi to the True Origin.
   Uses local constants:
	bF0  = semiminor axis, scaled.
	phi0 = latitude of true origin
   Caches values for later use.
*/

- (double) calcM: (double) phi
  {	double	dif,sum;

	dif = phi - xlate->phi0;
	sum = phi + xlate->phi0;

	M   = xlate->bF0 
		  * (xlate->M1 * dif - 
		     xlate->M2 * sin(    dif) * cos(    sum) +
		     xlate->M3 * sin(2.0*dif) * cos(2.0*sum) -
		     xlate->M4 * sin(3.0*dif) * cos(3.0*sum));

	return M;
  }


/*---------------------------------------------------------------------------*/
/* Calculate our estimated latitide, phi prime.
   Uses local constants:
	aF0  = semimajor axis, acled
	N0   = Grid north of True North
	phi0 = latitude of true origin
	convergence = difference at which we consider ourselves done
*/

- (double) calcPhiPrime: (double) N
  {	double	tmp;
	double	phiPrime = (N - xlate->N0)/xlate->aF0 + xlate->phi0;

	while(1)
	 {	[self calcM: phiPrime];
		tmp = N - xlate->N0 - M;
		if (fabs(tmp) < xlate->convergence) break;
		phiPrime += tmp/xlate->aF0;
	 }
	return phiPrime;
  }


/*===========================================================================*/
/* Conversion services */
/*===========================================================================*/
/* Calculate latitude and longitude given grid N and E. Accurate to 1 mm.
   Uses constants:
	lambda0 = longitude of grid origin
*/

- (void) utmToWorld
{	double	E,N;
	double	phiPrime;
	double	y,y2,y3,y4,y5,y6,y7;
	double	nu2,nu3,nu4,nu5,nu7;
	double	tanPhiPrime,tan2PhiPrime,tan4PhiPrime,tan6PhiPrime;
	double	secPhiPrime;
	double	VII, VIII,IX, X, XI, XII, XIIA;

	/* conversion is driven by the source constants */
	xlate = srcConstants;

	E = src[MISC_EASTINGS];
	N = src[MISC_NORTHINGS];

	phiPrime = [self calcPhiPrime: N];

	/* precalculate nu, rho and etaSqrd and powers of nu */
	[self blackboardCalc: phiPrime];
	nu2   = nu*nu;
	nu3   = nu2*nu;
	nu4   = nu2*nu2;
	nu5   = nu3*nu2;
	nu7   = nu4*nu3;

	/* precalculate all the trig values */

	tanPhiPrime  = tan(phiPrime);
	tan2PhiPrime = tanPhiPrime  * tanPhiPrime;
	tan4PhiPrime = tan2PhiPrime * tan2PhiPrime;
	tan6PhiPrime = tan4PhiPrime * tan2PhiPrime;
	secPhiPrime  = 1/cos(phiPrime);

	/* precalculate all the powers of delta E */
	y  = E - xlate->E0;
	y2 = y*y;
	y3 = y2*y;
	y4 = y2*y2;
	y5 = y3*y2;
	y6 = y3*y3;
	y7 = y4*y3;

	/* Calculate the terms used by the latitude and longitude equations */

	VII  = tanPhiPrime/(2.0*rho*nu);
	VIII = tanPhiPrime/(24.0*rho*nu3) *
		(5.0 + 3.0*tan2PhiPrime + etaSqrd - 9.0*etaSqrd*tan2PhiPrime);
	IX   = tanPhiPrime/(720.0*rho*nu5) *
		(61.0 + 90.0*tan2PhiPrime + 45.0*tan4PhiPrime);

	X    = secPhiPrime/nu;
	XI   = secPhiPrime/(6.0*nu3) * (nu/rho + 2.0*tan2PhiPrime);
	XII  = secPhiPrime/(120.0*nu5) *
		(5.0 + 28.0*tan2PhiPrime + 24.0*tan4PhiPrime);
	XIIA = secPhiPrime/(5040.0*nu7) *
		(61.0 + 662.0*tan2PhiPrime +
		 1320.0*tan4PhiPrime + 720.0*tan6PhiPrime);

	/* and finally, we give you the latitude and the longitude */	
	dst[MISC_LATITUDE]  = phiPrime     - y2*VII + y4*VIII + y6*IX;
	dst[MISC_LONGITUDE] = xlate->lambda0 + y*X - y3*XI + y5*XII - y7*XIIA;
	dst[MISC_ALTITUDE]  = src[MISC_ELEVATION];
}


/*---------------------------------------------------------------------------*/
/* Given World Coordinates, calculate the local grid coordinates in
   a UTM projection.
   Uses constants:
	N0,E0 = True origin offset from grid False origin
	lambda0 = longitude of True origin. 
*/

- (void) worldToUTM
{	double	phi,lambda;
	double	cosPhi,cos2Phi,cos3Phi,cos5Phi;
	double	tanPhi,tan2Phi,tan4Phi;
	double	P1, P2, P3, P4, P5, P6;
	double	I,II,III,IIIA,IV,V,VI;

	/* conversion is driven by the destination constants */
	xlate = dstConstants;

	/* assumes the internal storage format is double precision radians */
	phi = src[MISC_LATITUDE]; lambda = src[MISC_LONGITUDE];

	/* Calculate the Developed Arc of a Meridian from phi to the
	   True Origin. */
	[self calcM: phi];

	/* calculate nu, rho and etaSqrd */
	[self blackboardCalc: phi];

	/* precalculate all the trig values that are not on blackboard */
	cosPhi  = cos(phi);
	cos2Phi = cosPhi*cosPhi;
	cos3Phi = cos2Phi*cosPhi;
	cos5Phi = cos2Phi*cos3Phi;

	tanPhi  = tan(phi);
	tan2Phi = tanPhi  * tanPhi;
	tan4Phi = tan2Phi * tan2Phi;

	/* Precalculate all the powers of delta lambda */
	P1 = lambda - xlate->lambda0;
	P2 = P1 * P1;
	P3 = P2 * P1;
	P4 = P2 * P2;
	P5 = P3 * P2;
	P6 = P3 * P3;

	/* Calculate the terms used by the Easting and Northing equations */

	I   = M + xlate->N0;
	II  = nu/2.0   * sinPhi * cosPhi; 
	III = nu/24.0  * sinPhi * cos3Phi * (5.0 - tan2Phi + 9.0 * etaSqrd);
	IIIA= nu/720.0 * sinPhi * cos5Phi * (61.0 - 58.0 * tan2Phi + tan4Phi);
	IV  = nu       * cosPhi;
	V   = nu/6.0   * cos3Phi * (nu/rho - tan2Phi);

	VI  = nu/120.0 * cos5Phi * 
		(5.0 - 18.0*tan2Phi + tan4Phi + 14.0*etaSqrd - 
		 58.0 * etaSqrd * tan2Phi);

	/* And finally, we calculate the local grid coordinates */
	dst[MISC_NORTHINGS] = I         + P2*II + P4*III + P6*IIIA;
	dst[MISC_EASTINGS]  = xlate->E0 + P1*IV + P3*V   + P5*VI;
	dst[MISC_ELEVATION] = src[MISC_ALTITUDE];
 }


/*===========================================================================*/
/* Class methods */
/*===========================================================================*/
/* Initialize the class */

+ (void)initialize
  {[MiscPlanetCoordConverter setVersion:MISC_PLANET_CONVERT_VERSION_ID]; return;}

/*---------------------------------------------------------------------------*/
/* Only one converter of this type is every needed. Of course if we got
   into a really big multiprocess agora system there might be a case
   for multiple convertors. Since this object is shared by many, it
   doesn't particularly matter what zone it is allocated from.
*/

static id thePlanetCoordConverter = nil;

+ new
{
    if (!thePlanetCoordConverter)
	thePlanetCoordConverter = [[super alloc] init];
    return thePlanetCoordConverter;
}


/*---------------------------------------------------------------------------*/
/* Since there is only one object needed, alloc and allocFromZone: are 
   disabled and considered errors. free is simply overridden and made into a
   no op.

   Superalloc is included because we want to allow our subclasses to also 
   create single private instances.
*/

+alloc
{	MISC_ERROR_BAIL;
	return self;
 }

+ allocWithZone:(NSZone *)zone
{	MISC_ERROR_BAIL;
	return self;
 }

- (void)dealloc {
}

+superalloc {return [super alloc];}

/*===========================================================================*/
/* Initialization methods */
/*===========================================================================*/
/* We add a list of all the services which we are able to provide.
   Note that there is only one PlanetCoordConverter ever created. Once 
   initialized the same object is given to all comers and can not be destroyed.
*/

- init
  { Class	world, utm, ukutm, zutm;

    [super init];
    world  = [MiscWorldCoord   class];
    utm    = [MiscUTMCoord     class];
    ukutm  = [MiscUKUTMCoord   class];
    zutm   = [MiscZoneUTMCoord class];

    [self addService: @selector(utmToWorld)   convertsFrom: zutm  to: world]; 
    [self addService: @selector(worldToUTM)   convertsFrom: world to:  zutm]; 
    [self addService: @selector(utmToWorld)   convertsFrom: ukutm to: world]; 
    [self addService: @selector(worldToUTM)   convertsFrom: world to: ukutm]; 
    [self addService: @selector(utmToWorld)   convertsFrom: utm   to: world]; 
    [self addService: @selector(worldToUTM)   convertsFrom: world to: utm]; 
    [self addService: [self fastCopySelector] convertsFrom: ukutm to: ukutm]; 
    [self addService: [self fastCopySelector] convertsFrom: world to: world]; 

    /* This will method return NO if they are in different zones */
    [self addService: [self fastCopySelector] convertsFrom: zutm to: zutm]; 

    return self;
  }


/*===========================================================================*/
/* Archive methods */
/*===========================================================================*/

- awakeAfterUsingCoder:(NSCoder *)aDecoder
{
	[self release];
	return [MiscPlanetCoordConverter new];
}


@end
