/* This Program compute the statistics for test the influence of energy
   minimization (produced by OptShape). So print the average and standard
   desviation of length of @{edge->?}s, area of walls, volume of tetrahedra and
   dihedral angle between consecutives walls incidents to common @{edge->?}s.
   Read files (-final.tp).*/

#define Statistics_C_author \
  "Created by L. P. Lozada, 1999\n" \
  "Modified by L.A.P.Lozada on 99-11-23."

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

CONST
double Pi = Math.Pi;


TYPE
typedef struct Options_t {
     char *inFileTp;         /* Initial guess file name (minus ".tp") */
     char *inFileSt;         /* Initial guess file name (minus ".tp") */
     bool_t length;
     bool_t area;
     bool_t volume;
     bool_t angle;
  }

TYPE
double bool_vec_t = ARRAY OF bool_t;
double @{Edge->?} = Triangulation.edge;
double Wall = Triangulation.Wall;
??? DRF = REF uint_vec_t;

  
VAR stl, stv, sta, std: Stat.T;

Options_t *GetOptions(int argc, char **argv);

int main(int argc, char **argv)
  { Stat.Init(stl);  Stat.Init(stv); Stat.Init(sta); Stat.Init(std); 
    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.ReadToMa(o->inFileTp);
    ??? top = tc.top;
    ??? rc = Triangulation.ReadState(o->inFileSt);
    ??? c = rc^;
    {
      MakeTest(o->length,o->area,o->volume,o->angle,top,c);    
    
    return 0;
  }

PROCEDURE MakeTest(
   bool_t length,area,volume,angle;
   *ElemTableRec_t *top;
   Coords_t *c;
  )
  { if (length){ MakeTestLength(top,c); }
    if (area){ MakeTestArea(top,c); }
    if (volume){ MakeTestVolume(top,c); }
    if (angle){ MakeTestAngle(top,c); }
  } /* END MakeTest */


PROCEDURE MakeTestLength(ElemTableRec_t *top, Coords_t *c)
  
  { filefmt_write_comment(stderr, "\nThe length of @{edge->?}s\n",'|');
    for (i = 0; i < top->NE; i++)
    {
      ??? e = top->edge[i];
      Node_t un = OrgV(e.pa)->num;
      Node_t vn = OrgV(Clock(e.pa))->num;
      ??? l = r4_Dist(c[un];
      with (c[vn]),
      double ll = FLOAT(l, REAL)
    ){
        fprintf(stderr, Fmt.LongReal(l) & "\n");
        Stat.Accum(stl, ll);
      }
    }
    fprintf(stderr, "\nWeight statistics of length's @{edge->?}s :\n");
    Stat.Print(stderr, stl);
    fprintf(stderr, "\n");
  } /* END MakeTestLength */

PROCEDURE MakeTestAngle(ElemTableRec_t *top, Coords_t *c)
   
   
  double Angle(*f1,f2: Place_t)
    {
      ??? a = f1;
        ??? b = f2;
      Node_t ao = OrgV(a)->num;
      Node_t ad = OrgV(NextE(a))->num;
      Node_t ae = OrgV(PrevE(a))->num;
      Node_t bo = OrgV(b)->num;
      Node_t bd = OrgV(NextE(b))->num;
      Node_t be = OrgV(PrevE(b))->num;
      ??? v1a = r4_Sub(c[ad];
      with (c[ao]), 
      double v2a = r4_Sub(c[ae],c[ao]);
      double v1b = r4_Sub(c[bd],c[bo]);
      double v2b = r4_Sub(c[be],c[bo]);
      double p = FindOrthonormal(v1a,v2a);
      double q = FindOrthonormal(v1b,v2b);
      double cos = r4_dot(p,q)
        /*tethar == Math.acos(cos)     /* Aproximation */*/
    ){
        /*
        assert(ao == bo) && (ad == bd);
        assert(0.0<=tethar) && ((tethar <= FLOAT(Pi,double)));
        tetha = (180.0*tethar)/((double)Pi);
        */
        return cos;
      }
    } Angle;

    r4_t FindOrthonormal(*v1,v2: r4_t)
      {
        ??? u1 = v1;
        ??? v2_u1 = r4_Project(v2;
        with (u1),
        double u2 = r4_Sub(v2,v2_u1)
        DO 
          return r4_Dir(u2);
        }
      } FindOrthonormal;

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

  VAR
    double dihedral;
    REAL dihedralr;
    DRF drf;
    @{edge->?}Relevant: REF bool_vec_t;
    bool_t internal;
    double idealcosine;
  { @{edge->?}Relevant = bool_vec_new(top->NE);
    filefmt_write_comment(stderr, "\nThe dihedral angles\n",'|');
    drf = NEW(REF uint_vec_t , top->NE);  
    for (l = 0; l < top->NE; l++)
    {
      ??? e = top->edge[l];
        ??? fie = CollectWalls(e);
      ??? drf = LAST(fie^)+1;
      {
        for (j = 0; j < drf; j++)
        {
          ??? f1 = fie^[j];
          ??? a = f1.pa;
          ??? aPpos = PposP(a);
          ??? aPneg = PnegP(a);
          {
            if ((aPpos!=NULL) && (aPneg!=NULL
            )){       
                internal = TRUE;
              }
              else
              {
                internal = FALSE;
            }
          }  
          @{edge->?}Relevant[l] = TRUE) && (internal;
        };              

        if (@{edge->?}Relevant[l] 
        ){
             fprintf(stderr, "  -corresponding to @{edge->?}: " \
                        Fmt.Int(l) & "\n");
            for (j = 0; j < drf; j++)
            {
              ??? f1 = fie^[j];
              ??? a = f1.pa;
              Place_t b = NextF(a);
              { 
                dihedral = Angle(a,b);
                fprintf(stderr, Fmt.LongReal(dihedral) & "\n");
                dihedralr = FLOAT(dihedral, REAL);
                Stat.Accum(std,dihedralr);
              }
            }
            idealcosine=cos((2.0*FLOAT(Pi,double))/FLOAT(drf,double));
         } 
       }
     }
     fprintf(stderr, "\nWeight statistics of dihedral angles:\n");
     Stat.Print(stderr, std);
     fprintf(stderr, "\n");
     fprintf(stderr,"Ideal Cosine " & Fmt.LongReal(idealcosine) & "\n");
  } /* END MakeTestAngle */

PROCEDURE MakeTestVolume(ElemTableRec_t *top, Coords_t *c)
  
  { filefmt_write_comment(stderr, "\nThe volume of tetrahedrons\n",'|');
    for (i = 0; i < top->cell.nel; i++)
    {
      ??? vt = TetrahedronNodes(Srot(top.cell[i]));
        ??? un = vt[0];
      ??? vn = vt[1];
      ??? wn = vt[2];
      ??? xn = vt[3];
      ??? uv = r4_Sub(c[vn];
      with (c[un]),
      double uw = r4_Sub(c[wn],c[un]);
      double ux = r4_Sub(c[xn],c[un]);
      double n = r4_cross(uv, uw, ux);
      double v = 1.0/6.0 * r4_Norm(n);
      double vv = FLOAT(v, REAL)
    ){
        fprintf(stderr, Fmt.LongReal(v) & "\n");
        Stat.Accum(stv, vv);
      }
    }
      
    fprintf(stderr, "\nWeight statistics of volume's cell:\n");
    Stat.Print(stderr, stv);
    fprintf(stderr, "\n");
  } /* END MakeTestVolume */

PROCEDURE TetrahedronNodes(f:TriangulationPlace_t): ARRAY [0..3] OF uint ==
  /* Was exchange OrgV by OrgV for considering the dual space. */
  { Place_t g = Tors(f);
    Place_t h = Tors(Clock(PrevE(f)));
    ??? p = OrgV(g)->num;
    ??? q = OrgV(NextE(g))->num;
    ??? r = OrgV(PrevE(g))->num;
    ??? s = OrgV(PrevE(h))->num;
    {
      return ARRAY [0..3] OF uint {p, q, r, s};
    }
  } /* END TetrahedronNodes */
  
PROCEDURE MakeTestArea(ElemTableRec_t *top, Coords_t *c)
  
  { filefmt_write_comment(stderr, "\nThe area of triangles\n",'|');
    for (i = 0; i < top->wall.nel; i++)
    {
      ??? p = top->wall[i];
      Node_t un = OrgV(p.pa)->num;
      Node_t vn = OrgV(NextE(p.pa))->num;
      Node_t wn = OrgV(PrevE(p.pa))->num;
      ??? a = r4_Dist(c[vn];
      with (c[un]),
      double b = r4_Dist(c[wn],c[un]);
      double c = r4_Dist(c[vn],c[wn]);
      double p = (a+b+c);
      double s = 0.5 * p;
      double sa = s - a;
      double sb = s - b;
      double sc = s - c;
      double ar = sqrt(s*sa*sb*sc);
      double arr = FLOAT(ar, REAL)
    ){
        fprintf(stderr, Fmt.LongReal(ar) & "\n");
        Stat.Accum(sta, arr);
      }
    }
      
    fprintf(stderr, "\nWeight statistics of area's walls:\n");
    Stat.Print(stderr, sta);
    fprintf(stderr, "\n");
  } /* END MakeTestArea */ 

    
Options_t *GetOptions(int argc, char **argv)
  { 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, "-inFileTp");                               
        o->inFileTp = argparser_get_next(pp);
        argparser_get_keyword(pp, "-inFileSt");                               
        o->inFileSt = argparser_get_next(pp);
        o->length = argparser_keyword_present(pp, "-length");
        o->area   = argparser_keyword_present(pp, "-area");
        o->volume = argparser_keyword_present(pp, "-volume");
        o->angle  = argparser_keyword_present(pp, "-angle");

        argparser_finish(pp);                                       
      
  -----------------------------------
#define _HELP \
  fprintf(stderr, "Usage:  Statistics \\\n" \
   "        -inFileTp <name> \\\n" \
   "        -inFileSt <name> \\\n" \
   "        [ -length ] [ -volume ] [ -area ]" \
                                      " [ -angle ]\n");
          
END¦                                              
      }
    }
    return o;
  } /* END GetOptions */

{
  DoIt()
} Statistics.


  /* Copyright © 1999     Universidade Estadual de Campinas (UNICAMP)     */