void harp::spec_desisim::values ( vector_double & data ) const {

  data.resize ( nglobal_ );
  data.clear();

  fitsfile * fp;

  fits::open_read ( fp, path_ );

  // read the object flux

  fits::img_seek ( fp, objhdu_ );

  fits::img_read ( fp, data, false );

  // read the sky flux and sum

  vector_double skyflux ( data.size() );

  fits::img_seek ( fp, skyhdu_ );

  fits::img_read ( fp, skyflux, false );

  fits::close ( fp );

  for ( size_t i = 0; i < data.size(); ++i ) {
    data[i] += skyflux[i];
  }

  return;
}

void harp::spec_sim::inv_variance ( vector_double & data ) const {

  data.resize ( size_ );
  data.clear();

  return;
}

/*---------------------------------------------------------------
	Fills in the representative sector
		field in the gap structure:
  ---------------------------------------------------------------*/
void  CHolonomicND::calcRepresentativeSectorForGap(
	TGap					&gap,
	const poses::CPoint2D	&target,
	const vector_double		&obstacles)
{
    int     sector;
    int     sectors_to_be_wide = round( WIDE_GAP_SIZE_PERCENT * obstacles.size());
	int		TargetSector = direction2sector( atan2(target.y(),target.x()), obstacles.size() );

    if ( (gap.end-gap.ini) < sectors_to_be_wide )
    {
#if	1
		sector = round(0.5f*gap.ini+0.5f*gap.end);
#else
		double	min_dist_obs_near_ini=1, min_dist_obs_near_end=1;
		int		i;
		for ( i= gap.ini;i>=max(0,gap.ini-2);i--)
			min_dist_obs_near_ini = min(min_dist_obs_near_ini, obstacles[i]);
		for ( i= gap.end;i<=min((int)obstacles.size()-1,gap.end+2);i++)
			min_dist_obs_near_end = min(min_dist_obs_near_end, obstacles[i]);
		sector = round((min_dist_obs_near_ini*gap.ini+min_dist_obs_near_end*gap.end)/(min_dist_obs_near_ini+min_dist_obs_near_end));
#endif
    }
    else
    {
        // Para gaps anchos que NO contengan al target, cerca del borde
        //  mas cercano a este:
		//if ( TargetSector < gap.ini || TargetSector > gap.end )
  //      {
            int     dir;
            int     dist_ini = abs( TargetSector - gap.ini );
            int     dist_end = abs( TargetSector - gap.end );
            if (dist_ini<dist_end) {
                    sector = gap.ini;
                    dir = +1; }
            else {
                    sector = gap.end;
                    dir = -1; }
            sector = sector + dir * sectors_to_be_wide/2 ;
    //    }
    //    else
    //    {
    //        // Es un valle ancho con el Target dentro:
    //        // Buscar la maxima "distance" en un rango cerca del target:
    //        int     ini = max( gap.ini, TargetSector - sectors_to_be_wide / 2 );
    //        int     end = min( TargetSector + sectors_to_be_wide / 2, gap.end);

    //        sector = TargetSector;
    //        for (int i = ini;i<=end;i++)
				//if ( obstacles[i] > obstacles[sector] )
				//	sector = i;
    //    }
    }

	keep_max(sector, 0);
	keep_min(sector, (int)obstacles.size()-1);

    gap.representative_sector = sector;
}

/*---------------------------------------------------------------
		getAsVector
---------------------------------------------------------------*/
void CPose2D::getAsVector(vector_double &v) const
{
	v.resize(3);
	v[0]=m_coords[0];
	v[1]=m_coords[1];
	v[2]=m_phi;
}

void fixGPStimestamp(CObservationGPSPtr &obs, vector_double &time_changes, std::map<std::string,double> &DeltaTimes )
{
	if (!obs->has_GGA_datum && !obs->has_RMC_datum) return;

	CObservationGPS::TUTCTime	theTime;
	bool  hasTime=false;

	if (obs->has_GGA_datum && obs->GGA_datum.fix_quality>0)
	{
		theTime = obs->GGA_datum.UTCTime;
		hasTime = true;
	}
	else
	if (obs->has_RMC_datum && obs->RMC_datum.validity_char=='A' )
	{
		theTime = obs->RMC_datum.UTCTime;
		hasTime = true;
	}

	 // The last known delta_time for this sensor name
	if (DeltaTimes.find( obs->sensorLabel )==DeltaTimes.end())
		DeltaTimes[obs->sensorLabel] = 0;

	double &DeltaTime = DeltaTimes[obs->sensorLabel];

	if ( hasTime )
	{
		TTimeParts	timparts;
		mrpt::system::timestampToParts( obs->timestamp, timparts);

		DeltaTime  = 3600*theTime.hour + 60*theTime.minute  + theTime.sec;
		DeltaTime -= 3600*timparts.hour + 60*timparts.minute + timparts.second;

		if (theTime.hour < timparts.hour-2)
		{
			// The GPS time is one day ahead the "timestamp"
			DeltaTime += 3600*24;
		}
		else if (timparts.hour > theTime.hour+2)
		{
			// The "timstamp" is one day ahead the GPS time:
			DeltaTime -= 3600*24;
		}

		// Instead of delta, just replace:
		timparts.hour 		= theTime.hour;
		timparts.minute 	= theTime.minute;
		timparts.second 	= theTime.sec;

		obs->timestamp = buildTimestampFromParts(timparts);
	}
	else
	{
		// Use last delta
		obs->timestamp += mrpt::system::secondsToTimestamp(DeltaTime);
	}

	// Fix timestamp:
	time_changes.push_back( DeltaTime );
}

// The error function F(x):
void myFunction( const vector_double &x, const vector_double &y, vector_double &out_f)
{
	out_f.resize(1);

	// 1-cos(x+1) *cos(x*y+1)
	out_f[0] = 1 - cos(x[0]+1) * cos(x[0]*x[1]+1);
}

void ffff(const vector_double &x,const CQuaternionDouble &Q, vector_double &OUT)
{
	OUT.resize(3);
	CQuaternionDouble q(x[0],x[1],x[2],x[3]);
	q.normalize();
	q.rpy(OUT[2],OUT[1],OUT[0]);
}

void harp::eigen_decompose ( matrix_double const & invcov, vector_double & D, matrix_double & W, bool regularize ) {

  D.resize ( invcov.size1() );
  W.resize ( invcov.size1(), invcov.size2() );
  matrix_double temp ( invcov );

  int nfound;

  boost::numeric::ublas::vector < int > support ( 2 * invcov.size1() );

  boost::numeric::bindings::lapack::syevr ( 'V', 'A', boost::numeric::bindings::lower ( temp ), 0.0, 0.0, 0, 0, 0.0, nfound, D, W, support );

  if ( regularize ) {
    
    double min = 1.0e100;
    double max = -1.0e100;

    for ( size_t i = 0; i < D.size(); ++i ) {
      if ( D[i] < min ) {
        min = D[i];
      }
      if ( D[i] > max ) {
        max = D[i];
      }
    }

    double rcond = min / max;

    // pick some delta that is bigger than machine precision, but still tiny
    double epsilon = 10.0 * std::numeric_limits < double > :: epsilon();

    if ( rcond < epsilon ) {

      double reg = max * epsilon - min;
      //cerr << "REG offset = " << reg << " for min / max = " << min << " / " << max << endl;

      for ( size_t i = 0; i < D.size(); ++i ) {
        D[i] += reg;
      }

    }

  }

  return;
}

void harp::apply_inverse_norm ( vector_double const & S, vector_double & vec ) {

  for ( size_t i = 0; i < vec.size(); ++i ) {
    vec[i] /= S[i];
  }

  return;
}

void harp::spec_desisim::lambda ( vector_double & lambda_vals ) const {

  lambda_vals.resize ( nlambda_ );

  for ( size_t i = 0; i < nlambda_; ++i ) {
    lambda_vals[i] = crval + cdelt * (double)i;
  }

  return;
}

void harp::column_norm ( matrix_double const & mat, vector_double & S ) {

  S.resize( mat.size1() );
  S.clear();

  for ( size_t i = 0; i < mat.size2(); ++i ) {
    for ( size_t j = 0; j < mat.size1(); ++j ) {
      S[ j ] += mat( j, i );
    }
  }

  // Invert

  for ( size_t i = 0; i < S.size(); ++i ) {
    S[i] = 1.0 / S[i];
  }

  return;
}

/*---------------------------------------------------------------
					getHistogramNormalized
 ---------------------------------------------------------------*/
void CHistogram::getHistogramNormalized( vector_double &x, vector_double &hits ) const
{
	const size_t N = m_bins.size();
	linspace(m_min,m_max,N, x);

	hits.resize(N);
	const double K=m_binSizeInv/m_count;
	for (size_t i=0;i<N;i++)
		hits[i]=K*m_bins[i];
}

/** Returns a 1x7 vector with [x y z qr qx qy qz] */
void CPose3DQuat::getAsVector(vector_double &v) const
{
	v.resize(7);
	v[0] = m_coords[0];
	v[1] = m_coords[1];
	v[2] = m_coords[2];
	v[3] = m_quat[0];
	v[4] = m_quat[1];
	v[5] = m_quat[2];
	v[6] = m_quat[3];
}

void harp::sparse_mv_trans ( matrix_double_sparse const & AT, vector_double const & in, vector_double & out ) {

  // FIXME:  for now, we just use the (unthreaded) boost sparse matrix-vector product.  If this
  // operation dominates the cost in any way, we can add a threaded implementation here.

  size_t nrows = AT.size1();
  size_t ncols = AT.size2();

  if ( in.size() != nrows ) {
    std::ostringstream o;
    o << "length of input vector (" << in.size() << ") does not match number of rows in transposed matrix (" << nrows << ")";
    HARP_THROW( o.str().c_str() );
  }

  out.resize ( ncols );

  boost::numeric::ublas::axpy_prod ( in, AT, out, true );

  return;
}

/*---------------------------------------------------------------
	Fills in the representative sector
		field in the gap structure:
  ---------------------------------------------------------------*/
void  CHolonomicND::calcRepresentativeSectorForGap(
	TGap                        & gap,
	const mrpt::math::TPoint2D  & target,
	const vector_double         & obstacles)
{
	int sector;
	const unsigned int sectors_to_be_wide = round( options.WIDE_GAP_SIZE_PERCENT * obstacles.size());
	const unsigned int target_sector = direction2sector( atan2(target.y,target.x), obstacles.size() );

	if ( (gap.end-gap.ini) < sectors_to_be_wide )	//Select the intermediate sector
	{
#if	1
		sector = round(0.5f*gap.ini+0.5f*gap.end);
#else
		double	min_dist_obs_near_ini=1, min_dist_obs_near_end=1;
		int		i;
		for ( i= gap.ini;i>=max(0,gap.ini-2);i--)
			min_dist_obs_near_ini = min(min_dist_obs_near_ini, obstacles[i]);
		for ( i= gap.end;i<=min((int)obstacles.size()-1,gap.end+2);i++)
			min_dist_obs_near_end = min(min_dist_obs_near_end, obstacles[i]);
		sector = round((min_dist_obs_near_ini*gap.ini+min_dist_obs_near_end*gap.end)/(min_dist_obs_near_ini+min_dist_obs_near_end));
#endif
	}
	else	//Select a sector close to the target but spaced "sectors_to_be_wide/2" from it
	{
		unsigned int dist_ini = mrpt::utils::abs_diff(target_sector, gap.ini );
		unsigned int dist_end = mrpt::utils::abs_diff(target_sector, gap.end );

		if (dist_ini > 0.5*obstacles.size())
			dist_ini = obstacles.size() - dist_ini;
		if (dist_end > 0.5*obstacles.size())
			dist_end = obstacles.size() - dist_end;

		int dir;
		if (dist_ini<dist_end) {
				sector = gap.ini;
				dir = +1; }
		else {
				sector = gap.end;
				dir = -1; }

		sector = sector + dir*static_cast<int>(sectors_to_be_wide)/2;
	}

	keep_max(sector, 0);
	keep_min(sector, static_cast<int>(obstacles.size())-1 );

	gap.representative_sector = sector;
}

void harp::spec_sim::values ( vector_double & data ) const {

  double PI = std::atan2 ( 0.0, -1.0 );

  data.resize ( size_ );

  size_t bin = 0;

  size_t halfspace = (size_t)( atmspace_ / 2 );

  for ( size_t i = 0; i < nspec_; ++i ) {

    bool dosky = ( i % skymod_ == 0 ) ? true : false;

    size_t objoff = 2 * i;

    for ( size_t j = 0; j < nlambda_; ++j ) {

      double val = 0.0;

      val += background_ * sin ( 3.0 * (double)j / ( (double)atmspace_ * 2.0 * PI ) ) * sin ( 7.0 * (double)j / ( (double)atmspace_ * 2.0 * PI ) ) * sin ( 11.0 * (double)j / ( (double)atmspace_ * 2.0 * PI ) );

      val += 2.0 * background_;

      if ( ( ( j + halfspace ) % atmspace_ ) == 0 ) {
        val += atmpeak_;
      }

      if ( ! dosky ) {
        if ( ( ( j + objoff ) % objspace_ ) == 0 ) {
          val += objpeak_;
        }
      }

      data[ bin ] = val;

      ++bin;
    }
  }
  
  return;
}

/*---------------------------------------------------------------
						Search the best gap.
  ---------------------------------------------------------------*/
void  CHolonomicND::searchBestGap(
			vector_double		&obstacles,
			double				maxObsRange,
			TGapArray			&in_gaps,
			poses::CPoint2D		&target,
			int					&out_selDirection,
			double				&out_selEvaluation,
			TSituations			&out_situation,
			double				&out_riskEvaluation,
			CLogFileRecord_NDPtr log)
{
	// Para evaluar el risk:
	unsigned int min_risk_eval_sector	= 0;
	unsigned int max_risk_eval_sector	= obstacles.size()-1;
	unsigned int TargetSector		= direction2sector(atan2(target.y(),target.x()),obstacles.size());
	const double TargetDist	= std::max(0.01,target.norm());

    // El "risk" se calcula al final para todos los casos.

    // D1 : Camino directo?
    // --------------------------------------------------------
    const int freeSectorsNearTarget = 10;  // 3
    bool theyAreFree = true, caseD1 = false;
    if (TargetSector>(unsigned int)freeSectorsNearTarget &&
		TargetSector<(unsigned int)(obstacles.size()-freeSectorsNearTarget) )
    {
        for (int j=-freeSectorsNearTarget;j<=freeSectorsNearTarget;j++)
                if (obstacles[ TargetSector + j ]<0.95*TargetDist)
                        theyAreFree = false;
        caseD1 = theyAreFree;
    }

    if (caseD1)
    {
        // S1: Camino libre hacia target:
        out_selDirection	= TargetSector;

		// Si hay mas de una, la que llegue antes
		out_selEvaluation   =	1.0 + std::max( 0.0, (maxObsRange - TargetDist) / maxObsRange );
        out_situation		=	SITUATION_TARGET_DIRECTLY;
    }
    else
    {
        // Evaluar los GAPs (Si no hay ninguno, nada, claro):
        vector_double	gaps_evaluation;
        int				selected_gap		=-1;
        double			selected_gap_eval	= -100;

        evaluateGaps(
                obstacles,
				maxObsRange,
                in_gaps,
                TargetSector,
                TargetDist,
                gaps_evaluation );

		if (log) log->gaps_eval = gaps_evaluation;

        // D2: Hay algun gap que pase por detras del target?
        //   ( y no este demasiado lejos):
        // -------------------------------------------------
        for ( unsigned int i=0;i<in_gaps.size();i++ )
			if ( in_gaps[i].maxDistance >= TargetDist &&
				 abs((int)(in_gaps[i].representative_sector-(int)TargetSector)) <= (int)floor(MAX_SECTOR_DIST_FOR_D2_PERCENT * obstacles.size()) )
					if ( gaps_evaluation[i]>selected_gap_eval )
					{
						selected_gap_eval = gaps_evaluation[i];
						selected_gap = i;
					}


        // Coger el mejor GAP:
        //  (Esto ya solo si no se ha cogido antes)
        if ( selected_gap==-1 )
			for ( unsigned int i=0;i<in_gaps.size();i++ )
				if ( gaps_evaluation[i]>selected_gap_eval )
				{
						selected_gap_eval = gaps_evaluation[i];
						selected_gap = i;
				}
        //  D3:  No es suficientemente bueno? ( o no habia ninguno?)
        // ------------------------------------------------------------
        if ( selected_gap_eval <= 0 )
        {
            // S2: No way found
            // ------------------------------------------------------
            out_selDirection	= 0;
            out_selEvaluation	= 0.0f; // La peor
            out_situation		= SITUATION_NO_WAY_FOUND;
        }
        else
        {
            // El seleccionado:
            TGap    gap = in_gaps[selected_gap];

            int     sectors_to_be_wide = round( WIDE_GAP_SIZE_PERCENT * obstacles.size() );

//.........这里部分代码省略.........

/*---------------------------------------------------------------
						Navigate
  ---------------------------------------------------------------*/
void  CHolonomicND::navigate(
				poses::CPoint2D	&target,
				vector_double	&obstacles,
				double			maxRobotSpeed,
				double			&desiredDirection,
				double			&desiredSpeed,
				CHolonomicLogFileRecordPtr &logRecord)
{
	TGapArray			gaps;
	TSituations			situation;
	int					selectedSector;
	double				riskEvaluation;
	CLogFileRecord_NDPtr log;
	double				evaluation;

	// Create a log record for returning data.
	if (!logRecord.present())
	{
		log = CLogFileRecord_ND::Create();
		logRecord = log;
	}


	// Search gaps:
    gaps.clear();
	gapsEstimator(	obstacles,
					target,
					gaps );


	// Select best gap:
	searchBestGap(	obstacles,
					1.0f,
					gaps,
					target,
					selectedSector,
					evaluation,
					situation,
					riskEvaluation,
					log);

	if (situation == SITUATION_NO_WAY_FOUND)
	{
		// No way found!
		desiredDirection = 0;
		desiredSpeed = 0;
	}
	else
	{
		// A valid movement:
		desiredDirection = (double)(M_PI*(-1 + 2*(0.5f+selectedSector)/((double)obstacles.size())));

		// Speed control: Reduction factors
		// ---------------------------------------------
		double		targetNearnessFactor = max(0.20, min(1.0, 1.0-exp(-(target.norm()+0.01)/TARGET_SLOW_APPROACHING_DISTANCE)));
		//printf(" TARGET NEARNESS = %f\n",targetNearnessFactor);
		double		riskFactor = min(1.0, riskEvaluation / RISK_EVALUATION_DISTANCE );

		desiredSpeed = maxRobotSpeed * min(riskFactor,targetNearnessFactor);
	}

	last_selected_sector = selectedSector;

	// LOG --------------------------
	if (log)
	{
		// gaps:
		if (situation != SITUATION_TARGET_DIRECTLY )
		{
			int	i,n = gaps.size();
			log->gaps_ini.resize(n);
			log->gaps_end.resize(n);
			for (i=0;i<n;i++)
			{
				log->gaps_ini[i]  = gaps[i].ini;
				log->gaps_end[i]  = gaps[i].end;
			}
		}
        // Selection:
        log->selectedSector = selectedSector;
        log->evaluation = evaluation;
        log->situation = situation;
        log->riskEvaluation = riskEvaluation;
	}
}

/*---------------------------------------------------------------
						Find gaps in the obtacles.
  ---------------------------------------------------------------*/
void  CHolonomicND::gapsEstimator(
		vector_double		&obstacles,
		poses::CPoint2D		&target,
		TGapArray			&gaps_out )
{
	unsigned int	i,n;
	int				nMaximos=0;
    double			MaximoAbsoluto = -100;
	double			MinimoAbsoluto = 100;
    vector_int		MaximoIdx;
    vector_double	MaximoValor;

    // Hacer una lista con los maximos de las distancias a obs:
    // ----------------------------------------------------------
	MaximoIdx.resize(obstacles.size());
	MaximoValor.resize(obstacles.size());
	n = obstacles.size();

    for (i=1;i<(n-1);i++)
    {
		// Actualizar max. y min. absolutos:
		MaximoAbsoluto= max( MaximoAbsoluto, obstacles[i] );
		MinimoAbsoluto= min( MinimoAbsoluto, obstacles[i] );

		// Buscar maximos locales:
		if ( ( obstacles[i] >= obstacles[i+1] &&
			  obstacles[i] > obstacles[i-1] ) ||
			  ( obstacles[i] > obstacles[i+1] &&
			  obstacles[i] >= obstacles[i-1] ) )
		{
				MaximoIdx[nMaximos] = i;
				MaximoValor[nMaximos++] = obstacles[i];
		}
    }

    //  Crear GAPS:
    // --------------------------------------------------------
	TGapArray    gaps_temp;
   	gaps_temp.reserve( 150 );

	for (double factorUmbral = 0.975f;factorUmbral>=0.04f;factorUmbral-=0.05f)
	{
            double   umbral = factorUmbral* MaximoAbsoluto + (1.0f-factorUmbral)*MinimoAbsoluto;
			bool	dentro = false;
			int		sec_ini=0, sec_end;
			double	maxDist=0;

			for (i=0;i<n;i++)
			{
				if ( !dentro && (!i || obstacles[i]>=umbral) )
				{
					sec_ini = i;
					maxDist = obstacles[i];
					dentro = true;
				}
				else if (dentro && (i==(n-1) || obstacles[i]<umbral ))
				{
					sec_end = i;
					dentro = false;

					if ( (sec_end-sec_ini) > 2 )
					{
						// Add new gap:
						TGap	newGap;
						newGap.ini				= sec_ini;
						newGap.end				= sec_end;
						newGap.entranceDistance = min( obstacles[sec_ini], obstacles[sec_end] );
						newGap.maxDistance		= maxDist;

						gaps_temp.push_back(newGap);
					}
				}

				if (dentro) maxDist = max( maxDist, obstacles[i] );
			}
	}

    // Proceso de eliminacion de huecos redundantes:
    // -------------------------------------------------------------
	std::vector<bool>	borrar_gap;
	borrar_gap.resize( gaps_temp.size() );
        for (i=0;i<gaps_temp.size();i++)
			borrar_gap[i] = false;


    // Eliminar huecos con muy poca profundidad si estan dentro de otros:
	double	maxProfundidad = 0;
	for (i=0;i<gaps_temp.size();i++)
    {
		double profundidad =
				gaps_temp[i].maxDistance -
				gaps_temp[i].entranceDistance;
		maxProfundidad = max(maxProfundidad, profundidad);
	}

	for (i=0;i<gaps_temp.size();i++)
    {
//.........这里部分代码省略.........

// extracted from Particles3Dcomm.cpp
//
void Collective::read_particles_restart(
    const VCtopology3D* vct,
    int species_number,
    vector_double& u,
    vector_double& v,
    vector_double& w,
    vector_double& q,
    vector_double& x,
    vector_double& y,
    vector_double& z,
    vector_double& t)const
{
#ifdef NO_HDF5
  eprintf("Require HDF5 to read from restart file.");
#else
    if (vct->getCartesian_rank() == 0 && species_number == 0)
    {
      printf("LOADING PARTICLES FROM RESTART FILE in %s/restart.hdf\n",
        getRestartDirName().c_str());
    }
    stringstream ss;
    ss << vct->getCartesian_rank();
    string name_file = getRestartDirName() + "/restart" + ss.str() + ".hdf";
    // hdf stuff
    hid_t file_id, dataspace;
    hid_t datatype, dataset_id;
    herr_t status;
    size_t size;
    hsize_t dims_out[1];        /* dataset dimensions */
    int status_n;

    // open the hdf file
    file_id = H5Fopen(name_file.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    if (file_id < 0) {
      eprintf("couldn't open file: %s\n"
        "\tRESTART NOT POSSIBLE", name_file.c_str());
      //cout << "couldn't open file: " << name_file << endl;
      //cout << "RESTART NOT POSSIBLE" << endl;
    }


    //find the last cycle
    int lastcycle=0;
    dataset_id = H5Dopen2(file_id, "/last_cycle", H5P_DEFAULT); // HDF 1.8.8
    status = H5Dread(dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &lastcycle);
    status = H5Dclose(dataset_id);

    stringstream species_name;
    species_name << species_number;

    ss.str("");ss << "/particles/species_" << species_number << "/x/cycle_" << lastcycle;
    dataset_id = H5Dopen2(file_id, ss.str().c_str(), H5P_DEFAULT); // HDF 1.8.8
    datatype = H5Dget_type(dataset_id);
    size = H5Tget_size(datatype);
    dataspace = H5Dget_space(dataset_id); /* dataspace handle */
    status_n = H5Sget_simple_extent_dims(dataspace, dims_out, NULL);

    // get how many particles there are on this processor for this species
    status_n = H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
    const int nop = dims_out[0]; // number of particles in this process
    //Particles3Dcomm::resize_SoA(nop);
    {
      //
      // allocate space for particles including padding
      //
      const int padded_nop = roundup_to_multiple(nop,DVECWIDTH);
      u.reserve(padded_nop);
      v.reserve(padded_nop);
      w.reserve(padded_nop);
      q.reserve(padded_nop);
      x.reserve(padded_nop);
      y.reserve(padded_nop);
      z.reserve(padded_nop);
      t.reserve(padded_nop);
      //
      // define size of particle data
      //
      u.resize(nop);
      v.resize(nop);
      w.resize(nop);
      q.resize(nop);
      x.resize(nop);
      y.resize(nop);
      z.resize(nop);
      t.resize(nop);
    }
    // get x
    status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &x[0]);
    // close the data set
    status = H5Dclose(dataset_id);

    // get y
    ss.str("");ss << "/particles/species_" << species_number << "/y/cycle_" << lastcycle;
    dataset_id = H5Dopen2(file_id, ss.str().c_str(), H5P_DEFAULT); // HDF 1.8.8
    status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &y[0]);
    status = H5Dclose(dataset_id);

    // get z
//.........这里部分代码省略.........