direwolf/utm/LatLong-UTMconversion.c
WB2OSZ 8978f2de6c Version 1.0 - Initial commit
Changes to be committed:
	new file:   .gitattributes
	new file:   .gitignore
	new file:   APRStt-Implementation-Notes.pdf
	new file:   CHANGES.txt
	new file:   LICENSE-dire-wolf.txt
	new file:   LICENSE-other.txt
	new file:   Makefile.linux
	new file:   Makefile.win
	new file:   Quick-Start-Guide-Windows.pdf
	new file:   Raspberry-Pi-APRS.pdf
	new file:   User-Guide.pdf
	new file:   aclients.c
	new file:   aprs_tt.c
	new file:   aprs_tt.h
	new file:   atest.c
	new file:   audio.c
	new file:   audio.h
	new file:   audio_win.c
	new file:   ax25_pad.c
	new file:   ax25_pad.h
	new file:   beacon.c
	new file:   beacon.h
	new file:   config.c
	new file:   config.h
	new file:   decode_aprs.c
	new file:   decode_aprs.h
	new file:   dedupe.c
	new file:   dedupe.h
	new file:   demod.c
	new file:   demod.h
	new file:   demod_9600.c
	new file:   demod_9600.h
	new file:   demod_afsk.c
	new file:   demod_afsk.h
	new file:   digipeater.c
	new file:   digipeater.h
	new file:   direwolf.c
	new file:   direwolf.conf
	new file:   direwolf.desktop
	new file:   direwolf.h
	new file:   dsp.c
	new file:   dsp.h
	new file:   dtmf.c
	new file:   dtmf.h
	new file:   dw-icon.ico
	new file:   dw-icon.png
	new file:   dw-icon.rc
	new file:   dw-start.sh
	new file:   dwgps.c
	new file:   dwgps.h
	new file:   encode_aprs.c
	new file:   encode_aprs.h
	new file:   fcs_calc.c
	new file:   fcs_calc.h
	new file:   fsk_demod_agc.h
	new file:   fsk_demod_state.h
	new file:   fsk_filters.h
	new file:   fsk_gen_filter.h
	new file:   gen_packets.c
	new file:   gen_tone.c
	new file:   gen_tone.h
	new file:   hdlc_rec.c
	new file:   hdlc_rec.h
	new file:   hdlc_rec2.c
	new file:   hdlc_rec2.h
	new file:   hdlc_send.c
	new file:   hdlc_send.h
	new file:   igate.c
	new file:   igate.h
	new file:   kiss.c
	new file:   kiss.h
	new file:   kiss_frame.c
	new file:   kiss_frame.h
	new file:   kissnet.c
	new file:   kissnet.h
	new file:   latlong.c
	new file:   latlong.h
	new file:   ll2utm.c
	new file:   misc/README-dire-wolf.txt
	new file:   misc/strcasestr.c
	new file:   misc/strsep.c
	new file:   misc/strtok_r.c
	new file:   morse.c
	new file:   multi_modem.c
	new file:   multi_modem.h
	new file:   ptt.c
	new file:   ptt.h
	new file:   pttest.c
	new file:   rdq.c
	new file:   rdq.h
	new file:   redecode.c
	new file:   redecode.h
	new file:   regex/COPYING
	new file:   regex/INSTALL
	new file:   regex/LICENSES
	new file:   regex/NEWS
	new file:   regex/README
	new file:   regex/README-dire-wolf.txt
	new file:   regex/re_comp.h
	new file:   regex/regcomp.c
	new file:   regex/regex.c
	new file:   regex/regex.h
	new file:   regex/regex_internal.c
	new file:   regex/regex_internal.h
	new file:   regex/regexec.c
	new file:   rrbb.c
	new file:   rrbb.h
	new file:   server.c
	new file:   server.h
	new file:   symbols-new.txt
	new file:   symbols.c
	new file:   symbols.h
	new file:   symbolsX.txt
	new file:   textcolor.c
	new file:   textcolor.h
	new file:   tocalls.txt
	new file:   tq.c
	new file:   tq.h
	new file:   tt_text.c
	new file:   tt_text.h
	new file:   tt_user.c
	new file:   tt_user.h
	new file:   tune.h
	new file:   udp_test.c
	new file:   utm/LatLong-UTMconversion.c
	new file:   utm/LatLong-UTMconversion.h
	new file:   utm/README.txt
	new file:   utm/SwissGrid.cpp
	new file:   utm/UTMConversions.cpp
	new file:   utm/constants.h
	new file:   utm2ll.c
	new file:   version.h
	new file:   xmit.c
	new file:   xmit.h
2015-07-26 20:35:07 -04:00

191 lines
7.1 KiB
C

//LatLong- UTM conversion.c
//Lat Long - UTM, UTM - Lat Long conversions
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "constants.h"
#include "LatLong-UTMconversion.h"
/*Reference ellipsoids derived from Peter H. Dana's website-
http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
Department of Geography, University of Texas at Austin
Internet: pdana@mail.utexas.edu
3/22/95
Source
Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
*/
void LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long,
double *UTMNorthing, double *UTMEasting, char* UTMZone)
{
//converts lat/long to UTM coords. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees
//Written by Chuck Gantz- chuck.gantz@globalstar.com
double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
double k0 = 0.9996;
double LongOrigin;
double eccPrimeSquared;
double N, T, C, A, M;
//Make sure the longitude is between -180.00 .. 179.9
double LongTemp = (Long+180)-(int)((Long+180)/360)*360-180; // -180.00 .. 179.9;
double LatRad = Lat*deg2rad;
double LongRad = LongTemp*deg2rad;
double LongOriginRad;
int ZoneNumber;
ZoneNumber = (int)((LongTemp + 180)/6) + 1;
if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
ZoneNumber = 32;
// Special zones for Svalbard
if( Lat >= 72.0 && Lat < 84.0 )
{
if( LongTemp >= 0.0 && LongTemp < 9.0 ) ZoneNumber = 31;
else if( LongTemp >= 9.0 && LongTemp < 21.0 ) ZoneNumber = 33;
else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35;
else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37;
}
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
LongOriginRad = LongOrigin * deg2rad;
//compute the UTM Zone from the latitude and longitude
sprintf(UTMZone, "%d%c", ZoneNumber, UTMLetterDesignator(Lat));
eccPrimeSquared = (eccSquared)/(1-eccSquared);
N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
T = tan(LatRad)*tan(LatRad);
C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
A = cos(LatRad)*(LongRad-LongOriginRad);
M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
- (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
+ (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
- (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad));
*UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
+ (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
+ 500000.0);
*UTMNorthing = (double)(k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
+ (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
if(Lat < 0)
*UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
}
char UTMLetterDesignator(double Lat)
{
//This routine determines the correct UTM letter designator for the given latitude
//returns 'Z' if latitude is outside the UTM limits of 84N to 80S
//Written by Chuck Gantz- chuck.gantz@globalstar.com
char LetterDesignator;
if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X';
else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W';
else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V';
else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U';
else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T';
else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S';
else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R';
else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q';
else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P';
else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N';
else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M';
else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L';
else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K';
else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J';
else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H';
else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G';
else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F';
else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E';
else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D';
else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C';
else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
return LetterDesignator;
}
void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
double *Lat, double *Long )
{
//converts UTM coords to lat/long. Equations from USGS Bulletin 1532
//East Longitudes are positive, West longitudes are negative.
//North latitudes are positive, South latitudes are negative
//Lat and Long are in decimal degrees.
//Written by Chuck Gantz- chuck.gantz@globalstar.com
double k0 = 0.9996;
double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
double eccPrimeSquared;
double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
double N1, T1, C1, R1, D, M;
double LongOrigin;
double mu, phi1, phi1Rad;
double x, y;
int ZoneNumber;
char* ZoneLetter;
int NorthernHemisphere; //1 for northern hemispher, 0 for southern
x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
y = UTMNorthing;
ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
if (*ZoneLetter == '\0')
{
NorthernHemisphere = 1; //no letter - assume northern hemisphere
}
else if((*ZoneLetter >= 'N' && *ZoneLetter <= 'X') ||
(*ZoneLetter >= 'n' && *ZoneLetter <= 'x'))
{
NorthernHemisphere = 1; //point is in northern hemisphere
}
else
{
NorthernHemisphere = 0; //point is in southern hemisphere
y -= 10000000.0; //remove 10,000,000 meter offset used for southern hemisphere
}
LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
eccPrimeSquared = (eccSquared)/(1-eccSquared);
M = y / k0;
mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
+ (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
+(151*e1*e1*e1/96)*sin(6*mu);
phi1 = phi1Rad*rad2deg;
N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
T1 = tan(phi1Rad)*tan(phi1Rad);
C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
D = x/(N1*k0);
*Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
+(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
*Lat = *Lat * rad2deg;
*Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
*D*D*D*D*D/120)/cos(phi1Rad);
*Long = LongOrigin + *Long * rad2deg;
}