#Prvi geodetski zadatak
#Zadano: reducirane koordinate prve točke, duljina geodetske linije, smjerni kut geodetske linije
#Traži se: koordinate druge točke geodetske linije i drugi smjerni kut

#1 - unos reduciranih koordinata prve točke ( E,N )
#2 - računanje nereducirnih koordinata prve točke
#3 - parametri elipsoida GRS80 - a, b, e, n = (a - b )/(a + b ), e'
#4 - računanje veličina: A, ψ, c1, c2, c3, c4, c5
#5 - računanje ϕ(N1)
#6 - računanje R1 ( polumjer zakrivljenosti ), t1, η
#7 - računanje pomoćnih veličina - u, v

#8 - kuteve je potrebno pretvoriti u radiane - math.radians (x); math.degrees(x), delimiter za decimale je "."

import math #import math modula u kojem su trigonometrijske funkcije i pretvorba

#unos reduciranih kooridnata točke 1
print ("Unesite reducirane koordinate točke 1 u metrima ( N1, E1):")
N1 = float(input ("N1 = "))
E1 = float(input ("E1 = "))

#unos duljine geodetske linije
print ("Unesite duljinu geodetske linije:")
s12 = float(input ("s12 = "))

#unos zadanog smjernog kuta
print ("Unesite smjerni kut s točke 1 na točku 2 u D-M-S formatu:")
T12_dms = input ("T1-2 = ")

#izračun nereduciranih koordinata
N1_nered = N1 / 0.9999
E1_nered = (E1 - 500000) / 0.9999

#pretvaranje kuta D-M-S -> DEG -> RAD
def dms_deg_rad (kut):
    components = kut.split("-")
    st_deg = int ( components[0] )
    min_deg = int ( components[1] )
    sec_deg = float ( components[2])
    kut_deg = float (st_deg + min_deg/60.0 + sec_deg/3600.0)
    kut_rad = math.radians (kut_deg)
    return (kut_rad)

T12_rad = dms_deg_rad ( T12_dms )
T12_deg = math.degrees (T12_rad)

a = 6378137.00000000000000000  #velika poluos elipsoida GRS80 izražena u metrima
b = 6356752.31414034743838862  #mala poluos elipsoida GRS80 izražena u metrima
n = (a - b)/(a + b) #treća spljoštenost
e = math.sqrt((a**2 - b**2))/a #prvi numerički ekcentricitet
e_crtano = math.sqrt((a**2 - b**2))/b #drugi numerički ekcentricitet

#računanje kuta fi
A = a* (1 - n)*(1 - n**2)*(1 + 9/4*(n**2) + 225/64 * (n**4)) #pomoćni parametri
c1 = 3/2*n - 29/12*(n**3) + 553/80 * (n**5)
c2 = 21/8 * (n**2) - 1537/128 * (n**4) 
c3 = 151/24 * (n**3) - 32373/640 * (n**5)
c4 = 1097/64 * (n**4)
c5 = 8011 / 160 * (n**5)
psi = N1_nered / A
fi1 = psi + math.sin (2*psi)*(c1 +(c2 + (c3 + (c4 + c5*math.cos(2*psi))*math.cos(2*psi))*math.cos(2*psi))*math.cos(2*psi))
fi1_deg = math.degrees(fi1)

#računanje pomoćnih veličina
R1 = a / (math.sqrt(1 - (e**2) * (math.sin(fi1))**2 )) #polumjer zakrivljenosti po prvom vertikalu u točki T1
t1 = math.tan (fi1)
eta1 = e_crtano * math.cos (fi1)

#računanje u i v
u = s12 * math.cos ( T12_rad )
v = s12 * math.sin ( T12_rad )

#računanje nereduciranih koordinata točke 2
N2_nered = N1_nered + u + ((1 + eta1**2)/(2* R1**2) + (1 + 6* eta1**2)* E1_nered**2 / (24* R1**4) + E1_nered**4 / (720* R1**6))* E1_nered**2 * u + \
        t1* (-eta1**2 - eta1**4)* E1_nered**2 * u**2 / R1**3 + ((1 + eta1**2)/ R1**2 + (2*(1 + 3*eta1**2)* E1_nered**2 / (3*R1**4)) + 2*E1_nered**4 / (15 * R1**6))* E1_nered * u * v + \
        t1* (eta1**2 + eta1**4)* E1_nered**2 * v**2 / R1**3 + ((-1 - 3*eta1**2 + t1**2*eta1**2)*E1_nered**2/(3*R1**4) + (5* (-1)*E1_nered**4)/(18*R1**6))* u**3 + \
        2*t1*(-eta1**2 - eta1**4)*E1_nered*(u**2)*v/R1**3 + ((1 + eta1**2)/(3*R1**2)+ ((7 + 20*eta1**2 - 2*t1**2*eta1**2)*E1_nered**2/(6*R1**4) + 61*E1_nered**4/(72* R1**6)))*u*v**2 + \
        2*t1*(eta1**2 + eta1**4)*E1_nered*v**3/(3*R1**3) + ((-1 - 4*eta1**2 + 2*t1**2*eta1**2)*E1_nered/(3*R1**4) + 11* (-1)* E1_nered**3/(9*R1**6))*u**3*v + \
        5*t1*(-eta1**2 - eta1**4)*u**2*v**2/(6*R1**3) + ((2 + 6*eta1**2 - t1**2*eta1**2)*E1_nered/(3*R1**4) + 13*E1_nered**3/(9*R1**6))*u*v**3 + \
        t1*(eta1**2 + eta1**4)*v**4/(6*R1**3) + E1_nered**2*u**5/(15*R1**6) + ((-2 - 13*eta1**2+9*t1**2*eta1**2)/(30*R1**4) + 41*(-1)*E1_nered**2/(30*R1**6))*u**3*v**2 + \
        ((4 + 13*eta1**2 - 3*t1**2*eta1**2)/(30*R1**4) + 16*E1_nered**2/(15*R1**6))*u*v**4 + 2*E1_nered*u**5*v/(45*R1**6) + 26*(-1)*E1_nered*u**3*v**3/(45*R1**6) + 17*E1_nered*u*v**5/(45*R1**6)

E2_nered = E1_nered + v + ((1 + eta1**2)/(2*R1**2) + ((1+6*eta1**2)*E1_nered**2/(24*R1**4) + E1_nered**4/(720*R1**6)))*E1_nered**2*v + \
        ((-1 - eta1**2)/(2*R1**2) + ((-1 - 3*eta1**2)* E1_nered**2/(3*R1**4) + (-1)*E1_nered**4/(15*R1**6)))*E1_nered*u**2 + \
        2*t1*(-eta1**2 - eta1**4)*E1_nered**2*u*v/(R1**3) + ((1+ eta1**2)/(2*R1**2) + ((1+3*eta1**2)*E1_nered**2 /(3*R1**4 + E1_nered**4/(15*R1**6))))*E1_nered*v**2 + \
        2*t1*(eta1**2 + eta1**4)*E1_nered*u**3 / (3*R1**3)+ ((-1-eta1**2)/(6*R1**2) + ((-13 - 38*eta1**2 + 8*t1**2*eta1**2)*E1_nered**2/(12*R1**4)) + 121*(-1)*E1_nered**4/(144*R1**6))*u**2*v + \
        2*t1*(-eta1**2 - eta1**4)*E1_nered*u*v**2/R1**3 + ((1 + eta1**2)/(6*R1**2) + (5 + 14*eta1**2)*E1_nered**2/(12*R1**4) + 41*E1_nered**4/(144*R1**6))*v**3 + \
        ((1 + 6*eta1**2 - 4*t1**2*eta1**2)*E1_nered/(24*R1**4) + 5*E1_nered**3/(18*R1**6))*u**4 + t1*(eta1**2 + eta1**4)*u**3*v/(3*R1**3) + \
        ((-9-30*eta1**2 + 10*t1**2*eta1**2)*E1_nered/(12*R1**4) + 2*(-1)*E1_nered**3/R1**6)*u**2*v**2 + 2*t1*(-eta1**2-eta1**4)*u*v**3/(3*R1**3) + \
        ((5 + 14*eta1**2)*E1_nered/(24*R1**4) + 7*E1_nered**3/(18*R1**6))*v**4 + ((1 + 14*eta1**2 - 12*t1**2*eta1**2)/(120*R1**4) + 121*E1_nered**2/(240*R1**6))*u**4*v + \
        ((-9-38*eta1**2 + 18*t1**2*eta1**2)/(60*R1**4) + 209*(-1)*E1_nered**2/(120*R1**6))*u**2*v**3 + ((5 + 14*eta1**2)/(120*R1**4) + 61*E1_nered**2/(240*R1**6))*v**5 + \
        (-1)*E1_nered*u**6/(720*R1**6) + 179*E1_nered * u**4* v**2/(720*R1**6) + 479*(-1)*E1_nered*u**2*v**4/(720*R1**6) + 61*E1_nered*v**6/(720*R1**6)

#računanje reduciranih koordinata točke 2
E2 = round(E2_nered*0.9999 + 500000,3)
N2 = round(N2_nered*0.9999,3)

#računanje smjernog kuta T21 u radijanima (bez dodavanja 180 deg)
T21_rad = T12_rad + ((-1 - eta1**2)/(R1**2) + ((-1-6*eta1**2)*E1_nered**2/(6*R1**4) + (-1)*E1_nered**4/(120*R1**6)))*E1_nered*u + \
        2/R1**3*t1*(-eta1**2 - eta1**4)*E1_nered**2*v + 2/R1**3*t1*(eta1**2 + eta1**4)*E1_nered*u**2 + \
        (1/(2*R1**2)*(-1 - eta1**2) + (1/R1**4 * (-1 - 4*eta1**2 + t1**2*eta1**2)* E1_nered**2 + 1/(3*R1**6)*(-1)*E1_nered**4))*u*v + 2/R1**3*t1*(-eta1**2 - eta1**4)*E1_nered*v**2 + \
        (1/(6*R1**4)*(1 + 6*eta1**2 - 4*t1**2*eta1**2)*E1_nered + 13*E1_nered**3/(36*R1**6))*u**3 + 4/(3*R1**3)*t1*(eta1**2 + eta1**4)*u**2*v + \
        (1/(6*R1**4)*(-5-22*eta1**2 + 8*t1**2*eta1**2)*E1_nered + 41/(36*R1**6)*(-1)*E1_nered**3)*u*v**2 + 2/(3*R1**3)*t1*(-eta1**2 - eta1**4)*v**3 + \
        (1/(24*R1**4)*(1 + 14*eta1**2 - 12*t1**2*eta1**2) + 5*E1_nered**2/(6*R1**6))*u**3*v + ((-5-26*eta1**2 + 12*t1**2*eta1**2)/(24*R1**4) + 7*(-1)*E1_nered**2/(6*R1**6))*u*v**3 + \
        1/(120*R1**6)*(-1)*E1_nered*u**5 + 29/(60*R1**6)*E1_nered*u**3*v**2 + 61/(120*R1**6)*(-1)*E1_nered*u*v**4

#računanje smjernog kuta T21 u stupnjevima
if (T12_deg < 180):
    T21_deg = math.degrees (T12_rad) + 180
else:
    T21_deg = math.degrees (T12_rad) - 180
    
#pretvaranje smjernog kuta T21 DEG -> D-M-S format
def deg_dms (kut):
    st_dms = int (kut)
    min_dms = int ((kut - st_dms) * 60 )
    sec_dms = float (((( kut - st_dms) *60 )- min_dms) *60)
    sec_dms_zaokruzeno = round ( sec_dms, 5)
    kut_dms = str(st_dms)+"-"+str(min_dms)+"-"+str(sec_dms_zaokruzeno)
    return (kut_dms)

T21_dms = deg_dms (T21_deg)
fi1_dms = deg_dms (fi1_deg)
psi_dms = deg_dms (psi)

#ispis podataka
print ("Zadano: ")
print ("N1 = %.3f m, E1 = %.3f m, s12 = %.3f m, T1-2 = %s" %(N1, E1, s12, T12_dms))
print ("Izračunato: ")
print ("N2 = %.3f m, E2 = %.3f m, T2-1 = %s" %(N2, E2, T21_dms))

#ispis u .txt datoteku
prvi_geodetski_zadatak = open('prvi_geodetski_zadatak.txt', mode = 'w')
podaci = [
    'Program izračunava prvi geodetski zadatak!!'
    '\nZadano:',
    '  Točka 1',
    '   N1 = ' + str(N1),
    '   E1 = ' + str(E1),
    '   s12 = ' + str(s12),
    '   ϑT1-2 = ' + str(T12_dms),
    '\nIzračunato:',
    '  Točka 2',
    '   N2 = ' + str(N2),
    '   E2 = ' + str(E2),
    '   ϑT2-1 = ' + str(T21_dms),
    '\nPomoćne veličine:',
    ' R1 = ' + str(R1),
    ' u = ' + str(u),
    ' v = ' + str(v),
    ' A = ' + str(A),
    ' c1 = ' + str(c1),
    ' c2 = ' + str(c2),
    ' c3 = ' + str(c3),
    ' c4 = ' + str(c4),
    ' c5 = ' + str(c5),
    ' φ = ' + str(fi1_dms),
    ' ψ = ' + str(psi_dms),
    ' N1 nereducirano = ' + str(N1_nered),
    ' E1 nereducirano = ' + str(E1_nered),
    ' N2 nereducirano = ' + str(N2_nered),
    ' E2 nereducirano = ' + str(E2_nered),
    '\nParametri elipsoida GRS80:',
    ' a = ' + str(a),
    ' b = ' + str(b),
    ' n = ' + str(n),
    ' e = ' + str(e),
    ' e\' = ' + str(e_crtano)
]
prvi_geodetski_zadatak.write('\n'.join(podaci))
prvi_geodetski_zadatak.close()