/* -*-mode:c-*- */
/***************************************************************
*
*       Radiosity
*
*	Graphic driver for PostScript
*
*
***************************************************************/

#include <stdio.h>
#include <math.h>


#define SCREEN_WIDTH   (6.0*72)
#define SCREEN_HEIGHT  (4.8*72)
#define SCREEN_DEPTH   (65536)
#define ASPECT_RATIO  ((float)SCREEN_WIDTH/(float)SCREEN_HEIGHT)

#define PRE_CAT  (1)
#define POST_CAT (0)

#define DEFAULT_WINDOW_HEIGHT (2000.0)
#define DEFAULT_WINDOW_WIDTH  (DEFAULT_WINDOW_HEIGHT*ASPECT_RATIO)
#define DEFAULT_FRONT_PLANE_Z (2000.0)
#define DEFAULT_BACK_PLANE_Z  (-4000.0)
#define DEFAULT_PRP_Z         (10000.0)    /* Projection point Z coord. */


/*** Data structure used in radiosity algorithm ***/

typedef struct {
      float x, y, z ;
} Point ;

typedef struct {
      float r, g, b ;
} Rgb ;


typedef struct
{
      float v[4] ;                   /* x, y, z, and w */
} Vertex ;

typedef struct
{
      float m[4][4] ;                /* m[row][column], row vector assumed */
} Matrix ;


/**************************************************
*
*    Globals
*
***************************************************/

static Matrix  trans_mtx ;		/* WC -> DC */
static Vertex  prp ;			/* Projection point */
static Vertex  active_prp ;		/* Projection point in effect (WC) */
static float   view_rotx, view_roty ;	/* Viewing */
static float   view_zoom ;

static float   clip_right, clip_left ;	/* View volume (X) */
static float   clip_top, clip_bottom ;	/*             (Y) */
static float   clip_front, clip_back ;	/*             (Z) */


static FILE *ps_fd ;


void gset_unit_matrix(), gconcatenate_matrix(), gscale_matrix() ;
void gtranslate_matrix(), grotate_x_matrix(), grotate_y_matrix() ;
void grotate_z_matrix(), gtransform(), ginverse_matrix() ;


/************************************************************************
*
*    ps_open()
*    ps_close()
*
*************************************************************************/


int ps_open( file )
    char *file ;
{
      void init_transformation() ;
      void setup_transformation() ;

      if( (ps_fd = fopen( file, "w" )) == 0 )
      {
	    perror( file ) ;
	    return( 0 ) ;
      }

      /* Print out preamble */
      fprintf( ps_fd, "%%!PS-Adobe-1.0\n" ) ;
      fprintf( ps_fd, "%%%%EndComments\n" ) ;
      fprintf( ps_fd, "%%%%Pages: 1\n" ) ;
      fprintf( ps_fd, "%%%%EndProlog\n" ) ;
      fprintf( ps_fd, "%%%%Page: 1 1\n" ) ;
      fprintf( ps_fd, "\n" ) ;

      /* Default line cap/join */
      fprintf( ps_fd, "1 setlinecap 1 setlinejoin\n" ) ;

      /* Initialize transformation */
      init_transformation() ;
      setup_transformation() ;
}


void ps_close()
{
      if( ps_fd == 0 )
	    return ;


      fprintf( ps_fd, "showpage\n" ) ;
      fprintf( ps_fd, "%%%%Trailer\n" ) ;
      fclose( ps_fd ) ;

      ps_fd = 0 ;
}


/**************************************************
*
*    ps_linewidth()
*
***************************************************/

void ps_linewidth( w )
   
    float w ;
{
      if( ps_fd == 0 )
	    return ;

      fprintf( ps_fd, "%f setlinewidth\n", w ) ;
}


/**************************************************
*
*    ps_line()
*
***************************************************/

void ps_line( p1, p2 )

    Point *p1, *p2 ;
{
      Vertex  v1, v2 ;
      float x1, y1, x2, y2 ;

      if( ps_fd == 0 )
	    return ;

      v1.v[0] = p1->x ;  v1.v[1] = p1->y ;  v1.v[2] = p1->z ;  v1.v[3] = 1.0 ;
      v2.v[0] = p2->x ;  v2.v[1] = p2->y ;  v2.v[2] = p2->z ;  v2.v[3] = 1.0 ;
      gtransform( &v1, &v1, &trans_mtx ) ;
      gtransform( &v2, &v2, &trans_mtx ) ;
      x1 = v1.v[0] / v1.v[3] ;
      y1 = v1.v[1] / v1.v[3] ;
      x2 = v2.v[0] / v2.v[3] ;
      y2 = v2.v[1] / v2.v[3] ;


      fprintf( ps_fd, "newpath\n%f %f moveto\n", x1, y1 ) ;
      fprintf( ps_fd, "%f %f lineto\nstroke\n", x2, y2 ) ;
}


/**************************************************
*
*    ps_polygonedge()
*
***************************************************/

void ps_polygonedge( n, p_list )

    int n ;
    Point *p_list ;
{
      float dcx, dcy ;
      Vertex v ;
      int i ;

      if( ps_fd == 0 )
	    return ;

      /* Transform */
      v.v[0] = p_list[0].x ;
      v.v[1] = p_list[0].y ;
      v.v[2] = p_list[0].z ;
      v.v[3] = 1.0 ;
      gtransform( &v, &v, &trans_mtx ) ;
      dcx = v.v[0] / v.v[3] ; 
      dcy = v.v[1] / v.v[3] ; 
      fprintf( ps_fd, "newpath\n%f %f moveto\n", dcx, dcy ) ;

      for( i = 1 ; i < n ; i++ )
      {
	    /* Transform */
	    v.v[0] = p_list[i].x ;
	    v.v[1] = p_list[i].y ;
	    v.v[2] = p_list[i].z ;
	    v.v[3] = 1.0 ;
	    gtransform( &v, &v, &trans_mtx ) ;
	    dcx = v.v[0] / v.v[3] ; 
	    dcy = v.v[1] / v.v[3] ; 

	    fprintf( ps_fd, "%f %f lineto\n", dcx, dcy ) ;
      }
      
      fprintf( ps_fd, "closepath stroke\n" ) ;
}


/**************************************************
*
*    ps_polygon()
*
***************************************************/

void ps_polygon( n, p_list )

    int n ;
    Point *p_list ;
{
      float dcx, dcy ;
      Vertex v ;
      int i ;

      if( ps_fd == 0 )
	    return ;

      /* Transform */
      v.v[0] = p_list[0].x ;
      v.v[1] = p_list[0].y ;
      v.v[2] = p_list[0].z ;
      v.v[3] = 1.0 ;
      gtransform( &v, &v, &trans_mtx ) ;
      dcx = v.v[0] / v.v[3] ; 
      dcy = v.v[1] / v.v[3] ; 
      fprintf( ps_fd, "newpath\n%f %f moveto\n", dcx, dcy ) ;

      for( i = 1 ; i < n ; i++ )
      {
	    /* Transform */
	    v.v[0] = p_list[i].x ;
	    v.v[1] = p_list[i].y ;
	    v.v[2] = p_list[i].z ;
	    v.v[3] = 1.0 ;
	    gtransform( &v, &v, &trans_mtx ) ;
	    dcx = v.v[0] / v.v[3] ; 
	    dcy = v.v[1] / v.v[3] ; 

	    fprintf( ps_fd, "%f %f lineto\n", dcx, dcy ) ;
      }
      
      fprintf( ps_fd, "closepath fill\n" ) ;
}


/**************************************************
*
*    ps_spolygon()
*
***************************************************/

void ps_spolygon( n, p_list, c_list )

    int n ;
    Point *p_list ;
    Rgb   *c_list ;
{
      float dcx, dcy ;
      Vertex v ;
      int i ;
      float gray_scale ;

      if( ps_fd == 0 )
	    return ;

      /* Transform */
      v.v[0] = p_list[0].x ;
      v.v[1] = p_list[0].y ;
      v.v[2] = p_list[0].z ;
      v.v[3] = 1.0 ;
      gtransform( &v, &v, &trans_mtx ) ;
      dcx = v.v[0] / v.v[3] ; 
      dcy = v.v[1] / v.v[3] ; 
      fprintf( ps_fd, "newpath\n%f %f moveto\n", dcx, dcy ) ;

      for( i = 1 ; i < n ; i++ )
      {
	    /* Transform */
	    v.v[0] = p_list[i].x ;
	    v.v[1] = p_list[i].y ;
	    v.v[2] = p_list[i].z ;
	    v.v[3] = 1.0 ;
	    gtransform( &v, &v, &trans_mtx ) ;
	    dcx = v.v[0] / v.v[3] ; 
	    dcy = v.v[1] / v.v[3] ; 

	    fprintf( ps_fd, "%f %f lineto\n", dcx, dcy ) ;
      }
      
      gray_scale = c_list[0].g ;
      if( gray_scale > 1.0 )
	    gray_scale = 1.0 ;
      else if( gray_scale < 0.0 )
	    gray_scale = 0.0 ;
		  
      fprintf( ps_fd, "closepath %f setgray fill\n", gray_scale ) ;
}


/**************************************************
*
*    ps_clear()
*
***************************************************/

void ps_clear()
{
}


/**************************************************
*
*    ps_setup_view()
*
***************************************************/

void ps_setup_view( rot_x, rot_y, dist, zoom )

    float rot_x, rot_y, dist, zoom ;
{
      void setup_transformation() ;

      prp.v[0] = 0.0 ;
      prp.v[1] = 0.0 ;
      prp.v[2] = (float)dist ;
      prp.v[3] = 0.0 ;
      view_rotx = rot_x ;
      view_roty = rot_y ;
      view_zoom = zoom  ;

      setup_transformation() ;
}


/**************************************************
*
*    setup_transformation()
*
***************************************************/

void setup_transformation()
{
      float cf_z, cb_z ;
      int light ;
      Matrix pmat ;
      
      /* Set to unit matrix */
      gset_unit_matrix( &trans_mtx ) ;

      /* View orientation matrix */
      grotate_x_matrix( POST_CAT, &trans_mtx, view_rotx ) ;
      grotate_y_matrix( POST_CAT, &trans_mtx, view_roty ) ;

      /* Compute active (currently effective) projection point */
      ginverse_matrix( &pmat, &trans_mtx ) ;
      gtransform( &active_prp, &prp, &pmat ) ;

      /* Perspective projection */
      gset_unit_matrix( &pmat ) ;
      pmat.m[2][3] = - 1 / prp.v[2] ;
      gconcatenate_matrix( POST_CAT, &trans_mtx, &pmat ) ;

      cf_z = prp.v[2] * clip_front / ( prp.v[2] - clip_front ) ;
      cb_z = prp.v[2] * clip_back  / ( prp.v[2] - clip_back  ) ;

      /* Window-Viewport */
      gscale_matrix( POST_CAT, &trans_mtx,
		   (float)SCREEN_WIDTH  / (clip_right - clip_left),
		   (float)SCREEN_HEIGHT / (clip_top - clip_bottom),
		   (float)SCREEN_DEPTH  / (cf_z - cb_z) ) ;

      gtranslate_matrix( POST_CAT, &trans_mtx,
	      -(float)SCREEN_WIDTH * clip_left / (clip_right - clip_left),
	      -(float)SCREEN_HEIGHT* clip_top  / (clip_bottom - clip_top),
	      -(float)SCREEN_DEPTH * cb_z / (cf_z - cb_z) ) ;

      gtranslate_matrix( POST_CAT, &trans_mtx,
			(float)(1.0*72), (float)(0.5*72), 0 ) ;
}


/**************************************************
*
*    init_transformation()
*
***************************************************/

void init_transformation()
{
      /* Initialize matrix, just in case */
      gset_unit_matrix( &trans_mtx ) ;

      /* Initialize Projection point */
      prp.v[0] = 0.0 ;
      prp.v[1] = 0.0 ;
      prp.v[2] = DEFAULT_PRP_Z ;
      prp.v[3] = 0.0 ;
 
      /* Viewing */
      view_rotx = view_roty = 0.0 ;
      view_zoom = 1.0 ;
  
      /* Initialize view volume boundary */
      clip_right =  DEFAULT_WINDOW_WIDTH / 2.0 ;
      clip_left  = -DEFAULT_WINDOW_WIDTH / 2.0 ;
      clip_top   =  DEFAULT_WINDOW_HEIGHT / 2.0 ;
      clip_bottom= -DEFAULT_WINDOW_HEIGHT / 2.0 ; 
      clip_front =  DEFAULT_FRONT_PLANE_Z ;
      clip_back  =  DEFAULT_BACK_PLANE_Z ;
}


/********************************************
*
*    set_matrix()
*
*********************************************/

void  gset_matrix( mtx, m11, m12, m13, m14,  m21, m22, m23, m24,
		       m31, m32, m33, m34,  m41, m42, m43, m44 )

    Matrix *mtx ;
    float m11, m12, m13, m14 ;
    float m21, m22, m23, m24 ;
    float m31, m32, m33, m34 ;
    float m41, m42, m43, m44 ;
{
      mtx->m[0][0] = m11 ; mtx->m[0][1] = m12 ;
      mtx->m[0][2] = m13 ; mtx->m[0][3] = m14 ;

      mtx->m[1][0] = m21 ; mtx->m[1][1] = m22 ;
      mtx->m[1][2] = m23 ; mtx->m[1][3] = m24 ;

      mtx->m[2][0] = m31 ; mtx->m[2][1] = m32 ;
      mtx->m[2][2] = m33 ; mtx->m[2][3] = m34 ;

      mtx->m[3][0] = m41 ; mtx->m[3][1] = m42 ;
      mtx->m[3][2] = m43 ; mtx->m[3][3] = m44 ;
}    


/********************************************
*
*    set_unit_matrix()
*
*********************************************/

void  gset_unit_matrix( mtx )

    Matrix *mtx ;
{
      int  row, col ;
      
      /* Clear the matrix */
      for( row = 0 ; row < 4 ; row++ )
	    for( col = 0 ; col < 4 ; col++ )
		  mtx->m[row][col] = 0.0 ;

      /* Set 1.0s along diagonal line */
      for( row = 0 ; row < 4 ; row++ )
	    mtx->m[row][row] = 1.0 ;
}    


/********************************************
*
*    concatenate_matrix()
*
*          m1 <- m1 * m2 (precat = 1)
*          m1 <- m2 * m1 (precat = 0)
*
*********************************************/

void  gconcatenate_matrix( precat, m1, m2 )

    int precat ;
    Matrix *m1, *m2 ;
{
      int  row, col, scan ;
      Matrix *dest ;
      Matrix temp ;
      

      /* Swap pointer according to the concatenation mode */
      dest = m1 ;
      if( precat == 1 )
      {
	    m1 = m2 ;
	    m2 = dest ;
      }
      
      /* concatenate it */
      for( row = 0 ; row < 4 ; row++ )
	    for( col = 0 ; col < 4 ; col++ ) 
	    {
		  temp.m[row][col] = 0.0 ;
		  for( scan = 0 ; scan < 4 ; scan++ )
			temp.m[row][col] +=
			      m1->m[row][scan] * m2->m[scan][col];
	    }

      *dest = temp ;
}


/********************************************
*
*    scale_matrix()
*
*          m1 <- SCALE * m1 (precat = 1)
*          m1 <- m1 * SCALE (precat = 0)
*
*********************************************/

void  gscale_matrix( precat, m1, sx, sy, sz )

    int precat ;
    Matrix *m1 ;
    float sx, sy, sz ;
{
      Matrix smat ;
      
      /* Initialize to unit matrix */
      gset_unit_matrix( &smat ) ;

      /* Set scale values */
      smat.m[0][0] = sx ;
      smat.m[1][1] = sy ;
      smat.m[2][2] = sz ;

      /* concatenate */
      gconcatenate_matrix( precat, m1, &smat ) ;
}


/********************************************
*
*    translate_matrix()
*
*          m1 <- T * m1 (precat = 1)
*          m1 <- m1 * T (precat = 0)
*
*********************************************/

void  gtranslate_matrix( precat, m1, tx, ty, tz )

    int precat ;
    Matrix *m1 ;
    float tx, ty, tz ;
{
      Matrix tmat ;
      
      /* Initialize to unit matrix */
      gset_unit_matrix( &tmat ) ;

      /* Set scale values */
      tmat.m[3][0] = tx ;
      tmat.m[3][1] = ty ;
      tmat.m[3][2] = tz ;

      /* concatenate */
      gconcatenate_matrix( precat, m1, &tmat ) ;
}


/********************************************
*
*    rotate_x_matrix()
*    rotate_y_matrix()
*    rotate_z_matrix()
*
*          m1 <- ROT * m1 (precat = 1)
*          m1 <- m1 * ROT (precat = 0)
*
*********************************************/

void  grotate_x_matrix( precat, m1, rot )

    int precat ;
    Matrix *m1 ;
    float rot ;
{
      Matrix rmat ;
      float s_val, c_val ;
      
      /* Initialize to unit matrix */
      gset_unit_matrix( &rmat ) ;

      /* Set scale values */
      s_val = sin( rot * M_PI / 180.0 ) ;
      c_val = cos( rot * M_PI / 180.0 ) ;
      rmat.m[1][1] = c_val ;
      rmat.m[1][2] = s_val ;
      rmat.m[2][1] = -s_val ;
      rmat.m[2][2] = c_val ;

      /* concatenate */
      gconcatenate_matrix( precat, m1, &rmat ) ;
}


void  grotate_y_matrix( precat, m1, rot )

    int precat ;
    Matrix *m1 ;
    float rot ;
{
      Matrix rmat ;
      float s_val, c_val ;
      
      /* Initialize to unit matrix */
      gset_unit_matrix( &rmat ) ;

      /* Set scale values */
      s_val = sin( rot * M_PI / 180.0 ) ;
      c_val = cos( rot * M_PI / 180.0 ) ;
      rmat.m[0][0] = c_val ;
      rmat.m[0][2] = -s_val ;
      rmat.m[2][0] = s_val ;
      rmat.m[2][2] = c_val ;

      /* concatenate */
      gconcatenate_matrix( precat, m1, &rmat ) ;
}


void  grotate_z_matrix( precat, m1, rot )

    int precat ;
    Matrix *m1 ;
    float rot ;
{
      Matrix rmat ;
      float s_val, c_val ;
      
      /* Initialize to unit matrix */
      gset_unit_matrix( &rmat ) ;

      /* Set scale values */
      s_val = sin( rot * M_PI / 180.0 ) ;
      c_val = cos( rot * M_PI / 180.0 ) ;
      rmat.m[0][0] = c_val ;
      rmat.m[0][1] = s_val ;
      rmat.m[1][0] = -s_val ;
      rmat.m[1][1] = c_val ;

      /* concatenate */
      gconcatenate_matrix( precat, m1, &rmat ) ;
}


/********************************************
*
*    transform()
*
*          v1 <- v2 * mtx
*
*********************************************/

void gtransform( v1, v2, mtx )

    Vertex *v1, *v2 ;
    Matrix *mtx ;
{
      float x, y, z, w ;

      x  = v2->v[0] * mtx->m[0][0] ;
       y  = v2->v[0] * mtx->m[0][1] ;
        z  = v2->v[0] * mtx->m[0][2] ;
         w  = v2->v[0] * mtx->m[0][3] ;

      x += v2->v[1] * mtx->m[1][0] ;
       y += v2->v[1] * mtx->m[1][1] ;
        z += v2->v[1] * mtx->m[1][2] ;
         w += v2->v[1] * mtx->m[1][3] ;

      x += v2->v[2] * mtx->m[2][0] ;
       y += v2->v[2] * mtx->m[2][1] ;
        z += v2->v[2] * mtx->m[2][2] ;
         w += v2->v[2] * mtx->m[2][3] ;

      x += v2->v[3] * mtx->m[3][0] ;
       y += v2->v[3] * mtx->m[3][1] ;
        z += v2->v[3] * mtx->m[3][2] ;
         w += v2->v[3] * mtx->m[3][3] ;

      v1->v[0] = x ;
       v1->v[1] = y ;
        v1->v[2] = z ;
         v1->v[3] = w ;
}


/********************************************
*
*    inverse_matrix()
*
*          m1 <- inv(m2)
*
*********************************************/


void ginverse_matrix( m1, m2 )

    Matrix *m1, *m2 ;
{
      double detval ;
      double det(), cdet() ;
      
      /* det(m2) */
      detval = det( m2 ) ;

      /* Clamel's solution */
      m1->m[0][0] =  cdet( m2, 1,2,3, 1,2,3 ) / detval ;
      m1->m[0][1] = -cdet( m2, 0,2,3, 1,2,3 ) / detval ;
      m1->m[0][2] =  cdet( m2, 0,1,3, 1,2,3 ) / detval ;
      m1->m[0][3] = -cdet( m2, 0,1,2, 1,2,3 ) / detval ;

      m1->m[1][0] = -cdet( m2, 1,2,3, 0,2,3 ) / detval ;
      m1->m[1][1] =  cdet( m2, 0,2,3, 0,2,3 ) / detval ;
      m1->m[1][2] = -cdet( m2, 0,1,3, 0,2,3 ) / detval ;
      m1->m[1][3] =  cdet( m2, 0,1,2, 0,2,3 ) / detval ;

      m1->m[2][0] =  cdet( m2, 1,2,3, 0,1,3 ) / detval ;
      m1->m[2][1] = -cdet( m2, 0,2,3, 0,1,3 ) / detval ;
      m1->m[2][2] =  cdet( m2, 0,1,3, 0,1,3 ) / detval ;
      m1->m[2][3] = -cdet( m2, 0,1,2, 0,1,3 ) / detval ;

      m1->m[3][0] = -cdet( m2, 1,2,3, 0,1,2 ) / detval ;
      m1->m[3][1] =  cdet( m2, 0,2,3, 0,1,2 ) / detval ;
      m1->m[3][2] = -cdet( m2, 0,1,3, 0,1,2 ) / detval ;
      m1->m[3][3] =  cdet( m2, 0,1,2, 0,1,2 ) / detval ;
}


double det( m )

    Matrix *m ;
{
      double det_sum ;
      double cdet() ;
      
      /* Expand with respect to column 4 */
      det_sum = 0.0 ;
      if( m->m[0][3] != 0.0 )
	    det_sum -= m->m[0][3] * cdet( m, 1, 2, 3,  0, 1, 2 ) ;
      if( m->m[1][3] != 0.0 )
	    det_sum += m->m[1][3] * cdet( m, 0, 2, 3,  0, 1, 2 ) ;
      if( m->m[2][3] != 0.0 )
	    det_sum -= m->m[2][3] * cdet( m, 0, 1, 3,  0, 1, 2 ) ;
      if( m->m[3][3] != 0.0 )
	    det_sum += m->m[3][3] * cdet( m, 0, 1, 2,  0, 1, 2 ) ;

      return( det_sum ) ;
}


double cdet( m, r0, r1, r2, c0, c1, c2 )

    Matrix *m ;
    int r0, r1, r2, c0, c1, c2 ;
{
        double temp ;

        temp  = m->m[r0][c0] * m->m[r1][c1] * m->m[r2][c2] ;
        temp += m->m[r1][c0] * m->m[r2][c1] * m->m[r0][c2] ;
        temp += m->m[r2][c0] * m->m[r0][c1] * m->m[r1][c2] ;

        temp -= m->m[r2][c0] * m->m[r1][c1] * m->m[r0][c2] ;
        temp -= m->m[r1][c0] * m->m[r0][c1] * m->m[r2][c2] ;
        temp -= m->m[r0][c0] * m->m[r2][c1] * m->m[r1][c2] ;

        return( temp ) ;
}


/********************************************
*
*    normalize_vector()
*
*          v1 <- normalized( v2 )
*
*          W component is ignored.
*
*********************************************/


void gnormalize_vector( v1, v2 )

    Vertex *v1, *v2 ;
{
      float t0, t1, t2 ;

      t0 = v2->v[0] * v2->v[0] ;
       t1 = v2->v[1] * v2->v[1] ;
        t2 = v2->v[2] * v2->v[2] ;

      t0 = 1.0 / sqrt( t0 + t1 + t2 ) ;

      v1->v[0] = v2->v[0] * t0 ;
       v1->v[1] = v2->v[1] * t0 ;
        v1->v[2] = v2->v[2] * t0 ;
}


/********************************************
*
*    inner_product()
*
*          (v1.v2) <- inner_product( v1, v2 )
*
*          W component is ignored.
*
*********************************************/

float ginner_product( v1, v2 )

    Vertex *v1, *v2 ;
{
      float ip ;
      
      ip  = v1->v[0] * v2->v[0] ;
      ip += v1->v[1] * v2->v[1] ;
      ip += v1->v[2] * v2->v[2] ;

      return( ip ) ;
}


/********************************************
*
*    print_vector()
*
*
*********************************************/

void gprint_vector( v )

    Vertex *v ;
{
      printf( "(%g,%g,%g,%g)\n", v->v[0], v->v[1], v->v[2], v->v[3] ) ;
}