/* Last edited on 2024-12-21 11:35:45 by stolfi */

  /* Programa para testar o volume de um tetraedro no R4
   Mediante o uso do produto vetorial e mediante as
   coordenadas de seus nodes que o definem */

#include <stdio.h>

#include <math.h>
#include <r4.h>
 

double CalVol(*p1,p2,p3,p4: r4_t)
      d1,d2,d3,d4,vol : double;
  { with (
    double l1 = (r4_t){p1[0],p1[1],p1[2],1.0};
    double l2 = (r4_t){p2[0],p2[1],p2[2],1.0};
    double l3 = (r4_t){p3[0],p3[1],p3[2],1.0};
    double l4 = (r4_t){p4[0],p4[1],p4[2],1.0};

    double l5 = (r4_t){p1[0],p1[1],p1[3],1.0};
    double l6 = (r4_t){p2[0],p2[1],p2[3],1.0};
    double l7 = (r4_t){p3[0],p3[1],p3[3],1.0};
    double l8 = (r4_t){p4[0],p4[1],p4[3],1.0};

    double l9 =  (r4_t){p1[1],p1[2],p1[3],1.0};
    double l10 = (r4_t){p2[1],p2[2],p2[3],1.0};
    double l11 = (r4_t){p3[1],p3[2],p3[3],1.0};
    double l12 = (r4_t){p4[1],p4[2],p4[3],1.0};

    double l13 = (r4_t){p1[0],p1[2],p1[3],1.0};
    double l14 = (r4_t){p2[0],p2[2],p2[3],1.0};
    double l15 = (r4_t){p3[0],p3[2],p3[3],1.0};
    double l16 = (r4_t){p4[0],p4[2],p4[3],1.0}
  ){        
      d1 = r4_det(l1,l2,l3,l4);
      d2 = r4_det(l5,l6,l7,l8);
      d3 = r4_det(l9,l10,l11,l12);
      d4 = r4_det(l13,l14,l15,l16);

      vol = 1.0/6.0 * sqrt(d1*d1+d2*d2+d3*d3+d4*d4);
      return vol;
    }
 } CalVol;
    /*
    VAR
    p1,p2,p3,p4,A,B,C,n : r4_t;
    vol,vol1 : double;
    
    */

double *a,b,c,d;
{
  a =  79.05;
  b = 201.45;
  c =  38.25;
  d =  142.8;

  fprintf(Stdio.stdout, Fmt.Int(ROUND(a)) & Fmt.Int(ROUND(b)) \
                           Fmt.Int(ROUND(c)) & Fmt.Int(ROUND(d)) & "\n");

    
  /*
  p1 = (r4_t){ 1.000001,  1.000001,  1.00001,  1.00001};
  p2 = (r4_t){ 3.000001,  3.000001,  3.00000,  3.00000};
  p3 = (r4_t){-0.147001, -0.793001,  0.66301, -0.00002};
  p4 = (r4_t){-0.000005, -0.677001,  0.70002, -0.27301};
  vol = CalVol(p1,p2,p3,p4);
  fprintf(Stdio.stdout,
    "V. em funcao das coordenadas dos seus nodes extremos: " \
      Fmt.LongReal(vol, Fmt.Style.Fix,prec = 2) & "\n");

  A = r4_Sub(p2,p1);
  B = r4_Sub(p4,p1);
  C = r4_Sub(p3,p2);
  n = r4_cross(A,B,C);
  vol1 = 1.0/6.0 * r4_Norm(n);
  fprintf(Stdio.stdout,
    "V. em f. do produto vetorial dos nodes extremos: " \
      Fmt.LongReal(vol1, Fmt.Style.Fix,prec = 2) & "\n"); */
} Volume.