//---------------------------------------------------------------------------


#pragma hdrstop

#include "MT.h"
#include <math.h>

//---------------------------------------------------------------------------
void CreateFirstEdge(TProcess *p);


void FindIntersectionPatch(TProcess *p,
                           TPatch **Patch,
                           T3DPoint CutPlane,
                           double dCutPlane);



void FindNextPatch(TPatch *Patch,
                   TPatch **NewPatch);

T3DPoint getNormal(TProcess *p, TCubeEdge *Edge);

int IsEdgeInCube(TCubeEdge *Edge, TCube *Cube);

void LevelUpEdge(TCubeEdge **Edge);

void MostSubdividedEdge(TCubeEdge **Edge);

void FindFirstPatch(TProcess *p, TCube *Cube, TPatch **Patch);

void FindPatchForEdge(TProcess   *p,
                      TPatch     *StartingPatch,
                      T3DPoint    CutPlane,
                      double      dCutPlane,
                      T3DPoint    EndPlane,
                      double      dEndPlane,
                      TPatch      **Patch);

void MoveVertexToPatch(TProcess   *p,
                      T3DPoint    InPlane,
                      double      dInPlane,
                      T3DPoint    CutPlane,
                      double      dCutPlane,
                      T3DPoint   *Point,
                      T3DPoint   *Normal,
                      TPatch     **Patch);

void Normalize(T3DPoint *Vector);

void AddEdgeToHash(TProcess   *p,
                   TTriangleEdges *Edge);

void RemoveEdgeFromHash(TProcess   *p,
                   TTriangleEdges *Edge);

int CheckDistance(TProcess   *p,
                  T3DPoints  *A,
                  T3DPoints  *B,
                  T3DPoints  *C,
                  T3DPoints  *P1,
                  T3DPoints  *P2,
                  double      DiamMult,
                  T3DPoints **ClosestPoint,
                  T3DNormals **Normal,
                  TTriangleEdges **Edge);

void (*FindVertex)(TProcess *p,
                TPatch **Patch,
                T3DPoint Dir,
                T3DPoint StartPoint,
                T3DPoint StartNormal,
                T3DPoint CutPlane,
                double dCutPlane,
                T3DPoint *NewVertex,
                T3DPoint *NewNormal);

void FindVertex1(TProcess *p,
                TPatch **Patch,
                T3DPoint Dir,
                T3DPoint StartPoint,
                T3DPoint StartNormal,
                T3DPoint CutPlane,
                double dCutPlane,
                T3DPoint *NewVertex,
                T3DPoint *NewNormal);

void FindVertex2(TProcess *p,
                TPatch **Patch,
                T3DPoint Dir,                
                T3DPoint StartPoint,
                T3DPoint StartNormal,
                T3DPoint CutPlane,
                double dCutPlane,
                T3DPoint *NewVertex,
                T3DPoint *NewNormal);

void RefineVertex(TProcess *p, T3DPoint *Point, T3DPoint *Normal);

void CreateEdge(TProcess *p, TTriangleEdges **NewEdge, TTriangleEdges *Left, TTriangleEdges *Right);

void CreateTriangle(TProcess *p, T3DPoints *V1, T3DPoints *V2, T3DPoints *V3,
                             T3DNormals *N1, T3DNormals *N2, T3DNormals *N3);

int IsGoodDirection(TCubeEdge* Edge, TCubeEdge* PreviousEdge, TCubeEdge* NextEdge);

int isCornerEdge(TProcess *p, TCubeEdge *Edge, TCube *Cube);

TCubeEdge* FindNextCubeEdge(TCubeEdge *Edge, TCube *Cube, int *edge1, int *edge2);

TCube * FindNeighbour(TProcess *p, int edge1, int edge2, TCube *CurrentCube, TCubeEdge **NextEdge);

TCube * FindNeighbourInFace(TProcess *p,
                            TCubeEdge *PreviousEdge,
                            TCubeEdge *Edge,
                            TCube *CurrentCube,
                            TCubeEdge **NextEdge);

//---------------------------------------------------------------------------


void CreateTriangles(TProcess *p){
  T3DPoints *PNew;
  T3DNormals *NNew;
  T3DPoints *PNew1;
  T3DNormals *NNew1;
  TPatch *Patch;
  TTriangles *Triangle;
  TTriangleEdges *Edge1;
  TTriangleEdges *Edge2;
  TTriangleEdges *AE;
  TTriangleEdges *NextEdge;
  TTriangleEdges *ExistingEdge;

  TTriangleEdges *FirstRejected = NULL;

  T3DPoint Direction;
  T3DPoint DirLeft;
  T3DPoint DirRight;
  double AlphaNew;
  double AlphaLeft;
  double AlphaRight;
  int DistanceChecked;

  FindVertex = FindVertex2;

  p->Stats.Triangles.Start = clock();
//if no edges found, find one
  CreateFirstEdge(p);

//process all edges
  for(;p->ActiveEdge != NULL; p->ActiveEdge = NextEdge){
    if( p->Config.MAX_TR_COUNT && (p->Stats.Triangles.Count > p->Config.MAX_TR_COUNT)){
      MCALLBACK(p->Callback, p, S_MT, A_MT_MAX_COUNT)
      break;
    }
    if(p->quit) break;
    if(p->ActiveEdge == FirstRejected){
      MCALLBACK(p->Callback, p, S_MT, A_MT_CYCLING)
      break;
    }
    p->NewVertex = NULL;
    MCALLBACK(p->Callback, p, S_MT, A_MT_NEW_EDGE)
    NextEdge = p->ActiveEdge->NextInList;
    if(p->ActiveEdge->Inactive){
//      myfree(p->ActiveEdge);
      continue; // with another edge, this one is inactive
    }
    DistanceChecked = 0;

    AE = p->ActiveEdge;
    PNew = (T3DPoints *) mycalloc(1, sizeof(T3DPoints));
    NNew = (T3DNormals *) mycalloc(1, sizeof(T3DNormals));

//..think about angles check before finding new vertex?!?!?

    Direction = AE->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point;
    if(AE->Edge.Patch == NULL){
//....we need to find the patch
/**/     FindPatchForEdge(p,
                        AE->Edge.Vertices[1]->Patch,
                        (AE->Edge.Normals[0]->Normal + AE->Edge.Normals[1]->Normal)
                          % (Direction),
                       (-1.0) * AE->Edge.Vertices[1]->Point
                        * ((AE->Edge.Normals[0]->Normal + AE->Edge.Normals[1]->Normal)
                          % Direction),
                        Direction,
                        (-1.0) * (Direction)
                        * ((AE->Edge.Vertices[1]->Point + AE->Edge.Vertices[0]->Point) / 2.0),
                        &(AE->Edge.Patch)
                        );   /**/
/**      FindPatchForEdge(p,
                        AE->Edge.Vertices[1]->Patch,
                        (AE->Edge.Normals[1]->Normal)
                          % (Direction),
                       (-1.0) * AE->Edge.Vertices[1]->Point
                        * ((AE->Edge.Normals[1]->Normal)
                          % Direction),
                        Direction,
                        (-1.0) * (Direction)
                        * ((AE->Edge.Vertices[1]->Point + AE->Edge.Vertices[0]->Point) / 2.0),
                        &(AE->Edge.Patch)
                        );      /**/
      MCALLBACK(p->Callback, p, S_MT, A_MT_PATCH_FOR_EDGE)
      if(AE->Edge.Patch == NULL){
        if(FirstRejected == NULL) FirstRejected = AE;
        p->ActiveEdges->NextInList = AE;
        AE->NextInList = NULL;
        p->ActiveEdges = AE;
        continue;
      }
    }
    Patch = AE->Edge.Patch;
    FindVertex(p,
              &Patch,
              Direction,
              (AE->Edge.Vertices[0]->Point + AE->Edge.Vertices[1]->Point) / 2.0,
              AE->Edge.Normals[0]->Normal + AE->Edge.Normals[1]->Normal,
              Direction,
              (-1) * ((AE->Edge.Vertices[0]->Point + AE->Edge.Vertices[1]->Point) / 2.0)
                * (Direction),
              &(PNew->Point),
              &(NNew->Normal));

    if(Patch == NULL){
      if(FirstRejected == NULL) FirstRejected = AE;
      p->ActiveEdges->NextInList = AE;
      AE->NextInList = NULL;
      p->ActiveEdges = AE;
      continue;
    }

    p->NewVertex = PNew;
    MCALLBACK(p->Callback, p, S_MT, A_MT_PATCH_FOR_EDGE)
    p->NewVertex = NULL;
//..check angles
    AlphaNew = (Direction * (PNew->Point - AE->Edge.Vertices[0]->Point))
                / ( sqrt(Direction * Direction)
                    * sqrt((PNew->Point - AE->Edge.Vertices[0]->Point) * (PNew->Point - AE->Edge.Vertices[0]->Point)));

    AlphaLeft = -1.0;
//    if(AE->Edge.Normals[0]->Normal       <
    if(((PNew->Point - AE->Edge.Vertices[0]->Point) % Direction)
        * ((AE->Left->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point)
          % Direction) > 0.0){
//....the angle has good orientation
      AlphaLeft = ((AE->Left->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point) * Direction)
                  / ( sqrt(Direction * Direction)
                    * sqrt((AE->Left->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point) * (AE->Left->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point)));
    }

    AlphaRight = -1.0;
//    if(AE->Edge.Normals[1]->Normal     <
    if( ((PNew->Point - AE->Edge.Vertices[1]->Point) % Direction)
        * (((-1.0) * Direction)
          % (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point)) > 0.0){
//....the angle has good orientation
      AlphaRight = ((AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point) * ((-1.0) * Direction))
                / (sqrt(Direction * Direction)
                  * sqrt((AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point) * (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point)));
    }

    if((AlphaLeft > AlphaRight) && ( acos(AlphaLeft)
                                  < acos(AlphaNew) * p->Config.MAX_TN)){
      if(CheckDistance(p,
                     AE->Edge.Vertices[1],
                     AE->Edge.Vertices[0],
                     AE->Left->Edge.Vertices[0],
//                     NULL, NULL,
                     AE->Left->Left->Edge.Vertices[0],
                     AE->Right->Edge.Vertices[1],
                     1.0,
                     &PNew1,
                     &NNew1,
                     &ExistingEdge)){
//......use this vertex
        DistanceChecked = 1;
      }
      else {
        if( (AE->Left->Left != AE->Right) &&
          (AE->Edge.Normals[1]->Normal * AE->Left->Edge.Normals[0]->Normal < p->Config.ALPHA_DIVIDE)){
          Edge2 = AE->Left;

/** /      FindPatchForEdge(p,
                        AE->Edge.Vertices[1]->Patch,
                        (AE->Edge.Normals[1]->Normal + AE->Left->Edge.Normals[0]->Normal)
                          % (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point),
                       (-1.0) * AE->Edge.Vertices[1]->Point
                        * ((AE->Edge.Normals[1]->Normal + AE->Left->Edge.Normals[0]->Normal)
                          % (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point)),
                        AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point,
                        (-1.0) * (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point)
                        * ((AE->Left->Edge.Vertices[0]->Point + AE->Edge.Vertices[1]->Point) / 2.0),
                        &Patch
                        );

          if(Patch == NULL){
            p->ActiveEdges->NextInList = AE;
            AE->NextInList = NULL;
            p->ActiveEdges = AE;
            continue;
          }

          Edge1 = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));
          Edge2 = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));

          Edge1->Edge.Vertices[1] = (T3DPoints *) mycalloc(1, sizeof(T3DPoints));
          Edge1->Edge.Vertices[1]->Patch = Patch;
          Edge1->Edge.Vertices[1]->Point = (AE->Left->Edge.Vertices[0]->Point + AE->Edge.Vertices[1]->Point) / 2.0;

          Edge1->Edge.Normals[1] = (T3DNormals *) mycalloc(1, sizeof(T3DNormals));
          Edge1->Edge.Normals[1]->Normal = AE->Edge.Normals[1]->Normal + AE->Left->Edge.Normals[0]->Normal;
          Normalize(&(Edge1->Edge.Normals[1]->Normal));


          MoveVertexToPatch(p,
                            (AE->Edge.Normals[1]->Normal + AE->Left->Edge.Normals[0]->Normal)
                          % (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point),
                       (-1.0) * AE->Edge.Vertices[1]->Point
                        * ((AE->Edge.Normals[1]->Normal + AE->Left->Edge.Normals[0]->Normal)
                          % (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point)),
                        AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point,
                        (-1.0) * (AE->Edge.Vertices[1]->Point - AE->Left->Edge.Vertices[0]->Point)
                        * ((AE->Left->Edge.Vertices[0]->Point + AE->Edge.Vertices[1]->Point) / 2.0),
                            &(Edge1->Edge.Vertices[1]->Point),
                            &(Edge1->Edge.Normals[1]->Normal),
                            &Patch);
//........new vertex is too close to acive edge vertices, then reject edge
          if((Patch == NULL) ||
            ((Edge1->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point)
              * (Edge1->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point) < 0.001)
            || ((Edge1->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point)
              * (Edge1->Edge.Vertices[1]->Point - AE->Edge.Vertices[1]->Point) < 0.001)
            || ((Edge1->Edge.Vertices[1]->Point - AE->Left->Left->Edge.Vertices[0]->Point)
              * (Edge1->Edge.Vertices[1]->Point - AE->Left->Left->Edge.Vertices[0]->Point) < 0.001)   ){
            p->ActiveEdges->NextInList = AE;
            AE->NextInList = NULL;
            p->ActiveEdges = AE;
            continue;
          }

          Edge1->Left = AE->Left->Left; AE->Left->Left->Right = Edge1;

          Edge1->Edge.Vertices[1]->Index = p->Vertices->Index + 1;
          Edge1->Edge.Vertices[1]->Next = p->Vertices;
          p->Vertices = Edge1->Edge.Vertices[1];

          Edge1->Edge.Normals[1]->Index = p->Normals->Index + 1;
          Edge1->Edge.Normals[1]->Next = p->Normals;
          p->Normals = Edge1->Edge.Normals[1];

//          RefineVertex(p, &(Edge1->Edge.Vertices[1]->Point), &(Edge1->Edge.Normals[1]->Normal));

          Edge1->Edge.Vertices[0] = AE->Left->Edge.Vertices[0];
          Edge1->Edge.Normals[0] = AE->Left->Edge.Normals[0];
          AddEdgeToHash(p, Edge1);

          p->ActiveEdges->NextInList = Edge1;
          p->ActiveEdges = Edge1;

          Edge2->Left = Edge1;
//          Edge2 = NULL;
//          CreateEdge(p, &Edge2, Edge1, AE);


          CreateTriangle(p, AE->Edge.Vertices[0], Edge1->Edge.Vertices[0], Edge1->Edge.Vertices[1],
                            AE->Edge.Normals[0], Edge1->Edge.Normals[0], Edge1->Edge.Normals[1]);
          /**/
        }
        else{
          Edge2 = AE->Left;
        }

        CreateTriangle(p, AE->Edge.Vertices[1], AE->Edge.Vertices[0], Edge2->Left->Edge.Vertices[1],
                      AE->Edge.Normals[1], AE->Edge.Normals[0], Edge2->Left->Edge.Normals[1]);
        FirstRejected = NULL;
        Edge1 = NULL;
        CreateEdge(p, &Edge1, Edge2->Left, AE->Right);

        AddEdgeToHash(p, Edge1);
        if(Edge1->Left == Edge1->Right) {
//          RemoveEdgeFromHash(p, Edge1);
          Edge1->Inactive = 1;
          RemoveEdgeFromHash(p, Edge1->Left);
        }
        RemoveEdgeFromHash(p, AE->Left);
        RemoveEdgeFromHash(p, AE);

        MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)
        continue;

      }
    }

    if((AlphaRight > AlphaLeft) && ( (acos(AlphaRight) )
                                     < acos(AlphaNew) * p->Config.MAX_TN)){
      if(CheckDistance(p,
                     AE->Edge.Vertices[1],
                     AE->Edge.Vertices[0],
                     AE->Right->Edge.Vertices[1],
//                     NULL, NULL,
                     AE->Right->Right->Edge.Vertices[1],
                     AE->Left->Edge.Vertices[0], 
                     1.0,
                     &PNew1,
                     &NNew1,
                     &ExistingEdge)){
//......use this vertex
        DistanceChecked = 1;
      }
      else {
        if( (AE->Left->Left != AE->Right) &&
            (AE->Edge.Normals[0]->Normal * AE->Right->Edge.Normals[1]->Normal < p->Config.ALPHA_DIVIDE)){
          Edge2 = AE->Right;
/** /
          FindPatchForEdge(p,
                        AE->Right->Edge.Vertices[1]->Patch,
                        (AE->Edge.Normals[0]->Normal + AE->Right->Edge.Normals[1]->Normal)
                          % (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point),
                       (-1.0) * AE->Edge.Vertices[0]->Point
                        * ((AE->Edge.Normals[0]->Normal + AE->Right->Edge.Normals[1]->Normal)
                          % (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point)),
                        AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point,
                        (-1.0) * (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point)
                        * ((AE->Right->Edge.Vertices[1]->Point + AE->Edge.Vertices[0]->Point) / 2.0),
                        &Patch
                        );

          if(Patch == NULL){
            p->ActiveEdges->NextInList = AE;
            AE->NextInList = NULL;
            p->ActiveEdges = AE;
            continue;
          }

          Edge1 = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));
          Edge2 = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));

          Edge1->Edge.Vertices[0] = (T3DPoints *) mycalloc(1, sizeof(T3DPoints));
          Edge1->Edge.Vertices[0]->Patch = Patch;
//          Edge1->Edge.Vertices[0]->Point = (AE->Right->Edge.Vertices[1]->Point + AE->Edge.Vertices[0]->Point) / 2.0;

          Edge1->Edge.Normals[0] = (T3DNormals *) mycalloc(1, sizeof(T3DNormals));
//          Edge1->Edge.Normals[0]->Normal = AE->Edge.Normals[0]->Normal + AE->Right->Edge.Normals[1]->Normal;
//          Normalize(&(Edge1->Edge.Normals[0]->Normal));

          MoveVertexToPatch(p,
                            (AE->Edge.Normals[0]->Normal + AE->Right->Edge.Normals[1]->Normal)
                             % (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point),
                            (-1.0) * AE->Edge.Vertices[0]->Point
                            * ((AE->Edge.Normals[0]->Normal + AE->Right->Edge.Normals[1]->Normal)
                            % (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point)),
                            AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point,
                            (-1.0) * (AE->Right->Edge.Vertices[1]->Point - AE->Edge.Vertices[0]->Point)
                            * ((AE->Right->Edge.Vertices[1]->Point + AE->Edge.Vertices[0]->Point) / 2.0),
                            &(Edge1->Edge.Vertices[0]->Point),
                            &(Edge1->Edge.Normals[0]->Normal),
                            &Patch);

          if((Patch == NULL) ||
            ((Edge1->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point)
              * (Edge1->Edge.Vertices[0]->Point - AE->Edge.Vertices[0]->Point) < 0.001)
            || ((Edge1->Edge.Vertices[0]->Point - AE->Edge.Vertices[1]->Point)
              * (Edge1->Edge.Vertices[0]->Point - AE->Edge.Vertices[1]->Point) < 0.001)
            || ((Edge1->Edge.Vertices[0]->Point - AE->Right->Right->Edge.Vertices[1]->Point)
              * (Edge1->Edge.Vertices[0]->Point - AE->Right->Right->Edge.Vertices[1]->Point) < 0.001) ){
            p->ActiveEdges->NextInList = AE;
            AE->NextInList = NULL;
            p->ActiveEdges = AE;
            continue;
          }
//          RefineVertex(p, &(Edge1->Edge.Vertices[0]->Point), &(Edge1->Edge.Normals[0]->Normal));

          Edge1->Right = AE->Right->Right; AE->Right->Right->Left = Edge1;
          Edge1->Edge.Vertices[0]->Index = p->Vertices->Index + 1;
          Edge1->Edge.Vertices[0]->Next = p->Vertices;
          p->Vertices = Edge1->Edge.Vertices[0];

          Edge1->Edge.Normals[0]->Index = p->Normals->Index + 1;
          Edge1->Edge.Normals[0]->Next = p->Normals;
          p->Normals = Edge1->Edge.Normals[0];

          Edge1->Edge.Vertices[1] = AE->Right->Edge.Vertices[1];
          Edge1->Edge.Normals[1] = AE->Right->Edge.Normals[1];

          AddEdgeToHash(p, Edge1);

          p->ActiveEdges->NextInList = Edge1;
          p->ActiveEdges = Edge1;

          Edge2->Right = Edge1;

          CreateTriangle(p, Edge1->Edge.Vertices[1], AE->Edge.Vertices[1], Edge1->Edge.Vertices[0],
                      Edge1->Edge.Normals[1], AE->Edge.Normals[1], Edge1->Edge.Normals[0]);

          MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)    /**/
        }
        else {
          Edge2 = AE->Right;
        }

        CreateTriangle(p, AE->Edge.Vertices[1], AE->Edge.Vertices[0], Edge2->Right->Edge.Vertices[0],
                      AE->Edge.Normals[1], AE->Edge.Normals[0], Edge2->Right->Edge.Normals[0]);
        FirstRejected = NULL;
        Edge1 = NULL;
        CreateEdge(p, &Edge1, AE->Left, Edge2->Right);

        AddEdgeToHash(p, Edge1);
        
        if(Edge1->Left == Edge1->Right) {
          RemoveEdgeFromHash(p, Edge1);
          RemoveEdgeFromHash(p, Edge1->Left);
        }

        RemoveEdgeFromHash(p, AE->Right);
        RemoveEdgeFromHash(p, AE);

        MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)
        continue;
      }
    }

//..the triangle with the new vertex has to be created
    if(DistanceChecked || CheckDistance(p,
                     AE->Edge.Vertices[1],
                     AE->Edge.Vertices[0],
                     PNew,
                     AE->Right->Edge.Vertices[1],
                     AE->Left->Edge.Vertices[0],
//                     NULL, NULL,
                     1.5,
                     &PNew1,
                     &NNew1,
                     &ExistingEdge)){

      p->NewVertex = PNew1;
      MCALLBACK(p->Callback, p, S_MT, A_MT_PATCH_FOR_EDGE)
      p->NewVertex = NULL;

/**/      if((PNew1 == NULL) || CheckDistance(p,
                     AE->Edge.Vertices[1],
                     AE->Edge.Vertices[0],
                     PNew1,
                     NULL,
                     NULL,
                     1.0,
                     NULL, NULL, NULL)){     /**/
//......triangle from this edge cannot be created
//......add active edge at the end of the list
        if(FirstRejected == NULL) FirstRejected = AE;
        p->ActiveEdges->NextInList = AE;
        AE->NextInList = NULL;
        p->ActiveEdges = AE;
        continue;
/**/      }



//....use existing vertex


      if( (((ExistingEdge->Right->Edge.Vertices[1]->Point - ExistingEdge->Right->Edge.Vertices[0]->Point) *
          (AE->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point)
            > 0.0)
           && (
               ((ExistingEdge->Right->Edge.Vertices[1]->Point - ExistingEdge->Right->Edge.Vertices[0]->Point)
                % (AE->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point))
              * ((AE->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point)
                % (AE->Edge.Vertices[1]->Point - ExistingEdge->Edge.Vertices[1]->Point))
                > 0.0)
            )
          && !CheckDistance(p,
                     ExistingEdge->Right->Edge.Vertices[1],
                     ExistingEdge->Right->Edge.Vertices[0],
                     AE->Edge.Vertices[0],
                     AE->Edge.Vertices[0],
                     NULL,
                     1.0,
                     NULL, NULL, NULL)){
        if((ExistingEdge->Right != AE->Left)
            && (ExistingEdge != AE->Right)){
          Edge1 = NULL; Edge2 = NULL;
          RemoveEdgeFromHash(p, ExistingEdge->Right);
          CreateEdge(p, &Edge1, AE->Left, ExistingEdge->Right->Right);
          CreateEdge(p, &Edge2, ExistingEdge, AE);

          AddEdgeToHash(p, Edge1);
          AddEdgeToHash(p, Edge2);

          if(Edge1->Left == Edge1->Right){
            RemoveEdgeFromHash(p, Edge1);
            RemoveEdgeFromHash(p, Edge1->Left);
          }

          CreateTriangle(p, Edge2->Edge.Vertices[0], Edge1->Edge.Vertices[0], Edge1->Edge.Vertices[1],
                            Edge2->Edge.Normals[0], Edge1->Edge.Normals[0], Edge1->Edge.Normals[1]);
          FirstRejected = NULL;
          MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)
        }
      }
      else {
        if((((ExistingEdge->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point) *
          (AE->Edge.Vertices[1]->Point - ExistingEdge->Edge.Vertices[1]->Point)
            > 0.0)
            && (
               ((ExistingEdge->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point)
                % (AE->Edge.Vertices[1]->Point - ExistingEdge->Edge.Vertices[1]->Point))
              * ((AE->Edge.Vertices[1]->Point - ExistingEdge->Edge.Vertices[1]->Point)
                % (AE->Edge.Vertices[0]->Point - ExistingEdge->Edge.Vertices[1]->Point))
                > 0.0))
          && !CheckDistance(p,
                     ExistingEdge->Edge.Vertices[1],
                     ExistingEdge->Edge.Vertices[0],
                     AE->Edge.Vertices[1],
                     AE->Edge.Vertices[1],
                     NULL,
                     1.0,
                     NULL, NULL, NULL)){
          if(AE->Right != ExistingEdge){
            Edge1 = NULL; Edge2 = NULL;
            RemoveEdgeFromHash(p, ExistingEdge);
            CreateEdge(p, &Edge1, ExistingEdge->Left, AE->Right);
            CreateEdge(p, &Edge2, AE, ExistingEdge->Right);

            AddEdgeToHash(p, Edge1);
            AddEdgeToHash(p, Edge2);

            if(Edge1->Left == Edge1->Right){
              RemoveEdgeFromHash(p, Edge1);
              RemoveEdgeFromHash(p, Edge1->Left);
            }
            ExistingEdge = Edge2;
            CreateTriangle(p, Edge1->Edge.Vertices[1], Edge2->Edge.Vertices[1], Edge1->Edge.Vertices[0],
                            Edge1->Edge.Normals[1], Edge2->Edge.Normals[1], Edge1->Edge.Normals[0]);
            FirstRejected = NULL;

            MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)
          }
        }
        else {
          if(FirstRejected == NULL) FirstRejected = AE;
          p->ActiveEdges->NextInList = AE;
          AE->NextInList = NULL;
          p->ActiveEdges = AE;
          continue;
        }
      }



      Edge1 = NULL; Edge2 = NULL;
      CreateEdge(p, &Edge1, AE->Left, ExistingEdge->Right);
      CreateEdge(p, &Edge2, ExistingEdge, AE->Right);

      if(Edge1->Left == Edge1->Right){
        Edge1->Inactive = 1;
        RemoveEdgeFromHash(p, Edge1->Left);
//        RemoveEdgeFromHash(p, Edge1);
      }

      if(Edge2->Left == Edge2->Right){
        Edge2->Inactive = 1;
        RemoveEdgeFromHash(p, Edge2->Left);
//        RemoveEdgeFromHash(p, Edge2);
      }

      PNew = PNew1;
      NNew = NNew1; /**/
    }
    else {
      PNew->Index = p->Vertices->Index + 1;
      PNew->Patch = Patch;
      PNew->Next = p->Vertices;
      p->Vertices = PNew;

      NNew->Index = p->Normals->Index + 1;
      NNew->Next = p->Normals;
      p->Normals = NNew;

//      RefineVertex(p, &(PNew->Point), &(NNew->Normal));

      Edge1 = NULL; Edge2 = NULL;
      CreateEdge(p, &Edge1, AE->Left, NULL);
      CreateEdge(p, &Edge2, Edge1, AE->Right);
    }


    CreateTriangle(p, AE->Edge.Vertices[1], AE->Edge.Vertices[0], PNew,
                      AE->Edge.Normals[1], AE->Edge.Normals[0], NNew);
    FirstRejected = NULL;

    Edge1->Edge.Vertices[1] = PNew;  Edge1->Edge.Normals[1] = NNew;
    Edge2->Edge.Vertices[0] = PNew; Edge2->Edge.Normals[0] = NNew;

    AddEdgeToHash(p, Edge1);
    AddEdgeToHash(p, Edge2);
    RemoveEdgeFromHash(p, AE);
    MCALLBACK(p->Callback, p, S_MT, A_MT_TRIANGLE_CREATED)

//    if(p->Stats.Triangles.Count >= 9900) break;
  }

  p->Stats.Triangles.End = clock();
}
//---------------------------------------------------------------------------


void CreateFirstEdge(TProcess *p){
  TPatch *Patch1;
  TPatch *Patch2;
  int Direction2;
  TCubes * Cube = findCube(p,
          (int)(floor(p->Config.STARTING_POSITION.X / p->Config.BIG_STEP)) * p->Config.RC_RATIO,
          (int)(floor(p->Config.STARTING_POSITION.Y / p->Config.BIG_STEP)) * p->Config.RC_RATIO,
          (int)(floor(p->Config.STARTING_POSITION.Z / p->Config.BIG_STEP)) * p->Config.RC_RATIO);;
  //findCube(p, ((int) (3.0 / p->Config.BIG_STEP)) * p->Config.RC_RATIO , 0, 0);
  T3DPoints *P1 = (T3DPoints *) mycalloc(1, sizeof(T3DPoints));
  T3DNormals *N1 = (T3DNormals *) mycalloc(1, sizeof(T3DNormals));
  T3DPoints *P2 = (T3DPoints *) mycalloc(1, sizeof(T3DPoints));
  T3DNormals *N2 = (T3DNormals *) mycalloc(1, sizeof(T3DNormals));

  if(p->ActiveEdges == NULL ){
//..find active edge
    for(int i = 0; i < 12; i++){
      if( Cube->Cube.Edges[i] != NULL){
        FindFirstPatch(p, &(Cube->Cube), &Patch1);

        FindNextPatch(Patch1, &Patch2);
//......find first vertex
        FindVertex(p,
                    &Patch1,
                    Patch1->Edge->CubeEdges[1]->Intersection - Patch1->Edge->CubeEdges[0]->Intersection,
                    (Patch1->Edge->CubeEdges[0]->Intersection + Patch1->Edge->CubeEdges[1]->Intersection) / 2.0,
                    getNormal(p, Patch1->Edge->CubeEdges[0]) + getNormal(p, Patch1->Edge->CubeEdges[1]),
                    Patch1->Edge->CubeEdges[1]->Intersection - Patch1->Edge->CubeEdges[0]->Intersection,
                    (-1.0) * (Patch1->Edge->CubeEdges[1]->Intersection - Patch1->Edge->CubeEdges[0]->Intersection)
                      * ((Patch1->Edge->CubeEdges[0]->Intersection + Patch1->Edge->CubeEdges[1]->Intersection) / 2.0),
                    &(P1->Point),
                    &(N1->Normal));
        P1->Patch = Patch1;

//......find second vertex
        FindVertex(p,
                    &Patch2,
                    Patch2->Edge->CubeEdges[0]->Intersection - Patch2->Edge->CubeEdges[1]->Intersection,
                    (Patch2->Edge->CubeEdges[0]->Intersection + Patch2->Edge->CubeEdges[1]->Intersection) / 2.0,
                    getNormal(p, Patch2->Edge->CubeEdges[0]) + getNormal(p, Patch2->Edge->CubeEdges[1]),
                    Patch2->Edge->CubeEdges[0]->Intersection - Patch2->Edge->CubeEdges[1]->Intersection,
                    (-1.0)*(Patch2->Edge->CubeEdges[0]->Intersection - Patch2->Edge->CubeEdges[1]->Intersection)
                      * ((Patch2->Edge->CubeEdges[0]->Intersection + Patch2->Edge->CubeEdges[1]->Intersection) / 2.0),
                    &(P2->Point),
                    &(N2->Normal));

        P2->Patch = Patch2;
//......create triangle edge

        p->ActiveEdge = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));
        p->ActiveEdges = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));
//......set neighbourhood
        p->ActiveEdge->NextInList = p->ActiveEdges;
        p->ActiveEdge->Left = p->ActiveEdge->Right = p->ActiveEdges;
        p->ActiveEdges->Left = p->ActiveEdges->Left = p->ActiveEdges;
//......set first edge
        p->ActiveEdge->Edge.Vertices[0] = P1;
        p->ActiveEdge->Edge.Vertices[1] = P2;
        p->ActiveEdge->Edge.Normals[0] = N1;
        p->ActiveEdge->Edge.Normals[1] = N2;
//......set second edge
        p->ActiveEdges->Edge.Vertices[0] = P2;
        p->ActiveEdges->Edge.Vertices[1] = P1;
        p->ActiveEdges->Edge.Normals[0] = N2;
        p->ActiveEdges->Edge.Normals[1] = N1;

        P1->Index = 0;
        P2->Index = 1;
        P2->Next = P1;
        p->Vertices = P2;

        N1->Index = 0;
        N2->Index = 1;
        N2->Next = N1;
        p->Normals = N2;

        AddEdgeToHash(p, p->ActiveEdge);
        AddEdgeToHash(p, p->ActiveEdges);

        MCALLBACK(p->Callback, p, S_MT, A_MT_FIRST_EDGE)
        break;
      }
    }
  }
}


//---------------------------------------------------------------------------
void FindVertex1(TProcess *p,
                TPatch **Patch,
                T3DPoint Dir,
                T3DPoint StartPoint,
                T3DPoint StartNormal,
                T3DPoint CutPlane,
                double dCutPlane,
                T3DPoint *NewVertex,
                T3DPoint *NewNormal){

  T3DPoint TangentPlane;
  double dTangentPlane;
  T3DPoint Intersection1;
  T3DPoint Intersection2 = StartPoint;

  TPatch    *CurrentPatch = *Patch;
  TPatch    *NextPatch = *Patch;

  double    d1 = 0.0;
  double    d2 = 0.0;
  double    d1_old = 0;
  double    d2_old = 0;

  double    distance1;
  double    distance2 = 0.0;

  double    size1;
  double    size2;
  double    t;

//find triangle plane
  TangentPlane = StartNormal / sqrt(StartNormal * StartNormal);
  dTangentPlane = (-1.0) * (TangentPlane * StartPoint);

//loop
  while(1){
    Intersection1 = Intersection2;
    distance1 = distance2;
    d1_old = d1; d2_old = d2;

//..find next intersection
    FindIntersectionPatch(p, &NextPatch, CutPlane, dCutPlane);
    d1 = NextPatch->Edge->CubeEdges[0]->Intersection * CutPlane + dCutPlane;
    d2 = NextPatch->Edge->CubeEdges[1]->Intersection * CutPlane + dCutPlane;

    t = (d1 / (d1 - d2));
    Intersection2 = t * NextPatch->Edge->CubeEdges[1]->Intersection
                      + (1.0 - t) * NextPatch->Edge->CubeEdges[0]->Intersection;
//..compute distance from tangent plane
    distance2 = Intersection2 * TangentPlane + dTangentPlane;

    if(fabs(distance2) > p->Config.MAX_DISTANCE){
//....if intersection is far from tangent plane
      size1 = p->Config.MAX_DISTANCE;
      if(distance2 < 0.0){
        size1 *= -1.0;
      }
//....compute vertex
      if(fabs(distance1) > p->Config.MAX_DISTANCE){
        *NewVertex = Intersection1;
      }
      else {
        t = (distance2 - size1) / (distance2 - distance1);
        *NewVertex = t * Intersection1
                    + (1.0 - t) * Intersection2;
      }

      if( (size1 = ((*NewVertex) - StartPoint) * ((*NewVertex) - StartPoint))
           < p->Config.MIN_STEP3_QUAD){
//......if it is too close to start point
        if((size2 = (Intersection2 - StartPoint) * (Intersection2 - StartPoint))
             >= p->Config.MIN_STEP3_QUAD){
//........if Intersection2 is not too close
          size1 = sqrt((Intersection1 - StartPoint) * (Intersection1 - StartPoint));
          size2 = sqrt(size2);
          size1 -= p->Config.MIN_STEP3;
          size2 -= p->Config.MIN_STEP3;
//........compute vertex
          t = size2 / (size2 - size1);
          *NewVertex = t * Intersection1
                       + (1.0 - t) * Intersection2;
//........compute normal
          if(distance1 != 0.0){
//..........not first patch

            t = d2_old / (d2_old - d1_old);
            TangentPlane = t * getNormal(p, CurrentPatch->Edge->CubeEdges[0])
                        + (1.0 - t) * getNormal(p, CurrentPatch->Edge->CubeEdges[1]);
          }


          t = d2 / (d2 - d1);
          *NewNormal = t * getNormal(p, NextPatch->Edge->CubeEdges[0])
                          + (1.0 - t) * getNormal(p, NextPatch->Edge->CubeEdges[1]);
          t = size2 / (size2 - size1);
          *NewNormal = t * TangentPlane
                        + (1.0 - t) * (*NewNormal);
          break;
        }
      }
      else {
        if(size1 <= p->Config.MAX_STEP3_QUAD){  // size1 is the distance of vertex from startpoint
//........if it is NOT too far from start point
//........compute normal and break
          if(distance1 != 0.0){
//..........not first patch
//            TangentPlane = (d2_old * getNormal(CurrentPatch->Edge->CubeEdges[0])
//                          - d1_old * getNormal(CurrentPatch->Edge->CubeEdges[1]))
//                          / (d2_old - d1_old);
            TangentPlane = (d2_old / (d2_old - d1_old)) * getNormal(p, CurrentPatch->Edge->CubeEdges[0])
                        - (d1_old / (d2_old - d1_old)) * getNormal(p, CurrentPatch->Edge->CubeEdges[1]);
          }

          size1 = p->Config.MAX_DISTANCE;
          if(distance2 < 0.0){
            size1 *= -1.0;
          }

          distance1 -= size1;
          distance2 -= size1;
//          *NewNormal = (distance2 * TangentPlane - distance1
//                   * ((d2 * getNormal(NextPatch->Edge->CubeEdges[0])
//                       - d1 * getNormal(NextPatch->Edge->CubeEdges[1])))
//                    / (d2 - d1))
//                   / (distance2 - distance1);
          *NewNormal = (distance2 / (distance2 - distance1)) * TangentPlane
                      - (distance1 / (distance2 - distance1)) * (
                        (d2 / (d2 - d1)) * getNormal(p, NextPatch->Edge->CubeEdges[0])
                        - (d1 / (d2 - d1)) * getNormal(p, NextPatch->Edge->CubeEdges[1])
                      );
          break;
        } //endif(size1 <= p->Config.MAX_STEP3_QUAD)
      }
    }

    if( (size2 = (Intersection2 - StartPoint) * (Intersection2 - StartPoint))
        > p->Config.MAX_STEP3_QUAD){
//....if intersection is far from start point
      size1 = sqrt((Intersection1 - StartPoint) * (Intersection1 - StartPoint));
      size2 = sqrt(size2);
      size1 -= p->Config.MAX_STEP3;
      size2 -= p->Config.MAX_STEP3;
//....compute vertex
//      *NewVertex = (size2 * Intersection1
//                      - size1 * Intersection2)
//                    / (size2 - size1);
      *NewVertex = (size2 / (size2 - size1)) * Intersection1
                  - (size1 / (size2 - size1)) * Intersection2;
//....compute normal
      if(distance1 != 0.0){
//......not first patch
//        TangentPlane = (d2_old * getNormal(CurrentPatch->Edge->CubeEdges[0])
//                          - d1_old * getNormal(CurrentPatch->Edge->CubeEdges[1]))
//                          / (d2_old - d1_old);
          TangentPlane = (d2_old / (d2_old - d1_old)) * getNormal(p, CurrentPatch->Edge->CubeEdges[0])
                        - (d1_old / (d2_old - d1_old)) * getNormal(p, CurrentPatch->Edge->CubeEdges[1]);
      }
//      *NewNormal = (d2 * getNormal(NextPatch->Edge->CubeEdges[0])
//                       - d1 * getNormal(NextPatch->Edge->CubeEdges[1]))
//                       / (d2 - d1);
//      *NewNormal = (size2 * TangentPlane - size1
//                   * ( *NewNormal  ))
//                   / (size2 - size1);
      *NewNormal = (d2 / (d2 - d1)) * getNormal(p, NextPatch->Edge->CubeEdges[0])
                  - (d1 / (d2 - d1)) * getNormal(p, NextPatch->Edge->CubeEdges[1]);
      *NewNormal = (size2 / (size2 - size1)) * TangentPlane
                    - (size1 / (size2 - size1)) * (*NewNormal);
      break;
    }



//..continue with next Patch
    CurrentPatch = NextPatch;
    FindNextPatch(CurrentPatch, &NextPatch);
  }
//endloop
  Normalize(NewNormal);
  *Patch = NextPatch;
  return;
}

//---------------------------------------------------------------------------
void FindVertex2(TProcess *p,
                TPatch **Patch,
                T3DPoint Dir,
                T3DPoint StartPoint,
                T3DPoint StartNormal,
                T3DPoint CutPlane,
                double dCutPlane,
                T3DPoint *NewVertex,
                T3DPoint *NewNormal){


  T3DPoint Intersection1;
  T3DPoint Intersection2 = StartPoint;
  T3DPoint Intersection;


  T3DPoint Normal1;
  T3DPoint Normal2 = (StartNormal = (StartNormal / sqrt(StartNormal * StartNormal)));

  TPatch    *CurrentPatch = *Patch;
  TPatch    *NextPatch = *Patch;

  double    d1 = 0.0;
  double    d2 = 0.0;


  double    distance;
  double    distance1;
  double    distance2 = 0.0;

  double    Divergence1;
  double    Divergence2 = 1.0;
  double    t;
  int count = 0;

  int NormalFound = 0;


//loop
  while(1){
    Intersection1 = Intersection2;
    distance1 = distance2;
    Normal1 = Normal2;
    Divergence1 = Divergence2;

//..find next intersection
    FindIntersectionPatch(p, &NextPatch, CutPlane, dCutPlane);

    if(NextPatch == NULL) break;

    d1 = NextPatch->Edge->CubeEdges[0]->Intersection * CutPlane + dCutPlane;
    d2 = NextPatch->Edge->CubeEdges[1]->Intersection * CutPlane + dCutPlane;

    t = (d2 / (d2 - d1));
    Intersection = t * NextPatch->Edge->CubeEdges[0]->Intersection
                        + (1.0 - t) * NextPatch->Edge->CubeEdges[1]->Intersection;



    if( (Dir % (Intersection - StartPoint)) * StartNormal > 0.0 ){
      Intersection2 = Intersection;
      Normal2 = t * getNormal(p, CurrentPatch->Edge->CubeEdges[0])
            + (1.0 - t) * getNormal(p, CurrentPatch->Edge->CubeEdges[1]);

      *NewNormal = Normal2 / sqrt(Normal2 * Normal2);

/**/
      if(!NormalFound && ( Divergence2 = ((*NewNormal) * StartNormal)) < p->Config.ALPHA_STEP ){
        NormalFound = 1;
        t = (Divergence2 - p->Config.ALPHA_STEP) / (Divergence2 - Divergence1);

        *NewVertex = t * Intersection1 + (1.0 - t) * Intersection2;
        distance = ((*NewVertex) - StartPoint) * ((*NewVertex) - StartPoint);


        if(distance > p->Config.MAX_STEP3_QUAD){
          distance2 = sqrt((Intersection2 - StartPoint) * (Intersection2 - StartPoint));
          distance1 = sqrt(distance1);
          t = (distance2 - p->Config.MAX_STEP3) / (distance2 - distance1);
          *NewVertex = t * Intersection1 + (1.0 - t) * Intersection2;
          *NewNormal = t * Normal1 + (1.0 - t) * Normal2;
          break;
        }
        else {
          if(distance < p->Config.MIN_STEP3_QUAD){
            distance2 = (Intersection2 - StartPoint) * (Intersection2 - StartPoint);
            if(distance2 >= p->Config.MIN_STEP3_QUAD){
              distance2 = sqrt(distance2);
              distance1 = sqrt(distance1);
              t = (distance2 - p->Config.MIN_STEP3) / (distance2 - distance1);
              *NewVertex = t * Intersection1 + (1.0 - t) * Intersection2;
              *NewNormal = t * Normal1 + (1.0 - t) * Normal2;
              break;
            }
          }
          else {
            *NewNormal = t * Normal1 + (1.0 - t) * Normal2;
            break;
          }
        }
      }

      if(((distance2 = (Intersection2 - StartPoint) * (Intersection2 - StartPoint)) > p->Config.MIN_STEP3_QUAD)
        && NormalFound){
        distance2 = sqrt(distance2);
        distance1 = sqrt(distance1);
        t = (distance2 - p->Config.MIN_STEP3) / (distance2 - distance1);
        *NewVertex = t * Intersection1 + (1.0 - t) * Intersection2;
        *NewNormal = t * Normal1 + (1.0 - t) * Normal2;
        break;
      }

      if(distance2 > p->Config.MAX_STEP3_QUAD){
        distance2 = sqrt(distance2);
        distance1 = sqrt(distance1);
        t = (distance2 - p->Config.MAX_STEP3) / (distance2 - distance1);
        *NewVertex = t * Intersection1 + (1.0 - t) * Intersection2;
        *NewNormal = t * Normal1 + (1.0 - t) * Normal2;
        break;
      }

/**/
    }
//..continue with next Patch
    CurrentPatch = NextPatch;
    FindNextPatch(CurrentPatch, &NextPatch);

    if(count > 100){
      NextPatch = NULL;
      break;
    }
    count++;
  }
//endloop
  Normalize(NewNormal);
  *Patch = NextPatch;
  return;
}

//---------------------------------------------------------------------------
void FindIntersectionPatch(TProcess *p,
                           TPatch **Patch,
                           T3DPoint CutPlane,
                           double dCutPlane){
//cutplane: CutPlane.X * x + CutPlane.Y * y + CutPlane.Z * z + dCutPlane = 0
//compute the d from patch intersection
  double d1 = 1.0; // 1.0 becouse it's > then 0.0
  TPatch *NextPatch = NULL;// = *Patch;
  TPatch *FirstPatch = *Patch;

  double d2 = CutPlane * (*Patch)->CubeEdge->Intersection + dCutPlane;

  while(!((d1 <= 0.0) && (d2 > 0.0))){
    d1 = d2;
    if(NextPatch == FirstPatch){
      *Patch = NULL;
      break;
    }
    FindNextPatchEdge(p, *Patch, &NextPatch);
    *Patch = NextPatch;
    d2 = CutPlane * NextPatch->CubeEdge->Intersection + dCutPlane;

  }

  return;
}

//---------------------------------------------------------------------------
void FindNextPatchEdge(TProcess *p,
                       TPatch *Patch,
                       TPatch **NewPatch){
  TCube *CurrentCube;
  TCube *NeighbourCube = NULL;
  TCube *CurrentParent;
//  TCube *Parent2;
  TCubeEdge *Edge = NULL;
  TCubeEdge *NextEdge = NULL;
  int edge1, edge2;

// triangles.count == 209!!! breakpoint

  int index;
  int current = (Patch->Edge->Patches[0] == Patch)
                ? 0
                : 1;

  if( (*NewPatch = Patch->Next) == NULL){
//..next patch wasn't set yet
    CurrentCube = Patch->Edge->Cubes[current];

    for(index = 0; index < 4; index++){
      if(Patch->CubeEdge->PatchEdges[index] == NULL){
//......patch edge doesn't exist
        int i;
        i = 0;
        break;
      }

      if(Patch->CubeEdge->PatchEdges[index] == Patch->Edge){
//......find different then current one
        continue; // with another one
      }

      if(Patch->CubeEdge->PatchEdges[index]->Cubes[0] == CurrentCube){
//......the patch exists
        *NewPatch = Patch->CubeEdge->PatchEdges[index]->Patches[0];
        Patch->Next = *NewPatch;
        break;
      }
      if(Patch->CubeEdge->PatchEdges[index]->Cubes[1] == CurrentCube){
//......the patch edge exists
        *NewPatch = Patch->CubeEdge->PatchEdges[index]->Patches[1];
        Patch->Next = *NewPatch;
        break;
      }
    }

    if(index == 4){
      int i;
      i = 0;
    }

    if(*NewPatch == NULL){
//....still not found
//....create a new one
      *NewPatch = (TPatch *) mycalloc(1, sizeof(TPatch));
      Patch->Next = *NewPatch;
      if(Patch->CubeEdge->PatchEdges[index] == NULL){
//......need to create patchedge
        Patch->CubeEdge->PatchEdges[index] = (TPatchEdge *) mycalloc(1, sizeof(TPatchEdge));
        Patch->CubeEdge->PatchEdges[index]->Cubes[0] = CurrentCube;
        Patch->CubeEdge->PatchEdges[index]->Patches[0] = *NewPatch;
        Patch->CubeEdge->PatchEdges[index]->CubeEdges[0] = Patch->CubeEdge;

//......find the neighbour
        Edge = Patch->CubeEdge;
        CurrentParent = CurrentCube;
/**/
//......it's corner edge in currnet cube
        if(isCornerEdge(p, Edge, CurrentCube)){
//........find neighbour at lowest sbdlvl
          NextEdge = FindNextCubeEdge(Edge, CurrentCube, &edge1, &edge2);
          NeighbourCube = FindNeighbour(p, edge1, edge2, CurrentCube, &NextEdge);
        }
        else{
//          NeighbourCube = FindNeighbourInFace(p,
//                                             Patch->Edge->CubeEdges[current],
//                                             Edge, CurrentCube, &NextEdge);
//          int i;
//          i = 0;
//........else it's in the face
//........determine correct cube id


/**///wrong method
        while(NeighbourCube == NULL){
          for(int i = 0; i < 4; i++){
            if(Edge->Cubes[i] != NULL ){
              if(Edge->Cubes[i]->Cube.SubdivisionLevel > CurrentCube->SubdivisionLevel){
                if(&(Edge->Cubes[i]->Cube) == Patch->Edge->Cubes[current^1]) continue; // with next cube
//..............neighbours face is in currentcubes face
                for(int edge = 0; edge < 12; edge++){
                  if(Edge->Cubes[i]->Cube.Edges[edge] != Edge) continue; //with next edge

                  NeighbourCube = &(Edge->Cubes[i]->Cube);
                  NextEdge = NeighbourCube->Edges[
                          CROT[ CCALL[NeighbourCube->Configuration][C_ROT] ]
                          [ CCONF[CCALL[NeighbourCube->Configuration][C_CON] ][
                            CIROT[CCALL[NeighbourCube->Configuration][C_ROT] ][edge]
                          ][CRIGHT] ]
                          ];
//................next edge lies in currentcubes face
                  if(IsEdgeInCube(NextEdge, CurrentCube)) {
                    if( IsGoodDirection(Edge, Patch->Edge->CubeEdges[current], NextEdge) ) break;
                  }

                  NeighbourCube = NULL;
                  NextEdge = NULL;
                  break;
                }
                if(NeighbourCube != NULL) break;
              }
              else {
//..............currentcubes face is in neighbours face
                if(NextEdge == NULL){
//................find next edge in current cube
                  for(int edge = 0; edge < 12; edge++){
                    if(CurrentCube->Edges[edge] != Edge) continue; //with next edge
                    NextEdge = CurrentCube->Edges[
                          CROT[ CCALL[CurrentCube->Configuration][C_ROT] ]
                          [ CCONF[CCALL[CurrentCube->Configuration][C_CON] ][
                            CIROT[CCALL[CurrentCube->Configuration][C_ROT] ][edge]
                          ][CLEFT] ]
                          ];
                    break;
                  }
                }
//..............can't be the current one
                if(&(Edge->Cubes[i]->Cube) == CurrentParent) continue;
                NeighbourCube = &(Edge->Cubes[i]->Cube);
//..............next edge lies in currentcubes face
                if(IsEdgeInCube(NextEdge, NeighbourCube)) break;
                NeighbourCube = NULL;

              }
            } //endif(Edge->Cubes[i] != NULL)
          } // next i
//........NeighbourCube was found
          if(NeighbourCube != NULL) break;
//........move the edge one level up
          LevelUpEdge(&Edge);

//          CurrentParent = &(CurrentParent->Parent->Cube);
        } // endwhile(NeighbourCube == NULL)
        /**/
        }



        Patch->CubeEdge->PatchEdges[index]->Cubes[1] = NeighbourCube;
//......find the most subdivided edge NextEdge;
        MostSubdividedEdge(&NextEdge);
        Patch->CubeEdge->PatchEdges[index]->CubeEdges[1] = NextEdge;

//......add this patchedge to next cubedge
        for(int i = 0; i < 4; i++){
          if(NextEdge->PatchEdges[i] == NULL ){
            NextEdge->PatchEdges[i] = Patch->CubeEdge->PatchEdges[index];
            break;
          }
        }

        (*NewPatch)->CubeEdge = NextEdge;
      }
      else {// if NOT (Patch->CubeEdge->PatchEdges[index] == NULL)
        (*NewPatch)->CubeEdge = ( Patch->CubeEdge->PatchEdges[index]->CubeEdges[0] == Patch->CubeEdge
                                  ? Patch->CubeEdge->PatchEdges[index]->CubeEdges[1]
                                  : Patch->CubeEdge->PatchEdges[index]->CubeEdges[0]);
      } // endif(Patch->CubeEdge->PatchEdges[index] == NULL)
      (*NewPatch)->Edge = Patch->CubeEdge->PatchEdges[index];
      (*NewPatch)->Edge->Patches[( (*NewPatch)->Edge->Cubes[0] == CurrentCube
                                  ? 0
                                  : 1)] = (*NewPatch);
    } //endif(*NewPatch == NULL)
  } //endif( (*NewPatch = Patch->Next) == NULL)

  return;
}

//---------------------------------------------------------------------------
void FindNextPatch(TPatch *Patch,
                   TPatch **NewPatch){
  int next = (Patch == Patch->Edge->Patches[0] ? 1 : 0);
  *NewPatch = Patch->Edge->Patches[next];
  if(*NewPatch == NULL){
    (*NewPatch) = (TPatch *) mycalloc(1, sizeof(TPatch));
    (*NewPatch)->Edge = Patch->Edge;
    (*NewPatch)->CubeEdge = (Patch->CubeEdge == Patch->Edge->CubeEdges[0]
                              ? Patch->Edge->CubeEdges[1]
                              : Patch->Edge->CubeEdges[0]);
    Patch->Edge->Patches[next] = *NewPatch;
  }

  return;
}
//---------------------------------------------------------------------------
T3DPoint getNormal(TProcess *p, TCubeEdge *Edge){
  T3DPoint V[4];
  T3DPoint Normal;
  T3Index id, dir;
  TCubeVertex *v1;
  TCubeVertex *v2;
  double t;
  TCubes *Cube;
  int edge_right;

  if(Edge->NormalComputed){
    return Edge->Normal;
  }

  for(int i = 0; i < 4; i++){
//    V[i] = Edge->Intersection;
    if(Edge->Neighbours[i] == NULL){
      if(Edge->Neighbours[(i + 2) % 4] == NULL){
        if((Cube = Edge->Cubes[(i + 1) % 4]) == NULL){
          Cube = Edge->Cubes[(i + 3) % 4];  // not all can be NULL
        }
        for(int edge = 0; edge < 12; edge++){
          if(Cube->Cube.Edges[edge] != Edge) continue;
          edge_right = CROT[ CCALL[Cube->Cube.Configuration][C_ROT] ]
                              [ CCONF[CCALL[Cube->Cube.Configuration][C_CON] ][
                              CIROT[CCALL[Cube->Cube.Configuration][C_ROT]][edge]
                              ][CRIGHT] ];

          Edge->Neighbours[EDGE_NEIGHBOUR[ EDGE_BACK_NEIGHBOUR[edge][edge_right] ]] = Cube->Cube.Edges[edge_right];
//          Cube->Cube.Edges[edge]->Neighbours[ EDGE_NEIGHBOUR[ EDGE_BACK_NEIGHBOUR[edge][edge_right] ]]
//            = Cube->Cube.Edges[edge_right];
          break;
        }
        if(Edge->Neighbours[(i + 2) % 4] == NULL){
          V[i] = Edge->Neighbours[i]->Intersection;
          continue;
        }
      }

//      Edge->Neighbours[(i + 2) % 4]->Vertices[0]->Index
//      Edge->Neighbours[(i + 2) % 4]->Vertices[1]->Index

//      Edge->Vertices[0]->Index
//      Edge->Vertices[1]->Index
      if(!(Edge->Neighbours[(i + 2) % 4]->Vertices[0]->Index == Edge->Vertices[0]->Index)
        && !(Edge->Neighbours[(i + 2) % 4]->Vertices[0]->Index == Edge->Vertices[1]->Index)){
        id = Edge->Neighbours[(i + 2) % 4]->Vertices[0]->Index;
      } else {
        id = Edge->Neighbours[(i + 2) % 4]->Vertices[1]->Index;
      }

      dir.I = (Edge->Vertices[0]->Index.I - id.I)
              * ((Edge->Vertices[0]->Index.I == Edge->Vertices[1]->Index.I) ? 1 : 0);
      dir.J = (Edge->Vertices[0]->Index.J - id.J)
              * ((Edge->Vertices[0]->Index.J == Edge->Vertices[1]->Index.J) ? 1 : 0);
      dir.K = (Edge->Vertices[0]->Index.K - id.K)
              * ((Edge->Vertices[0]->Index.K == Edge->Vertices[1]->Index.K) ? 1 : 0);

      v1 = findCubeVertex(p,
                  Edge->Vertices[0]->Index.I + dir.I,
                  Edge->Vertices[0]->Index.J + dir.J,
                  Edge->Vertices[0]->Index.K + dir.K);

      v2 = findCubeVertex(p,
                  Edge->Vertices[1]->Index.I + dir.I,
                  Edge->Vertices[1]->Index.J + dir.J,
                  Edge->Vertices[1]->Index.K + dir.K);
//      if(Edge->Vertices[1])

      if(Edge->Vertices[0]->Value < 0.0){
        if(v1->Value >= 0.0){ // if v1 is outside
          v2 = Edge->Vertices[0];
        }
        else {// only v2 can be inside
          if(v2->Value < 0.0){
            v1 = Edge->Vertices[1];
          }
        }
      } else {
        if(v2->Value >= 0.0){
          v1 = Edge->Vertices[1];
        } else {
          if(v1->Value < 0.0){
            v2 = Edge->Vertices[0];
          }
        }
      }

      t = v1->Value / (v1->Value - v2->Value) ;
      V[i].X = ((t * (double) v2->Index.I )
               + ((1 - t) * (double) v1->Index.I)) * p->Config.STEP;
      V[i].Y = ((t * (double) v2->Index.J )
               + ((1 - t) * (double) v1->Index.J)) * p->Config.STEP;
      V[i].Z = ((t * (double) v2->Index.K )
               + ((1 - t) * (double) v1->Index.K)) * p->Config.STEP;

//      V[i] = Edge->Intersection;
    }
    else{
      V[i] = Edge->Neighbours[i]->Intersection;
//        NextEdge = Edge;
//        while(V[i] != Edge->Intersection){
//          for(int j = 0; j < 4; j++){
//            NextEdge = NextEdge->Neighbours[i];
//          }
//        }
    }
  }

  for(int i = 0; i < 4; i++){
    Normal = ((V[i] - Edge->Intersection) % (V[(i + 1) % 4] - Edge->Intersection));
    Normalize(&Normal);
    Edge->Normal = Edge->Normal + Normal;

  }


  Edge->Normal = Edge->Normal / 4.0;

  if(Edge->Vertices[0]->Value < 0.0){
    Edge->Normal = (-1.0) * Edge->Normal;
  }

  Edge->NormalComputed = 1;
  return Edge->Normal;
}

//---------------------------------------------------------------------------
void FindFirstPatch(TProcess *p, TCube *Cube, TPatch **Patch){
  TCube *FirstCube;
  TCube *NeighbourCube = NULL;
  TCubeEdge *FirstEdge;
  TCubeEdge *SecondEdge = NULL;
  TCubeEdge *Edge;
  int edge1, edge2;

//find not subdivided cube
  FirstCube = Cube;
  while(FirstCube->Subdivide){
    for(int i = 0; i < 8; i++){
      if((FirstCube->Children[i]->Cube.Configuration != 0)
          && (FirstCube->Children[i]->Cube.Configuration != 255)){
        FirstCube = &(FirstCube->Children[i]->Cube);
        break;
      }
    }
  }

//find first edge
  FirstEdge = FirstCube->Edges[CROT[ CCALL[FirstCube->Configuration][C_ROT] ][
                            C_TRAVERSAL_EDGES[CCALL[FirstCube->Configuration][C_CON]][0]
                              ]];

  MostSubdividedEdge(&FirstEdge);

  SecondEdge = FindNextCubeEdge(FirstEdge, FirstCube, &edge1, &edge2);
  NeighbourCube = FindNeighbour(p, edge1, edge2, FirstCube, &SecondEdge);

/** /
//find neighbour
//find the neighbour
  Edge = FirstEdge;
  while(NeighbourCube == NULL){
    for(int i = 0; i < 4; i++){
      if(Edge->Cubes[i] != NULL ){
        if(Edge->Cubes[i]->Cube.SubdivisionLevel > FirstCube->SubdivisionLevel){
//........neighbours face is in currentcubes face
          for(int edge = 0; edge < 12; edge++){
            if(Edge->Cubes[i]->Cube.Edges[edge] != Edge) continue; //with next edge

            NeighbourCube = &(Edge->Cubes[i]->Cube);
            SecondEdge = NeighbourCube->Edges[
                          CROT[ CCALL[NeighbourCube->Configuration][C_ROT] ]
                          [ CCONF[CCALL[NeighbourCube->Configuration][C_CON] ][
                            CIROT[CCALL[NeighbourCube->Configuration][C_ROT] ][edge]
                          ][CRIGHT] ]
                          ];
//..........next edge lies in currentcubes face
            if(IsEdgeInCube(SecondEdge, FirstCube)) break;

            NeighbourCube = NULL;
            SecondEdge = NULL;
            break;
          } //next edge
        }
        else {
//........currentcubes face is in neighbours face
          if(SecondEdge == NULL){
//..........find next edge in current cube
            for(int edge = 0; edge < 12; edge++){
              if(FirstCube->Edges[edge] != Edge) continue; //with next edge
              SecondEdge = FirstCube->Edges[
                      CROT[ CCALL[FirstCube->Configuration][C_ROT] ]
                        [ CCONF[CCALL[FirstCube->Configuration][C_CON] ][
                            CIROT[CCALL[FirstCube->Configuration][C_ROT] ][edge]
                          ][CLEFT] ]
                          ];
              break;
            }
          }
//........can't be the current one
          if(&(Edge->Cubes[i]->Cube) == FirstCube) continue;
          NeighbourCube = &(Edge->Cubes[i]->Cube);
//........next edge lies in currentcubes face
          if(IsEdgeInCube(SecondEdge, NeighbourCube)) break;
          NeighbourCube = NULL;

        }
      } //endif(Edge->Cubes[i] != NULL)
    } // next i
//..NeighbourCube was found
    if(NeighbourCube != NULL) break;
//..move the edge one level up
    LevelUpEdge(&Edge);
  } // endwhile(NeighbourCube == NULL)    /**/


  MostSubdividedEdge(&SecondEdge);

  (*Patch) = (TPatch *) mycalloc(1, sizeof(TPatch));
  (*Patch)->Edge = (TPatchEdge *) mycalloc(1, sizeof(TPatchEdge));
  (*Patch)->CubeEdge = SecondEdge;
  (*Patch)->Edge->CubeEdges[0] = FirstEdge;
  (*Patch)->Edge->CubeEdges[1] = SecondEdge;
  (*Patch)->Edge->Patches[0] = (*Patch);
  (*Patch)->Edge->Cubes[0] = FirstCube;
  (*Patch)->Edge->Cubes[1] = NeighbourCube;

  FirstEdge->PatchEdges[0] = (*Patch)->Edge;
  SecondEdge->PatchEdges[0] = (*Patch)->Edge;

}

//---------------------------------------------------------------------------
int IsEdgeInCube(TCubeEdge *Edge, TCube *Cube){
  T3Index id1 = Cube->Index;
  T3Index id2 = Cube->Index;
  int step = 16 / pow(2, Cube->SubdivisionLevel);

  id2.I += step;
  id2.J += step;
  id2.K += step;

  return ( (Edge->Vertices[0]->Index.I >= id1.I) && (Edge->Vertices[0]->Index.I <= id2.I)
          && (Edge->Vertices[0]->Index.J >= id1.J) && (Edge->Vertices[0]->Index.J <= id2.J)
          && (Edge->Vertices[0]->Index.K >= id1.K) && (Edge->Vertices[0]->Index.K <= id2.K))
         && ((Edge->Vertices[1]->Index.I >= id1.I) && (Edge->Vertices[1]->Index.I <= id2.I)
          && (Edge->Vertices[1]->Index.J >= id1.J) && (Edge->Vertices[1]->Index.J <= id2.J)
          && (Edge->Vertices[1]->Index.K >= id1.K) && (Edge->Vertices[1]->Index.K <= id2.K));
}
//---------------------------------------------------------------------------
void LevelUpEdge(TCubeEdge **Edge) {
  for(int i = 0; i < 4; i++){
    if((*Edge)->Cubes[i] != NULL){
      for(int edge = 0; edge < 12; edge++){
        if((*Edge)->Cubes[i]->Cube.Edges[edge] != (*Edge)) continue;
        *Edge = (*Edge)->Cubes[i]->Cube.Parent->Cube.Edges[edge];
        break;
      }

      break;
    }
  }
}

//---------------------------------------------------------------------------
void MostSubdividedEdge(TCubeEdge **Edge){
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  TCubeEdge *Edge1;
  TCubeEdge *Edge2;
  int end = 0;

  if(*Edge == NULL){
    int i;
    i =  0;
  }

  while(!end){
    end = 1;
//..while not found
    for(int cube = 0; cube < 4; cube++){
      if(((*Edge)->Cubes[cube] != NULL) && ((*Edge)->Cubes[cube]->Cube.Subdivide)){
//......found a cube at this edge, that was subdivided
        for(int edge = 0; edge < 12; edge++){
          if((*Edge)->Cubes[cube]->Cube.Edges[edge] != (*Edge)) continue;
//........found the edge index
          Edge1 = (*Edge)->Cubes[cube]->Cube.Children[EDGE_CHILDREN[edge][0]]->Cube.Edges[edge];
          if((Edge1 != NULL) && ((Edge1->Vertices[0]->Value < 0.0) != (Edge1->Vertices[1]->Value < 0.0))) {
//..........found the proper edge
            *Edge = Edge1;
            end = 0;
            break;
          }
          Edge2 = (*Edge)->Cubes[cube]->Cube.Children[EDGE_CHILDREN[edge][1]]->Cube.Edges[edge];
          if((Edge2 != NULL) && ((Edge2->Vertices[0]->Value < 0.0) != (Edge2->Vertices[1]->Value < 0.0))) {
//..........found the proper edge
            *Edge = Edge2;
            end = 0;
            break;
          }
        } // next edge
        break;
      } // endif(((*Edge)->Cubes[cube] != NULL) && ((*Edge)->Cubes[cube]->Cube.Subdivide))
    } // next cube
  } // endwhile(!end)

}


//---------------------------------------------------------------------------
void FindPatchForEdge(TProcess   *p,
                      TPatch     *StartingPatch,
                      T3DPoint    CutPlane,
                      double      dCutPlane,
                      T3DPoint    EndPlane,
                      double      dEndPlane,
                      TPatch      **Patch){

  double d = 1.0;
  double d1;
  double d2;
  int pass = 0;
  T3DPoint Intersection;
  TPatch *NextPatch = StartingPatch;

  while(d > 0.0){
    FindIntersectionPatch(p, &NextPatch, CutPlane, dCutPlane);
    if(NextPatch == NULL){
      *Patch = NULL;
      break;
    }
    d1 = NextPatch->Edge->CubeEdges[0]->Intersection * CutPlane + dCutPlane;
    d2 = NextPatch->Edge->CubeEdges[1]->Intersection * CutPlane + dCutPlane;
    Intersection = (d1 / (d1 - d2))* NextPatch->Edge->CubeEdges[1]->Intersection
                    - (d2 / (d1 - d2)) * NextPatch->Edge->CubeEdges[0]->Intersection;


    d = EndPlane * NextPatch->CubeEdge->Intersection + dEndPlane;
    if(d > 0.0){
      d = EndPlane * NextPatch->Edge->CubeEdges[0]->Intersection + dEndPlane;
      if(d < 0.0){
//        d = EndPlane * Intersection + dEndPlane;
      }
    }

//..continue with next Patch
    *Patch = NextPatch;
    FindNextPatch(*Patch, &NextPatch);
    if(pass > 100){
      *Patch = NULL;
      break;
    }
    pass++;
  }

  return;
}
//---------------------------------------------------------------------------
void Normalize(T3DPoint *Vector){
  double size = sqrt((*Vector) * (*Vector));

//  if(size > 0.000000000001){
    (*Vector) = (*Vector) / size;
//  }
//  else {
//    (*Vector).X = 0.0;
//    (*Vector).Y = 0.0;
//    (*Vector).Z = 0.0;
//  }
}

//---------------------------------------------------------------------------
void AddEdgeToHash(TProcess   *p,
                   TTriangleEdges *Edge){
  if(Edge->Inactive == 1) return;

  Edge->Index.I = floor(Edge->Edge.Vertices[1]->Point.X / p->Config.BIG_STEP);
  Edge->Index.J = floor(Edge->Edge.Vertices[1]->Point.Y / p->Config.BIG_STEP);
  Edge->Index.K = floor(Edge->Edge.Vertices[1]->Point.Z / p->Config.BIG_STEP);

  Edge->NextInHash = p->TriangleEdges[HASH(Edge->Index.I, Edge->Index.J, Edge->Index.K)];
  p->TriangleEdges[HASH(Edge->Index.I, Edge->Index.J, Edge->Index.K)] = Edge;
  if(Edge->NextInHash != NULL){
    Edge->NextInHash->PreviousInHash = Edge;
  }

  return;
}

//---------------------------------------------------------------------------
void RemoveEdgeFromHash(TProcess   *p,
                   TTriangleEdges *Edge){
  if(Edge->Inactive != 1){
    Edge->Inactive = 1;
    if(Edge->PreviousInHash == NULL){
      p->TriangleEdges[HASH(Edge->Index.I, Edge->Index.J, Edge->Index.K)] = Edge->NextInHash;
    }
    else {
      Edge->PreviousInHash->NextInHash = Edge->NextInHash;
    }

    if(Edge->NextInHash != NULL){
      Edge->NextInHash->PreviousInHash = Edge->PreviousInHash;
    }
  }
//?? free memory ???

}

//---------------------------------------------------------------------------
int CheckDistance1(TProcess   *p,
                  T3DPoints  *A,
                  T3DPoints  *B,
                  T3DPoints  *C,
                  T3DPoints  *P1,
                  T3DPoints  *P2,
                  double      DiamMult,
                  T3DPoints **ClosestPoint,
                  T3DNormals **Normal,
                  TTriangleEdges **ExistingEdge){
  T3Index Id1;
  T3Index Id2;
  T3DPoint Center;
  T3DPoint TriangleNormal;
  double Diameter;
  double Distance;
  double ClosestDistance;
  int SecondNormal;
  int Result = 0;
  double Angle = -1.0;

  Center = ( A->Point + B->Point + C->Point ) / 3.0;
  Diameter = ((C->Point - Center) * (C->Point - Center));
  Diameter *= DiamMult * DiamMult;
  ClosestDistance = 2 * Diameter;

  Id1.I = floor((Center.X - Diameter ) / p->Config.BIG_STEP) - 10;
  Id1.J = floor((Center.Y - Diameter ) / p->Config.BIG_STEP) - 10;
  Id1.K = floor((Center.Z - Diameter ) / p->Config.BIG_STEP) - 10;

  Id2.I = ceil((Center.X + Diameter ) / p->Config.BIG_STEP) + 10;
  Id2.J = ceil((Center.Y + Diameter ) / p->Config.BIG_STEP) + 10;
  Id2.K = ceil((Center.Z + Diameter ) / p->Config.BIG_STEP) + 10;

  TriangleNormal = (A->Point - B->Point) % (C->Point - A->Point);

  for(int i = Id1.I; i <= Id2.I; i++){
    for(int j = Id1.J; j <= Id2.J; j++){
      for(int k = Id1.K; k <= Id2.K; k++){
        for(TTriangleEdges *Edges = p->TriangleEdges[HASH(i,j,k)]; Edges != NULL; Edges = Edges->NextInHash){
          if( Edges->Inactive
              || (Edges->Index.I != i)
              || (Edges->Index.J != j)
              || (Edges->Index.K != k)) {
              continue; // with another edge
          }

          if( (Edges->Edge.Vertices[1] == A)
             || (Edges->Edge.Vertices[1] == B)
             || (Edges->Edge.Vertices[1] == C)
             || (Edges->Edge.Vertices[1] == P1)
             || (Edges->Edge.Vertices[1] == P2)){
             continue; // with another edge
          }

          if( (Distance = (Center - Edges->Edge.Vertices[1]->Point)
                         *(Center - Edges->Edge.Vertices[1]->Point)) <= Diameter){

//..........check normal orientation
            SecondNormal = 0;
            if(((Edges->Edge.Vertices[1]->Point - A->Point)
                % (Edges->Edge.Vertices[1]->Point - B->Point)) * Edges->Edge.Normals[1]->Normal > 0.0){
//............different orientation
              if(Edges->Right->Edge.Normals[0] != Edges->Edge.Normals[1]){
//..............vertex lying on a sharp edge
                if(((Edges->Edge.Vertices[1]->Point - A->Point)
                % ((Edges->Edge.Vertices[1]->Point - B->Point))) * Edges->Right->Edge.Normals[0]->Normal < 0.0) continue;  // with another edge
                SecondNormal = 1;
              }
              else{
                continue; // with another edge
              }
            }

//            if(((Edges->Edge.Vertices[1]->Point - A->Point)
//                % (Edges->Edge.Vertices[1]->Point - B->Point)) * TriangleNormal > 0.2){
//              continue;
//            }

//..........it's close to center
            if(ClosestPoint == NULL) return 1; // only need to know, there is such a vertex


//..........it is closer than the closest?
            if(Distance <= ClosestDistance){
/*              if(Edges->Edge.Vertices[1] == (*ClosestPoint)){
                if((A->Point - B->Point) * (Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point )
                  / sqrt((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point ) * (Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point ))
                  < Angle){
                  continue;
                }
              }
              Angle = (A->Point - B->Point) * (Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point )
                / sqrt((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point ) * (Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point ));*/
              Result = 1;
              (*ExistingEdge) = Edges;
              ClosestDistance = Distance;
              (*ClosestPoint) = Edges->Edge.Vertices[1];
              if( SecondNormal
                || ( TriangleNormal * Edges->Right->Edge.Normals[0]->Normal
                      > TriangleNormal * Edges->Edge.Normals[1]->Normal)) {
//..............use the other normal
                (*Normal) = Edges->Right->Edge.Normals[0];
              }
              else {
                (*Normal) = Edges->Edge.Normals[1];
              }
            }
          }
        }
      }
    }
  }

  return Result;
}

//---------------------------------------------------------------------------
int CheckDistance(TProcess   *p,
                  T3DPoints  *A,
                  T3DPoints  *B,
                  T3DPoints  *C,
                  T3DPoints  *P1,
                  T3DPoints  *P2,
                  double      DiamMult,
                  T3DPoints **ClosestPoint,
                  T3DNormals **Normal,
                  TTriangleEdges **ExistingEdge){
  T3Index Id1;
  T3Index Id2;
  T3DPoint Center;
  T3DPoint TriangleNormal;
  int isFirst;
  int isSecond;
  double Diameter;
  double Distance1;
  double Distance2;
  double ClosestDistance;
  int SecondNormal;
  int Result = 0;
  double Angle = -1.0;

  double a2 = (C->Point - B->Point) * (C->Point - B->Point);
  double b2 = (C->Point - A->Point) * (C->Point - A->Point);
  double c2 = (A->Point - B->Point) * (A->Point - B->Point);

//[A*(a^2*(b^2+c^2-a^2) +
//B*(b^2*(a^2+c^2-b^2) +
//C*(c^2*(a^2+b^2-c^2)]
//---------------------------------------------------------------
//2*(a^2*b^2+a^2*c^2+b^2*c^2)-(a^4+b^4+c^4)

  Center = (  (a2 *(b2 + c2 - a2)) * A->Point
            + (b2 *(a2 + c2 - b2)) * B->Point
            + (c2 *(a2 + b2 - c2)) * C->Point )
        / (2*(a2 * b2 + a2*c2 + b2*c2) - (a2*a2 + b2*b2 + c2*c2)) ;
//  Center = ( A->Point + B->Point + C->Point ) / 3.0;
  try{
    Diameter =  ((C->Point - Center) * (C->Point - Center));
    Diameter *=  DiamMult * DiamMult;
    ClosestDistance = sqrt(Diameter) + 2 * p->Config.MAX_STEP3;
    ClosestDistance *= ClosestDistance;
  }
  catch(...){
    return 1;
  }

  if(Diameter > 64 * p->Config.MAX_STEP3_QUAD) return 1;

  Id1.I = floor((Center.X - ClosestDistance ) / p->Config.BIG_STEP) - 1;
  Id1.J = floor((Center.Y - ClosestDistance ) / p->Config.BIG_STEP) - 1;
  Id1.K = floor((Center.Z - ClosestDistance ) / p->Config.BIG_STEP) - 1;

  Id2.I = ceil((Center.X + ClosestDistance ) / p->Config.BIG_STEP) + 1;
  Id2.J = ceil((Center.Y + ClosestDistance ) / p->Config.BIG_STEP) + 1;
  Id2.K = ceil((Center.Z + ClosestDistance ) / p->Config.BIG_STEP) + 1;

  TriangleNormal = (A->Point - B->Point) % (C->Point - B->Point);

  for(int i = Id1.I; i <= Id2.I; i++){
    for(int j = Id1.J; j <= Id2.J; j++){
      for(int k = Id1.K; k <= Id2.K; k++){
        for(TTriangleEdges *Edges = p->TriangleEdges[HASH(i,j,k)]; Edges != NULL; Edges = Edges->NextInHash){
          if( Edges->Inactive
              || (Edges->Index.I != i)
              || (Edges->Index.J != j)
              || (Edges->Index.K != k)) {
              continue; // with another edge
          }

/**/      if( (P1 == NULL) && ((Edges->Edge.Vertices[1] == A)
             || (Edges->Edge.Vertices[1] == B)
             || (Edges->Edge.Vertices[1] == C)
//             || (Edges->Edge.Vertices[1] == P1)
//             || (Edges->Edge.Vertices[1] == P2)
             || (Edges->Edge.Vertices[0] == A)
             || (Edges->Edge.Vertices[0] == B)
             || (Edges->Edge.Vertices[0] == C))){

             continue; // with another edge
          }     /**/

          if( (P2 != NULL) && ((A->Point - B->Point) * (Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point) > 0.0)){
            continue;
          }

          Distance1 = (Center - Edges->Edge.Vertices[1]->Point)
                      * (Center - Edges->Edge.Vertices[1]->Point);

          Distance2 = (Center - Edges->Edge.Vertices[0]->Point)
                         * (Center - Edges->Edge.Vertices[0]->Point);

          isFirst = (Distance1 <= ClosestDistance) &&
                    (((Edges->Edge.Vertices[1]->Point - A->Point)
                % (Edges->Edge.Vertices[1]->Point - B->Point)) * Edges->Edge.Normals[1]->Normal < 0.0);
          isFirst = isFirst && (Edges->Edge.Vertices[1] != A)
                            && (Edges->Edge.Vertices[1] != B)
                            && (Edges->Edge.Vertices[1] != C)
                            && (Edges->Edge.Vertices[1] != P1)
                            && (Edges->Edge.Vertices[1] != P2);

          isSecond = (Distance2 <= ClosestDistance) &&
                    (((Edges->Edge.Vertices[0]->Point - A->Point)
                % (Edges->Edge.Vertices[0]->Point - B->Point)) * Edges->Edge.Normals[0]->Normal < 0.0);
          isSecond = isSecond && (Edges->Edge.Vertices[0] != A)
                              && (Edges->Edge.Vertices[0] != B)
                              && (Edges->Edge.Vertices[0] != C)
                              && (Edges->Edge.Vertices[0] != P1)
                              && (Edges->Edge.Vertices[0] != P2);

          isSecond = isSecond && (!isFirst || (Distance2 < Distance1));
          isFirst = isFirst && !isSecond;                    

          if( isFirst || isSecond ){


/*
//..........check normal orientation
            if(((Edges->Edge.Vertices[1]->Point - A->Point)
                % (Edges->Edge.Vertices[1]->Point - B->Point)) * Edges->Edge.Normals[1]->Normal > 0.0){
//............different orientation
              continue; // with another edge
            }       */


            if( (((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                  % (Center - Edges->Edge.Vertices[0]->Point))
                  * ((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                  % (Center - Edges->Edge.Vertices[0]->Point))
                / ((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                  *(Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point))) > Diameter){
//............the edge is further than diameter
              continue;
            }

            if( ( (Distance1 < Diameter) || (Distance2 < Diameter))
              || (
                (((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                    * Edges->Edge.Vertices[1]->Point)
                  - ((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                    * Center))
                 * (((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                    * Edges->Edge.Vertices[0]->Point)
                  - ((Edges->Edge.Vertices[1]->Point - Edges->Edge.Vertices[0]->Point)
                    * Center)) < 0.0
                  )
                ){


//............it's close to center
              if(ClosestPoint == NULL) return 1; // only need to know, there is such a vertex

//............it is closer than the closest?
              Result = 1;
              if(isFirst){
                (*ExistingEdge) = Edges;
                ClosestDistance = Distance1;
                (*ClosestPoint) = Edges->Edge.Vertices[1];
                (*Normal) = Edges->Edge.Normals[1];
              } else {
                (*ExistingEdge) = Edges->Left;
                ClosestDistance = Distance2;
                (*ClosestPoint) = Edges->Edge.Vertices[0];
                (*Normal) = Edges->Edge.Normals[0];
              }
            }
          }
        }
      }
    }
  }

  return Result;
}

//---------------------------------------------------------------------------
void RefineVertex(TProcess *p, T3DPoint *Point, T3DPoint *Normal){
  T3DPoint P1;
  T3DPoint P2;
  double Value;

  return;
  P1 = (*Point);
  while(p->Function(P1.X, P1.Y, P1.Z) > 0.0){
    P1 = P1 + (0.01) * (*Normal);
  }
  P2 = P1 + (-0.01) * (*Normal);
  while(p->Function(P2.X, P2.Y, P2.Z) < 0.0){
    P2 = P2 + (-0.01) * (*Normal);
  }


  while(fabs(Value = p->Function(Point->X, Point->Y, Point->Z)) > 0.001){
    if(Value > 0.0){
      P2 = *Point;
    }
    else {
      P1 = *Point;
     }
    *Point = (P1 + P2) / 2.0;
  }

  Normal->X = (p->Function(Point->X + 0.00001, Point->Y, Point->Z)
              - p->Function(Point->X - 0.00001, Point->Y, Point->Z)) / 0.00002;
  Normal->Y = (p->Function(Point->X , Point->Y + 0.00001, Point->Z)
              - p->Function(Point->X , Point->Y - 0.00001, Point->Z)) / 0.00002;
  Normal->Z = (p->Function(Point->X , Point->Y, Point->Z + 0.00001)
              - p->Function(Point->X , Point->Y, Point->Z - 0.00001)) / 0.00002;

  *Normal = (*Normal) / sqrt((*Normal) * (*Normal));

  *Normal = (-1.0) * (*Normal);


}

//---------------------------------------------------------------------------
void MoveVertexToPatch(TProcess   *p,
                      T3DPoint    InPlane,
                      double      dInPlane,
                      T3DPoint    CutPlane,
                      double      dCutPlane,
                      T3DPoint   *Point,
                      T3DPoint   *Normal,
                      TPatch     **Patch){

  T3DPoint Intersection1;
  T3DPoint Intersection2;
  T3DPoint Normal1;
  T3DPoint Normal2;

  double d1, d2, t;


  FindIntersectionPatch(p, Patch, CutPlane, dCutPlane);

  if(*Patch == NULL){
    return;
  }

  d1 = (*Patch)->Edge->CubeEdges[0]->Intersection * CutPlane + dCutPlane;
  d2 = (*Patch)->Edge->CubeEdges[1]->Intersection * CutPlane + dCutPlane;

  t = (d1 / (d1 - d2));
  Intersection1 = t * (*Patch)->Edge->CubeEdges[1]->Intersection
                      + (1.0 - t) * (*Patch)->Edge->CubeEdges[0]->Intersection;
  Normal1 = t * getNormal(p, (*Patch)->Edge->CubeEdges[1])
            + (1.0 - t) * getNormal(p, (*Patch)->Edge->CubeEdges[0]);

  FindIntersectionPatch(p, Patch, (-1.0) * CutPlane, (-1.0) * dCutPlane);

  if(*Patch == NULL){
    return;
  }

  d1 = (*Patch)->Edge->CubeEdges[0]->Intersection * CutPlane + dCutPlane;
  d2 = (*Patch)->Edge->CubeEdges[1]->Intersection * CutPlane + dCutPlane;

  t = (d1 / (d1 - d2));
  Intersection2 = t * (*Patch)->Edge->CubeEdges[1]->Intersection
                      + (1.0 - t) * (*Patch)->Edge->CubeEdges[0]->Intersection;

  Normal2 = t * getNormal(p, (*Patch)->Edge->CubeEdges[1])
            + (1.0 - t) * getNormal(p, (*Patch)->Edge->CubeEdges[0]);

  d1 = Intersection1 * InPlane + dInPlane;
  d2 = Intersection2 * InPlane + dInPlane;
  t = (d1 / (d1 - d2));
  *Point = t * Intersection2
            + (1.0 - t) * Intersection1;

  *Normal = t * Normal2
            + (1.0 - t) * Normal1;

  Normalize(Normal);

}

//---------------------------------------------------------------------------
void CreateEdge(TProcess *p, TTriangleEdges **NewEdge, TTriangleEdges *Left, TTriangleEdges *Right){

  if(*NewEdge == NULL){
    *NewEdge = (TTriangleEdges *) mycalloc(1, sizeof(TTriangleEdges));
  }

  if(Left != NULL){
    (*NewEdge)->Left = Left; Left->Right = (*NewEdge);
    (*NewEdge)->Edge.Vertices[0] = Left->Edge.Vertices[1];
    (*NewEdge)->Edge.Normals[0] = Left->Edge.Normals[1];
  }

  if(Right != NULL){
    (*NewEdge)->Right = Right; Right->Left = (*NewEdge);
    (*NewEdge)->Edge.Vertices[1] = Right->Edge.Vertices[0];
    (*NewEdge)->Edge.Normals[1] = Right->Edge.Normals[0];
  }

  p->ActiveEdges->NextInList = (*NewEdge);
  p->ActiveEdges = (*NewEdge);
}


//---------------------------------------------------------------------------
void CreateTriangle(TProcess *p, T3DPoints *V1, T3DPoints *V2, T3DPoints *V3,
                             T3DNormals *N1, T3DNormals *N2, T3DNormals *N3){

  TTriangles * Triangle = (TTriangles *) mycalloc(1, sizeof(TTriangles));

  Triangle->Triangle.Vertices[0] = V1;
  Triangle->Triangle.Vertices[1] = V2;
  Triangle->Triangle.Vertices[2] = V3;

  Triangle->Triangle.Normals[0] = N1;
  Triangle->Triangle.Normals[1] = N2;
  Triangle->Triangle.Normals[2] = N3;

  Triangle->Next = p->FoundTriangles;
  p->FoundTriangles = Triangle;
  p->Stats.Triangles.Count++;
}

//---------------------------------------------------------------------------
int IsGoodDirection(TCubeEdge* Edge, TCubeEdge* PreviousEdge, TCubeEdge* NextEdge){
  T3Index delta;
  T3Index id;

  if(!(PreviousEdge->Vertices[0]->Index == Edge->Vertices[0]->Index)
    && !(PreviousEdge->Vertices[0]->Index == Edge->Vertices[1]->Index)){
    id = PreviousEdge->Vertices[0]->Index;
  } else {
    id = PreviousEdge->Vertices[1]->Index;
  }

  delta.I = (Edge->Vertices[0]->Index.I - id.I)
           * ((Edge->Vertices[0]->Index.I == Edge->Vertices[1]->Index.I) ? -1 : 0);
  delta.J = (Edge->Vertices[0]->Index.J - id.J)
           * ((Edge->Vertices[0]->Index.J == Edge->Vertices[1]->Index.J) ? -1 : 0);
  delta.K = (Edge->Vertices[0]->Index.K - id.K)
           * ((Edge->Vertices[0]->Index.K == Edge->Vertices[1]->Index.K) ? -1 : 0);

  if(!(NextEdge->Vertices[0]->Index == Edge->Vertices[0]->Index)
    && !(NextEdge->Vertices[0]->Index == Edge->Vertices[1]->Index)){
    id = NextEdge->Vertices[0]->Index;
  } else {
    id = NextEdge->Vertices[1]->Index;
  }

  if(delta.I != 0){
    if((Edge->Vertices[0]->Index.I - id.I) == 0){
      return 1;
    }
    return
      (delta.I < 0)
      !=
      (((Edge->Vertices[0]->Index.I - id.I)
           * ((Edge->Vertices[0]->Index.I == Edge->Vertices[1]->Index.I) ? -1 : 0)) < 0);
  }

  if(delta.J != 0){
    if((Edge->Vertices[0]->Index.J - id.J) == 0){
      return 1;
    }
    return
      (delta.J < 0)
      !=
      (((Edge->Vertices[0]->Index.J - id.J)
           * ((Edge->Vertices[0]->Index.J == Edge->Vertices[1]->Index.J) ? -1 : 0)) < 0);
  }

  if(delta.K != 0){
    if((Edge->Vertices[0]->Index.K - id.K) == 0){
      return 1;
    }
    return
      (delta.K < 0)
      !=
      (((Edge->Vertices[0]->Index.K - id.K)
           * ((Edge->Vertices[0]->Index.K == Edge->Vertices[1]->Index.K) ? -1 : 0)) < 0);
  }

  return 1;
}

int isCornerEdge(TProcess *p, TCubeEdge *Edge, TCube *Cube){
  T3Index IdCompare;
  int size = p->Config.RC_RATIO >> Cube->SubdivisionLevel;
  IdCompare.I = (Edge->Vertices[0]->Index.I == Edge->Vertices[1]->Index.I ? 1 : 0);
  IdCompare.J = (Edge->Vertices[0]->Index.J == Edge->Vertices[1]->Index.J ? 1 : 0);
  IdCompare.K = (Edge->Vertices[0]->Index.K == Edge->Vertices[1]->Index.K ? 1 : 0);

  return (!IdCompare.I || ((Cube->Index.I - Edge->Vertices[0]->Index.I) % size == 0) )
        && (!IdCompare.J || ((Cube->Index.J - Edge->Vertices[0]->Index.J) % size == 0) )
        && (!IdCompare.K || ((Cube->Index.K - Edge->Vertices[0]->Index.K) % size == 0) );
}

TCubeEdge* FindNextCubeEdge(TCubeEdge *Edge, TCube *Cube, int *edge1, int *edge2){
  TCubeEdge *NextEdge=NULL;
  T3Index IdCompare;
  IdCompare.I = (Edge->Vertices[0]->Index.I == Edge->Vertices[1]->Index.I ? 1 : 0);
  IdCompare.J = (Edge->Vertices[0]->Index.J == Edge->Vertices[1]->Index.J ? 1 : 0);
  IdCompare.K = (Edge->Vertices[0]->Index.K == Edge->Vertices[1]->Index.K ? 1 : 0);

  for(int edge = 0; edge < 12; edge++){
    if( (Cube->Edges[edge] != NULL) && (!IdCompare.I ||
          ((Cube->Edges[edge]->Vertices[0]->Index.I == Cube->Edges[edge]->Vertices[1]->Index.I)
          && (Cube->Edges[edge]->Vertices[0]->Index.I == Edge->Vertices[0]->Index.I)))
        && (!IdCompare.J ||
          ((Cube->Edges[edge]->Vertices[0]->Index.J == Cube->Edges[edge]->Vertices[1]->Index.J)
          && (Cube->Edges[edge]->Vertices[0]->Index.J == Edge->Vertices[0]->Index.J)))
        && (!IdCompare.K ||
          ((Cube->Edges[edge]->Vertices[0]->Index.K == Cube->Edges[edge]->Vertices[1]->Index.K)
          && (Cube->Edges[edge]->Vertices[0]->Index.K == Edge->Vertices[0]->Index.K)))
      ) {

      *edge1 = edge;
      *edge2 = CROT[ CCALL[Cube->Configuration][C_ROT] ]
                          [ CCONF[CCALL[Cube->Configuration][C_CON] ][
                            CIROT[CCALL[Cube->Configuration][C_ROT] ][edge]
                          ][CLEFT] ];
      NextEdge = Cube->Edges[ *edge2 ];
      break;
    }
  }

  if(NextEdge == NULL){
    int i;
    i = 0;
  }

  return NextEdge;
}

TCube * FindNeighbour(TProcess *p, int edge1, int edge2, TCube *CurrentCube, TCubeEdge **NextEdge){
  TCube *NeighbourCube = NULL;
  T3Index NeighbourIndex = CurrentCube->Index;
  int NeighbourEdge = EDGE_BACK_NEIGHBOUR[edge1][edge2];
  TCubeEdge *Edge = NULL;
  TCube *Parent = CurrentCube;
  int size = p->Config.RC_RATIO >> CurrentCube->SubdivisionLevel;
  int dir = C_EDGES_DIRS[edge1][edge2];


  NeighbourIndex.I += (dir == 0 ? (-1) * size : (dir == 1 ? size : 0));
  NeighbourIndex.J += (dir == 3 ? (-1) * size : (dir == 2 ? size : 0));
  NeighbourIndex.K += (dir == 5 ? (-1) * size : (dir == 4 ? size : 0));

  while(NeighbourCube == NULL){
    for(int i = 0; i < 4; i++){
      if((NeighbourCube == NULL) && (Parent->Edges[C_DIR_EDGES[dir][i]] != NULL)){
        for(int j = 0; j < 4; j++){
          if(Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j] != NULL){
            T3Index id = Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j]->Cube.Index;
//..........this cube is parent of neighbour cube
            if((id.I == ((NeighbourIndex.I / size) * size - ((NeighbourIndex.I % size) < 0 ? size : 0)))
              && (id.J == ((NeighbourIndex.J / size) * size - ((NeighbourIndex.J % size) < 0 ? size : 0)))
              && (id.K == ((NeighbourIndex.K / size) * size - ((NeighbourIndex.K % size) < 0 ? size : 0)))){

              NeighbourCube = &(Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j]->Cube);
              break;
            }
          }
        }
      }
    }
    if(NeighbourCube != NULL) break;
    size = size << 1;
    Parent = &(Parent->Parent->Cube);
  }

  while(NeighbourCube->Subdivide){
    Edge = NeighbourCube->Children[EDGE_CHILDREN[NeighbourEdge][0]]->Cube.Edges[NeighbourEdge];
    if((Edge != NULL) && ((Edge->Vertices[0]->Value < 0.0) != (Edge->Vertices[1]->Value < 0.0)) ){
      NeighbourCube = &(NeighbourCube->Children[EDGE_CHILDREN[NeighbourEdge][0]]->Cube);
    }
    else {
      NeighbourCube = &(NeighbourCube->Children[EDGE_CHILDREN[NeighbourEdge][1]]->Cube);
    }
  }

  if(CurrentCube->SubdivisionLevel < NeighbourCube->SubdivisionLevel){
//..if neighbour is more subdivided
//..find neightbour at lowest subdivision level


//..find next edge
    *NextEdge = NeighbourCube->Edges[ CROT[ CCALL[NeighbourCube->Configuration][C_ROT] ]
                          [ CCONF[CCALL[NeighbourCube->Configuration][C_CON] ][
                            CIROT[CCALL[NeighbourCube->Configuration][C_ROT] ][NeighbourEdge]
                          ][CRIGHT] ] ];
  }

  if(*NextEdge == NULL){
    int i;
    i = 0;
  }

  return NeighbourCube;
}

TCube * FindNeighbourInFace(TProcess *p,
                            TCubeEdge *PreviousEdge,
                            TCubeEdge *Edge,
                            TCube *CurrentCube,
                            TCubeEdge **NextEdge){
  TCube *NeighbourCube = NULL;
  TCube *Parent = CurrentCube;
  T3Index NeighbourIndex;
  int dir;
  int size = p->Config.RC_RATIO >> CurrentCube->SubdivisionLevel;

//determine neighbour index
//..determine cubesize from edge (4 posible ids)
//..must be out of current cube (2 possible ids)
//..in direction from previous edge (1 possible id)

//find direction in which the edge lies
//..=direction of neighbour index



//find neighbour in direction
  while(NeighbourCube == NULL){
    for(int i = 0; i < 4; i++){
      if((NeighbourCube == NULL) && (Parent->Edges[C_DIR_EDGES[dir][i]] != NULL)){
        for(int j = 0; j < 4; j++){
          if(Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j] != NULL){
            T3Index id = Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j]->Cube.Index;
//..........this cube is parent of neighbour cube
            if((id.I == ((NeighbourIndex.I / size) * size - ((NeighbourIndex.I % size) < 0 ? size : 0)))
              && (id.J == ((NeighbourIndex.J / size) * size - ((NeighbourIndex.J % size) < 0 ? size : 0)))
              && (id.K == ((NeighbourIndex.K / size) * size - ((NeighbourIndex.K % size) < 0 ? size : 0)))){

              NeighbourCube = &(Parent->Edges[C_DIR_EDGES[dir][i]]->Cubes[j]->Cube);
              break;
            }
          }
        }
      }
    }

    if(NeighbourCube != NULL) break;
    size = size << 1;
    Parent = &(Parent->Parent->Cube);
  }

//subdivide to lowest level
  while(NeighbourCube->Subdivide){

  }
//find correct next edge
  //*NextEdge = NeighbourCube->Edges[ CROT[ CCALL[NeighbourCube->Configuration][C_ROT] ]
  //                        [ CCONF[CCALL[NeighbourCube->Configuration][C_CON] ][
  //                         CIROT[CCALL[NeighbourCube->Configuration][C_ROT] ][NeighbourEdge]
  //                        ][CRIGHT] ] ];

}

#pragma package(smart_init)
