/* Programa que calcula os angulos diedrais de walls consecutivas incidentes
   a cada uma das arestas da triangulacao. */

#define AngleConse_C_author \
  "Created by L. P. Lozada, 1999-2000."

#include <math.h>
#include <Stat.h>
#include <Triangulation.h>
#define _GNU_SOURCE
#include <stdio.h>
#include <argparser.h>
#include <r4.h>

CONST
double Pi = Math.Pi;

TYPE
typedef struct Options_t {
      char *inFile;             /* Initial guess file name (minus ".tp") */
  }

double @{Edge->?} = Triangulation.edge;
double Wall = Triangulation.Wall;
double EndPoint = ARRAY [0..1] OF uint;
double DRF = REF uint_vec_t;
double Tri = RECORD u, v, w: uint; }  


void WriteCoord(FILE *wr, double x)
  
  { fprintf(wr, Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec = 2), 6) & " ");
  } /* END WriteCoord */

/* UNUSED */
void WritePoint(FILE *wr, r4_t *c)
  
  { WriteCoord(wr, c[0]);
    fprintf(wr, " ");
    WriteCoord(wr, c[1]);
    fprintf(wr, " ");
    WriteCoord(wr, c[2]);
    fprintf(wr, " ");
    WriteCoord(wr, c[3]);
    fprintf(wr, "  ");
  } /* END WritePoint */

/* UNUSED */
Tri ExtremusWall(*f: Wall)
  { ??? fun = f.node[0]->num;
    ??? fvn = f.node[1]->num;
    ??? fwn = f.node[2]->num;
    Tri tri = Tri{fun, fvn, fwn};
    {
      return tri;
    }
  } /* END ExtremusWall */


PROCEDURE Extremus@{Edge->?}(*e: @{Edge->?}) : EndPoint ==
   VAR uint stack[1+1];
  { ??? eun = e.node[0]->num;
    ??? evn = e.node[1]->num;
    {
      stack.el[0] = eun;
      stack.el[1] = evn;
      return stack;
    }
  } /* END Extremus@{Edge->?} */

/* UNUSED */
PROCEDURE @{Edge->?}Common(*f1,f2 : Wall) : EndPoint ==
      uint stack[1+1];
      n : uint = 0;
  { for (j = 0; j < 3; j++){
      for (i = 0; i < 3; i++)
      {
        ??? f1j = f1.node[j]->num;
          ??? f2i = f2.node[i]->num;
        {
          if (f1j == f2i)
          {
            stack.el[n] = f1j; n++;
          }
        }
      }
    }
    return stack;
  } /* END @{Edge->?}Common */    

void DoIt()
  
      sta: Stat.T;
      DRF drf;
      uint n;
  { Stat.Init(sta);
    Options_t *o = GetOptions(argc, argv);
    char *topo_cmt = jsprintf("Created by %s on %s", PROG_NAME, Today());
    /* Random_t coins = MakeRandomSource(4615); */
    ??? tc = Triangulation.Read(o->inFile);
    ??? top = tc.top;
    ??? c = tc.c^;
    {

      DRF ComputeDegreeRingWalls()
        VAR drf = NEW(REF uint_vec_t , top->NE);
        {
          for (l = 0; l < top->NE; l++)
          {
            ??? e = top->edge[l];
              ??? fie = CollectWalls(e);
            {
              drf^[l] =  ((fie^).nel);
            }
          }
        return drf;
      } ComputeDegreeRingWalls;

      /* UNUSED */
      PROCEDURE CollectIncidentWalls(*e: @{Edge->?}): REF ARRAY OF Tri ==
        uint NT = 0;
            VAR uint ct;
        {
          int NF = top->wall.nel;
          ??? t = NEW(REF ARRAY OF Tri, NF)^;
          {
            for (i = 0; i < NF; i++)
            {
              ct = 0;
              ??? f = top->wall[i];
                ??? fun = f.node[0]->num;
              ??? fvn = f.node[1]->num;
              ??? fwn = f.node[2]->num;
              ??? eun = e.node[0]->num;
              ??? evn = e.node[1]->num;
              {
                if (((fun==eun) || (fun==evn))){ ct++; }
                if (((fvn==eun) || (fvn==evn))){ ct++; }
                if (((fwn==eun) || (fwn==evn))){ ct++; }
                if (ct == 2)
                {
                  t[NT] = Tri{fun, fvn, fwn};
                  NT++;
                }
              } 
            }
            ??? r = NEW(REF ARRAY OF Tri, NT);
    {
              r^ = SUBARRAY(t,0,NT);
              return r;
            }
          }
        } CollectIncidentWalls;


        Wall_vec_t *CollectWalls(*e: @{Edge->?})==
        uint NT = 0;
            VAR uint ct;
        {
          int NF = top->wall.nel;
            ??? t = Wall_vec_new(NF)^;
          {
            for (i = 0; i < NF; i++)
            {
              ct = 0;
              ??? f = top->wall[i];
                ??? fun = f.node[0]->num;
              ??? fvn = f.node[1]->num;
              ??? fwn = f.node[2]->num;
              ??? eun = e.node[0]->num;
              ??? evn = e.node[1]->num;
              {
                if (((fun==eun) || (fun==evn))){ ct++; }
                if (((fvn==eun) || (fvn==evn))){ ct++; }
                if (((fwn==eun) || (fwn==evn))){ ct++; }
                if (ct == 2)
                {
                  t[NT] = f;
                  NT++;
                }
              } 
            }
            ??? r = Wall_vec_new(NT);
    {
              r^ = SUBARRAY(t,0,NT);
              return r;
            }
          }
        } CollectWalls;

 double ComputeAverageDiedralAngle(*fie: REF Wall_vec_t;
                                      EndPoint bi; 
                                      grau: uint)
     double tetham = 0.0;
     VAR { 
        ??? fie = fie;
          ??? de = grau;
        ??? bi = bi;
        {
          for (i = 0; i < de; i++)
          {
            ??? f1 = fie[i];
              ??? af1 = f1.pa^;
            ??? f2 = fie[(i+1) % de];
            ??? af2 = f2.pa^;
            ??? bf1 = WallBarycenter(af1,c);
            double bf2 = WallBarycenter(af2,c);
            double tetha = Angle(bf1,bf2,bi);
            {
              tetham = tetham + tetha;
            }
          }
          tetham = tetham/((double)de);
          return tetham;
        }            
      } ComputeAverageDiedralAngle;

      double Angle(*bf1,bf2 : r4_t; e : EndPoint)
       double *cos,tethar,tetha;
       {
         ??? e0 = e[0], e1 == e[1];
         double be = r4_Scale(1.0/2.0, r4_Add(c[e0],c[e1]));
         double a = r4_Sub(be,bf1);
         double b = r4_Sub(be,bf2);
         {
           cos = r4_Cos(a,b);
           tethar = sqrt( 2.0*(1.0-cos));
           assert(0.0 <= tethar) && ((tethar <= FLOAT(Pi,double)));
           tetha = (180.0*tethar)/((double)Pi);
           return tetha;
         } 
       } Angle;

      double *tethat, tetham, medio, dif;
          REAL diff;
      {
        drf = NEW(REF uint_vec_t , top->NE);  
        for (l = 0; l < top->NE; l++)
        {
          ??? e = top->edge[l];
            ??? fie = CollectWalls(e);
          ??? bi = Extremus@{Edge->?}(e);
          ??? drf = LAST(fie^)+1;
          {
            medio = ComputeAverageDiedralAngle(fie,bi,drf);
            /*
            fprintf(Stdio.stdout, "Aresta[" & Fmt.Int(e->num) & "]: " \
                       " Nodes Extremos:  " & Fmt.Int(bi[0]) & "-" \
                        Fmt.Int(bi[1]) & " , ");
            fprintf(Stdio.stdout,"DegreeRingWalls: " & Fmt.Int(drf) \
            ", Angulo ideal: " & Fmt.LongReal(360.0/FLOAT(drf,double))
            & "\n");
            fprintf(Stdio.stdout, "Walls Incidentes: \n");
            */
            dif = fabs(medio - 360.0/FLOAT(drf,double));
            diff = FLOAT(dif, REAL);
            Stat.Accum(sta,diff);
           /* 
           for(i = 0; i < fie.nel; i++)
           {
             ??? fun = fie^[i].node[0]->num;
               ??? fvn = fie^[i].node[1]->num;
             ??? fwn = fie^[i].node[2]->num;
             {
               fprintf(Stdio.stdout, Fmt.Int(fun) & " " & Fmt.Int(fvn) \
                          " " & Fmt.Int(fwn) & "\n");
             }
           }
          
          fprintf(Stdio.stdout, "Angles consecutives: ");
          tethat = 0.0;
          for (j = 0; j < drf; j++)
          {
            ??? f1 = fie[j];
            ??? af1 == f1.pa^;
            double f2 = fie[(j+1) % drf], af2 == f2.pa^;
            double bf1 = WallBarycenter(af1,c);
            double bf2 = WallBarycenter(af2,c);
            double tetha = Angle(bf1,bf2,bi);
            {
              WriteCoord(stdout,tetha);
              tethat = tethat + tetha;
            }
          }
          tetham = tethat/((double)drf);
           
           fprintf(Stdio.stdout, " ==>   Total: " \
              Fmt.Pad(Fmt.LongReal(tethat, Fmt.Style.Fix, prec = 1), 5) & " ");
           fprintf(Stdio.stdout, " ==>  Average: " \
               Fmt.Pad(Fmt.LongReal(tetham, Fmt.Style.Fix, prec = 1), 5)
               & "\n"); 
          }
          fprintf(Stdio.stdout, "\n");
          */
        }

        fprintf(stdout,
           "\nWeight statistics of diference between ideal and average diedral Angle:\n");
        Stat.Print(stdout, sta);
        fprintf(stdout, "\n");

        /*
        drf = ComputeDegreeRingWalls();
        n = ComputeDegreeTotal(drf);
        fprintf(Stdio.stdout, "Degree Ring Walls Total is : " & Fmt.Int(n)
                   & "\n");
        */
      }
    
    return 0;
  }

uint ComputeDegreeTotal(*drf: DRF)
      n : uint = 0;
  { for(l = 0; l < drf.nel; l++){
      ??? d = drf^[l];
      {
        n = n + d;
      }
    }
    return n;
  } /* END ComputeDegreeTotal */

  Options_t GetOptions ()
  { Options_t *o = (Options_t *)malloc(sizeof(Options_t));
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    
    argparser_get_keyword(pp, "-inFile");                               
        o->inFile = argparser_get_next(pp);                                
        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage: AngleConse" \
   "  -inFile <name>" \
   "\n");
          
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */


{
  DoIt()
} AngleConse.