#include <pz_shape.h>/* Last edited on 1999-08-06 22:17:06 by hcgl */

#include <pz_r3_chain.h>
#include <pz_double_chain.h>
#include <pz_fourier.h>
#include <r3.h>

#include <math.h>
#include <pz_types.h>


pz_double_chain_t *pz_shape_get(pz_r3_chain_t *c)
  {

    { /* with */
      ??? N = (c.ne);
      ??? rs = pz_double_chain_new(N), s = rs^;
      ??? L = Length(c);
      /* do */
      

      void Convert(
    unsigned i,
    unsigned j, 
    double f 
  )
        /*
          Converts the curve "c[i..j]" to signal "s[i+1..j-1]", 
          assuming "s[i]" and "s[j]" have been defined.
        */
        {
          if (( i < j-1 )){
            affirm((i-j) MOD 2 == 0 );

            { /* with */ 
              ??? k = (i + j) DIV 2;
              ??? h = - f * Angle(c[i], c[k], c[j]);
              /* do */
              s[k] = (s[i] + s[j] + h)/2.0e0;
              Convert(i,k,f/2.0e0); Convert(k,j,f/2.0e0);
            };
          }
        } /* Convert */
      
      {
        s[0] = 0.0e0;
        s[N-1] = 0.0e0;
        Convert(0,N-1,L/2.0e0);
      };
      return rs;
    }
  } /* Signal */
  

double  Length(
    pz_r3_chain_t *c 
  )
  /* Assumes the steps are all equal except for rounding error */
  /* and turning bias: */
  VAR s: double := 0.0e0;
  {
    for (i = 1; i < (c.ne ) ; i++){

      { /* with */
        ??? d2 = r3_dist_sqr(c[i], c[i-1]);
        /* do */ s := s + d2 END;;
    };
    return sqrt(s*FLOAT((c.ne)-1, double))
  } /* Length */
  

double  Angle(
    r3_t *a,
    r3_t *b,
    r3_t *c 
  )
  /*
    Computes the signed external angle between the vectors "b-a" and "c-b".
    Assumes the three points lie on the plane "Z==0".
  */
  {
    affirm(a[2] ==0.0e0 );
    affirm(b[2] ==0.0e0 );
    affirm(c[2] ==0.0e0 );

    { /* with */
      ??? u = r3_sub(b,a);
      ??? v = r3_sub(c,b);
      ??? s = u[0]*v[1] - u[1]*v[0];
      ??? c = u[0]*v[0] + u[1]*v[1];
      /* do */
      return atan2(s,c);
    }
  } /* Angle */


pz_fourier_t *pz_fourier_compute(pz_double_chain_t *p)
  {

    { /* with */ 
      ??? N = PowerOfTwoNotLessThan((p.ne)-1);
      /* do */
      if (( (p.ne) != N )){
        affirm((p.ne) == N+1 );
        affirm(p[0] == p[(p.ne - 1)] );;
      };

      { /* with */
        ??? q = pz_double_chain_new(N)^;
        ??? ra = NEW(REF pz_fourier_t, N), a = ra^;
        /* do */
        q = SUBARRAY(p, 0, N);
        pz_fourier.FFT(q,a);
        return ra;
      };
    }
  } /* pz_fourier_compute */


pz_power_spectrum_t *pz_power_spectrum_compute(pz_fourier_t *a)
  VAR w: double;
  {

    { /* with */
      ??? N = (a.ne);
      ??? fMax = N DIV 2;
      ??? rpwr = pz_double_chain_new(fMax+1), pwr = rpwr^;
      /* do */
      affirm(N MOD 2 == 0 );
      for (f = 0; f <= fMax ; f++){

        { /* with */
          ??? t = a[f];
          /* do */ w := t*t END;
        if (( f > 0 ) ANDAND ( f < fMax )){

          { /* with */
            ??? t = a[N-f];
            /* do */ w := w + t*t END;
        };
        pwr[f] = w;;
      };
      return rpwr;;
    }
  } /* pz_power_spectrum_compute */



unsigned  PowerOfTwoNotLessThan(
    unsigned n 
  ) 
  /* Smallest power of two not less than "n" */
  VAR m: unsigned := 1;
  {
    affirm(n > 0 );
    while ( m < n ){ m = m+m ;};
    return m
  } /* PowerOfTwoNotLessThan */

{