trans_seq.c

Go to the documentation of this file.

#include "trans_seq.h"

#include "geom.h"
#include "basic_algebra.h"

#include <math.h>
#include <string.h>

/*************************************************************************/
/*************************************************************************/
/*************************************************************************/

void InitCtTrans(THTransform *ct,TTrans *t)
{
  t->t=CT_TRANS;
  t->s=1;
  t->u=NO_UINT;
  t->v=NO_UINT;
  NEW(t->ct,1,THTransform);
  HTransformCopy(t->ct,ct);
  t->vect=NULL;
  t->p=NULL;
}

void InitVarTrans(unsigned int tp,int s,unsigned int v,TTrans *t)
{
  if ((tp==CT_TRANS)||(IS_PATCH_TRANS(tp))||(tp==TV))
    Error("InitVarTrans should not be used for CT_TRANS, TV, or PA transforms");
  t->t=tp;
  t->s=s;
  t->v=v;
  t->u=NO_UINT;
  t->ct=NULL;
  t->vect=NULL;
  t->p=NULL;
}

void InitTVTrans(int s,unsigned int v,double *vect,TTrans *t)
{
  t->t=TV;
  t->s=s;
  t->v=v;
  t->u=NO_UINT;
  t->ct=NULL;
  NEW(t->vect,3,double);
  memcpy(t->vect,vect,3*sizeof(double));
  t->p=NULL;
}

void InitPatchTrans(unsigned int tp,int s,unsigned int u,unsigned int v,
                    double **p,TTrans *t)
{
  unsigned int i;

  if (!IS_PATCH_TRANS(tp))
    Error("Using InitPatchTrans to define a non-PA transform");

  if (u==v)
    Error("Repeated parameter in InitPatchTrans");
  
  t->t=tp;
  t->s=s;
  t->v=v;
  t->u=u;
  t->ct=NULL;
  t->vect=NULL;
  NEW(t->p,7,double*);
  for(i=0;i<7;i++)
    { NEW(t->p[i],3,double); }
  if (tp==PA)
    {
      memcpy(t->p[0],p[0],3*sizeof(double)); /* (p_0)  t->p_0=p_0 */
      DifferenceVector(3,p[2],p[0],t->p[1]); /* (p_02) t->p_1=p_2-p_0 */
      DifferenceVector(3,p[1],p[0],t->p[2]); /* (p_01) t->p_2=p_1-p_0 */
      DifferenceVector(3,p[0],p[1],t->p[3]); /* (d)    t->p_3=p_0-p_1-p_2+p_3 */
      DifferenceVector(3,t->p[3],p[2],t->p[3]);
      AccumulateVector(3,p[3],t->p[3]); /* t->p[3]+=p[3] */
    }
  else
    {
      /* For derived patches, we already have the four points computed */
      for(i=0;i<4;i++)
        memcpy(t->p[i],p[i],3*sizeof(double));
    }

  CrossProduct(t->p[1],t->p[2],t->p[4]); /* (n0) t_p_4=p_0_2 x p_0_1 */
  CrossProduct(t->p[1],t->p[3],t->p[5]); /* (n1) t_p_5=p_0_2 x d */
  CrossProduct(t->p[3],t->p[2],t->p[6]); /* (n2) t_p_6=d x p_0_1 */
}

void CopyTrans(TTrans *t_dst,TTrans *t_src)
{
  unsigned int i;

  t_dst->t=t_src->t;
  t_dst->s=t_src->s;
  t_dst->v=t_src->v;
  t_dst->u=t_src->u;
  if (t_src->ct!=NULL)
    {
      NEW(t_dst->ct,1,THTransform);
      HTransformCopy(t_dst->ct,t_src->ct);
    }
  else
    t_dst->ct=NULL;

  if (t_src->vect!=NULL)
    {
      NEW(t_dst->vect,3,double);
      memcpy(t_dst->vect,t_src->vect,3*sizeof(double));
    }
  else
    t_dst->vect=NULL;

  if (t_src->p!=NULL)
    {
       NEW(t_dst->p,7,double*);
       for(i=0;i<7;i++)
         {
           NEW(t_dst->p[i],3,double);
           memcpy(t_dst->p[i],t_src->p[i],3*sizeof(double));
         }
    }
  else
    t_dst->p=NULL;
}

void TransInvert(TTrans *ti,TTrans *t)
{
  CopyTrans(ti,t);
  if (ti->t==CT_TRANS)
    HTransformInverse(ti->ct,ti->ct);
  else
    ti->s=-ti->s;
}

boolean TransHasVar(unsigned int v,TTrans *t)
{
  return((v==t->v)||(v==t->u));
}

void EvaluateVectorsPATrans(double u,double v,double *x,double *y,double *h,
                            TTrans *t)
{
 if (!IS_PATCH_TRANS(t->t))
    Error("Using EvaluatePATrans to evaluate a non-PA transform");

  /* x=n0+u*n1+v*n2 */
  SumVectorScale(3,t->p[4],u,t->p[5],x);
  SumVectorScale(3,x,v,t->p[6],x);

  /* y=p_02+v*d */
  SumVectorScale(3,t->p[1],v,t->p[3],y);
                
  /* t=p_0+u*p_02+v*p_01+u*v*d */
  SumVectorScale(3,t->p[0],u,t->p[1],h);
  SumVectorScale(3,h,v,t->p[2],h);
  SumVectorScale(3,h,u*v,t->p[3],h);
}

void EvaluatePATrans(double u,double v,THTransform *a,TTrans *t)
{
  double x[3],y[3],z[3],h[3];

  EvaluateVectorsPATrans(u,v,x,y,h,t);
  Normalize(3,x);
  Normalize(3,y);                
  /* z=x X y  */
  CrossProduct(x,y,z);

  HTransformFromVectors(x,y,z,h,a);
  if (t->s<0)
    HTransformInverse(a,a);
}

void DeleteTrans(TTrans *t)
{
  if (t->ct!=NULL)
    {
      HTransformDelete(t->ct);
      free(t->ct);
      t->ct=NULL;
    }
  if (t->vect!=NULL)
    {
      free(t->vect);
      t->vect=NULL;
    }
  if (t->p!=NULL)
    {
      unsigned int i;
      
      for(i=0;i<7;i++)
        free(t->p[i]);
      free(t->p);
      t->p=NULL;
    }
}


/*************************************************************************/
/*************************************************************************/
/*************************************************************************/

void InitTransSeq(TTransSeq *ts)
{
  ts->m=INIT_NUM_TERMS_TS;
  ts->n=0;
  NEW(ts->t,ts->m,TTrans*);
}

void CopyTransSeq(TTransSeq *ts_dst,TTransSeq *ts_src)
{
  unsigned int i;

  ts_dst->m=ts_src->m;
  ts_dst->n=ts_src->n;
  
  NEW(ts_dst->t,ts_dst->m,TTrans *);
  for(i=0;i<ts_dst->n;i++)
    {
      NEW(ts_dst->t[i],1,TTrans);
      CopyTrans(ts_dst->t[i],ts_src->t[i]);
    }
}

void ResetTransSeq(TTransSeq *ts)
{ 
  unsigned int i;

  for(i=0;i<ts->n;i++)
    {
      DeleteTrans(ts->t[i]);
      free(ts->t[i]);
    }
  ts->n=0;
}

boolean IsEmptyTransSeq(TTransSeq *ts)
{
  return(ts->n==0);
}

boolean HasCtRotTransSeq(TTransSeq *ts)
{
  unsigned int i;
  boolean found;

  found=FALSE;
  for(i=0;((!found)&&(i<ts->n));i++)
    found=((ts->t[i]->t==CT_TRANS)&&(!HTransformIsTranslation(ts->t[i]->ct)));

  return(found);
}

unsigned int TransSeqSize(TTransSeq *ts)
{
  return(ts->n);
}

TTrans *GetElementFromTransSeq(unsigned int i,TTransSeq *ts)
{
  if (i>=ts->n)
    Error("Index out of range in GetVarElementFromTransSeq");

  return(ts->t[i]);
}

void AddTrans2TransSeq(TTrans *t,TTransSeq *ts)
{
  /* For CT_TRANS we try to compact, if possible */
  if (t->t==CT_TRANS)
    {
      if (t->s>0)
        AddCtTrans2TransSeq(t->ct,ts);
      else
        {
          THTransform ict;

          HTransformInverse(t->ct,&ict);
          AddCtTrans2TransSeq(&ict,ts);
          HTransformDelete(&ict);
        }
    }
  else
    {
      /* For variable transforms compact only if the previous trasform is
         the identity (this happens for just initialized transform sequences) */
      if ((ts->n>0)&&(ts->t[ts->n-1]->t==CT_TRANS)&&(HTransformIsIdentity(ts->t[ts->n-1]->ct)))
        {
          DeleteTrans(ts->t[ts->n-1]);
          CopyTrans(ts->t[ts->n-1],t);
        }
      else
        {
          if (ts->n==ts->m)
            { MEM_DUP(ts->t,ts->m,TTrans *); }
          NEW(ts->t[ts->n],1,TTrans);
          CopyTrans(ts->t[ts->n],t);
          ts->n++;
        }
    }
}

void AddCtTrans2TransSeq(THTransform *t,TTransSeq *ts)
{
  if ((ts->n>0)&&(ts->t[ts->n-1]->t==CT_TRANS))
    {
      HTransformProduct(ts->t[ts->n-1]->ct,t,ts->t[ts->n-1]->ct);
      /* Identity matrix is removed (execept the base one) */
      if ((HTransformIsIdentity(ts->t[ts->n-1]->ct))&&(ts->n>1))
        {
          ts->n--;
          DeleteTrans(ts->t[ts->n]);
          free(ts->t[ts->n]);
        }
    }
  else
    {
      if (ts->n==ts->m)
        {
          MEM_DUP(ts->t,ts->m,TTrans *);
        }
      NEW(ts->t[ts->n],1,TTrans);
      InitCtTrans(t,ts->t[ts->n]);
      ts->n++;
    }
}

void AddVarTrans2TransSeq(unsigned int t,int s,unsigned int v,TTransSeq *ts)
{
  if (t==CT_TRANS)
    Error("Adding a constant varaible transform in AddVarTrans2TransSeq");

  if (IS_PATCH_TRANS(t))
    Error("Do not use AddVarTrans2TransSeq to add a patch transform to a Trans Seq");

  if (t==TV)
    Error("Do not use AddVarTrans2TransSeq to add a displacement transform to a Trans Seq");

  if ((ts->n>0)&&(ts->t[ts->n-1]->t==CT_TRANS)&&(HTransformIsIdentity(ts->t[ts->n-1]->ct)))
    InitVarTrans(t,s,v,ts->t[ts->n-1]);
  else
    {
      if (ts->n==ts->m)
        { MEM_DUP(ts->t,ts->m,TTrans *); }
      NEW(ts->t[ts->n],1,TTrans);
      InitVarTrans(t,s,v,ts->t[ts->n]);
      ts->n++;
    }
}

void AddDispTrans2TransSeq(int s,unsigned int v,double *vect,TTransSeq *ts)
{
  if ((ts->n>0)&&(ts->t[ts->n-1]->t==CT_TRANS)&&(HTransformIsIdentity(ts->t[ts->n-1]->ct)))
    InitTVTrans(s,v,vect,ts->t[ts->n-1]);
  else
    {
      if (ts->n==ts->m)
        { MEM_DUP(ts->t,ts->m,TTrans *); }
      NEW(ts->t[ts->n],1,TTrans);
      InitTVTrans(s,v,vect,ts->t[ts->n]);
      ts->n++;
    }
}

void AddPatchTrans2TransSeq(unsigned int t,int s,unsigned int u,unsigned int v,
                            double **p,TTransSeq *ts)
{
  if ((ts->n>0)&&(ts->t[ts->n-1]->t==CT_TRANS)&&(HTransformIsIdentity(ts->t[ts->n-1]->ct)))
    InitPatchTrans(t,s,u,v,p,ts->t[ts->n-1]);
  else
    {
      if (ts->n==ts->m)
        { MEM_DUP(ts->t,ts->m,TTrans *); }
      NEW(ts->t[ts->n],1,TTrans);
      InitPatchTrans(t,s,u,v,p,ts->t[ts->n]);
      ts->n++;
    }
}

boolean VarIncludedinTransSeq(unsigned int v,TTransSeq *ts)
{
  boolean found;
  unsigned int i;

  found=FALSE;
  i=0;
  while ((!found)&&(i<ts->n))
    {
      found=TransHasVar(v,ts->t[i]);
      i++;
    }
  return(found);
}

void UpdateUsedDOF(unsigned int *dof,TTransSeq *ts)
{
  unsigned int i;

  for(i=0;i<ts->n;i++)
    {
      switch (ts->t[i]->t)
        {
        case CT_TRANS:
          break;
        case TX:
          dof[TX]++;
          break;
        case TY:
          dof[TY]++;
          break;
        case TZ:
          dof[TZ]++;
          break;
        case TV:
          /* assuming a general displacement. */
          dof[TX]++;
          dof[TY]++;
          dof[TZ]++;
          break;
        case PA:
        case dPA_U:
        case dPA_V:
        case ddPA_UU:
        case ddPA_UV:
        case ddPA_VV:
          /* Assuming a general patch */
          dof[TX]++;
          dof[TY]++;
          dof[TZ]++;
          dof[RX]++;
          dof[RY]++;
          dof[RZ]++;
          break;
        case RX:
        case dRX:
        case ddRX:
          dof[RX]++;
          break;
        case RY:
        case dRY:
        case ddRY:
          dof[RY]++;
          break;
        case RZ:
        case dRZ:
        case ddRZ:
          dof[RZ]++;
          break;
        default:
          Error("Unkown transform type in UpdateUsedDOF");
        }
    }
}
 
void ShiftVariablesInTransSeq(unsigned int nv,TTransSeq *ts)
{
  unsigned int i;

   for(i=0;i<ts->n;i++)
    {
      if (TransHasVar(nv,ts->t[i]))
        Error("Removing a variable used in a transform sequence");
      else
        {
          if ((ts->t[i]->v!=NO_UINT)&&(ts->t[i]->v>nv))
            ts->t[i]->v--;
          if ((ts->t[i]->u!=NO_UINT)&&(ts->t[i]->u>nv))
            ts->t[i]->u--;
        }
    }
}

void SimplifyTransSeq(TTransSeq *ts)
{
  if ((ts->n>2)&&(ts->t[0]->t==CT_TRANS)&&(ts->t[ts->n-1]->t==CT_TRANS))
    {
      ts->n--;
      HTransformProduct(ts->t[ts->n]->ct,ts->t[0]->ct,ts->t[0]->ct);
      DeleteTrans(ts->t[ts->n]);
      free(ts->t[ts->n]);
    }
}

void DeriveTransSeq(unsigned int v,unsigned int *n,TTransSeq ***tsd,TTransSeq *ts)
{
  unsigned int i;
  
  *n=0;
  for(i=0;i<ts->n;i++)
    {
      if (TransHasVar(v,ts->t[i]))
        (*n)++;
    }
  
  /* *n will be typically 0 or 1 but it can be larger if a variable is repeated
     in a loop equation.  */
  if (*n==0)
    tsd=NULL;
  else
    {
      THTransform a;
      unsigned int r,t;
      double sgn;

      NEW(*tsd,*n,TTransSeq*);
      for(r=0;r<*n;r++)
        {
          /* We derive in turns for each time varible 'v' appears in the 
             sequence. '*n' is the maximum times and 'r' is the time with
             respect to which we are deriving now. 't' is the numbers
             of times the variable has been encounterd so far in the
             scan over the sequence. */
          t=0;
          NEW((*tsd)[r],1,TTransSeq);
          InitTransSeq((*tsd)[r]);
          for(i=0;i<ts->n;i++)
            {
              if ((!TransHasVar(v,ts->t[i]))||(t!=r))
                AddTrans2TransSeq(ts->t[i],(*tsd)[r]);
              else
                {
                  /* The r-th time variable 'v' appears in the sequence */
                  switch (ts->t[i]->t)
                    {
                    case TX:
                      HTransformZero(&a);
                      HTransformSetElement(AXIS_X,AXIS_H,(double)ts->t[i]->s,&a);
                      AddCtTrans2TransSeq(&a,(*tsd)[r]);
                      HTransformDelete(&a);
                      break;
                    case TY:
                      HTransformZero(&a);
                      HTransformSetElement(AXIS_Y,AXIS_H,(double)ts->t[i]->s,&a);
                      AddCtTrans2TransSeq(&a,(*tsd)[r]);
                      HTransformDelete(&a);
                      break;
                    case TZ:
                      HTransformZero(&a);
                      HTransformSetElement(AXIS_Z,AXIS_H,(double)ts->t[i]->s,&a);
                      AddCtTrans2TransSeq(&a,(*tsd)[r]);
                      HTransformDelete(&a);
                      break;
                    case TV:
                      sgn=(double)ts->t[i]->s;
                      HTransformZero(&a);
                      HTransformSetElement(AXIS_X,AXIS_H,sgn*ts->t[i]->vect[0],&a);
                      HTransformSetElement(AXIS_Y,AXIS_H,sgn*ts->t[i]->vect[1],&a);
                      HTransformSetElement(AXIS_Z,AXIS_H,sgn*ts->t[i]->vect[2],&a);
                      AddCtTrans2TransSeq(&a,(*tsd)[r]);
                      HTransformDelete(&a);
                      break;
                    case RX:
                      AddVarTrans2TransSeq(dRX,ts->t[i]->s,ts->t[i]->v,
                                           (*tsd)[r]);
                      break;
                    case RY:
                      AddVarTrans2TransSeq(dRY,ts->t[i]->s,ts->t[i]->v,
                                           (*tsd)[r]);
                      break;
                    case RZ:
                      AddVarTrans2TransSeq(dRZ,ts->t[i]->s,ts->t[i]->v,
                                           (*tsd)[r]);
                      break;
                    case PA:
                      /* For parametrized transforms, we keep all the 
                         information changing the type of transform */
                      if (v==ts->t[i]->u) /* Derivative w.r.t. the 1st param */
                        AddPatchTrans2TransSeq(dPA_U,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      else                /* Derivative w.r.t. the 2nd param */
                        AddPatchTrans2TransSeq(dPA_V,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      break;
                    case dRX:
                      AddVarTrans2TransSeq(ddRX,ts->t[i]->s,
                                           ts->t[i]->v,(*tsd)[r]);
                      break;
                    case dRY:
                      AddVarTrans2TransSeq(ddRY,ts->t[i]->s,
                                           ts->t[i]->v,(*tsd)[r]);
                      break;
                    case dRZ:
                      AddVarTrans2TransSeq(ddRZ,ts->t[i]->s,
                                           ts->t[i]->v,(*tsd)[r]);
                      break;
                    case dPA_U: 
                      if (v==ts->t[i]->u)  /* du du */
                        AddPatchTrans2TransSeq(ddPA_UU,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      else                 /* du dv */
                        AddPatchTrans2TransSeq(ddPA_UV,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      break;
                    case dPA_V: 
                      if (v==ts->t[i]->u)  /* dv du */
                        AddPatchTrans2TransSeq(ddPA_UV,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      else                 /* dv dv */
                        AddPatchTrans2TransSeq(ddPA_VV,ts->t[i]->s,ts->t[i]->u,
                                               ts->t[i]->v,ts->t[i]->p,(*tsd)[r]);
                      break;
                    default:
                      Error("Unknown transform type in DeriveTransSeq");
                    }
                }
              if (TransHasVar(v,ts->t[i]))
                t++;
            }
        }
    }
}

void EvaluateTransSeq(double *v,THTransform *r,TTransSeq *ts)
{
  unsigned int i;
  THTransform a;

  HTransformIdentity(r);
  for(i=0;i<ts->n;i++)
    {
       if (ts->t[i]->t==CT_TRANS)
         HTransformProduct(r,ts->t[i]->ct,r);
       else
         {
           double val;
           double d1,d2,d3;

           switch (ts->t[i]->t)
             {
             case TX:
               d1=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumTrans(d1,0,0,r);
               break;
             case TY:
               d1=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumTrans(0,d1,0,r);
               break;
             case TZ:
               d1=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumTrans(0,0,d1,r);
               break;
             case TV:
               val=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               d1=ts->t[i]->vect[0];
               d2=ts->t[i]->vect[1];
               d3=ts->t[i]->vect[2];
               HTransformAcumTrans(val*d1,val*d2,val*d3,r);
               break;      
             case RX:      
               val=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumRot(RX,sin(val),cos(val),r);
               break;
             case RY:
               val=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumRot(RY,sin(val),cos(val),r);
               break;
             case RZ:
               val=((double)(ts->t[i]->s))*v[ts->t[i]->v];
               HTransformAcumRot(RZ,sin(val),cos(val),r);
               break;
             case PA:
               EvaluatePATrans(v[ts->t[i]->u],v[ts->t[i]->v],&a,ts->t[i]);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case dRX:
               /* With negative sign we have to evaluate the derivative of the
                  inverse. However, for rotations this derivative can be
                  computed directly. */
               val=v[ts->t[i]->v];
               HTransformRx2(((double)(ts->t[i]->s))*cos(val),-sin(val),&a);
               HTransformSetElement(AXIS_X,AXIS_X,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case dRY:
               val=v[ts->t[i]->v];
               HTransformRy2(((double)(ts->t[i]->s))*cos(val),-sin(val),&a);
               HTransformSetElement(AXIS_Y,AXIS_Y,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case dRZ:
               val=v[ts->t[i]->v];
               HTransformRz2(((double)(ts->t[i]->s))*cos(val),-sin(val),&a);
               HTransformSetElement(AXIS_Z,AXIS_Z,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case dPA_U:
               {
                 double v1,v2;
                 double nv[3],dnv[3],nx,ny;
                 double x[3],y[3],t[3];
                 double dx[3],dy[3]={0,0,0},dz[3];
                 
                 v1=v[ts->t[i]->u];
                 v2=v[ts->t[i]->v];

                 /* Get the non-normalized x,y,t vectors (t is not used) */
                 EvaluateVectorsPATrans(v1,v2,x,y,t,ts->t[i]);

                 /* dx = diff(x_normalized,u) */
                 /*   x_norm x = norm(x)
                      For a generic vector'v' depending on a parameter 'u',
                      and with n_v=norm(v) (i.e., n_f is also funciton of 'u')
                        diff(v/n_v,u)=(1/n_v^2)*[diff(v,u)*n_v-v*diff(n_v,u)]
                      where
                        diff(n_v,u)=1/n_v*diff(v\tr*v,u)
                                   =1/n_v*v\tr*diff(v,u)
                      and, thus
                         diff(v/n_v,u)= (1/n_v^2)*[diff(v,u)*n_v
                                                   -v*1/n_v*v\tr*diff(v,u)]
                                          = diff(v,u)/n_v
                                            -(v/n_v)*(v\tr/n_v)*(diff(v,u)/n_v)

                      Thus, if v_norm=v/nv and dv_norm=diff(v,u)/norm(v)
                      then     diff(v/n_v,u)=dv_norm-(v_norm\tr*dv_norm)*v_norm

                      This is used every time we have to evaluate the 
                      differential of a normalized vector w.r.t. a given variable 
                      where we have the formula giving the  unormalized vector.
                  */
                 nx=Norm(3,x);
                 ScaleVector2(1/nx,3,x,nv);
                 /* ts->t[i]->p[5] = n1 = diff(x,u) */
                 ScaleVector2(1/nx,3,ts->t[i]->p[5],dnv);  
                 SumVectorScale(3,dnv,-DotProduct(nv,dnv),nv,dx);
                 
                 /* dy = diff(y_norm,u) is zero */

                 /* dz = diff(z_norm,u) with z_norm = x_norm X y_norm */
                 ny=Norm(3,y);
                 CrossProduct(dx,y,dz);
                 ScaleVector2(1/ny,3,dz,dz);

                 /* dt = diff(p,u) = p_02+v*d = y */

                 /* Set up the matrix */
                 if (ts->t[i]->s>0)
                   {
                     HTransformFromVectors(dx,dy,dz,y,&a);
                     HTransformSetElement(AXIS_H,AXIS_H,0,&a);
                   }
                 else
                   {
                     /* trans=T*R -> trans\inv=R\tr*(-T)=[R\tr -R\tr*t]
                        diff(trans\inv)=[diff(R\tr) -[diff(R\tr)*t+R\tr*diff(t)]]
                      */
                     THTransform Rt,dR;
                     double h[3],z[3],y1[3];
                     double zero[3]={0,0,0};

                     /* R\tr */
                     ScaleVector2(1/nx,3,x,x);
                     ScaleVector2(1/ny,3,y,y1);
                     CrossProduct(x,y1,z);
                     HTransformFromVectors(x,y1,z,zero,&Rt);
                     HTransformTranspose(&Rt,&Rt);
                     
                     /* a=diff(R\tr) = diff(R)\tr*/
                     HTransformFromVectors(dx,dy,dz,zero,&dR);
                     HTransformSetElement(AXIS_H,AXIS_H,0,&dR);
                     HTransformTranspose(&dR,&a);

                     /* Now set the homogeneous part of 'a' */
                     HTransformApplyRot(t,h,&a);  /* h=diff(R\tr)*t */
                     HTransformApplyRot(y,z,&Rt); /* z=R\tr*diff(t) */
                     AccumulateVector(3,z,h);     /* h=h+z */

                     HTransformSetElement(AXIS_X,AXIS_H,-h[0],&a);
                     HTransformSetElement(AXIS_Y,AXIS_H,-h[1],&a);
                     HTransformSetElement(AXIS_Z,AXIS_H,-h[2],&a);
                   }
                 HTransformProduct(r,&a,r);
                 HTransformDelete(&a);
               } 
               break;
             case dPA_V:
               {
                 double v1,v2;
                 double nv[3],dnv[3],nx,ny;
                 double x[3],y[3],t[3];
                 double dx[3],dy[3],dz[3],dt[3];
                 
                 v1=v[ts->t[i]->u];
                 v2=v[ts->t[i]->v];

                 /* Get the non-normalized x,y,t vectors (only x and y are used) */
                 EvaluateVectorsPATrans(v1,v2,x,y,t,ts->t[i]);

                 /* dx = diff(x_normalized,v) */
                 nx=Norm(3,x);
                 ScaleVector2(1/nx,3,x,nv);
                 /* ts->t[i]->p[6] = n2 = diff(x,v) */
                 ScaleVector2(1/nx,3,ts->t[i]->p[6],dnv);  
                 SumVectorScale(3,dnv,-DotProduct(nv,dnv),nv,dx);
                 
                 /* dy = diff(y_normalized,v) */
                 ny=Norm(3,y);
                 ScaleVector2(1/ny,3,y,nv);
                 /* ts->t[i]->p[3] = d = diff(y,v) */
                 ScaleVector2(1/ny,3,ts->t[i]->p[3],dnv); 
                 SumVectorScale(3,dnv,-DotProduct(nv,dnv),nv,dy);

                 /* dz = diff(z_norm,u) with z_norm = x_norm X y_norm */
                 /*  diff(z_norm,v) = diff(x_norm,v) X y/norm(y) 
                                      + x/norm(x) X diff(y_norm,v) */
                 CrossProduct(dx,y,dz);
                 ScaleVector2(1/ny,3,dz,dz);
                 CrossProduct(x,dy,nv);
                 SumVectorScale(3,dz,1/nx,nv,dz);

                 /* dt = diff(t,v) = p_01+u*d */
                 SumVectorScale(3,ts->t[i]->p[2],v1,ts->t[i]->p[3],dt);

                 if (ts->t[i]->s>0)
                   {
                     /* Set up the matrix */
                     HTransformFromVectors(dx,dy,dz,dt,&a);
                     HTransformSetElement(AXIS_H,AXIS_H,0,&a);
                   }
                 else
                   {
                     /* trans=T*R -> trans\inv=R\tr*(-T)=[R\tr -R\tr*t]
                        diff(trans\inv)=[diff(R\tr) -[diff(R\tr)*t+R\tr*diff(t)]]
                      */
                     THTransform Rt,dR;
                     double h[3],z[3];
                     double zero[3]={0,0,0};

                     /* R\tr */
                     ScaleVector2(1/nx,3,x,x);
                     ScaleVector2(1/ny,3,y,y);
                     CrossProduct(x,y,z);
                     HTransformFromVectors(x,y,z,zero,&Rt);
                     HTransformTranspose(&Rt,&Rt);
                     
                     /* a=diff(R\tr) = diff(R)\tr*/
                     HTransformFromVectors(dx,dy,dz,zero,&dR);
                     HTransformSetElement(AXIS_H,AXIS_H,0,&dR);
                     HTransformTranspose(&dR,&a);

                     /* Now set the homogeneous part of 'a' */
                     HTransformApplyRot(t,h,&a);   /* h=diff(R\tr)*t */
                     HTransformApplyRot(dt,z,&Rt); /* z=R\tr*diff(t) */
                     AccumulateVector(3,z,h);      /* h=h+z */

                     HTransformSetElement(AXIS_X,AXIS_H,-h[0],&a);
                     HTransformSetElement(AXIS_Y,AXIS_H,-h[1],&a);
                     HTransformSetElement(AXIS_Z,AXIS_H,-h[2],&a);
                   }
                 HTransformProduct(r,&a,r);
                 HTransformDelete(&a);
               }
               break;
             case ddRX:
               /* With negative sign we have to evaluate the second derivative of the
                  inverse. However, for rotations this derivative can be
                  computed directly. */
               val=v[ts->t[i]->v];
               HTransformRx2(-((double)(ts->t[i]->s))*sin(val),-cos(val),&a);
               HTransformSetElement(AXIS_X,AXIS_X,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case ddRY:
               val=v[ts->t[i]->v];
               HTransformRy2(-((double)(ts->t[i]->s))*sin(val),-cos(val),&a);
               HTransformSetElement(AXIS_Y,AXIS_Y,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case ddRZ:
               val=v[ts->t[i]->v];
               HTransformRz2(-((double)(ts->t[i]->s))*sin(val),-cos(val),&a);
               HTransformSetElement(AXIS_Z,AXIS_Z,0,&a);
               HTransformSetElement(AXIS_H,AXIS_H,0,&a);
               HTransformProduct(r,&a,r);
               HTransformDelete(&a);
               break;
             case ddPA_UU:
             case ddPA_UV:
             case ddPA_VV:
               Error("The evaluation of the second derivative of parametrized-patch transforms is not implemented yet");
               break;
             default:
               Error("Unknown transform type in EvaluateTransSeq");
             }
         }
    }
}

void PrintTransSeq(FILE *f,char **varNames,TTransSeq *ts)
{
  unsigned int i,j,k,t;
  double val;

  for(i=0;i<ts->n;i++)
    {
      t=ts->t[i]->t;
      if (t==CT_TRANS)
        {
          if (HTransformIsIdentity(ts->t[i]->ct))
            fprintf(f,"Id");
          else
            {
              //HTransformPrettyPrint(f,ts->t[i]->ct);
              
              fprintf(f,"T(");
              for(j=0;j<DIM_SP;j++)
                {
                  for(k=0;k<(DIM_SP+1);k++)
                    {
                      val=HTransformGetElement(j,k,ts->t[i]->ct);
                      PrintReal(f,val);
                      if (k<DIM_SP)
                        fprintf(f,",");
                    }
                  if (j<DIM_SP-1)
                    fprintf(f,";");
                }
              fprintf(f,")");
            }
        }
      else
        {
          switch (t)
            {
            case TX:
              fprintf(f,"Tx(");
              break;
            case TY:
              fprintf(f,"Ty(");
              break;
            case TZ:
              fprintf(f,"Tz(");
              break;
            case TV:
              fprintf(f,"Tv(");
              PrintReal(f,ts->t[i]->vect[0]);fprintf(f,",");
              PrintReal(f,ts->t[i]->vect[1]);fprintf(f,",");
              PrintReal(f,ts->t[i]->vect[2]);fprintf(f,",");
              break;
            case RX:          
              fprintf(f,"Rx(");
              break;
            case RY:
              fprintf(f,"Ry(");
              break;
            case RZ:
              fprintf(f,"Rz(");
              break;
            case PA:
            case dPA_U:
            case dPA_V:
            case ddPA_UU:
            case ddPA_UV:
            case ddPA_VV:
              {
                unsigned int k;
                double v[3];

                switch (t)
                  { 
                  case PA:
                    fprintf(f,"Pa(");
                    break;
                  case dPA_U:
                    fprintf(f,"dPa_u(");
                    break;
                  case dPA_V:
                    fprintf(f,"dPa_v(");
                    break;
                  case ddPA_UU:
                    fprintf(f,"ddPa_uu(");
                    break;
                  case ddPA_UV:
                    fprintf(f,"ddPa_uv(");
                    break;
                  case ddPA_VV:
                    fprintf(f,"ddPa_vv(");
                    break;
                  }     
                /* Print the 4 vectors defining the patch */
                if (t!=PA)
                  {
                    /* For derived patches we print the already
                       processed patch extremes +
                       Derived patches typically do not appear in
                       equation files */
                    for(j=0;j<4;j++)
                      {
                        for(k=0;k<3;k++)
                          {
                            PrintReal(f,ts->t[i]->p[j][k]);
                            if (k<2) fprintf(f,",");
                          }
                        fprintf(f,";");
                      }
                  }
                else
                  {
                    /* For oritinal patches we print the
                       orignal extremes. This facilitates the
                       unsderstanding of equation files. */
                    /* p_0 */
                    for(k=0;k<3;k++)
                      {
                        PrintReal(f,ts->t[i]->p[0][k]);
                        if (k<2) fprintf(f,",");
                      }
                    fprintf(f,";");
                    /* p_1=p_01+p_0 */
                    SumVector(3,ts->t[i]->p[2],ts->t[i]->p[0],v);
                    for(k=0;k<3;k++)
                      {
                        PrintReal(f,v[k]);
                        if (k<2) fprintf(f,",");
                      }
                    fprintf(f,";");
                    /* p_2=p_02+p_0 */
                    SumVector(3,ts->t[i]->p[1],ts->t[i]->p[0],v);
                    for(k=0;k<3;k++)
                      {
                        PrintReal(f,v[k]);
                        if (k<2) fprintf(f,",");
                      }
                  }
                fprintf(f,";");
                /* p_3=d+p_01+p_02+p_0 */
                SumVector(3,ts->t[i]->p[3],ts->t[i]->p[2],v);
                AccumulateVector(3,ts->t[i]->p[1],v);
                AccumulateVector(3,ts->t[i]->p[0],v);
                for(k=0;k<3;k++)
                  {
                    PrintReal(f,v[k]);
                    if (k<2) fprintf(f,",");
                  }
                fprintf(f,";");

                /* Now the variables */
                if (varNames!=NULL)
                  { 
                    char *n;
                    
                    n=varNames[ts->t[i]->u];
                    PRINT_VARIABLE_NAME(f,n);
                  }
                else
                  fprintf(f,"v_%u",ts->t[i]->u);
                fprintf(f,",");
              }
              break;
            case dRX:
              fprintf(f,"dRx(");
              break;
            case dRY:
              fprintf(f,"dRy(");
              break;
            case dRZ:
              fprintf(f,"dRz(");
              break;
            case ddRX:
              fprintf(f,"ddRx(");
              break;
            case ddRY:
              fprintf(f,"ddRy(");
              break;
            case ddRZ:
              fprintf(f,"ddRz(");
              break;
            default:
              Error("Unknown transform type in EvaluateTransSeq");
            }
          if (varNames!=NULL)
            { 
              char *n;

              n=varNames[ts->t[i]->v];
              PRINT_VARIABLE_NAME(f,n);
            }
          else
            fprintf(f,"v_%u",ts->t[i]->v);
          fprintf(f,")");
          if (ts->t[i]->s<0)
            fprintf(f,"^-1");
        }
      if (i<(ts->n-1))
        fprintf(f,"*");
    }
}

void DeleteTransSeq(TTransSeq *ts)
{
  ResetTransSeq(ts);
  free(ts->t);
}