C++ VSET函数代码示例

double checkReactions(gsl_vector *haz, gsl_vector *residuals, gsl_vector *fast_params, 
                      double deltat,  gsl_vector *sps, gsl_matrix *PostPre)
{
  int i;
  double small_residual = -1000000.0;

  double total_hazard = 0.0;
  for(i=0; i < fast_params->size; i++) {
    if(checkHazard(haz, PostPre, sps, i, deltat) && 
       checkSpecies(PostPre, sps, i))  {
      VSET(residuals, i, small_residual);
    } else if(VGET(residuals, i) < (small_residual + 10000))  {
      total_hazard += VGET(haz, i);
      VSET(residuals, i, log(gsl_ran_flat(r, 0.0, 1.0)));
      VSET(fast_params, i, 0.0);
      //      running_totals1[i]++;
    } else {
      total_hazard += VGET(haz, i);
      VSET(fast_params, i, 0.0);
      //      running_totals1[i]++;
    }
    //printf("%d ", running_totals1[i]);
  }
  //  printf("\n");
  return total_hazard;        
}

void
render_flos_work(render_t *render, struct tie_s *tie, struct tie_ray_s *ray, vect_t *pixel) {
    struct tie_id_s id, tid;
    vect_t vec;
    fastf_t angle;
    struct render_flos_s *rd;

    rd = (struct render_flos_s *)render->data;

    if (tie_work(tie, ray, &id, render_hit, NULL) != NULL) {
	VSET(*pixel, 0.0, 0.5, 0.0);
    } else
	return;

    VSUB2(vec, ray->pos, id.pos);
    VUNITIZE(vec);
    angle = VDOT(vec, id.norm);

    /* Determine if direct line of sight to fragment */
    VMOVE(ray->pos, rd->frag_pos);
    VSUB2(ray->dir, id.pos, rd->frag_pos);
    VUNITIZE(ray->dir);

    if (tie_work(tie, ray, &tid, render_hit, NULL)) {
	if (fabs (id.pos[0] - tid.pos[0]) < TIE_PREC
	    && fabs (id.pos[1] - tid.pos[1]) < TIE_PREC
	    && fabs (id.pos[2] - tid.pos[2]) < TIE_PREC)
	{
	    VSET(*pixel, 1.0, 0.0, 0.0);
	}
    }

    VSCALE(*pixel, *pixel, (0.5+angle*0.5));
}

double Atmosphere::GADD( lwpp::VolumetricAccess &va, double pp[ 3 ], double y, double falloff,
   int use_lighting, int atmos_type, double li[ 3 ] )
{
   LWItemID light = iteminfo->first( LWI_LIGHT, NULL );
   double ilu[ 3 ], ldir[ 3 ];

   /* accumulate contributions from each light */

   if ( use_lighting && ( va.getFlags() & LWVEF_COLOR )) {
      VSET( li, 0.0 );
      while ( light ) {
         va.illuminate( light, pp, ldir, ilu );
         VADD( li, ilu );
         light = iteminfo->next( light );
      }
   }
   else
      VSET( li, 1.0 );

   /* return GADD (geometric atmospheric density distribution) */

   if ( atmos_type == 1 )
      return exp( -falloff * y );
   else
      return 1.0;
}

/*
  Hybrid method
*/
void updateHazard(gsl_vector *haz, gsl_vector *params, gsl_vector *sps)
{
  VSET(haz, 0, VGET(params, 0));
  VSET(haz, 1, VGET(params, 1));
  VSET(haz, 2, VGET(params, 2)*VGET(sps, 0));        
  VSET(haz, 3, VGET(params, 3)*VGET(sps, 1));        
  VSET(haz, 4, VGET(params, 4)*VGET(sps, 0)*VGET(sps, 1)); 
}           

void
PL_FORTRAN(f3vect, F3VECT)(FILE **fp, float *fx, float *fy, float *fz, float *tx, float *ty, float *tz, float *fl , float *tl)
{
    point_t from, to;

    VSET(from, *fx, *fy, *fz);
    VSET(to, *tx, *ty, *tz);
    tp_3vector(*fp, from, to, *fl, *tl);
}

void mvn_sample(gsl_vector *mean_cand, gsl_matrix *var)
{
  /* Takes a mean vec, mean and var matrix, 
   * var and gives vector of MVN(mean,var) realisations, x 
   */
  int i, j;
  int dimen = var -> size1;
  double value;
  gsl_matrix *disp;
  gsl_vector *ran;
  gsl_matrix *fast_species;
  
  fast_species = gsl_matrix_alloc(2, 2);
  gsl_matrix_set_identity(fast_species);
  
  for(i=0;i<dimen; i++) {
    if(MGET(var, i, i) <0.00000000001) {
      MSET(var, i, i, 1.0);
      MSET(fast_species, i, i, 0.0);
    }
  }
  
  disp = gsl_matrix_alloc(2, 2);
  ran = gsl_vector_alloc(2);
  gsl_matrix_memcpy(disp, var);
  if(postive_definite == 1) {
    gsl_linalg_cholesky_decomp(disp);
    for(i=0;i<dimen;i++) {
      for (j=i+1;j<dimen;j++) {
        MSET(disp,i,j,0.0);
      }
    }
  }else{
    value = pow(MGET(disp, 0 ,0), 0.5);
    gsl_matrix_set_identity(disp);
    MSET(disp, 0,0, value);
    MSET(disp, 1,1, value);       
  }

  for (j=0;j<dimen;j++) {
    VSET(ran,j,gsl_ran_gaussian(r,1.0));
  }

  /*remove update from slow species*/
  gsl_matrix_mul_elements(disp, fast_species);
    
  /*Add noise to mean cand*/
  gsl_blas_dgemv(CblasNoTrans,1.0, disp, ran, 1.0, mean_cand);
  for(i=0; i<2; i++)  {
    if(VGET(mean_cand,i)<=0.0001 && MGET(fast_species, i, i) > 0.000001)
      VSET(mean_cand,i,0.0001);
  }
  gsl_vector_free(ran);
  gsl_matrix_free(disp);
  gsl_matrix_free(fast_species);
}

static void SuperQ_Event( sqData *tool, int code )
{
   switch ( code )
   {
      case LWT_EVENT_DROP:

      /* The drop action is caused when the user clicks in the blank
         area of the display or uses the keyboard equivalent.  If the
         tool is active, we force a rejection of any interactive
         action partly complete.  For inactive tools we drop through
         to the reset action. */

      if ( tool->active ) {
         tool->update = LWT_TEST_REJECT;
         break;
      }

      case LWT_EVENT_RESET:

      /* The reset action corresponds to the reset command on the
         numeric panel, and causes us to snap the SuperQ factor back
         to zero.  Resets are also implicit when the user drops an
         inactive tool, thus the passthru above. */

         tool->nsides = 24;
         tool->nsegments = 12;
         VSET( tool->org, 0.0 );
         VSET( tool->center, 0.0f );
         VSET( tool->holex, 0.0f );
         VSET( tool->holez, 0.0f );
         VSET( tool->rad, 0.5 );
         tool->shape = 0;
         tool->bf1 = 2.0;
         tool->bf2 = 2.0;
         tool->axis = 1;
         tool->diam = 0.50;
         tool->holex[ 0 ] = ( float )( tool->diam * 0.5 );
         tool->holez[ 2 ] = tool->holex[ 0 ];
         tool->update = LWT_TEST_UPDATE;
         tool->dirty = 1;
         setBoxCorners( tool );
         break;

      case LWT_EVENT_ACTIVATE:

      /* Activation can be triggered from the numeric window or with a
         keystroke, and it should restart the edit operation with its
         current settings. */

         tool->update = LWT_TEST_UPDATE;
         tool->active = 1;
         tool->dirty = 1;
         break;
   }
}

void
bn_mat_lookat(mat_t rot, const vect_t dir, int yflip)
{
    mat_t first;
    mat_t second;
    mat_t prod12;
    mat_t third;
    vect_t x;
    vect_t z;
    vect_t t1;
    fastf_t hypot_xy;
    vect_t xproj;
    vect_t zproj;

    /* First, rotate D around Z axis to match +X axis (azimuth) */
    hypot_xy = hypot(dir[X], dir[Y]);
    bn_mat_zrot(first, -dir[Y] / hypot_xy, dir[X] / hypot_xy);

    /* Next, rotate D around Y axis to match -Z axis (elevation) */
    bn_mat_yrot(second, -hypot_xy, -dir[Z]);
    bn_mat_mul(prod12, second, first);

    /* Produce twist correction, by re-orienting projection of X axis */
    VSET(x, 1, 0, 0);
    MAT4X3VEC(xproj, prod12, x);
    hypot_xy = hypot(xproj[X], xproj[Y]);
    if (hypot_xy < 1.0e-10) {
	bu_log("Warning: bn_mat_lookat:  unable to twist correct, hypot=%g\n", hypot_xy);
	VPRINT("xproj", xproj);
	MAT_COPY(rot, prod12);
	return;
    }
    bn_mat_zrot(third, -xproj[Y] / hypot_xy, xproj[X] / hypot_xy);
    bn_mat_mul(rot, third, prod12);

    if (yflip) {
	VSET(z, 0, 0, 1);
	MAT4X3VEC(zproj, rot, z);
	/* If original Z inverts sign, flip sign on resulting Y */
	if (zproj[Y] < 0.0) {
	    MAT_COPY(prod12, rot);
	    MAT_IDN(third);
	    third[5] = -1;
	    bn_mat_mul(rot, third, prod12);
	}
    }

    /* Check the final results */
    MAT4X3VEC(t1, rot, dir);
    if (t1[Z] > -0.98) {
	bu_log("Error:  bn_mat_lookat final= (%g, %g, %g)\n", V3ARGS(t1));
    }
}

/* Make a leaf node out of an ARB4 */
void
do_leaf(char *name)
{
    point_t pt[4];

    VSET(pt[0], 0, 0, 0);
    VSET(pt[1], 100, 0, 0);
    VSET(pt[2], 50, 100*sin60, 0);
    VSET(pt[3], 50, 100*sin60/3, 100*sin60);

    mk_arb4(outfp, name, &pt[0][X]);
}

st_part_at * timingsInitPrior(double X, double Y)
{
  st_part_at *part;
  part = (st_part_at *) malloc(sizeof(st_part_at));
  part->params = gsl_vector_alloc(5);
  part->sps = gsl_vector_alloc(2);
  part->res = gsl_vector_alloc(5);

  VSET(part->sps, 0, X);
  VSET(part->sps, 1, Y);
  return(part);
}    

/**
 * Calculate a bounding RPP around a metaball
 */
int
rt_metaball_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol))
{
    struct rt_metaball_internal *mb;
    point_t center;
    fastf_t radius;
    mb = (struct rt_metaball_internal *)ip->idb_ptr;
    RT_METABALL_CK_MAGIC(mb);
    VSETALL(center, 0);

    radius = rt_metaball_get_bounding_sphere(&center, mb->threshold, mb);

    VSET((*min), center[X] - radius, center[Y] - radius, center[Z] - radius);
    VSET((*max), center[X] + radius, center[Y] + radius, center[Z] + radius);
    return 0;
}

/**
 * Returns a vector of the orbital element el for planet pl, built from
 * the list
 * @param kl
 * @param pl
 * @param el
 * @return 
 */
gsl_vector* KL_getElements(const ok_list* kl, const int pl, const int el) {
    gsl_vector* v = gsl_vector_alloc(kl->size);
    
    for (int i = 0; i < kl->size; i++)
        VSET(v, i, MGET(kl->kernels[i]->elements, pl, el));
    return v;
}

void
bn_wrt_point_direc(mat_t out, const mat_t change, const mat_t in, const point_t point, const vect_t direc, const struct bn_tol *tol)
{
    static mat_t t1;
    static mat_t pt_to_origin, origin_to_pt;
    static mat_t d_to_zaxis, zaxis_to_d;
    static vect_t zaxis;

    /* build "point to origin" matrix */
    MAT_IDN(pt_to_origin);
    MAT_DELTAS_VEC_NEG(pt_to_origin, point);

    /* build "origin to point" matrix */
    MAT_IDN(origin_to_pt);
    MAT_DELTAS_VEC_NEG(origin_to_pt, point);

    /* build "direc to zaxis" matrix */
    VSET(zaxis, 0.0, 0.0, 1.0);
    bn_mat_fromto(d_to_zaxis, direc, zaxis, tol);

    /* build "zaxis to direc" matrix */
    bn_mat_inv(zaxis_to_d, d_to_zaxis);

    /* apply change matrix...
     * t1 = change * d_to_zaxis * pt_to_origin * in
     */
    bn_mat_mul4(t1, change, d_to_zaxis, pt_to_origin, in);

    /* apply origin_to_pt matrix:
     * out = origin_to_pt * zaxis_to_d *
     * change * d_to_zaxis * pt_to_origin * in
     */
    bn_mat_mul3(out, origin_to_pt, zaxis_to_d, t1);
}

/*
 * Replica of STEP function:
 * FUNCTION build_axes() :
 LIST [3:3] OF direction;
 LOCAL
 d1, d2 : direction;
 END_LOCAL;
 d1 := NVL(normalise(axis), dummy_gri || direction([0.0,0.0,1.0]));
 d2 := first_proj_axis(d1, ref_direction);
 RETURN [d2, normalise(cross_product(d1,d2)).orientation, d1];

 END_FUNCTION;
/////////
*/
void
Axis2Placement3D::BuildAxis() {
    double d1[3] = VINIT_ZERO;
    double d2[3] = VINIT_ZERO;
    double d1Xd2[3] = VINIT_ZERO;

    if (axis == NULL) {
	VSET(d1,0.0,0.0,1.0);
    } else {
	VMOVE(d1,axis->DirectionRatios());
	VUNITIZE(d1);
    }
    if (ref_direction == NULL) {
	FirstProjAxis(d2,d1,NULL);
    } else {
	FirstProjAxis(d2,d1,ref_direction->DirectionRatios());
    }
    VCROSS(d1Xd2,d1,d2);
    VUNITIZE(d1Xd2);
    VMOVE(p[0],d2);
    VMOVE(p[1],d1Xd2);
    VMOVE(p[2],d1);

    return;
}

Create( void *priv, void *context, LWError *err )
{
   AtmosphereData *dat;

   dat = calloc( 1, sizeof( AtmosphereData ));
   if ( !dat ) {
      *err = "Couldn't allocate instance data.";
      return NULL;
   }

   dat->lo      = -4.0;
   dat->hi      = 0.0;
   dat->fa      = 0.5;
   dat->lum     = 0.5;
   dat->opa     = 1.0;
   dat->den     = 1.0;
   dat->res     = 2;
   dat->type    = 1;
   dat->bck     = 1;
   dat->useTxtr = 0;
   dat->march   = 0;
   VSET( dat->col, 1.0 );

   return dat;
}