#Utilizando el método de Coticchia - Surace, para emisferio norte
#Datos de coordenadas geográficas
# longitud (73° 7'24.97"O) 
# Latitud (  7° 6'18.42"N)
#Datos de semieje mayor y semieje menor
# WGS84

#semiejemayor(a) = 6378137m
#semiejemenor(b) = 6356752.3142m
a <- (6378137)
a
## [1] 6378137
b <- (6356752.31424)
b
## [1] 6356752
#Hallando la excentricidad
excentr <- (((a)^2-(b)^2)^0.5/b)
excentr
## [1] 0.08209444
#Hallando excentricidad al cuadrado
e <-(excentr)^2
e
## [1] 0.006739497
#Hallando radio polar de la  curvatura
c <- (((a)^2)/b)
c
## [1] 6399594
#Hallando el aplanamiento
aplan <- (a-b)/(a)
aplan
## [1] 0.003352811
#Conversion de latitud y longitud en grados decimales
long <- (73+(7/60)+((24.97)/3600))
long
## [1] 73.1236
lat <- (7+(6/60)+(18.42/3600))
lat
## [1] 7.105117
#Conversión de latitud y longitud en radianes
longr <- long*(pi/180)
longr
## [1] 1.276248
latr <- lat*(pi/180)
latr
## [1] 0.1240077
# Cálculo de los signos latitud y longitud.
#valores positivos longitud cuando se escribe como este "E" "(+)"
#Latitud norte, tiene valor positivo, "N" es "(+)"

# Valores negativos longitud oeste "O" es (-)
#                      Latitud Sur "S" es (-)

signlong <-(long*(-1))
signlong
## [1] -73.1236
signlat <-(lat*(1))
signlat
## [1] 7.105117
# Cálculo del huso
huso <- ((signlong/6)+31)
huso
## [1] 18.81273
# Hallando el meridiano central del huso, se usa el valor entero del uso.
# poner el valor entrro del huso.
lamda1 <- ((18*6)-183)
lamda1
## [1] -75
# Distancia angular entre la longitud de un punto y el meridiano.
# Se pone signo negativo por la ubicación de la longitud en radianes
distang <- (-longr - (lamda1*(pi/180)))
distang
## [1] 0.03274931
# Cálculo de los parametros de la ecuación de Coticchia-Surace
#

A <- (cos(latr)*sin(distang))
A
## [1] 0.03249202
#
Xi <- (0.5*log((1+A)/(1-A)))
Xi
## [1] 0.03250346
#
eta <- (atan((tan(latr))/(cos(distang)))-latr)
eta
## [1] 6.584941e-05
#
nu <-((c)/(1+e*(cos(latr)^2))^0.5*0.9996)
nu
## [1] 6375912
# Donde 0.9996 es un valor de factor.
zeta <- (e/2)*Xi^2*(cos(latr))^2
zeta
## [1] 3.505588e-06
#
A1 <-sin(2*latr)
A1
## [1] 0.2454805
#
A2 <- A1*(cos(latr))^2
A2
## [1] 0.2417249
#
j2 <- (latr + (A1/2))
j2
## [1] 0.2467479
#
j4 <- ((3*j2+A2)/4)
j4
## [1] 0.2454922
#
j6 <- ((5*j4+A2*(cos(latr))^2)/3)
j6
## [1] 0.4884958
#
alfa <-((3/4)*e)
alfa
## [1] 0.005054623
#
beta <-((5/3)*(alfa)^2)
beta
## [1] 4.258202e-05
#
gama <- ((35/27)*(alfa)^3)
gama
## [1] 1.674058e-07
#
Bo <-(0.9996*c*(latr-(alfa*j2)+(beta*j4)-(gama*j6)))
Bo
## [1] 785369.2
# Hallando  los valores de las coordendas
X <- ((Xi*nu*(1+(zeta/3))) +500000)
X
## [1] 707239.4
Y <- (eta*nu*(1+zeta))+Bo
Y
## [1] 785789
# Los valores de las coordenadas en Google Earth
Xg <- 707239.43 
Yg <- 785789.03 
#Hallando la diferencia de valores con los valores de Google Earth
dX <- abs(X-Xg)
dX
## [1] 0.002195957
#
dy <- abs(Y-Yg)
dy
## [1] 0.007561922
# Hallando el error de la distancia absoluta a un nivel de confianza del 95% en metros.
Error <- (((dX)^2+(dy)^2)^0.5)*1.7308
Error
## [1] 0.01362887