import arcpy
import numpy
import math
import os

#inShapefile1 = arcpy.GetParameterAsText(0)
#inShapefile1 = "c:/gis_toolbox/input_data/pisek/sada1/XYpf1.shp"
#inShapefile2 = arcpy.GetParameterAsText(1)
#inShapefile2 = "c:/gis_toolbox/input_data/pisek/sada1/XYpf2.shp"
nazev = arcpy.GetParameterAsText(1)
path = arcpy.GetParameterAsText(2)
outSHPpoints =nazev+".shp"
inList = arcpy.GetParameterAsText(0)
LIST=[]
for fc in inList.split(';'):
    LIST.append(fc)



#path="c:/gis_toolbox/input_data/pisek/python_test"
#cesta ke slozce
#outSHPpoints="bodovkahk.shp"
#nazev *.shp
BODU_V_PF = arcpy.GetParameterAsText(3)
BODU_V_PF=int(BODU_V_PF)
#BODU_V_PF = 25
BODU_MEZI_PF = arcpy.GetParameterAsText(4)
BODU_MEZI_PF=int(BODU_MEZI_PF)
#BODU_MEZI_PF = 30

print "start"


#nacteni profilu zacatek--------------------------------------
def nacti_PF(inShapefile):
    cursor = arcpy.SearchCursor(inShapefile)
    a=0
    for row in cursor:
        a=a+1
    field = "y"
    y = "POINT_Y"
    x = "POINT_X"
    z = "zetko"
    #pole je indexovano od nuly jako v C++
    #matice je radek, sloupec
    matrix = [[0]*3 for i in range(a)]
    i=0
    cursor = arcpy.SearchCursor(inShapefile)
    for row in cursor:
        matrix [i][0]=row.getValue(y)
        matrix [i][1]=row.getValue(x)
        matrix [i][2]=row.getValue(z)
        matrix [i][2]=float(matrix[i][2])# konverze string na float
        i=i+1
    return matrix
#nacteni profilu konec--------------------------------------

#uhel zacatek--------------------------------------
def uhel(bod1,bod2):
    lon1=bod1[1]
    lat1=bod1[0]
    lon2=bod2[1]
    lat2=bod2[0]
    lon1RAD = lon1*math.pi/180
    lat1RAD = lat1*math.pi/180
    lon2RAD = lon2*math.pi/180
    lat2RAD = lat2*math.pi/180
    uhel = math.atan2(math.sin(lon2RAD-lon1RAD)*math.cos(lat2RAD),math.cos(lat1RAD)*math.sin(lat2RAD)-math.sin(lat1RAD)*math.cos(lat2RAD)*math.cos(lon2RAD-lon1RAD))
    uhel=uhel*180/math.pi
    
    return uhel
#uhel konec--------------------------------------

#vzdalenost zacatek--------------------------------------
def vzdalenost(bod1,bod2):
    lon1=bod1[1]
    lat1=bod1[0]
    lon2=bod2[1]
    lat2=bod2[0]
    lon1RAD = lon1*math.pi/180
    lat1RAD = lat1*math.pi/180
    lon2RAD = lon2*math.pi/180
    lat2RAD = lat2*math.pi/180
    delka = math.acos(math.sin(lat1RAD)*math.sin(lat2RAD)+math.cos(lat1RAD)*math.cos(lat2RAD)*math.cos(lon2RAD-lon1RAD))*6378.137
    return delka #vysledek je v kilometrech
#vzdalenost konec--------------------------------------

#dopocet bodu zacatek--------------------------------------
def dopocet_bodu(bod,uhel,vzdalenost_KM):
    vzdalenost=vzdalenost_KM
    lon0=bod[1]
    lat0=bod[0]
    R=6378.137
    uhelRAD = uhel*math.pi/180;
    lon0RAD = lon0*math.pi/180;
    lat0RAD = lat0*math.pi/180;
    lat1RAD=math.asin(math.sin(lat0RAD)*math.cos(vzdalenost/R) + math.cos(lat0RAD)*math.sin(vzdalenost/R)*math.cos(uhelRAD))
    lon1RAD=lon0RAD + math.atan2(math.sin(uhelRAD)*math.sin(vzdalenost/R)*math.cos(lat0RAD),math.cos(vzdalenost/R)-math.sin(lat0RAD)*math.sin(lat1RAD))
    lon1=lon1RAD*180/math.pi
    lat1=lat1RAD*180/math.pi
    novy=[lat1,lon1]
    return novy
#dopocet bodu konec--------------------------------------

#nalezeni vysky bodu zacatek--------------------------------------
def vyska_nejblizsi(stary_profil,bod):
    vzdal_nejbliz = 100#vzdalenost nejblizsiho stareho od hledaneho
    zacatek = [0,0]
    zacatek[0]=stary_profil[0][0]
    zacatek[1]=stary_profil[0][1]
    for i in range(int(len(stary_profil))):
        bod_stary=[0,0]
        bod_stary[0]=stary_profil[i][0]
        bod_stary[1]=stary_profil[i][1]
        vzdal=vzdalenost(bod_stary,bod)*1000
        if vzdal<vzdal_nejbliz:
            vzdal_nejbliz=vzdal
            nejblizsi=bod_stary
            x1=i#index nejblizsiho bodu
    vzdal_bod=vzdalenost(zacatek,bod)*1000#vzdalenost hledaneho bodu od zacatku
    vzdal_nejbliz_poc=vzdalenost(zacatek,nejblizsi)*1000#vzdalenost nejblizsiho bodu od zacatku
    if vzdal_bod<vzdal_nejbliz_poc:#index souseda od nejblizsiho, mezi nimiz lezi bovy bod
        x2=x1
        x1=x1-1
    else:
        x2=x1+1
    levy=[0,0,0]
    levy[0]=stary_profil[x1][0]
    levy[1]=stary_profil[x1][1]
    levy[2]=stary_profil[x1][2]
    pravy=[0,0,0]
    pravy[0]=stary_profil[x2][0]
    pravy[1]=stary_profil[x2][1]
    pravy[2]=stary_profil[x2][2]
    vzdal_soused=vzdalenost(levy,pravy)*1000
    vzdal_levy_bod=vzdalenost(levy,bod)*1000
    prevyseni_soused=abs(stary_profil[x1][2]-stary_profil[x2][2])
    if levy[2]>pravy[2]:
        vyska=levy[2]-vzdal_levy_bod*(prevyseni_soused/vzdal_soused)
    else:
        vyska=levy[2]+vzdal_levy_bod*(prevyseni_soused/vzdal_soused)
    return vyska
#nalezeni vysky bodu konec--------------------------------------

#predelani profilu zacatek--------------------------------------
def resample_PF(stary,bodu_v_PF):
    stary_bodu = len(stary)
    zacatek = [0,0,0]
    zacatek[0]=stary[0][0]
    zacatek[1]=stary[0][1]
    zacatek[2]=stary[0][2]
    konec= [0,0,0]
    konec [0]=stary[stary_bodu-1][0]
    konec [1]=stary[stary_bodu-1][1]
    konec [2]=stary[0][2]
    sirka = vzdalenost(zacatek,konec)*1000
    rozestup = sirka/(bodu_v_PF-1)
    profil = [[0]*3 for i in range(int(bodu_v_PF))]
    profil[0][0]=zacatek[0]
    profil[0][1]=zacatek[1]
    profil[0][2]=zacatek[2]
    profil[bodu_v_PF-1][0]=konec[0]
    profil[bodu_v_PF-1][1]=konec[1]
    profil[bodu_v_PF-1][2]=konec[2]
    for i in range(int(bodu_v_PF-2)):#dopocet souradnic pro vnitrni body
        smer=uhel(zacatek,konec)
        pomocna=dopocet_bodu(zacatek,smer,((i+1)*rozestup/1000))
        profil[i+1][0]=pomocna[0]
        profil[i+1][1]=pomocna[1]
        profil[i+1][2]=vyska_nejblizsi(stary,pomocna)
    return profil
#predelani profilu konec--------------------------------------

#prevraceni profilu zacatek--------------------------------------
def prevrat_PF(profil1, profil2):
    a = len(profil2)
    b=3
    c = len(profil1)
    X1 = profil1[0][0]
    Y1=profil1[0][1]
    X2 = profil1[c-1][0]
    Y2=profil1[c-1][1]
    X3 = profil2[0][0]
    Y3=profil2[0][1]
    X4 = profil2[a-1][0]
    Y4=profil2[a-1][1]
    y = ((Y1-Y3)*(Y2-Y4)*(X1-X2)+(X2-X4)*(Y1-Y3)*Y2-(Y2-Y4)*(X1-X3)*Y1)/((Y1-Y3)*(X2-X4)-(Y2-Y4)*(X1-X3))
    x=X1+((X1-X3)/(Y1-Y3))*(y-Y1)
    if (((X1<=x and x<=X3)or(X3<=x and x<=X1))and((X2<=x and x<=X4)or(X4<=x and x<=X2))and((Y1<=y and y<=Y3)or(Y3<=y and y<=Y1))and((Y2<=y and y<=Y4)or(Y4<=y and y<=Y2))):
        A=math.fmod(a,2)
        B=(a/2)-A
        if A==0:
            B=B-1
        for i in range(int(B)):
            for j in range(int(b)):
                BOX=profil2[i][j]
                profil2[i][j]= profil2[a-i-1][j]
                profil2[a-i-1][j]= BOX      
    return profil2
#prevraceni profilu konec--------------------------------------

#cast segmentu zacatek--------------------------------------
def cast_segmentu(bod1,bod2,mezi_profily):
    sirka = vzdalenost(bod1,bod2)
    uhell = uhel(bod1,bod2)
    rozestup = sirka/(mezi_profily+1)
    a=mezi_profily +2
    segment = [[0]*3 for i in range(int(a))]
    segment[0][0]=bod1[0]
    segment[0][1]=bod1[1]
    segment[0][2]=bod1[2]
    segment[a-1][0]=bod2[0]
    segment[a-1][1]=bod2[1]
    segment[a-1][2]=bod2[2]
    pomocna=[0,0]
    for i in range(int(mezi_profily)):
        pomocna=dopocet_bodu(bod1,uhell,(rozestup*(i+1))); 
        segment[i+1][0]=pomocna[0]
        segment[i+1][1]=pomocna[1]
        segment[i+1][2]=bod1[2]-((rozestup*i)*((bod1[2]-bod2[2])/sirka))
    return segment
#cast segmentu konec--------------------------------------

#cely segment zacatek--------------------------------------
def cely_segment(profil1,profil2,bodu_v_PF,mezi_profily):
    profil1=resample_PF(profil1,bodu_v_PF)
    profil2=resample_PF(profil2,bodu_v_PF)
    a = (mezi_profily+2)*(bodu_v_PF);
    segment = [[0]*3 for l in range(int(a))]
    i=0
    for b in range(int(bodu_v_PF)):
        bod1=[0,0,0]
        bod1[0]=profil1[b][0]
        bod1[1]=profil1[b][1]
        bod1[2]=profil1[b][2]
        bod2=[0,0,0]
        bod2[0]=profil2[b][0]
        bod2[1]=profil2[b][1]
        bod2[2]=profil2[b][2]
        pomocny_profil=[[0]*3 for m in range(int(mezi_profily+2))]
        pomocny_profil=cast_segmentu(bod1,bod2,mezi_profily)
        j=0;
        for k in range(int(mezi_profily+2)):
            segment[i][0]=pomocny_profil[k][0]
            segment[i][1]=pomocny_profil[k][1]
            segment[i][2]=pomocny_profil[k][2]
            #segment[i][0]=pomocny_profil[j][0]
            #segment[i][1]=pomocny_profil[j][1]
            #segment[i][2]=pomocny_profil[j][2]
            j=j+1
            i=i+1
    return segment
#cely segment konec--------------------------------------


items=len(LIST)
segmentu = items-1
eee=[[0]*3 for m in range(int(BODU_V_PF*(BODU_MEZI_PF+2)*segmentu))]
k=0
for i in range(int(items-1)):
    profil1=nacti_PF(LIST[i])
    profil2=nacti_PF(LIST[i+1])
    profil2=prevrat_PF(profil1, profil2)
    ddd=cely_segment(profil1,profil2,BODU_V_PF,BODU_MEZI_PF)
    for j in range(int(BODU_V_PF*(BODU_MEZI_PF+2))):
    #for j in range(int(BODU_V_PF*(BODU_MEZI_PF+2)-1)):
        eee[k][0]=ddd[j][0]
        eee[k][1]=ddd[j][1]
        eee[k][2]=ddd[j][2]
        k=k+1
        
    
    



#vytvoreni a ulozeni *.SHP zacatek--------------------------------------
arcpy.env.workspace = path
#coordinate system WGS_84 (zadan kodem)
sr = arcpy.SpatialReference(4326)
arcpy.CreateFeatureclass_management(path,outSHPpoints,"POINT", has_z="ENABLED",spatial_reference=sr)
arcpy.AddField_management(outSHPpoints, "vyska", "float", "", "", "", "", "", "")

cur = arcpy.InsertCursor(path+"\\"+outSHPpoints) 
pnt = arcpy.Point()
      
f=len(eee)

for d in range(int(f)):
    pnt.X = eee[d][1]
    pnt.Y = eee[d][0]
    pnt.Z = eee[d][2]
    feat = cur.newRow()
    feat.setValue("SHAPE", arcpy.Point(pnt.X, pnt.Y, pnt.Z))
    feat.setValue("vyska",eee[d][2])
    cur.insertRow(feat)
    cur.updateRow  
#vytvoreni a ulozeni *.SHP konec--------------------------------------



arcpy.AddMessage("kral je mrtev, at zije kral")
#arcpy.AddError("kral je mrtev, at zije kral")


print "hotovo"



