# -*- coding: utf-8 -*-
"""
Created on Mon january 28 2021

@author: Adrian Ringenbach
#slightly updated by Marlene Kronenberg, July 2022

"""

# Stem outlines with given WindDirection 


import math
import numpy as np
import fileinput


import rasterio
from rasterio import rio
from rasterio.rio import sample
from rasterio.rio.sample import sample
from rasterio.enums import Resampling
from rasterio.plot import show

import pandas as pd 
from random import *


#Detection of the PivotPoint
def PivotPoint(x0, y0, z0, dbh, windDir, rootPlate):              #future adaption --> z_out world 2pixel(x_out, y_out)
    if rootPlate == 'y' :
        rpFactor = 3
    else:
        rpFactor = 1.3
        
    x_out = x0+rpFactor*(dbh/(2*100))*-math.sin(np.radians(windDir))
    y_out = y0+rpFactor*(dbh/(2*100))*-math.cos(np.radians(windDir))
    z_out = z0

    return  (x_out, y_out, z_out, rpFactor)   


def StemBase(x_out, y_out, z_out, winDir, angles, rpFactor, dbh): #input: XYZ from Pivot Point, etc...
   length = len(angles)  #number Steps will be Calculated (from Vertical in 5° steps till -90° = 36)
   #StemBaseOut=np.zeros((length, 3))
   StemBaseX=np.zeros((length,1))
   StemBaseY=np.zeros((length,1))
   StemBaseZ=np.zeros((length,1))
   i=0
   
   for theta in angles:#Neigungswinkel theta von 0 - 180 Grad in 5° schritten erhöhen.
        theta=theta+270   #adaption for horizontal start (vs. vertical start of TreeTop fct)
        if theta>=360:
            theta=theta-360     
        xtemp=x_out+(rpFactor*dbh/100*math.sin(np.radians(theta))*math.sin(np.radians(winDir)))
        ytemp=y_out+(rpFactor*dbh/100*math.sin(np.radians(theta))*math.cos(np.radians(winDir)))
        ztemp= z_out+rpFactor*dbh/100*math.cos(np.radians(theta))
        StemBaseX[i]= xtemp
        StemBaseY[i]= ytemp
        StemBaseZ[i]= ztemp
        #StemBaseOut[i, 0]= xtemp
        #StemBaseOut[i, 1]= ytemp
        #StemBaseOut[i, 2]= ztemp
        i=i+1
            
   return (StemBaseX, StemBaseY, StemBaseZ)

          

#Detection of the Treetop pOSITION  during Fall (e.g. every 5° )
def TreeTop(x0, y0, z0, windDir, H, angles):
   length = len(angles)
   #number Steps will be Calculated (from Vertical in 5° steps till -90° = 36)
   out=np.zeros((length, 3))
   X=np.zeros((length,1))
   Y=np.zeros((length,1))
   Z=np.zeros((length,1))
   i=0
   
   for theta in angles:#Neigungswinkel theta von 0 - 180 Grad in 5° schritten erhöhen.
        while windDir>=360:
            windDir=windDir-360
        while windDir<0:
            windDir=windDir+360
        x_out=x0[i]+(H*math.sin(np.radians(theta))*-math.sin(np.radians(windDir)))
        y_out=y0[i]+(H*math.sin(np.radians(theta))*-math.cos(np.radians(windDir)))
        z_out= z0[i]+H*math.cos(np.radians(theta))
        X[i]= x_out
        Y[i]= y_out
        Z[i]= z_out
        out[i, 0]= x_out
        out[i, 1]= y_out
        out[i, 2]= z_out
        i=i+1
            
   return (out, X, Y, Z)

###Calc location of Lower Stem Base betweem StemBase and PivotPoint in case of Rootplate=y
def lowerStBase(x_out, y_out, z_out, rpFactor, StemBaseX, StemBaseY, StemBaseZ, dbh, angles, windDir):
    length = len(angles)
       #number Steps will be Calculated (from Vertical in 5° steps till -90° = 36)
       #out=np.zeros((length, 3))
    lowerStBaseX=np.zeros((length,1))
    lowerStBaseY=np.zeros((length,1))
    lowerStBaseZ=np.zeros((length,1))
    i=0
      
    for theta in angles:#Neigungswinkel theta von 0 - 180 Grad in 5° schritten erhöhen.
         theta=theta+270   #adaption for horizontal start (vs. vertical start of TreeTop fct)
         if theta>=360:
            theta=theta-360     
         lowerStBaseX_out=x_out+(1.7*dbh/100*math.sin(np.radians(theta))*-math.sin(np.radians(windDir)))
         lowerStBaseY_out=y_out+(1.7*dbh/100*math.sin(np.radians(theta))*-math.cos(np.radians(windDir)))
         lowerStBaseZ_out= z_out+1.7*dbh/100*math.cos(np.radians(theta))
         lowerStBaseX[i]= lowerStBaseX_out
         lowerStBaseY[i]= lowerStBaseY_out
         lowerStBaseZ[i]= lowerStBaseZ_out
         i=i+1
                
    return (lowerStBaseX, lowerStBaseY, lowerStBaseZ)





#Stem-outline  From PivotPoint To TreeTop (solved as a 2D Problem: lower edge (unterkante Baum))
def SO_Pivo_Top(Ttop_X, Ttop_Y, Ttop_Z, Pivot_X, Pivot_Y, Pivot_Z, NrPts):
    length=len(Ttop_X)
    SO_x=np.zeros((length,NrPts))
    SO_y=np.zeros((length,NrPts))
    SO_z=np.zeros((length, NrPts))
    for i in range(0,length):
        SO_x[i]= np.transpose(np.linspace(Pivot_X[i], Ttop_X[i], NrPts))
        SO_y[i]= np.transpose(np.linspace(Pivot_Y[i], Ttop_Y[i], NrPts))
        SO_z[i]= np.transpose(np.linspace(Pivot_Z[i], Ttop_Z[i], NrPts))
    
    return (SO_x,SO_y, SO_z)




#https://stackoverflow.com/questions/33596416/sample-raster-cell-value-at-coordinates:
def world2Pixel_old(geoT, x, y, dem_res): #gt 
  ulX = geoT[2] #gt[0]
  ulY = geoT[5]#gt[3]
  xDist = geoT[0]#gt[1]
  yDist = geoT[4]#gt[5]
  rtnX = geoT[1]#gt[2]
  rtnY = geoT[3]#gt[4]
  line = int((x - ulX) / xDist)                 # not as the weblink..  line and pixel inverted, as wrong resultes (Dw. in the air or below ground)
  pixel = int((ulY - y) / (-yDist))             # not as the weblink..  line and pixel inverted, as wrong resultes (Dw. in the air or below ground)
  
  return (pixel, line)   
 


def world2Pixel_new(geoT, x, y, dem_res): #gt          Test_ fehler suche
  #src = gdal.Open(input_dhm); 
  #ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform() #ulx, xres, xskew, uly, yskew, yres 
  ulX = geoT[2]     #gt[0]
  ulY = geoT[5]#gt[3]
  xDist = geoT[0]#gt[1]
  yDist = geoT[4]#gt[5]
  rtnX = geoT[1]#gt[2]
  rtnY = geoT[3]#gt[4]
  line =int((x - ulX) / xDist)         # not as the weblink..  line and pixel inverted, as wrong resultes (Dw. in the air or below ground)
  pixel = int(( ulY-y) / (-yDist))     # not as the weblink..  line and pixel inverted, as wrong resultes (Dw. in the air or below ground)
  value=dem_res[pixel+1, line+1]      
  return (value)



def getMaxAngle(angles, Npoints, dem, geoT, SO_x, SO_y, SO_z, min_height):
    positive=0
    negative=0
    maxAngle=180 
  
    for i in range(0,len(angles)): #[18, 20]: #range(0,len(X)):
        if negative == 0:
            for j in range(int(Npoints/10),int(8*Npoints/10) ):                         #lowest 10% are not looked at, as the tree may overturn, and top 20 are neglected as it is too flexible)
                #pixel, line = world2Pixel_new(geoT, SO_x[i,j], SO_y[i,j])
                value= world2Pixel_new(geoT, SO_x[i,j], SO_y[i,j], dem)
                #band = dhm_gdal.GetRasterBand(1)
                #value=dem[pixel+1, line+1]
                #value = band.ReadAsArray(pixel, line, 1, 1)[0, 0]
                #print(z_valueArry,value )
                if value <  SO_z[i,j]-min_height:
                    positive=positive+1
                else:
                    if maxAngle>angles[i]:
                        maxAngle=angles[i]
                    negative=negative+1
                    #print('max. Angle is ', angles[i])
    
    return(maxAngle) #angles[i],


def getAnglePivotStem(a, b, c):                 
    a_np = np.array([a[0], a[1],a[2]])
    b_np = np.array([b[0], b[1],b[2]])
    c_np = np.array([c[0], c[1],c[2]])

    ba = a_np - b_np
    bc = c_np - b_np

    cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))
    angle = np.degrees(np.arccos(cosine_angle))
    return (angle)         



def getPivotOpp(maxAngle, winDir, dbh, rpFactor, x_pivot,y_pivot,z_pivot):
    #if maxAngle>90:
      #  Dir=180-winDir
    #else:
    Dir=winDir
    P_opp_x = x_pivot+((dbh/100.0*rpFactor) * math.cos(np.radians(maxAngle))*math.sin(np.radians(Dir)))
    P_opp_y = y_pivot+((dbh/100.0*rpFactor) * math.cos(np.radians(maxAngle))*-math.cos(np.radians(Dir)))
    P_opp_z = z_pivot+(dbh/100.0*rpFactor) * math.sin(np.radians(maxAngle))
    
    return(P_opp_x, P_opp_y, P_opp_z)     
#upper Stem-outline  From opposite PivotPoint To TreeTop (solved as a 2D Problem: upper  edge (unterkante Baum))

def finTtop_cord(maxAngle, angles, X, Y, Z):
    idx=0
    for i in range (0, len(angles)):
        if angles[i] == maxAngle:
            idx=i
    finTtopX = X[idx]
    finTtopY = Y[idx]
    finTtopZ = Z[idx]
    return(finTtopX, finTtopY, finTtopZ, idx)
   
def upSO_PiOp_Top(finTtopX, finTtopY, finTtopZ, P_opp_x, P_opp_y, P_opp_z, Npoints):     
    length=len(finTtopX)
    uSO_x= np.transpose(np.linspace(P_opp_x, finTtopX, Npoints))
    uSO_y= np.transpose(np.linspace(P_opp_y, finTtopY, Npoints))
    uSO_z= np.transpose(np.linspace(P_opp_z, finTtopZ, Npoints))
    
    return (uSO_x,uSO_y, uSO_z)


def resample(input_dhm, upscale_factor):
    #upscale_factor=10 #20cm resolution

    with rasterio.open(input_dhm) as dataset:
    
        # resample data to target shape
        data = dataset.read(
            out_shape=(
                dataset.count,
                int(dataset.height * upscale_factor),
                int(dataset.width * upscale_factor)
            ),
            resampling=Resampling.bilinear
        )
        # scale image transform
        transform = dataset.transform * dataset.transform.scale(
            (dataset.width / data.shape[-1]),
            (dataset.height / data.shape[-2])
        )
        data=data[0, :, :]                #necessary to get the correct dem_res shape: (x, y) instead of (1,x, y)
    return(data, transform)



def calcTreeAxis(SO_x,SO_y, SO_z, uSO_x, uSO_y, uSO_z, winDir, idx):   #middle between upper and lower Stem-outline
    zAx=np.zeros(len(uSO_z[0,:]))
    xAx=np.zeros(len(uSO_z[0,:]))
    yAx=np.zeros(len(uSO_z[0,:]))

    for i in range(0, len(zAx)):
        zAx[i]=(SO_z[idx,:][i]+uSO_z[0,i])/2
        xAx[i]=(SO_x[idx,:][i]+uSO_x[0,i])/2
        yAx[i]=(SO_y[idx,:][i]+uSO_y[0,i])/2
    
    
    
    return(xAx, yAx, zAx)
    
def normalize(v):
    norm = np.linalg.norm(v)
    if norm == 0: 
       return v
    return v / norm   
    
#CALCULATE STEM DISCS based on the normal vector (= stem Axis, v3) and the DBH), 
#based on a new local coordinate system with V1, v2, v3, 
#(https://stackoverflow.com/questions/27714014/3d-point-on-circumference-of-a-circle-with-a-center-radius-and-normal-vector)
    
#Step1: calculate normal vector and normalize its length: 
def VecSys(xAx, yAx, zAx, finTtopX, finTtopY, finTtopZ):    
    v3=[xAx[0]-finTtopX, yAx[0]-finTtopY, zAx[0]-finTtopZ]
    v3_n=normalize(v3)
    #new "X-vector"
    v1=[zAx[0]-finTtopZ, [0],(-xAx[0]+finTtopX) ]
    v1_n=normalize(v1)
    #new "y-vector"
    v2=np.cross(v3_n[:,0], v1_n[:,0])
    v2_n= normalize(v2)
    
    return(v1_n, v2_n, v3_n)
    
    
def dist3D(x1, y1, z1, x2, y2, z2):
    dist=((x1-x2)**2+(y1-y2)**2+(z1-z2)**2)**0.5
    return(dist)
    
def StemDisc(v1_n, v2_n, v3_n, xAx, yAx, zAx, r):    
    a=range(0,360,30) #Nr. of points per StemDisc --> 360/40 = 9 pts
    
    circ_x=np.zeros(len(a))
    circ_y=np.zeros(len(a))
    circ_z=np.zeros(len(a))
    
    for i in range(0,len(a)):
        circ_x[i] = xAx[0] + r * (v1_n[0][0] * math.cos(np.radians(a[i])) + v2_n[0] * math.sin(np.radians(a[i])))
        circ_y[i] = yAx[0] + r * (v1_n[1][0] * math.cos(np.radians(a[i])) + v2_n[1] * math.sin(np.radians(a[i])))
        circ_z[i] = zAx[0] + r * (v1_n[2][0] * math.cos(np.radians(a[i])) + v2_n[2] * math.sin(np.radians(a[i])))

    return(circ_x, circ_y, circ_z)


#    test_radius=np.zeros(9)
#    for i in range(0, 9):
#        test_radius[i]=((xAx[0]-circ_x[i])**2+(yAx[0]-circ_y[i])**2+(zAx[0]-circ_z[i])**2)**0.5
#    
#    circ_x.mean()
#    circ_y.mean()
#    circ_z.mean()
    
    
def energyAbs(dbh, age):     #diameter and 0/1 integer if Fomitopsis pinicola (=fungi) is present
    biw_fresh=238890 # [j /m2] break impact work without for fresh Spruce logs
    biw = 131890 # [j /m2] break impact work without fomitopsis pinicola for 8-10 year old spruce logs (M. Ammann, 2006)
    biw_fp = 25590  # [j /m2] break impact work with fomitopsis pinicola for 8-10 year old spruce logs(M. Ammann, 2006)    
    
    if age==0:
        eA=     round(math.pi*((dbh/100)/2)**2*biw_fresh,0) 
    if age==1:
        eA=     round(math.pi*((dbh/100)/2)**2*biw,0)           #energy Absobed [j]
    if age==2:                              
        eA=     round(math.pi*((dbh/100)/2)**2*biw_fp,0)        #energy Absobed [j]
    return (eA)    



def fp_prob(distrib,mu):              #Fomitopsis Pinicola: wood Rotting -fungi #gaus verteilung mit Mu =0  0.1 0.2 0.3 ....   simga < 0.5
    if distrib=='random':
        fp_prob=randint(1, 2)
    if distrib == 'binominal':        
        fp_prob=np.random.binomial(1,mu)
    return(fp_prob)


#____________________________________________   

def DWG(input_dem, treefile, out_folder, state, mainWindDir, sigmaWinDir, rp_percentage, upscaleFact, WindThorwPercentage, DBH_thrs):
    
    tf=pd.DataFrame()
    tf = pd.read_csv(treefile, header=None,  sep=r"\s+", usecols=[0,1,2,3], names=['x', 'y', 'H', 'BHD'])
    
    dem_res, geoT= resample(input_dem, upscaleFact) # here, inputDEM has alread 20 cm res #Upscalefactor = 10: 2m Res --> 20 cm Res
    
    
    angles=range(0,180,5)
    Npoints=100
    
    
    winDir=np.zeros([len(tf),1])
    rp_state=np.zeros([len(tf),1])
    rp_counter = 0
    wt_state=np.zeros([len(tf),1]) #Windthrow Status: 0-> Standing 1 lying
    
    
    np.random.seed(1) #make the random binomial distribution static, so that the same Trees resist a windthrow between different scenrios  
    wt_state=np.random.binomial(1,WindThorwPercentage,size=len(tf))
    
    np.random.seed(1) #np.random.seed(1) # to make the wind direction static within the fresh and 10 yr results
    winDir=np.random.normal(mainWindDir, sigmaWinDir, size=(len(tf)))



    for i in range(0,len(tf)):   
    #winDir[i]=np.random.normal(mainWindDir, sigmaWinDir, size=(len(tf)))[i] 
        while winDir[i]>=360:                               #adapt all fututre calulation (cos, rand, etc.) to the range of 0-360° --> Fewer errors
            winDir[i]=winDir[i]-360
        while winDir[i]<0:
            winDir[i]=winDir[i]+360
            #np.random.randint(-20,high=21, )[i] 
          
        rp_state[i]=np.random.binomial(1,rp_percentage,size=len(tf))[i] 
    
    
    if state== 'fresh':  
        min_height = 2.15/3.33
        rp_dia_fact = 3
    if state== 'old':
        min_height = 0.86/3.33 
        rp_dia_fact = 2.2
    
    dw_list_path=np.zeros(int(len(tf)+sum(rp_state)),dtype = object)
    dw_list_eA= np.zeros([int(len(tf)+sum(rp_state)),1])
    
    #dw_complete= np.zeros([len(tf),2])
    
    dw_complete = [[0] * 2 for i in range(int(len(tf)+sum(rp_state)))]
    #dw_complete_10yr = [[0] * 2 for i in range(2*len(tf))]
            
    dw_complete_cmd_line = [[0] * 4 for i in range(int(len(tf)+sum(rp_state)))]        
            
            
    for i in range(0,len(wt_state)):
        if i == 0 or i%50==0:
            print('processing dw log ' +str(i)) # print for every 50th log a status update
   
        if wt_state[i]== 0:
    #print(wt_state[i])
            continue
        elif tf.iloc[i,3]< DBH_thrs: #skip tree if DBH is not > 8 cm
            continue
        else:
        
            tree=tf.iloc[i,:]
            H=tf.iloc[i,2]
            dbh=tf.iloc[i,3]*100
        
            x0=tf.iloc[i,0]
            y0=tf.iloc[i,1]
            z0=world2Pixel_new(geoT, x0, y0, dem_res)
            
            
        
            if rp_state[i]==0:
                rp='n'
            else:
                rp='y'
            #winDir=round(uniform(mainWindDir-20, mainWindDir+20),0)   #old aproach, windDir just in the for loop determined. Disadvantage: no seeds possible
            
            #defining the x, y, z, Coordinates of the Pivot Point
            x_pivot, y_pivot, z_pivot, rpFactor = PivotPoint(x0, y0, z0, dbh,winDir[i], rp)
            z_pivot=world2Pixel_new(geoT, x_pivot, y_pivot, dem_res)+0.1
             
            StemBaseX, StemBaseY, StemBaseZ=StemBase(x_pivot, y_pivot, z_pivot, winDir[i], angles, rpFactor, dbh)
            lowerStBaseX, lowerStBaseY, lowerStBaseZ= lowerStBase(x_pivot, y_pivot, z_pivot, rpFactor, StemBaseX, StemBaseY, StemBaseZ,  dbh, angles, winDir[i])
                  
    #        fig = plt.figure()
    #        ax = plt.axes(projection="3d")
    #        ax.set_aspect('equal', 'box')
    #        #ax.scatter3D( circ_x,circ_y,circ_z)
    #        #ax.scatter3D( rp_x,rp_y,rp_z)
    #        ax.scatter3D( x0,y0,z0)
    #        ax.scatter3D( x_pivot, y_pivot, z_pivot)
    #        ax.scatter3D( StemBaseX, StemBaseY, StemBaseZ)
    #        ax.scatter3D( lowerStBaseX, lowerStBaseY, lowerStBaseZ)  
    #        ax.scatter3D( X, Y, Z)
            
            ##Save circ_x,circ_y,circ_z and finTtopX, finTtopY, finTtopZ into pts file and export it
            # new array for pts file 
                
                    
            out, X, Y, Z= TreeTop(StemBaseX, StemBaseY, StemBaseZ, winDir[i], H, angles)
          
            #defining the coordinates of the lower Stem Outline 
            SO_x,SO_y, SO_z =SO_Pivo_Top(X, Y, Z, lowerStBaseX, lowerStBaseY, lowerStBaseZ,Npoints)
            
            #checking the max Angle 
            maxAngle=getMaxAngle(angles, Npoints, dem_res, geoT, SO_x, SO_y, SO_z, min_height)
            #print(maxAngle)
            #claculating the angle of the TreeConus
            beta=getAnglePivotStem([x0, y0, z0], [SO_x[0,0],SO_y[0,0], SO_z[0,0]],[SO_x[0,5],SO_y[0,5], SO_z[0,5]] )
            P_opp_x, P_opp_y, P_opp_z = getPivotOpp(maxAngle, winDir[i], dbh, rpFactor, x_pivot,y_pivot,z_pivot)
            finTtopX, finTtopY, finTtopZ, idx=finTtop_cord(maxAngle, angles, X, Y, Z)
            uSO_x,uSO_y, uSO_z= upSO_PiOp_Top(finTtopX, finTtopY, finTtopZ, P_opp_x, P_opp_y, P_opp_z, Npoints)
        #    fig = plt.figure()
        #    ax = plt.axes(projection="3d")
        #    ax.scatter3D(SO_x[idx,:], SO_y[idx,:], SO_z[idx,:]);
        #    ax.scatter3D(uSO_x[0,:],uSO_y[0,:], uSO_z[0,:])
        
            for j in range(0,Npoints):
                data_x, data_y=world2Pixel_old(geoT, uSO_x[0,j], uSO_y[0,j], dem_res)
                dem_res[data_x,data_y ] = uSO_z[0,j]+min_height    #+0.5 due to branches etc...   
          
            
            xAx,yAx,zAx=calcTreeAxis(SO_x,SO_y, SO_z, uSO_x, uSO_y, uSO_z, winDir[i],idx)
            v1, v2, v3=VecSys(xAx, yAx, zAx, finTtopX, finTtopY, finTtopZ)
            
            circ_x, circ_y, circ_z=StemDisc(v1, v2, v3, xAx, yAx, zAx, (dbh*1.3/(2*100)))
            
            
            #Create Rootplate out of two disc. Disc1 = efold of stemfeet, disc2 in distance of 1 Dbh away   
            if rp_state[i]==1:
                rp_x, rp_y, rp_z=StemDisc(v1, v2, v3, xAx, yAx, zAx, (dbh*3/(2*100)))
                rp_x2, rp_y2, rp_z2=StemDisc(v1, v2, -v3, xAx[0]+(dbh/100)*v3[0], yAx+(dbh/100)*v3[1], zAx+(dbh/100)*v3[2], (dbh*rp_dia_fact/(2*100)))
                
                for j in range(0,len(rp_x)):
                    data_x, data_y=world2Pixel_old(geoT, rp_x[j], rp_y[j], dem_res)
                    dem_res[data_x,data_y ] = rp_z[j]
                    data_x, data_y=world2Pixel_old(geoT, rp_x2[j], rp_y2[j], dem_res)
                    dem_res[data_x,data_y ] = rp_z2[j]
            
            
        #    fig = plt.figure()
        #    ax = plt.axes(projection="3d")
        #    ax.set_aspect('equal', 'box')
        #    ax.scatter3D( circ_x,circ_y,circ_z)
        #    ax.scatter3D( rp_x,rp_y,rp_z)
        #    ax.scatter3D( rp_x2,rp_y2,rp_z2)
        
            ##Save circ_x,circ_y,circ_z and finTtopX, finTtopY, finTtopZ into pts file and export it
            # new array for pts file 
            shape_export=(circ_x.shape[0]+1, 3) #length circ + treetop by 3 (x,y,z coordinates)
            pts_dw_exp=np.zeros(shape_export)

            for k in range(0, len(circ_x)):
                pts_dw_exp[k,0]=circ_x[k]-finTtopX
                pts_dw_exp[k,1]=circ_y[k]-finTtopY
                pts_dw_exp[k,2]=circ_z[k]-finTtopZ
                
            pts_dw_exp[circ_x.shape[0],0]=0
            pts_dw_exp[circ_x.shape[0],1]=0
            pts_dw_exp[circ_x.shape[0],2]=0
            
            #write pts coodinates 
            dw_i=str(i)
            filename_pts='dw_'+dw_i+'.pts'
            np.savetxt(out_folder+'\\'+filename_pts, pts_dw_exp, fmt='%f', delimiter=' ', newline='\n')
            #print("5pts file exported")
            
            
            
            #Save rootplates (rp1 and rp2), as relative coordinates from Roottop
            #As 2 discs are saved, length of shape_rp_export = 2 *len(shape)
            if rp_state[i]==1:
                shape_rp_export=(2*circ_x.shape[0], 3) #length circ + treetop by 3 (x,y,z coordinates)
                pts_rp_exp=np.zeros(shape_rp_export)
            
                for k in range(0, len(circ_x)):
                    pts_rp_exp[k,0]=rp_x[k]-finTtopX
                    pts_rp_exp[k,1]=rp_y[k]-finTtopY
                    pts_rp_exp[k,2]=rp_z[k]-finTtopZ
                    pts_rp_exp[len(circ_x)+k,0]=rp_x2[k]-finTtopX
                    pts_rp_exp[len(circ_x)+k,1]=rp_y2[k]-finTtopY
                    pts_rp_exp[len(circ_x)+k,2]=rp_z2[k]-finTtopZ
            
                
                #calculate pts coodinates 
                filename_rp_pts='rp_'+dw_i+'.pts'
                np.savetxt(out_folder+'\\'+filename_rp_pts, pts_rp_exp, fmt='%f', delimiter=' ', newline='\n')
            
       
            #calculate centroid point (--> Treetop agiert als "Centroid")
            txt_exp=np.zeros([1,3])
            
            txt_exp[0,0]=finTtopX
            txt_exp[0,1]=finTtopY
            txt_exp[0,2]=finTtopZ   
                
 #           #Export routine: centroid point
 #           #dw_i=str(i)
 #           filename_txt='dw_'+dw_i+'.txt'
 #           np.savetxt(out_folder+'\\'+filename_txt, txt_exp, fmt='%f', delimiter=' ', newline='\n')
            
 #           if rp_state[i]==1:
 #               filename_rp_txt='rp_'+dw_i+'.txt'
 #               np.savetxt(out_folder+'\\'+filename_rp_txt, txt_exp, fmt='%f', delimiter=' ', newline='\n')
        
            #print("txt file exported")
            
            #entry in dw_list
            dw_list_path[i] = out_folder+'\\'+filename_pts
            
            if rp_state[i]==1:
                dw_list_path[int(len(tf)+rp_counter)] =out_folder+'\\'+filename_rp_pts
            
            
            
            if state=='fresh': 
                dw_list_eA[i] = energyAbs(dbh,0) # 0 --> fresh 
                if rp_state[i]==1:
                    dw_list_eA[int(len(tf)+rp_counter)] = energyAbs(dbh,0)
            
            if state=='old': 
                dw_list_eA[i] = energyAbs(dbh,(fp_prob('binominal', 0.333)+1))    #binominal distribution between "old" =1 and old with fungi ="
                if rp_state[i]==1:
                    dw_list_eA[int(len(tf)+rp_counter)] = energyAbs(dbh,0)
            
              #fp_prob('random',0,0))   --> First run: fresh Deadwood. Afterwards, activte comment..  
              #--> No Decay of rootplate, Energyabsorbtion of Fresh deadwood
            
            dw_complete[i][0] = dw_list_path[i] 
            dw_complete[i][1] = dw_list_eA[i]
            if rp_state[i]==1:
                dw_complete[int(len(tf)+rp_counter)][0] = dw_list_path[int(len(tf)+rp_counter)] 
                dw_complete[int(len(tf)+rp_counter)][1] = dw_list_eA[int(len(tf)+rp_counter)]
            
            
            dw_complete_cmd_line[i][0] = dw_list_eA[i]
            dw_complete_cmd_line[i][1] = finTtopX
            dw_complete_cmd_line[i][2] = finTtopY
            dw_complete_cmd_line[i][3] = finTtopZ
            
            if rp_state[i]==1:
                dw_complete_cmd_line[int(len(tf)+rp_counter)][0] = dw_list_eA[int(len(tf)+rp_counter)]
                dw_complete_cmd_line[int(len(tf)+rp_counter)][1] = finTtopX
                dw_complete_cmd_line[int(len(tf)+rp_counter)][2] = finTtopY
                dw_complete_cmd_line[int(len(tf)+rp_counter)][3] = finTtopZ
                
            rp_counter=rp_counter+rp_state[i]
        
            
        
        
        # nonZero_dw_complete=np.delete(dw_complete, idx)
        
        # np.nonzero(wt_state==0)[0]                         
        # dw_complete[int(idx[0])]
        
        # df.drop(df.index[df['myvar'] == 'specific_name'], inplace = True)
        
          
    #Export routine: 
    np.savetxt(out_folder+'\\'+'dw_list_'+state+'paths.txt', dw_list_path, fmt='%s', delimiter=' ', newline='\n')
    for line in fileinput.FileInput(out_folder+'\\'+'dw_list_'+state+'paths.txt',inplace=1):
        if line.strip() == '0': continue # remove 0 Values due to wt-state<1 not all lines are processed  https://stackoverflow.com/questions/25684209/remove-line-breaks-from-txt-file
        print(line.replace('None',"").rstrip()) 
    np.savetxt(out_folder+'\\'+'dw_list_'+state+'energyThreshhold.txt', dw_list_eA, fmt='%f', delimiter=' ', newline='\n')
    for line in fileinput.FileInput(out_folder+'\\'+'dw_list_'+state+'energyThreshhold.txt',inplace=1):
        if line.strip() == '0.000000': continue # remove 0 Values due to wt-state<1 not all lines are processed  https://stackoverflow.com/questions/25684209/remove-line-breaks-from-txt-file
        print(line.replace('None',"").rstrip())
    #print("txt file exported")
    
  
    
    #dw_complete=zip(dw_list_path,dw_list_eA)
    np.savetxt(out_folder+'\\'+'dw_list_complete_'+state+'.txt', dw_complete, fmt='%s, %s', delimiter=' ', newline='\n')
    #find and replace [ ] and ,  (According https://stackoverflow.com/questions/4746190/find-and-replace-within-a-text-file-using-python 
    #and https://stackoverflow.com/questions/25940101/removing-lines-from-a-csv-with-python-also-adds-an-extra-line?noredirect=1#comment40610727_25940101) 
    for line in fileinput.FileInput(out_folder+'\\'+'dw_list_complete_'+state+'.txt',inplace=1):
        if line.strip() == '0, 0': continue # remove 0 Values due to wt-state<1 not all lines are processed  https://stackoverflow.com/questions/25684209/remove-line-breaks-from-txt-file
        line = line.replace("[","")
        line = line.replace("]","0")
        print(line.replace('None',"").rstrip()) 
        
           
    #print('dw_list_complete_'+ state + '.txt file exported')

    #dw_complete=zip(dw_list_path,dw_list_eA)
    np.savetxt(out_folder+'\\'+'dw_list_cmd_line_values_'+state+'.txt', dw_complete_cmd_line, fmt='%1.3f, %1.2f,%1.2f, %1.2f', delimiter=',', newline='\n')
    #find and replace [ ] and ,  (According https://stackoverflow.com/questions/4746190/find-and-replace-within-a-text-file-using-python 
    #and https://stackoverflow.com/questions/25940101/removing-lines-from-a-csv-with-python-also-adds-an-extra-line?noredirect=1#comment40610727_25940101) 
    for line in fileinput.FileInput(out_folder+'\\'+'dw_list_cmd_line_values_'+state+'.txt',inplace=1):
        if line.strip() == '0.000, 0.00,0.00, 0.00': continue # remove 0 Values due to wt-state<1 not all lines are processed  https://stackoverflow.com/questions/25684209/remove-line-breaks-from-txt-file
        line = line.replace(","," ")
        print(line.replace('None',"").rstrip()) 
        
    #########################################
    #########################################
    #########################################
    #Export new treefile with remaining standing trees (MK 20220726)       
    
    tfraw=pd.DataFrame()
    tfraw = pd.read_csv(treefile, header=None,  sep=r"\s+")#, usecols=[0,1,2,3], names=['x', 'y', 'H', 'BHD'])
    tfraw = pd.read_csv(treefile, header=None,  sep=r"\s+", usecols=[0,1,2,3,4,5,6,7,8,9], names=['x', 'y', 'H', 'BHD','PHD','z1','EN','z2','z3','z4'])

    #print(tfraw)

    #one, header=None, index=Non 
    index1 = np.where(wt_state == 0)     
    index1n = pd.DataFrame(index1)     
    index1ns = index1n.iloc[0,:]
    index2 = np.where(tfraw.iloc[:,3] < DBH_thrs) #include small trees below dbh-threshold
    
    print(index1ns)
    print('index 2')
    print(index2)
               
    tf_new1=pd.DataFrame()
    tf_new1 = [tfraw.loc[j]for j in index1ns]    
    tf_new1 = np.reshape(tf_new1, (len(pd.DataFrame(index1).columns), len(tfraw.columns)))
    tf_new1_df = pd.DataFrame(tf_new1, columns=['x', 'y', 'H', 'BHD','PHD','z1','EN','z2','z3','z4'])

    tf_new2=pd.DataFrame()
    tf_new2 = [tfraw.loc[j] for j in index2]   
    tf_new2 = np.reshape(tf_new2, (len(pd.DataFrame(index2).columns), len(tfraw.columns)))
    tf_new2_df = pd.DataFrame(tf_new2, columns=['x', 'y', 'H', 'BHD','PHD','z1','EN','z2','z3','z4'])

    #combine both dataframes, omit duplicate
    tf_new_df = pd.concat([tf_new1_df,tf_new2_df]).drop_duplicates().reset_index(drop=True)
    #print(tf_new_df)
    
    #safe tupdated tree file!  
    #tf_new_df.to_csv(out_folder+'\\tree_file_standing_trees.txt',sep=',', na_rep='', float_format=None, columns=None, index_label=None, mode='w', encoding=None, compression='infer', quoting=None, quotechar='"', line_terminator=None, chunksize=None, date_format=None, doublequote=True, escapechar=None, decimal='.', errors='strict', storage_options=None)
    #tf_new_df.to_csv(out_folder+'\\tree_file_standing_trees.txt', header=False, index=False)
    tf_new_df.to_csv(out_folder+'\\tree_file_standing_trees_tab.txt', sep='\t',  header=False, index=False)
    tf_new_df.to_csv(out_folder+'\\tree_file_standing_trees.txt', sep=' ',  header=False, index=False)


     
    return(dem_res, winDir, rp_state, wt_state)

    

    
