# -*- coding: utf-8 -*-
"""
Created on Tue Nov 27 16:14:14 2018

@author: noordam
"""
import helpFunctions
import math
import shutil
import DGeoStabParser

class NWO:
    
    def __init__(self):
        self.localXCoords=[]
        self.localYCoords=[]
        
        self.globalXCoords=[]
        self.globalYCoords=[]
        
        self.NWOXBoundaryType=0
        self.NWOYBoundaryType=0
        self.epsilon=1e-1
        
#==============================================================================
#  local coordinates of NWO
#==============================================================================
    def CreateCustomNWO(self, h1,n1,h2,n2,b):
        x1Loc =0;
        y1Loc =0;
        
        x2Loc = x1Loc+h1 * n1;
        y2Loc = y1Loc-h1;
        
        x3Loc = x2Loc + math.sqrt(b**2-(h2-h1)**2);
        y3Loc = y1Loc-h2;
        
        x4Loc = x3Loc + h2*n2;
        y4Loc=y1Loc;
        
        #bTot = x4Loc-x1Loc;
        
        self.localXCoords = [x1Loc, x2Loc, x3Loc, x4Loc];
        self.localYCoords = [y1Loc, y2Loc, y3Loc, y4Loc];

        
    def CalculateGlobalCoordinatesCustomNWO(self,newX1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine):
        x1Glob = newX1Glob
        xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine)
        y1Glob =  helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1)
        #==============================================================================
        # calculate global coordinates
        #==============================================================================
        x1Loc=self.localXCoords[0]
        x2Loc=self.localXCoords[1]
        x3Loc=self.localXCoords[2]
        x4Loc=self.localXCoords[3]

        y2Loc=self.localYCoords[1]
        y3Loc=self.localYCoords[2]
        
        bTot = x4Loc-x1Loc;
        
        
        intersectionIsFound = False
        ii=0
        while not intersectionIsFound and ii<len(xCoordsOriginalSurfaceLine):
            try:
                xCoordLine4Start = xCoordsOriginalSurfaceLine[coordIndexLineStart1+ii]
                yCoordLine4Start = yCoordsOriginalSurfaceLine[coordIndexLineStart1+ii]
                
                xCoordLine4End = xCoordsOriginalSurfaceLine[coordIndexLineStart1+1+ii]
                yCoordLine4End = yCoordsOriginalSurfaceLine[coordIndexLineStart1+1+ii]
                
                [x4Glob,y4Glob] = helpFunctions.computeCircleLineIntersection(xCoordLine4Start,yCoordLine4Start,xCoordLine4End, yCoordLine4End, x1Glob, y1Glob, bTot)
                intersectionIsFound = True;
            except:
                ii=ii+1
                if ii==xCoordsOriginalSurfaceLine:
                    print('No intersection is found for point 4')
                pass
        
        dy = y4Glob - y1Glob;
        alpha = math.asin(dy/bTot)
        
        x2Glob = math.cos(alpha)*x2Loc -math.sin(alpha)*y2Loc + x1Glob;   
        y2Glob = math.sin(alpha)*x2Loc +math.cos(alpha)*y2Loc + y1Glob;
        
        x3Glob = math.cos(alpha)*x3Loc -math.sin(alpha)*y3Loc + x1Glob;
        y3Glob = math.sin(alpha)*x3Loc +math.cos(alpha)*y3Loc + y1Glob;
        
        
            
        # limit the NWO coordinates such that there is no backwards propegating surfaceline
        if x2Glob<x1Glob+self.epsilon:
            x2Glob=x1Glob+self.epsilon
        if x3Glob>x4Glob-self.epsilon:
            x3Glob=x4Glob-self.epsilon
        
        self.globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob];
        self.globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob];
  

    def CreateBoundaryAtPleistocene(self,newX1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine,fileName):
        
        x1Glob = newX1Glob
        xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine)
        y1Glob =  helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1)
        
       # file = r".\R032+20_Bishop_freatische lijn_aangepaste schematisatie.sti"
        dgs = DGeoStabParser.DGeoStability(fileName)
        
        x2Glob= newX1Glob+self.epsilon
        y2Glob = dgs.getLayerTop(0, x2Glob)
        
        x3Glob = xCoordsOriginalSurfaceLine[-1]-self.epsilon
        y3Glob = dgs.getLayerTop(0, x3Glob)

        x4Glob = xCoordsOriginalSurfaceLine[-1]
        y4Glob = yCoordsOriginalSurfaceLine[-1]
        
        self.globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob];
        self.globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob];
        
# =============================================================================
# 
# 
# 
#  x1Glob = newX1Glob
#         xCoordsLine1, yCoordsLine1, coordIndexLineStart1 = helpFunctions.GetNearestLineAtX(x1Glob,xCoordsOriginalSurfaceLine,yCoordsOriginalSurfaceLine)
#         y1Glob =  helpFunctions.GetYatX(x1Glob,xCoordsLine1,yCoordsLine1)
#         #==============================================================================
#         # calculate global coordinates
#         #==============================================================================
#         
#         # =============================================================================
#         #     Add forbidden line to NWO
#         # =============================================================================
#         old_name='D:\\noordam\\Documents\\DAMNWO\\calculations\\test2\\input\\Geometries\\DWP_3.sti'
#         new_name= workingDirectory+'Geometries\\DWP_3_'+str(x1Glob)+'.sti'
#         shutil.copy(old_name, new_name)
#         
#        # readAndWriteFiles.AddForbiddenLine(x1Glob,x2Glob,y1Glob,y2Glob,xCoordsNewSurfaceLine,yCoordsNewSurfaceLine,new_name)
#         
#         #file = r".\R032+20_Bishop_freatische lijn_aangepaste schematisatie.sti"
#         dgs = DGeoStabParser.DGeoStability(new_name)
#         
#         
#         
#         intersectionIsFound = False
#         ii=0
#         while not intersectionIsFound and ii<len(xCoordsOriginalSurfaceLine):
#             try:
#                 xCoordLine4Start = xCoordsOriginalSurfaceLine[coordIndexLineStart1+ii]
#                 yCoordLine4Start = yCoordsOriginalSurfaceLine[coordIndexLineStart1+ii]
#                 
#                 xCoordLine4End = xCoordsOriginalSurfaceLine[coordIndexLineStart1+1+ii]
#                 yCoordLine4End = yCoordsOriginalSurfaceLine[coordIndexLineStart1+1+ii]
#                 
#                 [x4Glob,y4Glob] = helpFunctions.computeCircleLineIntersection(xCoordLine4Start,yCoordLine4Start,xCoordLine4End, yCoordLine4End, x1Glob, y1Glob, bTot)
#                 intersectionIsFound = True;
#             except:
#                 ii=ii+1
#                 if ii==xCoordsOriginalSurfaceLine:
#                     print('No intersection is found for point 4')
#                 pass
#             
#         dy = y4Glob - y1Glob;
#         alpha = math.asin(dy/bTot)
#         
#         x2Glob = math.cos(alpha)*x2Loc -math.sin(alpha)*y2Loc + x1Glob;   
#         y2Glob = dgs.getLayerTop(0, x2Glob)
#         
#         #y2Glob = math.sin(alpha)*x2Loc +math.cos(alpha)*y2Loc + y1Glob;
#         
#         x3Glob = math.cos(alpha)*x3Loc -math.sin(alpha)*y3Loc + x1Glob;
#         y3Glob = dgs.getLayerTop(0, x3Glob)
#        # y3Glob = math.sin(alpha)*x3Loc +math.cos(alpha)*y3Loc + y1Glob;
#         
#         
#             
#         # limit the NWO coordinates such that there is no backwards propegating surfaceline
#         if x2Glob<x1Glob+epsilon:
#             x2Glob=x1Glob+epsilon
#         if x3Glob>x4Glob-epsilon:
#             x3Glob=x4Glob-epsilon
#         
#         globalXCoords = [x1Glob, x2Glob, x3Glob, x4Glob];
#         globalYCoords = [y1Glob, y2Glob, y3Glob, y4Glob];      
# =============================================================================