This is MiscPlanetCoordConverter.m in view mode; [Download] [Up]
/*====================== 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 <misckit/miscgiskit.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;
return;
}
/*---------------------------------------------------------------------------*/
/* 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(TRUE)
{ [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 */
+ initialize
{[MiscPlanetCoordConverter setVersion:MISC_PLANET_CONVERT_VERSION_ID]; return self;}
/*---------------------------------------------------------------------------*/
/* 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
{ [self error:" %s class should not be sent '%s' messages\n",
[[self class] name], sel_getName(_cmd)];
return self;
}
+allocFromZone: (NXZone *)zone
{ [self error:" %s class should not be sent '%s' messages\n",
[[self class] name], sel_getName(_cmd)];
return self;
}
- free {return self;}
+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 */
/*===========================================================================*/
- finishUnarchiving
{
[self free];
return [MiscPlanetCoordConverter new];
}
@end
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.